lin_evol_opt_rb_operators_separate_bases.m

function [LL_I, LL_E, bb, M_E, M_b, M_EE, M_Eb, M_bb, M_I, M_II, M_IE, M_Ib, ...
     varargout] = lin_evol_opt_rb_operators_separate_bases(model, detailed_data)
%function [LL_I, LL_E, bb, M_E, M_b, M_EE, M_Eb, M_bb, M_I, M_II, M_IE, M_Ib, ...
%    {M_EdEd, M_IdId, M_bdbd, M_Ed, M_Id, M_bd, M_IEd, M_IId, M_Ibd, M_EEd, ...
%    M_EId, M_Ebd, M_EdId, M_Edbd, M_Idbd]} = ...
%    lin_evol_opt_rb_operators_separate_bases(model, detailed_data)
%
%
% function computing the time-dependent reduced basis operators and
% vectors (Galerkin Projection) including derivative operators. 
% Splittet into different decomp modes for returning operators, components
% or coefficients.
% We assume different reduced bases for the original solution and the
% derivative solutions.
% The derivative reduced bases are stores in detailed_data.RB_der{i}.
%
% In contrast to lin_evol_rb_operators the time step (I+dt*L) is not
% included in the operators.
%
% When choosing decomp mode = 1, it gives back projections on ALL available
% bases in RB_der. A selection of parameters to optimize is not done.
% When choosing decomp_mode = 2, it gives back the coefficients to ALL
% parameters in the model. A selection of parameters to optimize is not
% done.
%
% 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
% compute_derivative_info:  1 if derivative information is wanted
%                           0 if derivative information is not wanted
%
%
% 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  
%
% required fields of detailed_data:
%   RB: the reduced basis for the original solution
%   RB_der: cell array of derivative reduced bases
%
% generated fields of varargout:
%   L_E_dd:    galerkin projection on derivative space
%   L_I_dd:
%   dL_E_sd:    galerkin projection of operator applied on solution on
%               derivative space
%   dL_I_sd:
%   db:         b projected on derivative space
%   K...:       K error matrices (formulas in detail in Paper)
%
% In 'coefficients' mode the detailed data is empty.

% Markus Dihlmann 24.03.2011

% determine affine_decomposition_mode as integer  
decomp_mode = model.decomp_mode;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% decomp_mode 0                                       %%%
%%% complete calculation matrices                       %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if decomp_mode == 0 % complete: simple projection on RB

  %eval('u0 = ',params.init_value_algorithm,'(grid,params)']);
  [L_I, L_E, b] = model.operators_ptr(model,detailed_data);
    
  % 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;
  
  %reduced Operators
  LL_I = detailed_data.RB' * A * L_I * detailed_data.RB;
  
  %    LL_E    :  matrix 
  LL_E = detailed_data.RB' * A * L_E * detailed_data.RB;
  
  %    bb      :  vector 
  bb = detailed_data.RB' * A * b;
  
    
  %Error estimation operators 
  
  % often used /Phi^T L_E^T A
  RB_L_E = detailed_data.RB' * L_E' * A;
    
  M_E = RB_L_E * detailed_data.RB;
  M_b = b' * A * detailed_data.RB;
  M_EE = RB_L_E * L_E * detailed_data.RB;
  M_Eb = RB_L_E * b;
  M_bb = b' * A * b;
   

  % Implicit operators
  % often used /Phi^T L_I^T A
  RB_L_I = detailed_data.RB' * L_I' * A ;
  
  M_I = RB_L_I * detailed_data.RB;
  M_II = RB_L_I * L_I * detailed_data.RB;
  M_IE = RB_L_I * L_E * detailed_data.RB;
  M_Ib = RB_L_I * b;
   

  % creation of derivative matrices
  % get components and coefficients and do linear combination
