dictionary_extend_reduced_data.m

function reduced_data =  dictionary_extend_reduced_data(model,reduced_data,detailed_data)
%function reduced_data =  dictionary_extend_reduced_data(model,reduced_data,detailed_data)
%
% function extending reduced_data if single new RB has been added
% in detailed_data.RB
%
%% The following also works, but results in slow linear combination
%and repetitive linear combinations over N extension loop.
%    reduced_data= lin_stat_gen_reduced_data(model,detailed_data);

% B. Haasdonk 13.3.2015

A_comp = detailed_data.A_comp;
f_comp = detailed_data.f_comp;
Qa = length(A_comp);
%reduced_data.Qa = Qa;
Qf = length(f_comp);
%reduced_data.Qf = Qf;

%reduced_data.AN_comp = cell(1,length(A_comp));
for q = 1:Qa
  AN_comp_vec1 = detailed_data.RB'*A_comp{q}*detailed_data.RB(:,end);
  AN_comp_vec2 = detailed_data.RB(:,end)'*A_comp{q}*detailed_data.RB;
  reduced_data.AN_comp{q}  = ...
      [reduced_data.AN_comp{q}, AN_comp_vec1(1:end-1);
       AN_comp_vec2(1:end)];
end;

%reduced_data.fN_comp = cell(1,length(f_comp));
for q = 1:Qf
  reduced_data.fN_comp{q} = [reduced_data.fN_comp{q} ;
                    detailed_data.RB(:,end)' * f_comp{q}];
end;

if model.compute_output_functional
  % assumption: nonparametic output functional, then simple RB
  % evaluation possible  
  l_comp = ...
      model.operators_output(model,detailed_data);
  Q_l = length(l_comp);
  reduced_data.lN_comp = cell(1,Q_l);
  for q = 1:Q_l
    reduced_data.lN_comp{q} = detailed_data.RB' * l_comp{q}; 
  end;
end;

N = model.get_rb_size(model,detailed_data);
reduced_data.N = N;

% plus error estimation quantities

% G = (v_r^q, v_r^q) = v_r^q' * K * v_r^q
% with {v_r^q}_q = (v_f^q, v_a^qq')_{qq'}
% G = [Gff, Gfa; Gfa', Gaa];
%
% (v_f^q,v_f^q') = v_f^q' * K * v_f^q

% matrices of coefficient vectors of Riesz-representers:
% K * v_f^q = f^q      (coefficient vector equations)
K = model.get_inner_product_matrix(detailed_data);
% search solution in H10, i.e. set dirichlet DOFs

%v_f = zeros(size(detailed_data.RB,1),Qf);
K_v_f = reduced_data.K_v_f;
v_f = reduced_data.v_f;
%K_v_f = zeros(size(detailed_data.RB,1),Qf);
%for q = 1:Qf
%  K_v_f(:,q) = f_comp{q};
%  v_f = K \ K_v_f;
%  %  v_f(:,q) = K \ f_comp{q};
%end;

K_v_a_N = zeros(size(K,1),Qa);
v_a_N = zeros(size(K,1),Qa);
%v_a = cell(N,1);
%K_v_a = cell(N,1);
%for n = 1:N
%  K_v_a{n} = zeros(size(K,1),Qa);
%  v_a{n} = zeros(size(K,1),Qa);
%  for q = 1:Qa
%    K_v_a{n}(:,q) = (A_comp{q}*detailed_data.RB(:,n));
%    v_a{n} = K \ K_v_a{n};
%    %    v_a{n}(:,q) = K \ K_v_a{n}(:,q);
%  end;
%end;
for q = 1:Qa
  K_v_a_N(:,q) = A_comp{q}*detailed_data.RB(:,end);
end;
v_a_N = K \ K_v_a_N;

reduced_data.K_v_a{N} = K_v_a_N; 
reduced_data.v_a{N} = v_a_N; 

v_a = reduced_data.v_a;
K_v_a = reduced_data.K_v_a;

% finally assemble G = [G_ff, G_fa{n}; G_fa{n}', G_aa{n}];
%Q_r = Qf + N * Qa;
%G = zeros(Q_r,Q_r);
%%G(1:Qf,1:Qf) = reduced_data.Gff;
%G(1:Qf,1:Qf) = v_f' * K_v_f;
%for n = 1:N
%  G(1:Qf,Qf+(n-1)*Qa +(1:Qa)) = ...
%      v_f' * K_v_a{n};
%  %      reduced_data.Gfa{n};
%  G(Qf+(n-1)*Qa +(1:Qa),1:Qf) = ...
%      transpose(v_f' * K_v_a{n});
%  %      reduced_data.Gfa{n}';
%end;
%for n1 = 1:N
%  for n2 = 1:N
%    G(Qf+(n1-1)*Qa+(1:Qa),Qf+(n2-1)*Qa +(1:Qa)) = ...
%       v_a{n1}' * K_v_a{n2};
%    %  reduced_data.Gaa{n1,n2};
%  end;
%end;

% assemble Qf x Qf matrix Hff = (<vf1,vf1>, ... ,<vf1,vfQf>)
%                                      .........
%                               (<vfQf,vf1>, ... ,<vfQf,vfQf>)

%reduced_data.bar_Hff = v_f' * K_v_f;

% assemble N x QaQf  matrix 
%          Haf = (<vf1,va11>...<vf1,va1Qa> | <vf2,va11>  ... <vf2,va1Qa> | ... <vfQf,va1Qa>)
%                 <vf1,va21>...<vf1,va2Qa> | <vf2,va21>  ... <vf2,va2Qa> | ... <vfQf,va2Qa>)
%                      ............
%                 <vf1,vaN1>...<vf1,vaNQa> | <vf2,vaN1>  ... <vf2,vaNQa> | ... <vfQf,vaNQa>)

%Haf = zeros(N,Qa * Qf);
for q = 1:Qf
  %  for n = 1:N
  %    % insert row vector in block matrix:
  %    vaf = v_f(:,q)' * K_v_a{n};
  %    Haf(n,(1+(q-1)*Qa):(q*Qa)) = vaf;
  % insert row vector in block matrix:
  vaf = v_f(:,q)' * K_v_a_N;
  reduced_data.bar_Haf(N,(1+(q-1)*Qa):(q*Qa)) = vaf;
end;
Haf = reduced_data.bar_Haf;

% assemble N^2 x Qa^2  matrix  
%          Haa = (<va11,va11>...<va11,va1Qa> | <va12,va11>  ... <va12,va1Qa> | <va1Qa,va11>... <va1Qa,va1Qa>)
%                      ............
%                 <vaN1,va11>...<vaN1,va1Qa> | <vaN2,va11>  ... <vaN2,va1Qa> | <vaNQa,va11>... <vaNQa,va1Qa>)
%                 -------------------------------------------------------------------------------------------
%                 <va11,va21>...<va11,va2Qa> | <va12,va21>  ... <va12,va2Qa> | <va1Qa,va21>... <va1Qa,va2Qa>
%                      ............
%                 <vaN1,va21>...<vaN1,va2Qa> | <vaN2,va21>  ... <vaN2,va2Qa> | <vaNQa,va21>... <vaNQa,va2Qa>)
%                 -------------------------------------------------------------------------------------------
%                                                       ............
%                 -------------------------------------------------------------------------------------------
%                 <va11,vaN1>...<va11,vaNQa> | <va12,vaN1>  ... <va12,vaNQa> | <va1Qa,vaN1>... <va1Qa,vaNQa>
%                      ............
%                 <vaN1,vaN1>...<vaN1,vaNQa> | <vaN2,vaN1>  ... <vaN2,vaNQa> | <vaNQa,vaN1>... <vaNQa,vaNQa> )

