lin_evol_rb_operators_primal_dual.m

Go to the documentation of this file.

function [LL_I, LL_E, bb, K_II, K_IE, K_EE, m_I, m_E, m, K_IdId, K_IdE, LL_I_correct, LL_E_correct, bb_correct, scm_offline_data] = ...
    lin_evol_rb_operators_primal_dual(model, detailed_data)
%[LL_I, LL_E, bb, K_II, K_IE, K_EE, m_I, m_E, m, K_IdId, K_IdE, LL_I_correct, LL_E_correct, bb_correct, scm_offline_data] = ...
%    lin_evol_rb_operators_primal_dual(model, detailed_data)
%
% function computing the time-dependent reduced basis operators and
% vectors. 
%
% Function supports affine decomposition, i.e. different operation modes
% guided by optional field decomp_mode in params. See also the
% contents.txt for general explanation
%
% Required fields of model
% operators_algorithm: name of function for computing the
%                 L_E,L_I,b-operators with arguments (grid, params)
%                 example fv_operators_implicit
%
% 
% Optional fields of model:
%   mu_names : names of fields to be regarded as parameters in vector mu
%   decomp_mode: operation mode of the function
%     -0= 'complete' (default): no parameter dependence or decomposition is 
%                performed. output is as described above.
%     -1= 'components': For each output argument a cell array of output
%                 arguments is returned representing the q-th component
%                 independent of the parameters given in mu_names  
%     -2= 'coefficients': For each output argument a cell array of output
%                 arguments is returned representing the q-th coefficient
%                 dependent of the parameters given in mu_names  
%
% In 'coefficients' mode the detailed data is empty.
%
% Bernard Haasdonk 23.7.2006
%
% extended to primal/dual formulaton by Dominik Garmatter:
%
% see lin_evol_gen_reduced_data_primal_dual.m for the various cases that
% are included in this function. scm_offline_data are only generated if
% decomp_mode = 1, because you want to produce those data in your
% reduced_data and in gen_reduced_data the
% rb_operators are always called with decomp_mode = 1.
%
% special Note: in the european_option_pricing model a functional named
% 'theta' includes a partial timederivative and therefore produces a
% slightly different dual problem. Therefore the 'theta'-case is sometimes
% treated differently. Don't care about this case if you don't use the
% european_option_pricing model.


% Dominik Garmatter 20.07 2012

bb = [];
m_I = [];
m_E = [];
m = [];
bb_correct = [];
LL_I_correct = [];
LL_E_correct = [];
K_IdId = [];
K_IdE = [];

decomp_mode = model.decomp_mode;

if decomp_mode == 0 % complete: simple projection on RB
    
% if you want the improved output you need both detailed_data, which
% should be connected in the current 'detailed_data'. This function splits
% them, so the operators can work.
if model.want_improved_output && model.want_dual == 0
    [detailed_data, detailed_data_dual] = lin_evol_split_detailed_data(detailed_data);
end

if isfield(detailed_data, 'L_I_comp') && isfield(detailed_data, 'L_E_comp')
% use the components given in detailed_data to assemble the dual or primal
% components of the respective evolution scheme. The required coefficients
% are now computed.
  old_mode = model.decomp_mode;
  model.decomp_mode = 2;
  if model.want_dual % dual
    [~, ~, ~, L_I_coeff, L_E_coeff] = model.operators_ptr(model, detailed_data);
    if strcmp(model.name_output_functional, 'theta') == 1
      b_coeff = model.operators_output(model, detailed_data);
      b = lincomb_sequence(detailed_data.b_comp, b_coeff);
    end
  else % primal
    [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model,detailed_data);
    b = lincomb_sequence(detailed_data.b_comp, b_coeff);
  end
  model.decomp_mode = old_mode;
  L_I = lincomb_sequence(detailed_data.L_I_comp, L_I_coeff);
  L_E = lincomb_sequence(detailed_data.L_E_comp, L_E_coeff);
  
else % there are no components in detailed_data compute the quantities (decomp_mode = 0)
  if model.want_dual % dual
    [~, ~, ~, L_I, L_E] = model.operators_ptr(model, detailed_data);
    if strcmp(model.name_output_functional, 'theta') == 1
        b = model.operators_output(model, detailed_data);
    end
  else % primal
    [L_I, L_E, b] = model.operators_ptr(model,detailed_data);
  end