%   if model.compute_derivative_info
%     if isfield(model.optimization,'params_to_optimize')
%         indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
%     else
%         indx_mu_i = 1:length(get_mu(model));
%     end
%     nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
% 
%     % create dummy files for derivative matrices --> memory
%     % preallocation
%     if ~isfield(model,'t')
%         model.t=1;
%     end
% 
%     model.decomp_mode = 1;
%     [dummy_LL_I, dummy_LL_E, dummy_bb, dummy_M_E, dummy_M_b, dummy_M_EE, ...
%     dummy_M_Eb, dummy_M_bb, dummy_M_I, dummy_M_II, dummy_M_IE, ...
%     dummy_M_Ib, dummy_M_EdEd, dummy_M_IdId, dummy_M_bdbd, ...
%     dummy_M_Ed, dummy_M_Id, dummy_M_bd, dummy_M_IEd, dummy_M_IId, ...
%     dummy_M_Ibd, dummy_M_EEd, dummy_M_EId, dummy_M_Ebd, dummy_M_EdId,...
%     ]=rb_operators(model,detailed_data); %dummy_M_Edbd, dummy_M_Idbd] = rb_operators(model,detailed_data); --> MEdbd, M_Idbd not needed
% 
%     if isempty(dummy_M_EdEd)
%         M_EdEd = zeros(1,1,nr_indx_mu_i);
%     else
%         M_EdEd = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_IdId)
%         M_IdId = zeros(1,1,nr_indx_mu_i);
%     else
%         M_IdId = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_Ed)
%         M_Ed = zeros(1,1,nr_indx_mu_i);
%     else
%         M_Ed = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_Id)
%         M_Id = zeros(1,1,nr_indx_mu_i);
%     else
%         M_Id = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_IEd)
%         M_IEd = zeros(1,1,nr_indx_mu_i);
%     else
%         M_IEd = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_IId)
%         M_IId = zeros(1,1,nr_indx_mu_i);
%     else
%         M_Id = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_EEd)
%         M_EEd = zeros(1,1,nr_indx_mu_i);
%     else
%         M_EEd = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_EId)
%         M_EId = zeros(1,1,nr_indx_mu_i);
%     else
%         M_EId = zeros(nRB,nRB,nr_indx_mu_i);
%     end
% 
%     if isempty(dummy_M_EdId)
%         M_EdId = zeros(1,1,nr_indx_mu_i);
%     else
%         M_EdId = zeros(nRB,nRB,nr_indx_mu_i);
%     end
%     
%     M_bdbd = zeros(1  ,1  ,nr_indx_mu_i);
%     M_bd   = zeros(1  ,nRB,nr_indx_mu_i);
%     M_Ibd  = zeros(nRB,1  ,nr_indx_mu_i);
%     M_Ebd  = zeros(nRB,1  ,nr_indx_mu_i);
%     M_Edbd = zeros(nRB,1  ,nr_indx_mu_i);
%     M_Idbd = zeros(nRB,1  ,nr_indx_mu_i);
%     % end of: create dummy files for derivative matrices --> memory
%     % preallocation
% 
% 
%     %get components
%     model.decomp_mode = 1;
%     [compLL_I, compLL_E, compbb, compM_E, compM_b, compM_EE, compM_Eb,...
%       compM_bb, compM_I, compM_II, compM_IE, compM_Ib, compM_EdEd,...
%       compM_IdId, compM_bdbd, compM_Ed, compM_Id, compM_bd, compM_IEd,...
%       compM_IId, compM_Ibd, compM_EEd, compM_EId, compM_Ebd,...
%       compM_EdId, compM_Edbd, compM_Idbd] = rb_operators(model,detailed_data);
% 
%   
%     %get coefficients
%     model.decomp_mode = 2;
%     [sLL_I, sLL_E, sbb, sM_E, sM_b, sM_EE, sM_Eb, sM_bb, sM_I, sM_II,...
%       sM_IE, sM_Ib, sM_EdEd, sM_IdId, sM_bdbd, sM_Ed, sM_Id, sM_bd,...
%       sM_IEd, sM_IId, sM_Ibd, sM_EEd, sM_EId, sM_Ebd, sM_EdId,...
%       sM_Edbd, sM_Idbd] = rb_operators(model,detailed_data);
% 
%   
%     %do linear combination of components and coefficients
%     for i=1:nr_indx_mu_i
%         M_EdEd(:,:,i) = lincomb_sequence(compM_EdEd,sM_EdEd(:,i));
%         M_IdId(:,:,i) = lincomb_sequence(compM_IdId,sM_IdId(:,i));
%         M_bdbd(:,:,i) = lincomb_sequence(compM_bdbd,sM_bdbd(:,i));
%         M_Ed(:,:,i)   = lincomb_sequence(compM_Ed  ,sM_Ed(:,i));
%         M_Id(:,:,i)   = lincomb_sequence(compM_Id  ,sM_Id(:,i));
%         M_bd(:,:,i)   = lincomb_sequence(compM_bd  ,sM_bd(:,i));
%         M_IEd(:,:,i)  = lincomb_sequence(compM_IEd ,sM_IEd(:,i));
%         M_IId(:,:,i)  = lincomb_sequence(compM_IId ,sM_IId(:,i));
%         M_Ibd(:,:,i)  = lincomb_sequence(compM_Ibd ,sM_Ibd(:,i));
%         M_EEd(:,:,i)  = lincomb_sequence(compM_EEd ,sM_EEd(:,i));
%         M_EId(:,:,i)  = lincomb_sequence(compM_EId ,sM_EId(:,i));
%         M_Ebd(:,:,i)  = lincomb_sequence(compM_Ebd ,sM_Ebd(:,i));
%         M_EdId(:,:,i) = lincomb_sequence(compM_EdId,sM_EdId(:,i));
%         M_Edbd(:,:,i) = lincomb_sequence(compM_Edbd,sM_Edbd(:,i));
%         M_Idbd(:,:,i) = lincomb_sequence(compM_Idbd,sM_Idbd(:,i));
%     end
%     
%       varargout( 1) = {M_EdEd};
%       varargout( 2) = {M_IdId};
%       varargout( 3) = {M_bdbd};
%       varargout( 4) = {M_Ed};
%       varargout( 5) = {M_Id};
%       varargout( 6) = {M_bd};
%       varargout( 7) = {M_IEd};
%       varargout( 8) = {M_IId};
%       varargout( 9) = {M_Ibd};
%       varargout(10) = {M_EEd};
%       varargout(11) = {M_EId};
%       varargout(12) = {M_Ebd};
%       varargout(13) = {M_EdId};
%       varargout(14) = {M_Edbd};
%       varargout(15) = {M_Idbd};
%   
%   end %if model.compute_derivative_info  