%% old expensive rule of computation:
%Haa = zeros(N*N,Qa*Qa);
%% insert row vectors in block matrix:
%for q = 1:Qa % loop over column-blocks
%  for n1 = 1:N % loop over block-rows
%    for n2 = 1:N % loop over rows within one block;
%      vaa = v_a{n2}(:,q)' * K_v_a{n1};  
%      Haa((n1-1)*N + n2, (1+(q-1)*Qa):q*Qa) = vaa;      
%    end;
%  end;
%end;

% new rule
bar_Haa = zeros(N*N,Qa*Qa);
%distribute old matrix into new one:
dummy = ones(N,N);
dummy(N,:) = 0; 
dummy(:,N) = 0; 
new_row_ind = find(dummy);
bar_Haa(new_row_ind,:) = reduced_data.bar_Haa;

% insert last row in each block-row:
n2 = N;
for n1 = 1:N % loop over block-rows
  vaa = transpose(v_a_N' * K_v_a{n1});  
  bar_Haa((n1-1)*N + n2, :) = vaa(:);      
end;
% insert last block-row:
n1 = N; % loop over block-rows
for n2 = 1:N-1 % loop over rows within last block row;
  for q = 1:Qa % loop over column-block
    vaa = transpose(v_a{n2}' * K_v_a_N);  
    bar_Haa((n1-1)*N + n2, :) = vaa(:);      
  end;
end;

%if ~isequal(bar_Haa,Haa)
%  disp('bar_Haa not equal Haa, please check!');
%  keyboard;
%end;

reduced_data.bar_Haa = bar_Haa;

reduced_data.used_dictionary_gen_reduced_data = 1;