end
   
% First compute Quantities that are required in every case. But remember
% the included matrices are different (primal or dual)

  % apply all operator contributions to RB:
  L_E_RB = L_E * detailed_data.RB;
  L_I_RB = L_I * detailed_data.RB;
  
  % fill output quantities
  %    LL_I    :  matrix 
  A = detailed_data.W;
  
  LL_I = detailed_data.RB' * A * L_I_RB;
  
  %    LL_E    :  matrix 
  LL_E = detailed_data.RB' * A * L_E_RB;
  
  %    K_II   : matrix 
  K_II = L_I_RB' * A * L_I_RB;
 
  %    K_IE   : matrix 
  K_IE =  L_I_RB' * A *  L_E_RB;
    
  %    K_EE   :  matrix 
  K_EE =  L_E_RB' * A * L_E_RB;
  
  if model.want_dual % additional requirements for the dual residual in the 0 -> 1 timestep
    L_Id_RB = speye(size(L_I)) * detailed_data.RB;
    K_IdId = L_Id_RB' * A * L_Id_RB;
    K_IdE = L_Id_RB' * A * L_E_RB;
  end
  
% the Quantities that come into play when an inhomogeneity is involved (b
% in the primal case or v_l in the dual theta case in the second timestep)
  if ~model.want_dual  || strcmp(model.name_output_functional, 'theta') == 1
    %    bb      :  vector 
    bb = detailed_data.RB' * A * b;
    
    %    m_I    :  vector 
    m_I  =  L_I_RB' * A * b;
  
    %    m_E    :  vector 
    m_E  = L_E_RB' * A * b;

    %    m      :  scalar
    m  =        b' * A * b;
    if model.want_improved_output && model.want_dual == 0
      % These Quantities are required for the correction term used in the
      % improved output
      bb_correct = b' * A * detailed_data_dual.RB;
      LL_I_correct = L_I_RB' * A * detailed_data_dual.RB;
      LL_E_correct = L_E_RB' * A * detailed_data_dual.RB;
    end
  end
  
elseif decomp_mode == 1
    
% if you want the improved output you need both detailed_data, which
% should be connected in the current 'detailed_data'. This function splits
% them, so the operators can work.
if model.want_improved_output && model.want_dual == 0
    [detailed_data, detailed_data_dual] = lin_evol_split_detailed_data(detailed_data);
end

if isfield(detailed_data, 'L_I_comp') && isfield(detailed_data, 'L_E_comp')
% simply extract the components out of the detailed_data
  if model.want_dual
    if strcmp(model.name_output_functional, 'theta') == 1
      b = detailed_data.b_comp;
    end
  else
    b = detailed_data.b_comp;
  end
  L_I = detailed_data.L_I_comp;
  L_E = detailed_data.L_E_comp;
  
else % there are no components in detailed_data compute them (decomp_mode = 1)
  if model.want_dual % dual
    [~, ~, ~, L_I, L_E] = model.operators_ptr(model, detailed_data);
    if strcmp(model.name_output_functional, 'theta') == 1
        b = model.operators_output(model, detailed_data);
    end
  else % primal
    [L_I, L_E, b] = model.operators_ptr(model,detailed_data);
  end
end

% in gen_reduced_data the rb_operators are always called with
% decomp_mode = 1. Therefore only in this case the computation
% of scm_offline_data can be required. Generate them here if desired.
if ~isfield(model, 'use_scm')
    model.use_scm = 0;
end
if model.use_scm
    % these 2 fields are required in detailed_data for the scm to work!
    if ~isfield(detailed_data, 'L_I_comp')
        detailed_data.L_I_comp = L_I;
    end
    if ~isfield(detailed_data, 'L_E_comp')
        detailed_data.L_E_comp = L_E;
    end
  scm_offline_data = scm_offline(model, detailed_data);
end
  
  Q_L_I = size(L_I,1);
  Q_L_E = size(L_E,1);

  Nmax = size(detailed_data.RB,2);
  H = size(detailed_data.RB,1);
  A = detailed_data.W;