%set decomp_mode back to 0
model.decomp_mode = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% decomp_mode 1                                      %%%
%%% calculation of components                          %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif decomp_mode == 1
  
  %eval('u0 = ',params.init_value_algorithm,'(grid,params)']);
  [L_I, L_E, b] = model.operators_ptr(model, detailed_data);

  Q_L_I = length(L_I);
  Q_L_E = length(L_E);
  Q_b = length(b);
  
  
  Nmax = size(detailed_data.RB,2);
  H = size(detailed_data.RB,1);
  A = detailed_data.W;
  
  
  %reduced operator components:
  %    LL_I    :  matrices 
  LL_I = cell(Q_L_I,1);
  if ~isempty(L_I)
    LL_I(:) = {zeros(Nmax,Nmax)};
  for q = 1:Q_L_I;
    LL_I{q} = detailed_data.RB' * A * L_I{q} * detailed_data.RB;
  end;
  end
  %    LL_E    :  matrices
  LL_E = cell(Q_L_E,1);
  LL_E(:) = {zeros(Nmax,Nmax)};
  for q = 1:Q_L_E;
    LL_E{q} = detailed_data.RB' * A * L_E{q} * detailed_data.RB;
  end;
  
  %    bb      :  vectors 
  bb = cell(Q_b,1);
  bb(:) = {zeros(Nmax,1)};
  for q = 1:Q_b;
    bb{q} = detailed_data.RB' * A * b{q};
  end; 
  
   
  % init output quantities  
  M_E = cell(Q_L_E,1);          M_E(:) = {zeros(Nmax,Nmax)};
  M_b = cell(Q_b,1);            M_b(:) = {zeros(Nmax,1)};
  M_EE = cell(Q_L_E*Q_L_E,1);   M_EE(:) = {zeros(Nmax,Nmax)};
  M_Eb = cell(Q_L_E*Q_b,1);     M_Eb(:) = {zeros(Nmax,1)};
  M_bb = cell(Q_b*Q_b,1);       M_bb(:) = {zeros(1,1)};
  M_I = cell(Q_L_I,1);          M_I(:) = {zeros(Nmax,Nmax)};
  M_II = cell(Q_L_I*Q_L_I,1);   M_II(:) = {zeros(Nmax,Nmax)};  
  M_IE = cell(Q_L_I*Q_L_E,1);   M_IE(:) = {zeros(Nmax,Nmax)};
  M_Ib = cell(Q_L_I*Q_b,1);     M_Ib(:) = {zeros(Nmax,1)}; 

  
  % often used matrices:
  RB_L_E = cell(Q_L_E,1);
  RB_L_E(:) = {zeros(Nmax,H)};
  for q=1:Q_L_E
    RB_L_E{q}(:,:)= detailed_data.RB' * L_E{q}' * A;
  end;
  
  RB_L_I = cell(Q_L_I,1);
  RB_L_I(:) = {zeros(Nmax,H)};
  for q=1:Q_L_I
    RB_L_I{q}(:,:)= detailed_data.RB' * L_I{q}' * A;
  end;
  
  L_E_RB = cell(Q_L_E,1);
  L_E_RB(:) = {zeros(H,Nmax)};
  for q=1:Q_L_E
      L_E_RB{q} = L_E{q}*detailed_data.RB;
  end
  
  L_I_RB = cell(Q_L_I,1);
  L_I_RB(:) = {zeros(H,Nmax)};
  for q=1:Q_L_I
      L_I_RB{q} = L_I{q}*detailed_data.RB;
  end
  
  
  % fill output quantities
  for q=1:Q_L_E
      M_E{q} = RB_L_E{q} * detailed_data.RB;
  end
    
  for q=1:Q_b
      M_b{q} = b{q}' * A * detailed_data.RB;
  end
  
  for q1=1:Q_L_E
      for q2= 1:Q_L_E
          M_EE{q2+(q1-1)*Q_L_E} = RB_L_E{q1} * L_E_RB{q2};
      end
  end
  
  for q1=1:Q_L_E
      for q2=1:Q_b
          M_Eb{q2+(q1-1)*Q_b} = RB_L_E{q1} * b{q2};
      end;
  end
  
  for q1=1:Q_b
      for q2=1:Q_b
          M_bb{q2+(q1-1)*Q_b} = b{q1}' * A * b{q2};
      end
  end
  
  for q=1:Q_L_I
      M_I{q} = RB_L_I{q} * detailed_data.RB;
  end

  for q1=1:Q_L_I
      for q2= 1:Q_L_I
          M_II{q2+(q1-1)*Q_L_I} = RB_L_I{q1} * L_I_RB{q2};
      end
  end  
  
  for q1=1:Q_L_I
      for q2= 1:Q_L_E
          M_IE{q2+(q1-1)*Q_L_E} = RB_L_I{q1} * L_E_RB{q2}; 
      end
  end  
  
  for q1=1:Q_L_I
      for q2=1:Q_b
          M_Ib{q2+(q1-1)*Q_b} = RB_L_I{q1} * b{q2};
      end;
  end  
  
  
  if model.compute_derivative_info    
      
     if ~isfield(detailed_data,'RB_der')
         error('No derivative bases in detailed_data')
     end
     