% First compute Quantities that are required in every case. But remember
% the included matrices are different (primal or dual)

  % matrices:
  K_II = cell(Q_L_I*Q_L_I,1);
  K_II(:) = {zeros(Nmax,Nmax)};
  K_IE = cell(Q_L_I*Q_L_E,1);
  K_IE(:) = {zeros(Nmax,Nmax)};
  K_EE = cell(Q_L_E*Q_L_E,1);
  K_EE(:) = {zeros(Nmax,Nmax)};
  LL_I = cell(Q_L_I,1);
  LL_I(:) = {zeros(Nmax,Nmax)};
  LL_E = cell(Q_L_E,1);
  LL_E(:) = {zeros(Nmax,Nmax)};
  
  % apply all operator contributions to RB:
  L_E_RB = cell(Q_L_E,1);
  L_E_RB(:) = {zeros(H,Nmax)};
  L_I_RB = cell(Q_L_I,1);
  L_I_RB(:) = {zeros(H,Nmax)};
  for q=1:Q_L_E
    L_E_RB{q}(:,:)= L_E{q} * detailed_data.RB;
  end;
  for q=1:Q_L_I
    L_I_RB{q}(:,:)= L_I{q} * detailed_data.RB;
  end;  
  
  if model.want_dual % additional requirements for the dual residual in the 0 -> 1 timestep
    L_Id_RB = cell(1,1);
    L_Id_RB(:) = {zeros(H,Nmax)};
    K_IdId = cell(1,1);
    K_IdId(:) = {zeros(Nmax,Nmax)};
    K_IdE = cell(1*Q_L_E,1);
    K_IdE(:) = {zeros(Nmax,Nmax)};
  
    L_Id_RB{1}(:,:) = speye(H) * detailed_data.RB;
    K_IdId{1}(:,:) = L_Id_RB{1}' * A * L_Id_RB{1};
    for q2 = 1:Q_L_E
      K_IdE{q2} = L_Id_RB{1}' * A * L_E_RB{q2};
    end
  end 
  
  % fill output quantities
  %    LL_I    :  matrices 
  for q = 1:Q_L_I;
    LL_I{q} = detailed_data.RB' * A * L_I_RB{q};
  end
  
  %    LL_E    :  matrices 
  for q = 1:Q_L_E;
    LL_E{q} =  detailed_data.RB' * A * L_E_RB{q};
  end


  %    K_II   : matrices 
  for q1 = 1:Q_L_I;
    for q2 = 1:Q_L_I;
      K_II{ (q1-1)*Q_L_I + q2}(:,:) = ... 
          L_I_RB{q1}' * A * L_I_RB{q2};
    end;
  end;
  
  %    K_IE   : time sequence of matrices 
  % K_IE{1} = Gram matrix of L_I_RB{1} and L_E_RB{1}
  % K_IE{2} = Gram matrix of L_I_RB{1} and L_E_RB{2} 
  % so L_E is "inner loop". This must fit to the sigma-ordering
  for q1 = 1:Q_L_I;
    for q2 = 1:Q_L_E;
      K_IE{ (q1-1)*Q_L_E + q2} = ...
          L_I_RB{q1}' * A * L_E_RB{q2};
    end;
  end;
  
  %    K_EE   :  matrices 
  for q1 = 1:Q_L_E;
    for q2 = 1:Q_L_E;
      K_EE{ (q1-1)*Q_L_E + q2} = ...
          L_E_RB{q1}' * A * L_E_RB{q2};
    end;
  end;
  
% the Quantities that come into play when an inhomogeneity is involved (b
% in the primal case or v_l in the dual theta case in the second timestep)
  if ~model.want_dual || strcmp(model.name_output_functional, 'theta') == 1
    Q_b = length(b); 
    % init output quantities
    % scalars:
    m = cell(Q_b*Q_b,1);
    m(:) = {zeros(1,1)};
    % vectors:
    m_I = cell(Q_L_I*Q_b,1);
    m_I(:) = {zeros(Nmax,1)};
    m_E = cell(Q_L_E*Q_b,1);
    m_E(:) = {zeros(Nmax,1)};
    bb = cell(Q_b,1);
    bb(:) = {zeros(Nmax,1)};
    
    %    bb      :  vectors 
    for q = 1:Q_b;
        bb{q} = detailed_data.RB' * A * b{q};
    end;
  
    %    m_I    :  vectors 
    for q1 = 1:Q_L_I;
      for q2 = 1:Q_b;
        m_I{ (q1-1)*Q_b + q2} = ...
        L_I_RB{q1}' * A *  b{q2};
      end;
    end;
  
    %    m_E    :  vectors 
    for q1 = 1:Q_L_E;
      for q2 = 1:Q_b;
        m_E{ (q1-1)*Q_b + q2} = ...
        L_E_RB{q1}' * A * b{q2};
      end;
    end;
  
    %    m      :  scalars 
    for q1 = 1:Q_b;
      for q2 = 1:Q_b;
        m{ (q1-1)*Q_b + q2} = ...
        b{q1}' * A * b{q2};
      end;
    end
  end
  
  if model.want_improved_output && model.want_dual == 0
      % These Quantities are required for the correction term used in the
      % improved output (only in the primal case available)
      LL_I_correct = cell(Q_L_I,1);
      LL_I_correct(:) = {zeros(Nmax,Nmax)};
      LL_E_correct = cell(Q_L_E,1);
      LL_E_correct(:) = {zeros(Nmax,Nmax)};
      bb_correct = cell(Q_b,1);
      bb_correct(:) = {zeros(Nmax,1)};
      
      for q = 1:Q_L_I;
        LL_I_correct{q} = L_I_RB{q}' * A * detailed_data_dual.RB;
      end
      
      for q = 1:Q_L_E;
        LL_E_correct{q} =  L_E_RB{q}' * A * detailed_data_dual.RB;
      end
      
      for q = 1:Q_b;
        bb_correct{q} = b{q}' * A * detailed_data_dual.RB;
      end
  end

else % decomp_mode== 2 -> coefficients

% computing the coeffs dependant on the controlfield 'model.want_dual'
% 1 = dual Operators are computet, 0 = primal ones are.
  if model.want_dual
    [~, ~, ~, L_I, L_E] = model.operators_ptr(model,[]);
    if strcmp(model.name_output_functional, 'theta') == 1
      b = model.operators_output(model, detailed_data);
    end
  else
    [L_I, L_E, b] = model.operators_ptr(model,[]);
  end
  
% First compute Quantities that are required in every case. But remember
% the included coefficients can be different (primal or dual)

  Q_L_I = length(L_I(:));
  Q_L_E = length(L_E(:));
  
  % rb-operator coefficients can simply be forwarded
  LL_E = L_E;
  LL_I = L_I;
  
  % pairs: products of sigmas:
  K_II = repmatrows(L_I(:), Q_L_I) .*  repmat(L_I(:), Q_L_I, 1); 
  K_IE = repmatrows(L_I(:), Q_L_E) .*  repmat(L_E(:), Q_L_I, 1); 
  K_EE = repmatrows(L_E(:), Q_L_E) .*  repmat(L_E(:), Q_L_E, 1);
  
  if model.want_dual % additional requirements for the dual residual in the 0 -> 1 timestep: (remember L_I ist simply the Identity there)
    L_Id   = 1;
    Q_L_Id = length(L_Id(:));
    K_IdId = repmatrows(L_Id(:), Q_L_Id) .*  repmat(L_Id(:), Q_L_Id, 1); 
    K_IdE  = repmatrows(L_Id(:), Q_L_E) .*  repmat(L_E(:), Q_L_Id, 1); 
  end

% the Quantities that come into play when an inhomogeneity is involved (b
% in the primal case or v_l in the dual theta case in the second timestep)
  if ~model.want_dual || strcmp(model.name_output_functional, 'theta') == 1
    Q_b = length(b(:));
    bb = b;
    m_I  = repmatrows(L_I(:), Q_b)   .*  repmat(b(:)  , Q_L_I, 1); 
    m_E  = repmatrows(L_E(:), Q_b)   .*  repmat(b(:)  , Q_L_E, 1); 
    m    = repmatrows(b(:)  , Q_b)   .*  repmat(b(:)  , Q_b  , 1);
    if model.want_improved_output && model.want_dual == 0
      % These Quantities are required for the correction term used in the
      % improved output (only in the primal case available)
      LL_I_correct = L_I;
      LL_E_correct = L_E;
      bb_correct = b;
    end
  end

end %end if