%      %Projection on which bases is actually needed?
%     if isfield(model,'compute_derivative_indices')
%         indx_mu_i = find(model.compute_derivative_indices); %indexes i of \mu_i for which derivative information is needed
%     else
%         indx_mu_i = 1:length(get_mu(model));
%     end
%     nr_indx_mu_i= length(indx_mu_i); 
%     RB_der=detailed_data.RB_der(indx_mu_i); 

     RB_der=detailed_data.RB_der;
     nr_indx_mu_i = length(RB_der);
     
     %initialize cells
     L_E_dd = cell(nr_indx_mu_i,1);
     L_I_dd = cell(nr_indx_mu_i,1);
     dL_E_sd = cell(nr_indx_mu_i,1);
     dL_I_sd = cell(nr_indx_mu_i,1);
     db = cell(nr_indx_mu_i,1);
     K_E = cell(nr_indx_mu_i,1);
     K_I = cell(nr_indx_mu_i,1);
     K_EE = cell(nr_indx_mu_i,1);
     %K_Eb = cell(nr_indx_mu,1);
     %K_bb = cell(nr_indx_mu,1);
     K_II =  cell(nr_indx_mu_i,1);
     %K_Ib = cell(nr_indx_mu,1);
     K_IE = cell(nr_indx_mu_i,1);
     K_EdEd = cell(nr_indx_mu_i,1);
     K_IdId = cell(nr_indx_mu_i,1);
     K_bdbd = cell(nr_indx_mu_i,1);
     K_Ed = cell(nr_indx_mu_i,1);
     K_Id = cell(nr_indx_mu_i,1);
     K_bd = cell(nr_indx_mu_i,1);
     K_IEd = cell(nr_indx_mu_i,1);
     K_IId = cell(nr_indx_mu_i,1);
     K_Ibd = cell(nr_indx_mu_i,1);
     K_EEd = cell(nr_indx_mu_i,1);
     K_EId = cell(nr_indx_mu_i,1);
     K_Ebd = cell(nr_indx_mu_i,1);
     K_EdId = cell(nr_indx_mu_i,1);
     K_Edbd = cell(nr_indx_mu_i,1);
     K_Idbd =  cell(nr_indx_mu_i,1);
     
     
     
     
     for i=1:length(RB_der)
       
      NRB = size(RB_der{i},2);
      %make Galerkin projection for every reduced basis
      
      %    dL_E_dd    :  components
      L_E_dd{i} = cell(Q_L_E,1);
      L_E_dd{i}(:) = {zeros(NRB,NRB)};
      for q = 1:Q_L_E;
        L_E_dd{i}{q} = RB_der{i}' * A * L_E{q} * RB_der{i};
      end;
      
      %    dL_I_dd    :  components
      L_I_dd{i} = cell(Q_L_I,1);
      if ~isempty(L_I)
      L_I_dd{i}(:) = {zeros(NRB,NRB)};
      for q = 1:Q_L_I;
        L_I_dd{i}{q} = RB_der{i}' * A * L_I{q} * RB_der{i};
      end;
      end
      
      %    dL_E_sd    :  components
      dL_E_sd{i} = cell(Q_L_E,1);
      dL_E_sd{i}(:) = {zeros(NRB,Nmax)};
      for q = 1:Q_L_E;
        dL_E_sd{i}{q} = RB_der{i}' * A * L_E{q} * detailed_data.RB;
      end;
      
      %    dL_I_sd    :  components
      dL_I_sd{i} = cell(Q_L_I,1);
      if ~isempty(L_I)
      dL_I_sd{i}(:) = {zeros(NRB,Nmax)};
      for q = 1:Q_L_I;
        dL_I_sd{i}{q} = RB_der{i}' * A * L_I{q} * detailed_data.RB;
      end;
      end
      
      %    db      :  vectors 
      db{i} = cell(Q_b,1);
      db{i}(:) = {zeros(NRB,1)};
      for q = 1:Q_b;
        db{i}{q} = RB_der{i}' * A * b{q};
      end; 
      
     
      
     %      K_E
     K_E{i} = cell(Q_L_E,1);          
     K_E{i}(:) = {zeros(NRB,NRB)};
     for q=1:Q_L_E
         K_E{i}{q} = (RB_der{i}'*L_E{q}'*A)*RB_der{i};
     end
     
     %      K_b
     %K_b{i} = cell(Q_b,1);          
     %K_b(:) = {zeros(1,NRB)};
     %for q=1:Q_b
     %    K_b{i}{q} = b{q}'*A*RB_der{i}
     %end
     
     %      K_I
     K_I{i} = cell(Q_L_I,1);          
     K_I{i}(:) = {zeros(NRB,NRB)};
     for q=1:Q_L_I
         K_I{i}{q} = (RB_der{i}'*L_I{q}'*A)*RB_der{i};
     end
     
     % K_EE
     K_EE{i} = cell(Q_L_E*Q_L_E,1);
     K_EE{i}(:) = {zeros(NRB,NRB)};
      for q1=1:Q_L_E
          for q2= 1:Q_L_E
              K_EE{i}{q2+(q1-1)*Q_L_E} = (RB_der{i}'*L_E{q1}'*A)*(L_E{q2}*RB_der{i});
          end
      end
      
      %K_II
      K_II{i} = cell(Q_L_I*Q_L_I,1);
      K_II{i}(:) = {zeros(NRB,NRB)};
      for q1=1:Q_L_I
          for q2= 1:Q_L_I
              K_II{i}{q2+(q1-1)*Q_L_I} = (RB_der{i}' *L_I{q1}'*A)*(L_I{q2}*RB_der{i});
          end
      end  

      % K_IE
      K_IE{i} = cell(Q_L_I*Q_L_E,1);
      K_IE{i}(:) = {zeros(NRB,NRB)};
      for q1=1:Q_L_I
          for q2= 1:Q_L_E
              K_IE{i}{q2+(q1-1)*Q_L_E} = (RB_der{i}'*L_I{q1}'* A) * (L_E{q2} * RB_der{i}); 
          end
      end  

      %erster buchstabe---aeussere schleife
     %K_EdEd
     K_EdEd{i} = cell(Q_L_E*Q_L_E,1);
     K_EdEd{i}(:) = {zeros(Nmax,Nmax)};
     for q1=1:Q_L_E
         for q2=1:Q_L_E
             K_EdEd{i}{q2+(q1-1)*Q_L_E} = (detailed_data.RB'*L_E{q1}'*A)*(L_E{q2}*detailed_data.RB);
         end
     end
     
     % K_IdId
     K_IdId{i} = cell(Q_L_I*Q_L_I,1);
     K_IdId{i}(:) = {zeros(Nmax,Nmax)};
     for q1=1:Q_L_I
         for q2=1:Q_L_I
             K_IdId{i}{q2+(q1-1)*Q_L_I} = (detailed_data.RB'*L_I{q1}'*A)*(L_I{q2}*detailed_data.RB);
         end
     end
     
     %K_bdbd
     K_bdbd{i} = cell(Q_b*Q_b,1);
     K_bdbd{i} = {0};
     for q1 = 1:Q_b
         for q2 = 1:Q_b
             K_bdbd{i}{q2+(q1-1)*Q_b} = b{q1}'*A*b{q2};
         end
     end
     
     % K_Ed
     K_Ed{i} = cell(Q_L_E,1);
     K_Ed{i}(:) = {zeros(NRB,Nmax)};
     for q = 1:Q_L_E
         K_Ed{i}{q} = (detailed_data.RB'*L_E{q}'*A)*(RB_der{i});
     end
     
     %K_Id
     K_Id{i} = cell(Q_L_I,1);
     K_Id{i}(:) = {zeros(NRB,Nmax)};
     for q = 1:Q_L_I
         K_Id{i}{q} = RB_der{i}'*A*L_I{q}*detailed_data.RB;
     end
     
     %K_bd
     K_bd{i} = cell(Q_b,1);
     K_bd{i}(:) = {zeros(1,NRB)};
     for q= 1:Q_b
         K_bd{i}{q} = b{q}'*A*RB_der{i};
     end
     
     %K_IEd
     K_IEd{i} = cell(Q_L_I * Q_L_E,1);
     K_IEd{i}(:) = {zeros(NRB,Nmax)};
     for q1=1:Q_L_I
         for q2=1:Q_L_E
             K_IEd{i}{q2+(q1-1)*Q_L_E} = (RB_der{i}' * L_I{q1}'*A)*(L_E{q2}*detailed_data.RB);
         end
     end
     
     %K_IId
     K_IId{i} = cell(Q_L_I*Q_L_I,1);
     K_IId{i}(:) = {zeros(NRB,Nmax)};
     for q1=1:Q_L_I
         for q2=1:Q_L_I
             K_IId{i}{q2+(q1-1)*Q_L_I} = (RB_der{i}' * L_I{q1}'*A)*(L_I{q2}*detailed_data.RB);
         end
     end
     
     %K_Ibd
     K_Ibd{i} = cell(Q_L_I*Q_b,1);
     K_Ibd{i}(:) = {zeros(NRB,1)};
     for q1 = 1:Q_L_I
         for q2 = 1: Q_b
             K_Ibd{i}{q2+(q1-1)*Q_b} = (RB_der{i}'*L_I{q1}') * A * b{q2};
         end
     end
     
     %K_EEd
     K_EEd{i} = cell(Q_L_E*Q_L_E,1);
     K_EEd{i}(:) = {zeros(NRB,Nmax)};
     for q1=1:Q_L_E
         for q2=1:Q_L_E
             K_EEd{i}{q2+(q1-1)*Q_L_E} = (RB_der{i}' * L_E{q1}'*A)*(L_E{q2}*detailed_data.RB);
         end
     end
     
     %K_EId
     K_EId{i} = cell(Q_L_E*Q_L_I,1);
     K_EId{i}(:) = {zeros(NRB,Nmax)};
     for q1=1:Q_L_E
         for q2=1:Q_L_I
             K_EId{i}{q2+(q1-1)*Q_L_I} = (RB_der{i}' * L_E{q1}'*A)*(L_I{q2}*detailed_data.RB);
         end
     end
     
     %K_Ebd
     K_Ebd{i} = cell(Q_L_E*Q_b,1);
     K_Ebd{i}(:) = {zeros(NRB,1)};
     for q1 = 1:Q_L_E
         for q2 = 1:Q_b
             K_Ebd{i}{q2+(q1-1)*Q_b} = RB_der{i}' * L_E{q1}' * A * b{q2};
         end
     end
     
     %K_EdId
     K_EdId{i} = cell(Q_L_E*Q_L_I,1);
     K_EdId{i}(:) = {zeros(Nmax,Nmax)};
     for q1 = 1:Q_L_E
         for q2 = 1:Q_L_I
             K_EdId{i}{q2+(q1-1)*Q_L_I} = (detailed_data.RB'*L_E{q1}' * A) * (L_I{q2}'* detailed_data.RB);
         end
     end
     
     %K_Edbd
     K_Edbd{i} = cell(Q_L_E*Q_b,1);
     K_Edbd{i}(:) = {zeros(Nmax,1)};
     for q1 = 1:Q_L_E
         for q2 = 1:Q_b
             K_Edbd{i}{q2+(q1-1)*Q_b} = detailed_data.RB'*L_E{q1}'*A*b{q2};
         end
     end
    
     %K_Idbd
     K_Idbd{i} = cell(Q_L_I*Q_b,1);
     K_Idbd{i}(:) = {zeros(Nmax,1)};
     for q1 = 1:Q_L_I
         for q2 = 1:Q_b
             K_Idbd{i}{q2+(q1-1)*Q_b} = detailed_data.RB'*L_I{q1}'*A*b{q2};
         end
     end
     
     end %for over all RB_der     
     
     
     varargout(1) = {L_E_dd};
     varargout(2) = {L_I_dd};
     varargout(3) = {dL_E_sd};
     varargout(4) = {dL_I_sd};
     varargout(5) = {db};
     varargout(6) = {K_E}; 
     varargout(7) = {K_I}; 
     varargout(8) = {K_EE};
     varargout(9) = {K_II};
     varargout(10) = {K_IE};
     varargout(11) = {K_EdEd}; 
     varargout(12) = {K_IdId};
     varargout(13) = {K_bdbd};
     varargout(14) = {K_Ed};
     varargout(15) = {K_Id};
     varargout(16) = {K_bd};
     varargout(17) = {K_IEd};
     varargout(18) = {K_IId};
     varargout(19) = {K_Ibd};
     varargout(20) = {K_EEd};
     varargout(21) = {K_EId};
     varargout(22) = {K_Ebd};
     varargout(23) = {K_EdId};
     varargout(24) = {K_Edbd};
     varargout(25) = {K_Idbd};
     
     
%       varargout( 1) = {M_EE};  % M_EdEd
%       varargout( 2) = {M_II};  % M_IdId
%       varargout( 3) = {M_bb};  % M_bdbd
%       varargout( 4) = {M_E} ;  % M_Ed
%       varargout( 5) = {M_I} ;  % M_Id
%       varargout( 6) = {M_b} ;  % M_bd
%       varargout( 7) = {M_IE};  % M_IEd 
%       varargout( 8) = {M_II};  % M_IId
%       varargout( 9) = {M_Ib};  % M_Ibd
%       varargout(10) = {M_EE};  % M_EEd
%       varargout(11) = {M_IE'}; % M_EId
%       varargout(12) = {M_Eb};  % M_Ebd
%       varargout(13) = {M_IE'}; % M_EdId
%       varargout(14) = {M_Eb};  % M_Edbd
%       varargout(15) = {M_Ib};  % M_Idbd
  end
  
  
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% decomp_mode 2               
%%% calculation of coefficients 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else % decomp_mode== 2 -> coefficients 
  [L_I, L_E, b] = model.operators_ptr(model,[]);
  [sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model); %gives coefficients to all the parameters in the model!
    
  Q_L_I = length(L_I(:));
  Q_L_E = length(L_E(:));
  Q_b = length(b(:));
  
  % rb-operator coefficients can simply be forwarded
  LL_E = L_E;
  LL_I = L_I;
  bb = b;
  
  
  % rb error operators -coefficients
  M_E = L_E;
  M_b = b;
  M_I = L_I;
  M_EE = repmatrows(L_E(:),Q_L_E) .* repmat(L_E(:),Q_L_E,1);
  M_Eb = repmatrows(L_E, Q_b) .* repmat(b,Q_L_E,1);
  M_bb = repmatrows(b, Q_b) .* repmat(b,Q_b,1);
  M_II = repmatrows(L_I(:),Q_L_I) .* repmat(L_I(:),Q_L_I,1);
  M_Ib = repmatrows(L_I, Q_b) .* repmat(b,Q_L_I,1);
  M_IE = repmatrows(L_I(:),Q_L_E) .* repmat(L_E(:),Q_L_I,1); 
  
  % computation of derivative coefficients
  if model.compute_derivative_info
%       
%     if isfield(model.optimization,'params_to_optimize')
%           indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
%     else
%           indx_mu_i = length(get_mu(model));
%     end
%     nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
      
%      if isfield(model,'compute_derivative_indices')
%         nr_indx_mu_i = length(model.compute_derivative_indices);
%      else
%         nr_indx_mu_i = length(get_mu(model));
%      end
        
      nr_indx_mu_i=length(get_mu(model));
      %obsolete: use only those coefficients of parameters which are to be optimimzed
      %sLe = sLe';%sLe(indx_mu_i,:)';
      %sb = sb';%(indx_mu_i,:)';
      % in case of no implicit operator (e.g. sLi=[]) repmat does not work!
      if ~isempty(sLi)
          %sLi = sLi(indx_mu_i,:)';
      else
          sLi=repmat(repmatrows(L_I(:),Q_L_I),1,nr_indx_mu_i); %sLi ist a 0xn empty matrix, not only a 0x1 empty matrix
      end
      
      varargout(1) = {L_E};
      varargout(2) = {L_I};
      varargout(3) = {sLe'};
      varargout(4) = {sLi'};
      varargout(5) = {sb'};
      
      %error matrices
      varargout(6)= {L_E}; %K_E
      %varargout(7) = b;%K_b
      varargout(7) = {L_I}; %K_I
      varargout(8) = {repmatrows(L_E(:),Q_L_E) .* repmat(L_E(:),Q_L_E,1)}; %K_EE
      %varargout(9) = repmatrows(L_E, Q_b) .* repmat(b,Q_L_E,1); %K_Eb
      %varargout(10) = repmatrows(b, Q_b) .* repmat(b,Q_b,1); %K_bb
      varargout(9) = {repmatrows(L_I(:),Q_L_I) .* repmat(L_I(:),Q_L_I,1)}; %K_II
      %varargout(12) = repmatrows(L_I, Q_b) .* repmat(b,Q_L_I,1); %K_Ib
      varargout(10) = {repmatrows(L_I(:),Q_L_E) .* repmat(L_E(:),Q_L_I,1)};  %K_IE
      varargout(11) = {repmatrows(sLe',Q_L_E) .* repmat(sLe',Q_L_E,1)};  % K_EdEd 
      varargout(12) = {repmatrows(sLi',Q_L_I) .* repmat(sLi',Q_L_I,1)};  % K_IdId
      varargout(13) = {repmatrows(sb',Q_b) .* repmat(sb',Q_b,1)};        % K_bdbd
      varargout(14) = {sLe'};                                           % K_Ed
      varargout(15) = {sLi'};                                           % K_Id
      varargout(16) = {sb'};                                            % K_bd
      varargout(17) = {repmat(repmatrows(L_I(:),Q_L_E),1,nr_indx_mu_i) .* repmat(sLe',Q_L_I,1)}; % K_IEd
      varargout(18) = {repmat(repmatrows(L_I(:),Q_L_I),1,nr_indx_mu_i) .* repmat(sLi',Q_L_I,1)}; % K_IId
      varargout(19) = {repmat(repmatrows(L_I(:),Q_b  ),1,nr_indx_mu_i) .* repmat(sb',Q_L_I,1)};  % K_Ibd
      varargout(20) = {repmat(repmatrows(L_E(:),Q_L_E),1,nr_indx_mu_i) .* repmat(sLe',Q_L_E,1)}; % K_EEd 
      varargout(21) = {repmat(repmatrows(L_E(:),Q_L_I),1,nr_indx_mu_i) .* repmat(sLi',Q_L_E,1)}; % K_EId
      varargout(22) = {repmat(repmatrows(L_E(:),Q_b  ),1,nr_indx_mu_i) .* repmat(sb',Q_L_E,1)};  % K_Ebd
      varargout(23) = {repmatrows(sLe',Q_L_I) .* repmat(sLi',Q_L_E,1)};                           % K_EdId
      varargout(24) = {repmatrows(sLe',Q_b) .* repmat(sb',Q_L_E,1)};                              % K_Edbd
      varargout(25) = {repmatrows(sLi',Q_b) .* repmat(sb',Q_L_I,1)};                              % K_Idbd
      
  end %if model.compute_derivative_info

end