function [t_m,u_m,stats,varargout] = rowmap(rhs,tspan,u0,opts,varargin)
% ROWMAP ... A ROW-code with Krylov techniques for large stiff ODEs.
%
%    ROWMAP solves the inital value problem for stiff systems of
% first order ODEs
% 
%            du                                 u = (  u(1),...,  u(n))
%           ---- = rhs(t,u) , u(t0) = u0      rhs = (rhs(1),...,rhs(n)) 
%            dt
%
%    Here n is the number of equations of the ODE system.
%           
%    ROWMAP is based on the ROW-methods of order 4 of the codes ROS4 and
% RODAS of Hairer and Wanner [4] and RODASP of Steinebach [5].
% The appearing large linear equation systems are solved with Krylov
% techniques. By a special multiple Arnoldi process the order of the basic
% method is preserved with small Krylov dimensions. Step size control is 
% done by embedding with a method of order 3. To reduce rounding errors 
% in the Arnoldi process, the selective reorthogonalization of Giraud and 
% Langou is used [6].
%
%
% Version of May 25, 2009.
%
% Authors:
%   St. Beck, H. Podhaisky, R. Weiner,
%   Institut für Mathematik, Universitaet Halle
%   06099 Halle, Germany
%   email: steffen.beck@mathematik.uni-halle.de
%          helmut.?podhaisky@mathematik.uni-halle.de
%          ruediger.?weiner@mail.mathematik.uni-halle.de,
% 
%   B.A. Schmitt,
%   Fachbereich Mathematik und Informatik, Universitaet Marburg
%   35032 Marburg, Germany
%   email: schmitt@mathematik.uni-marburg.de
% 
%--------------------------------------------------------------------------
% REFERENCES
%--------------------------------------------------------------------------
%
% [1]  R. Weiner, B.A. Schmitt an H. Podhaisky: 
%         ROWMAP - a ROW-code with Krylov techniques for large stiff ODEs. 
%         Applied Numerical Mathematics 25 (1997)
% [2]  B.A. Schmitt and R. Weiner:
%         Matrix-free W-methods using a multiple Arnoldi Iteration,
%         APNUM 18(1995), 307-320
% [3]  R.Weiner and B.A.Schmitt:
%         Order results for Krylov-W-methods, Computing 61, 69 -- 89 (1998).
% [4]  E. Hairer and G. Wanner:
%         Solving Ordinary Differential Equations II, Springer-Verlag,
%         Second Edition, 1996
% [5]  Steinebach, G.: 
%         Order-reduction of ROW-methods for DAEs and Method of Lines 
%         Applications. Preprint-Nr. 1741, 1995.
% [6]  L. Giraud and J. Langou
%         A robust criterion for the modified Gram-Schmidt algorithm with 
%         selective reorthogonalization.
%         SIAM J. Sci. Comput. 25 (2003), no. 2, 417--441 (electronic).
%
%
%--------------------------------------------------------------------------
% SYNOPSIS
%--------------------------------------------------------------------------
%
% [tout,yout] = rowmap(rhs,tspan,u0)
% [tout,yout] = rowmap(rhs,tspan,u0,opts)
% [tout,yout] = rowmap(rhs,tspan,u0,[],P1,P2,...)
% [tout,yout] = rowmap(rhs,tspan,u0,opts,P1,P2,...)
% [tout,yout,stats] = rowmap(rhs,tspan,u0)
% ...
%
% Details for the input and output variables you can find below.
%
% NOTE: It makes a difference whether the differential equation is 
%       autonomous or not, use the option field opts.Autonomous
%       to specify this. This influences the whole process of integration
%       with rowmap. See the Options section below.
%
%--------------------------------------------------------------------------
% INPUT PARAMETERS
%--------------------------------------------------------------------------
% 
% Mandatory arguments
%=====================
%
% rhs       rhs = @name_of_odefun
%           Is a function_handle(@) to the function that evaluates the 
%           right-hand side of the differential equations. 
%
%           possibilities for creating your ODE function:
%
%              1. If you only need the values of t and u use: 
%                   function du = name_of_odefun(t,u)
%                   ...
%              2. If you want to use some parameters P1,P2,... which have 
%                 the same values over the whole integration use:
%                   function du = name_of_odefun(t,u,P1,P2,...)
%                   ...
%                 or in combination with function_handle(@) you can
%                 use 
%                           rhs = @(t,u)name_of_odefun(t,u,P1,P2,...)
%                 
%                 this has a better performance.
%
%              3. If you want to use some parameters P1,P2,... which have
%                 changing values during the integration use:
%                   function [du,P1,P2,...] = name_of_odefun(t,u,P1,P2,...)
%                   ...
%                 or you use nested functions in combination with
%                 function_handle(@)
%           
%              In all cases u and du are column vectors.
%
%           NOTE: 
%            If you want to use parameters in your right side function you
%            must use the following call to ROWMAP:
%               
%                  [t,u,stats] = rowmap(rhs,tspan,u0,opts,P1,P2...) 
%            
%            This call passes the additional parameters P1,P2,... to the 
%            ODE function as rhs(t,u,P1,P2...), and to all functions 
%            specified in OPTIONS. Use opts = [] as a place holder if 
%            no options are set.
% 
% tspan     tspan = [t0 tfinal] 
%           A vector specifying the interval of integration, [t0 tfinal]. 
%           The solver imposes the initial condition u0 at tspan(1), and 
%           integrates from t0 = tspan(1) to tfinal = tspan(end). 
%           The solver returns the numerical solution at the final time 
%           tfinal.
%           If you want to have a numerial solution for each integration 
%           step use the structure field opts.Solout. See the options
%           section below.
%           
% u0        u0 = u(t0) = [u1(t0),u2(t0),...,un(t0)]^T
%           The column vector of initial conditions.
% 
% Options
%=========
% 
% opts      when it is given to ROWMAP it must be a 1x1 structure array
%           it can contain the following fields, in brackets are the 
%           default values:
%
%       opts.Scheme         is a String of the underlying ROW-method
%        ['GRK4T']          one of:   'GRK4T'
%                                     'LSTAB'
%                                     'RODAS'
%                                     'RODASP'
% 
%       opts.AbsTol         absolute error tolerance
%        [1e-03]            scalar or vector of length n
%                           
%       opts.RelTol         relative error tolerance
%        [1e-03]            scalar or vector of length n
%
%       opts.Autonomous     is of class logical
%        [false]            describes if the given right-hand side function
%                           is independent of t (autonomous) or may depend 
%                           on t (nonautonomous)
%
%                               aut = true      rhs is autonomous
%                               aut = false     rhs is nonautonomous
%      
%       opts.KryTol         scalar
%        [0.1]              tolerance for the iterative solution of the
%                           linear equations
%
%                           refTol = h_m*sqrt(1/n*sum((res./...
%                                              (atol+rtol.*abs(u_m))).^2));
%
%                           if refTol < KryTol the iteration process is
%                                              stoped
% 
%       opts.ReorthTol      scalar
%        [1e+04]            tolerance for reorthogonalization by the so
%                           called L-criterion see [5]
%
% 
%       opts.MaxIter        scalar
%        [63]               the iterative solution of the linear equations
%                           in the first stage of integration will be 
%                           aborted if MaxIter is reached
%                           
% 
%       opts.InitialStep    scalar
%        [h_init]           first step size, default:
%
%                           h_init = 0.1*(sqrt(n)/norm(1./...
%                                             (atol+rtol.*abs(u_m))))^0.25;
% 
%       opts.StepControl    parameters for the step size control
%        [0.8,0.25,2]       [a,amin,amax] vector of the length 3
%         
%                           h_new = h_old * min( amax , ...
%                                           max( a * err^-(p+1) , amin ) );
% 
%       opts.MinStep        scalar
%        [1e-14]            if the step size is smaller then MinStep the
%                           integration is aborted
% 
%       opts.MaxStepNum     scalar
%        [1e+04]            if the number of steps(accepted + rejected)
%                           reaches MaxNumStep the integration is aborted
% 
%       opts.Jacv           function_handle(@) for the
%        [@AvDQ]            user-supplied function for calculating the
%                           jacobian-vektor product
%                           If this option is not provided Jv is computed  
%                           by finite differences in the internal function 
%                           AvDQ. 
%
%                           Function Jacv must have the following form:
%
%                           function [Jv,P1,P2,...,Pi] = ...
%                             Jacv(rhs,t_m,u_m,norm_u_m,f_m,v,P1,P2,...,Pj)
%
%                           rhs,t_m,u_m,norm_u_m,f_m,v 
%                              are input arguments which must exist                           
%                           P1,P2,...,Pj
%                              are optional input arguments, j<=the number
%                              of arguments, which are given to rowmap
%                           Jv 
%                              is the computed Jacobian vector product%
%                           P1,P2,...,Pi
%                              are optional return arguments, i<=j
% 
%       opts.Fdt            function_handle(@) for the
%        [@FdtDQ]           user-supplied function for calculating the
%                           time derivative ft of f(t,u)
%                           If this option is not provided ft is computed
%                           by finite differences internally.
%
%                           Function Fdt must have the following form:
%
%                           function [Ft,P1,P2,...,Pi] = ...
%                                         Fdt(rhs,t_m,u_m,f_m,P1,P2,...,Pj)
%
%                           rhs,t_m,u_m,f_m
%                              are input arguments which must exist
%                           Ft 
%                              is the computed time derivative of f(t,u) %                              
%
%                           P1,P2,...,Pj and P1,P2,...,Pi
%                              as mentioned in opts.Jacv
% 
%       opts.Solout         function_handle(@) for the
%        []                 user-supplied function exporting the numerical
%                           solution during integration
%
%                           Function solout must have the following form:
%                           
%                           function [P1,P2,...,Pi] = ...
%                            solout(told,tnew,uold,unew,fold,fnew,...
%                                                            ,P1,P2,...,Pj)
%                           told,tnew,uold,unew,fold,fnew
%                                  are input arguments which must be exist
%
%                           no output argumends are needed but you can pass
%                           some parameters P1,P2,... back if you want, 
%                           they will be stored on your declared variabels 
%                           P1,P2,...
%
%                           P1,P2,...,Pj and P1,P2,...,Pi
%                                  as mentioned in opts.Jacv
%
%--------------------------------------------------------------------------
% OUTPUT PARAMETERS
%--------------------------------------------------------------------------
%
% t_m       the last point in time for which a numerical solution has been 
%           computed, normaly tfinal
%
% u_m       the numerical solution at the point t_m
%
% stats     structure array with some statistics
%
%   stats.Scheme     underlying ROW-Method
%   stats.AbsTol     absolute error tolerance
%   stats.RelTol     relative error tolerance
%   stats.KryTol     tolerance for the Krylov method
%   stats.ReorthTol  tolerance for reorthogonalization
%   stats.nsteps     number of total(accepted + rejected)t ime steps
%   stats.nacc       number of accepted time steps
%   stats.nrej       number of rejected time steps
%   stats.Time       CPU time for time integration
%   stats.FCN        number of calls of the ode function
%   stats.nJacv      number of calls of the function Jacv (internally AvDQ)
%   stats.nFdt       number of calls of the function Fdt
%   stats.nOrth      number of orthogonalizations
%   stats.nReorth    number of reorthogonalization
%   stats.Info       Reports on success upon return
%   stats.Ave_iter   Average of Krylov dimension of the first stage over 
%                    all timesteps(stats.nsteps)
%


% -------------------------------------------------------------------------
% ---------------------------- Initialisation -----------------------------
% -------------------------------------------------------------------------

RM_vers = mfilename; % filename of the current ROWMAP version

% struct for the number of input and output arguments, the third entry of
% each field is the number of the necessary arguments of the corresponding 
% function
ninout = struct('funcName',{'RM'   ,'rhs','Jacv','Fdt','Solout'},...
                'nin'     ,{nargin , []  , []   , []  , []     },...
                'ninnec'  ,{ 3     , 2   , 6    , 4   , 6      },...
                'ninpar'  ,{ []    , []  , []   , []  , []     },...
                'nout'    ,{nargout, []  , []   , []  , []     },...
                'noutnec' ,{ 3     , 1   , 1    , 1   , 0      },...
                'noutpar' ,{ []    , []  , []   , []  , []     });
npar = struct('RM',{[],[]},'rhs',{[],[]},'Jacv',{[],[]},...
              'Fdt',{[],[]},'Solout',{[],[]});
fnames = {ninout.funcName}';

% test if the number of input arguments is correct and if options are
% chosen
if ninout(1).nin < ninout(1).ninnec + 1
    opts = [];   % all default options are used
    if ninout(1).nin < ninout(1).ninnec
        error(['Not enough input arguments. See ',RM_vers,'.'])        
    end
else
    ninout(1).ninnec = ninout(1).ninnec + 1;
end

if any(size(tspan) > 1); sztp = size(tspan);else sztp = [2 2];end
szu0 = size(u0);

% additional parameters P1,P2,... are saved in par and will be passed to
% the rhs file, and to all other user-supplied functions
par = varargin;
szpar = length(par);

% collecting the necessary variables in a cell array together with their
% properties they should have
varin = {rhs   ,'function_handle',[1 1]                ,'first'
         tspan ,'double'         ,[sztp(1),1;1,sztp(2)],'second'
         u0    ,'double'         ,[szu0(1),1;1,szu0(2)],'third'};

% checking the class an the size of the necessary variables
for i = 1:size(varin,1)
    if ~any(strcmp(class(varin{i,1}),varin{i,2})) % checking class
        error('ROWMAP_Initialisation:Varin:class',...
          ['The ',varin{i,4},' input argument is of the wrong class.\n',...
           'See ',RM_vers,'.'])        
    end
    if ~any(all(varin{i,3}' == ...     % checking size
        repmat(size(varin{i,1})',1,size(varin{i,3},1))))
        error('ROWMAP_Initialisation:Varin:class',...
            ['The ',varin{i,4},' input argument has the wrong size.\n',...
             'See ',RM_vers,'.'])
    end
end

u0    = u0(:);          % making column vectors
tspan = tspan(:);

n = length(u0);         % reading the number of equations

szm = [1 1];
if isa(opts,'struct')
    if any(size(opts)>1)
        error('ROWMAP_Initialisation:opts:size',...
              ['Your option structure has the wrong size. It must be',...
               'of the size 1x1. See ',RM_vers,'.'])
    end
    fnsopts = fieldnames(opts);   % fieldnames of the selected options            
    ifield = strcmpi(fnsopts,'Scheme');
    if any(ifield); szm = size(opts.(fnsopts{ifield}));
    end
else
    fnsopts = {};   % fieldnames of the selected options            
end

szv = length(RM_vers);

% default options and the properties they should have
% fieldname        fieldvalue    fieldclass           fieldsize
defopts = struct(...
  'Scheme'        ,{'GRK4T'     ,'char'        ,[1 szm(2);szm(1) 1]},...
  'AbsTol'        ,{1e-03       ,'double'            ,[1 1;n 1;1 n]},...
  'RelTol'        ,{1e-03       ,'double'            ,[1 1;n 1;1 n]},...
  'KryTol'        ,{0.1         ,'double'            ,[1 1]}        ,...
  'ReorthTol'     ,{1e+04       ,{'double','logical'},[1 1]}        ,...
  'MaxIter'       ,{63          ,'double'            ,[1 1]}        ,...
  'InitialStep'   ,{[]          ,'double'            ,[1 1]}        ,...
  'StepControl'   ,{[0.8,0.25,2],'double'            ,[3 1;1 3]}    ,...
  'MinStep'       ,{1e-14       ,'double'            ,[1 1]}        ,...
  'MaxStepNum'    ,{1e+04       ,'double'            ,[1 1]}        ,...
  'Jacv'          ,{@AvDQ       ,'function_handle'   ,[1 1]}        ,...
  'Fdt'           ,{@FdtDQ      ,'function_handle'   ,[1 1]}        ,...
  'Solout'        ,{[]          ,'function_handle'   ,[1 1]}        ,...
  'Autonomous'    ,{false       ,'logical'           ,[1 1]}        ,...
  'IterStat'      ,{false       ,'logical'           ,[1 1]}        ,...
  's'             ,{4           ,'double'            ,[1 1]}        ,...
  'Vers'          ,{RM_vers     ,'char'              ,[1 szv;szv 1]},...
  'LinSysErrPlot' ,{false       ,'logical'           ,[1 1]}             );

doprops = defopts(2:3); % defopts properties
defopts = defopts(1);   % defopts

% reading the fieldnames of the default options
fnsallopts = fieldnames(defopts);

% changing default options by user-supplied
if ~isempty(opts)                     % test if any option is changed
    if isa(opts,'struct')             % test if opts is a structure array        
        for i = 1:length(fnsopts)     % defining the chosen options
            ifield = strcmpi(fnsallopts,fnsopts{i});
            % test if the chosen option has the rigth fieldname
            if any(ifield)
                optval = opts.(fnsopts{i});
                fprop = {doprops(1:2).(fnsallopts{ifield})};
                if ~any(strcmp(class(optval),fprop{1})) % checking class
                    error('ROWMAP_Initialisation:opts:class',...
                        ['The option field "',(fnsallopts{ifield}),'" ',...
                         'is of the wrong class. See ',RM_vers,'.'])
                end
                if ~any(all(fprop{2}' == ...   % checking size
                    repmat(size(optval)',1,size(fprop{2},1))))
                    error('ROWMAP_Initialisation:opts:size',...
                        ['The option field "',(fnsallopts{ifield}),'" ',...
                         'has the wrong size. See ',RM_vers,'.'])
                end
                % change default option with the user defined option
                defopts.(fnsallopts{ifield}) = opts.(fnsopts{i});
            else
                error('ROWMAP_Initialisation:WrongOptionNames',['Wrong '...
                      'option field.\n The option "',fnsopts{i},'" isn`'...
                      't supported in ',RM_vers,'. See ',RM_vers,'.'])
            end
        end
    else
        error(['The argument opts musst be a structure array. See ',...
               RM_vers])
    end
end

opts = defopts; % replacing the argument opts by defopts, maybe some 
                % default options are changed if the user has defined some

% reading the number of in- and output arguments of the user-supplied
% functions
ninout(2).nin  =  nargin(rhs);
ninout(2).nout = nargout(rhs);
for i = 3:length(fnames)
    if ~isempty(opts.(fnames{i}))
        ninout(i).nin  =  nargin(opts.(fnames{i}));
        ninout(i).nout = nargout(opts.(fnames{i}));
    else
        ninout(i).nin  = ninout(i).ninnec;
        ninout(i).nout = ninout(i).noutnec;
    end
end

% checking the number of in- and output arguments of the user-supplied
% functions and compute the number of the passing variables of each
% function
for i = 1:length(fnames)
    % checking the number of input variables, compute number of variables 
    % which are forwarded to the corresponding function
    errninpar = ['The number of parameters, which are defined ',...
        'as input variables in your \n"',fnames{i},...
        '" function is greater as the number of parameters',...
        ' which are given to ',RM_vers,'.'];
    if ninout(i).nin >=0
        if ninout(i).ninnec > ninout(i).nin
            error('ROWMAP_Init:ninoutnec',...
                ['The number of your defined input variables in your \n"',...
                fnames{i},'" function is to small.']);
        else                
            ninout(i).ninpar  =  ninout(i).nin  - ninout(i).ninnec;
            if ninout(i).ninpar > ninout(1).ninpar
                error('ROWMAP_init:NinparToLarge1',errninpar);
            end
        end
    else
        if abs(ninout(i).nin) -1 - ninout(i).ninnec > ninout(1).ninpar
            error('ROWMAP_init:NinparToLarge2',errninpar);
        end
        ninout(i).ninpar = szpar;
    end
    % checking the number of output variables, compute number of variables 
    % which are sended back from the corresponding function
    errnoutpar = ['The number of parameters, which are defined ',...
        'as output variables in your \n"',fnames{i},...
        '" function is greater as the number of parameters',...
        ' which are given to ',RM_vers,'.'];
    if ninout(i).nout >=0
        if ninout(i).noutnec > ninout(i).nout && i>1
            error('ROWMAP_Init:ninoutnec',...
                ['The number of your defined output variables in your \n"',...
                fnames{i},'" function is to small.']);
        elseif i == 1 && ninout(i).noutnec > ninout(i).nout
            ninout(i).noutpar = 0;
        else          
            ninout(i).noutpar =  ninout(i).nout - ninout(i).noutnec;
            if ninout(i).noutpar > ninout(1).ninpar
                error('ROWMAP_init:NoutparToLarge1',errnoutpar);
            end
        end
    else
        if abs(ninout(i).nout) -1 - ninout(i).noutnec > ninout(1).ninpar
            error('ROWMAP_init:NoutparToLarge2',errnoutpar);
        end
        ninout(i).noutpar = szpar;
    end
    [npar.(fnames{i})] = deal(ninout(i).ninpar,ninout(i).noutpar);
end

if ~any(strcmpi(fnsopts,'jacv'))
    [npar.Jacv] = npar.rhs;
end

if any(strcmpi(fnsopts,'solout'))
    iout = true;
else
    iout = false;
end

stats = struct('Scheme',opts.Scheme,...
               'AbsTol',opts.AbsTol,...
               'RelTol',opts.RelTol,...
               'KryTol',opts.KryTol,...
               'ReorthTol',opts.ReorthTol,...
               'nsteps',[],...
               'nacc',[],...
               'nrej',[],...
               'Time',[],...
               'FCN',[],...
               'nJacv',[],...
               'nFdt',[],...
               'nOrth',[],...
               'nReorth',[],...
               'Info',[],...
               'Ave_iter',[],...
               'K',[]);

fprintf('--------------------------------------------------------------\n')
fprintf('Starting ROWMAP with the following options:\n')
disp(opts)
fprintf('--------------------------------------------------------------\n')
               
% Transfering the chosen options to the ROWMAP variables                
% -------------------------------------------------------------------------

method = opts.Scheme;
kd_max = opts.MaxIter;
aut    = opts.Autonomous;

% Abrufen der Verfahrensparameter
[a_ij,b,be,c_ij,c_i,d_i,g_ij,kd_inc,b_par] = rowcoe(method,aut,kd_max);
gamma   = g_ij(1,1);
d_i = d_i';

% Krylovparameter
kd      = cumsum(kd_inc); % Krylovdimensionen der einzelnen Stufen
kd_init = [0 kd];         % Ausgangskrylovdimensionen
s       = length(b);      % number of stages
s_K     = kd(s);          % maximale Krylovdimension
kd      = [0,kd];         % kd erweitern fï¿½r Stufenï¿½bergang
kd1_ave = 0;              % Durchschnitt der Krylovdim der 1. Stufe

% Matrizen initialisieren
sqrt_n = sqrt(n);          % Wurzel aus n zum skalieren
K      = zeros(n,s);       % Steigungwertmatrix
Q      = zeros(n,s_K + s); % Q auf maximale Groesse vordefinieren
H      = zeros(s_K);       % H auf maximale Groesse vordefinieren
Hs     = H;                % Matrix zur Lsg. des GS in jeder Stufe
Giv    = cell(s,1);        % Cellarry fuer GD der jew. Stufe
dimre  = {NaN(1,s_K),0};   % Vektor der Reotho. und num_b_in speichert


% Werte fuer ODE
u_m   = u0;               % Anfangswert  
t_m   = tspan(1);         % Anfangszeit
tend  = tspan(end);       % Endzeit

% Schrittzaehler
nmax     = opts.MaxStepNum; % maximale Anzahl von Schritten
nstep    = 0;               % Gesamtschrittzahl
nacc     = 0;               % Anzahl akzeptierter Schritte
fcn      = 0;               % Anzahl Funktionsaufrufe
% Anzahl Fkt.aufrufe fuer w_i in jedem Integrationsschritt
fcn_inc  = sum(b_par{1,3}) + not(aut)*length(b_par{1,3}); 
% Anzahl der Orthogonalisierungen der rechten Seiten je Schritt 
n_orth_b = sum(b_par{1,1})-1;  

% Toleranzen
atol  = opts.AbsTol(:); % Toleranzen fï¿½r Schrittweitensteuerung bzw. 
rtol  = opts.RelTol(:); % Abbruch im 1ten Arnoldi-Prozess
theta = opts.KryTol;    % Toleranz Residuum 1ste Stufe
reTOL = opts.ReorthTol; % Reorthogonalisierungstoleranz (L-Kriterium)

% Schrittweitenparameter
% erste Schrittweite
if isempty(opts.InitialStep)
    h_new = 0.1*(sqrt(n)/norm(1./(atol+rtol.*abs(u_m))))^0.25;  
else
    h_new = opts.InitialStep;
end
h_min = opts.MinStep;     % wenn h < h_min Abbruch
h_m   = h_new;            % Schrittweite fuer Eingangstest (h_m > h_min)
alpha = opts.StepControl; % Schrittweitensteurungsparameter

% defining the the function for Jacobian*v
funcAv = opts.Jacv;

% defining the computation of the time derivativ of rhs(t,y)
Fdt = opts.Fdt;

% Parameter fuer grafischen Ausgabe des Fehlers in Q
grafQ   = opts.LinSysErrPlot; % grafische Ausgabe des Fehlers in Q aktivieren?
gQx     = 0;                  % x-Koordinate des naechsten Plotteils initialisieren
n_orth  = [0,0];              % Anzahl der Nachorthos und Gesamtorthos zaehlen 
Wi      = zeros(n,s);         % Matrix fuer rechte Seiten, fuer Nachorthoplot umgeformt
res     = zeros(kd(2),1);     % Vektor fuer residuen in erste Stufe fuer Nachorthoplot
Kerr    = zeros(n,s);         % Steigungwertmatrix fuer Nachortoplot
indHerr = cell(s,1);          % cellarray fuer Indices in H
if grafQ 
    colgQ = colormap(hsv(s)); % Farben fuer grafQ-Plot
end
gQlines = {false , false       % ob Linien nach jeder Krylovdim gezogen werden sollen [dim(1:s-1) dim(s)]
           '--k' , 'k' };     % und welchen Style sie haben sollen

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------

% Spaltenidices fï¿½r mit H berechnete Q erzeugen
posQ_init = colQ(kd,s,s_K,b_par{1,1});

% Stufenuebergangsparameter, st ... stage
st = struct('kd',    kd,...
            'K',     ones(s_K,1),...
            'posQ',  posQ_init,...
            'b_in',  b_par{1,1},...
            'b_exp', b_par{1,2},...
            'cf',    b_par{1,3},...
            's',     s               );

% Schleife ueber tspan = [t0,te]
rmintern = tic;
while (t_m < tend) && (h_new > h_min) && (nstep < nmax)
    err = 2;
    rejected = false;

    norm_u_m   = norm(u_m)/sqrt_n;
    [f_m,par{1:npar(2).rhs}] = rhs(t_m,u_m,par{1:npar(1).rhs});
    fcn  = fcn  + 1;          % Anzahl Fkt.Aufrufe
    
    % bei nichtautonomer rechter Seite Ableitung der rechten Seite
    % nach t durch Differenzenquotient berechnen oder durch vom
    % Benutzer definierter funktion Fdt
    if ~aut
        [ft,par{1:npar(2).Fdt}] = Fdt(rhs,t_m,u_m,f_m,par{1:npar(1).Fdt});
        Fdti = ft*d_i;
    end
    
    % Schleife fï¿½r Fehler
    while (err > 1 || isnan(err)) && (h_new > h_min) && (nstep < nmax)
        
        st.b_in    = b_par{1,1};
        st.b_exp   = b_par{1,2};
        if ~rejected
            dimre  = {NaN(1,s_K),0};% Vektor der Reotho. und num_b_in speichert
%             Giv(1:s,1) = {[]};
            st.kd      = kd_init;
            st.posQ    = posQ_init;            
            res        = zeros(kd(2),1);
        else
            st.kd(1,2:s+1) = cumsum([st.kd(2),kd_inc(2:s)]);
            st.posQ    = colQ( st.kd, st.s, s_K, st.b_in);
            dimre{1,1}(st.kd(2)+1:s_K) = NaN;% Vektor der Reotho. und num_b_in speichert
            dimre{1,2} = 0;
        end
        h_m        = h_new;
        delta      = h_m*gamma;
               
        % Schleife ueber die Stufen 1 bis s        
        for i=1:s
            
            % Berechnung f_i wenn nicht identisch mit f_i aus letzter Stufe
            if st.cf(1,i)
                if i > 1
                    % Stufenwert
                    u_mi = u_m + h_m * K(:,1:i-1) * a_ij(1:i-1,i);
                    % Funktionswert vom Stufenwert
                    [f_i,par{1:npar(2).rhs}] = rhs(t_m + c_i(i)*h_m ,...
                        u_mi,par{1:npar(1).rhs});
                end
                                        
            end % Ende neues w_i kommt in Krylovr.
            
            if i > 1
                w_i = f_i + K(:,1:i-1) * c_ij(1:i-1,i);
            else
                w_i = f_m; % rechte Seite in erster Stufe
            end

            % bei nichtautonomer DGL, Term mit Ableitung nach der
            % Zeit multipliziert mit d_i zur rechten Seite addieren
            if ~aut;  w_i = w_i + h_m*Fdti(:,i);  end
            
            Wi(:,i) = w_i;
            
            % GS dieser Stufe mit Hilfe von Stage loesen
            [Q,H,Hs,K_i,st,Giv,dimre,par,res,Qrej] = ...
                Stage(u_m,norm_u_m,t_m,f_m,h_m,n,Q,H,Hs,w_i,delta,Giv,...
                st,i,reTOL,atol,rtol,theta,s_K,dimre,rhs,...
                par,npar,funcAv,res,rejected);
            
            % fuer Fehlerberechnung Indices speichern
            if grafQ; 
                [indH,indHerr{i,1}] = ...
                         ind_pos( [st.b_exp(:,i)';st.kd(i:i+1)] , st.posQ);
            end
            
            Kerr(:,i) = K_i; % fuer Nachorthoplot
                     
            % tatsaechenlichen Steigungswert dieser Stufe berechnen, kd und
            % entsprechen posQ an erste Stufe anpassen            
            if i > 1
                K(:,i) = K_i - K(:,1:i-1)*c_ij(1:i-1,i);
            elseif ~rejected
                Qrejsave = Qrej;
                K(:,i)   = K_i;
                if  st.kd(2) ~= kd_init(2)                    
                    % Krylovdimensionen an erste Stufe anpassen
                    st.kd(1,2:s+1) = cumsum([st.kd(2),kd_inc(2:s)]);                    
                    % Spaltenindices in Q fuer H an kd anpassen
                    st.posQ      = colQ(st.kd,s,s_K,st.b_in);                    
                end                
            else
                K(:,i) = K_i;
                Q(:,st.posQ(st.kd(2))) = Qrejsave;
            end
        end         % Ende der Schleife ueber die Stufen
        
        % Orhogonalisierung, Fkt.aufrufe zaehlen
        if rejected; n_orthi = st.kd(1,s+1) + sum(st.b_in) - 1 - st.kd(2);
        else         n_orthi = st.kd(1,s+1) + sum(st.b_in) - 1 ;
        end
        nreorth = sum(dimre{1,1}>0);
        n_orth = n_orth + [nreorth , n_orthi];
        
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>        
        % grafische Ausgabe des Fehlers in Q mit evtl. NO

        if grafQ
            gQx = errQplot(Q,res,delta,Kerr,K,Wi,c_ij,dimre{1,1},...
                st.kd,indHerr,gQx,nstep,colgQ,gQlines,par,H);
            % gQx = errQplot(Q,res,delta,Kerr,Ki,Wi,c_ij,dimre,kd,iH,...
            %     gQx,nstep,colgQ,gQlines,par,H_map)
        end

% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
        
        u_new  = u_m + h_m * K * b;       % u_{m+1} wird berechnet
        u_diff = h_m * K * be;            % Differenz zwischen Loesung des Verfahrens und                                            
       
        % Fehlerschaetzung --> [1] S.311
        % err  = sqrt( 1/n * sum( ( (u_diff) ./ (atol+rtol.*abs(u_m)) ).^2) );       
        err  = norm ( (u_diff) ./ (atol+rtol.*abs(u_m)))/sqrt(n);
        
        % neue Schrittweite --> [1] S.311
        h_new   = h_m * min( alpha(3),max(alpha(1)*err^(-0.25),alpha(2)) );  
        nstep   = nstep + 1; % zaehlen der Gesamtschritte
        
        kd1_ave = kd1_ave + st.kd(1,2);
        
        fprintf(['\n n = %3.0f  rej = %3.0f  t = %4.3f',...
                '  h_m = %6.2e  err = %6.2e  1st_kd = %2.0f',...
                '  NumRO = %3.0f numb_in = %2.0f b_in = %s'],...
                nstep,nstep-nacc-1,t_m,h_m,err,st.kd(2),nreorth,...
                dimre{1,2},sprintf('%1.0f',st.b_in))             
        if err > 1 || isnan(err)           
            fprintf('   Step rejected')
            if aut;
                rejected = true;
            end
        end
    end                       % Ende der Err-Schleife    
    
    u_old = u_m;
    u_m   = u_new;
    t_old = t_m;
    t_m   = t_m + h_m;
    f_old = f_m;
    if (t_m + h_new) > tend;  
        h_new = tend - t_m;
    end
    nacc = nacc + 1;          % zaehlen der akzeptierten Schritte
    
    % call of the user supplied function Solout
    if iout
        [par{1:npar(2).Solout}] ...
            = opts.Solout(t_old , t_m , u_old , u_m , f_old , f_m,...
                                            par{1:npar(1).Solout});
    end
    
end                           % Ende der Hauptschleife
Zeit = toc(rmintern);         % benoetigte Zeit fuer Integration

nrej = nstep-nacc;
nJacv  = n_orth(1,2) - n_orth_b*nstep;
fcn  = nstep*fcn_inc - nrej;
nFdt = nacc;

% stats.profile = profile('info');
% names = {stats.profile.FunctionTable.FunctionName}';
% calls = {stats.profile.FunctionTable.NumCalls}';
% 
% fnrhs  = func2str(rhs);
% fnrhs(strfind(fnrhs,'/')) = '';
% fnJacv = func2str(funcAv);
% fnJacv(strfind(fnJacv,'/')) = '';
% 
% if ~isempty(Fdt);    
%     fnFdt = func2str(Fdt);
%     fnFdt(strfind(fnFdt,'/')) = '';
% end
% 
% j = strfind(names,'>');
% 
% for i = 1:length(names)
%     names{i}(j{i}) = '';
%     if length(names{i}) >= length(fnJacv)
%         if strcmp(names{i}(end-length(fnJacv)+1:end),fnJacv)
%             stats.Jacv = calls{i};
%         end
%     end
%     if ~isempty(Fdt)
%         if length(names{i}) >= length(fnFdt)
%             if strcmp(names{i}(end-length(fnFdt)+1:end),fnFdt)
%                 stats.nFdt = calls{i};
%             end
%         end
%     end
%     if length(names{i}) >= length(fnrhs)
%         if strcmp(names{i}(end-length(fnrhs)+1:end),fnrhs)
%             stats.FCN = calls{i};
%         end
%     end
% end
% stats.FCN = stats.FCN - stats.Jacv;

stats.nFdt   = nFdt;
stats.FCN    = fcn;
stats.nJacv  = nJacv;
stats.nsteps = nstep;
stats.nacc   = nacc;
stats.nrej   = nrej;
stats.Time   = Zeit;
stats.nOrth  = n_orth(2);
stats.nReorth      = n_orth(1);
stats.Ave_iter = kd1_ave/nstep;

if h_new <= h_min && t_m < tend-eps
    str = 'stepsize to small';
    stats.Info = str;
    fprintf('\n%s\n\n',str)
elseif nstep >= nmax
    str = 'number of steps exceeded the given maximum';
    stats.Info = str;
    fprintf('\n%s\n\n',str)
elseif t_m >= tend-eps  
    str = 'Computation successful';
    stats.Info = str;
    fprintf('\n%s\n\n',str)
else
    str = 'Computation interrupted by a unkown reason';
    stats.Info = str;
    fprintf('\n%s\n\n',str)
end

varargout = {par{1:ninout(1).noutpar}};

if grafQ
    % figure an Bildschirm anpassen
    % set(0,'Units','normalized')
%     set(0,'Units','pixels')
%     scrsz = get(0,'ScreenSize');
%     set(gcf,'Position',[0 scrsz(4)*8/10 scrsz(3) scrsz(4)*8/10])
    
    % Anzahl NO anzeigen
    t = [num2str(n_orth(1,1)),'(',num2str(n_orth(1,2)),')'];
    text(-0.6,10^-16.5,t,'HorizontalAlignment','right','FontSize',20)
end

% st.K = st.K(1:st.kd(s+1));
% posb = st.kd(1:s)+1;
% posb = posb(st.b_in);
% indb = 1:s; indb = indb(st.b_in);
% st.K(posb) = indb;
% 
% for i = 1:s
%     for j = st.kd(i) + 1 + st.b_in(i) : st.kd(i+1)
%         st.K(j) = st.K(j - i + ~st.b_in(i));
%     end
% end
% 
% % disp([(1:st.kd(s+1))',st.K]);
% stats.K = st.K;

% *************************************************************************
% =========================================================================
% --------------------------------- Stage ---------------------------------
% =========================================================================
% *************************************************************************

function [Q,H,Hs,xs,st,Giv,dimre,par,res_out,Qrej] = ...
              Stage(u_m,norm_u_m,t_m,f_m,h_m,n,Q,H,Hs,b,delta,Giv,...
                    st,i,reTOL,atol,rtol,theta,s_K,dimre,rhs,...
                    par,npar,funcAv,res_out,rejected)
                      
Qrej = [];
% rechte Seite initialisieren
bs = zeros(st.kd(i+1),1);
% bs1 = zeros(st.kd(i+1),1);

% neue rechte Seite in Q einfuegen, wenn fuer diese Stufe nicht schon in K
if st.b_in(i)
    
    posb = st.kd(i) + 1;   % Position des b_i in Q
    [ Q(:,posb), dimre, bs(1:posb,1), b_in] = ...
        orthov( Q, posb, true, b, [], [1 1], [posb,st.kd(i)],...
                                                          reTOL, dimre, 0);
%     [v,dimre,h,b_in] = orthov(Q,posQ,genQ,v,h_in,m,n,TOL,dimre,cn)

    if b_in && i > 1
        
        % 1. st anpassen
        st.b_in(i)        = false;
        st.b_exp(2*i:end) = st.b_exp(2*i:end) - 1;        

        % Krylovdimensionen anpassen
        st.kd(i+1:end) = st.kd(i+1:end)-1;

        % 2. Q anpassen
        posQnew              = colQ( st.kd, st.s, s_K, st.b_in);
        ind                  = st.kd(i) + ( -st.b_exp(1,i) + 1 : 0);                        
        Q(:,posQnew( ind ) ) = Q(:,st.posQ( ind ) );        
        st.posQ              = posQnew;
        
    end
else
    posb = st.kd(i) + 1;   % Position des b_i in Q
    [ v, dimre, bs(1:posb,1)] = ...
        orthov( Q, posb, true, b, [], [1 1], [posb,st.kd(i)],...
        reTOL, dimre, 0);
end

if rejected && i == 1
    Hs(1:st.kd(2),1:st.kd(2)) = eye(st.kd(2)) - ...
                                            delta*H(1:st.kd(2),1:st.kd(2));
    
    Giv_i = cell(1,st.kd(2)-1); % Givensdrehungsmatrix fuer diese Stufe
    for j = 1:st.kd(2)-1
        Giv_i{1,j} = [ zeros(2,1); j ; j+1];
    end
    
    % Berechnung der neuen Givensrotationen fuer diese Stufe
    [Hs(1:st.kd(2),1:st.kd(2)),Giv_i] = ...
        givrot(Hs(1:st.kd(2),1:st.kd(2)),1:st.kd(2)-1,Giv_i);
    
    Giv{i,1} = [Giv_i{1,:}]; % Givensdrehungen fuer diese Stufe speichern
    
    cs_pq = [Giv{1:i,1}]; % Drehungen fuer bsolve aus Giv auslesen
    for j = 1:size(cs_pq,2);
        bs(cs_pq(3:4,j),1) = [ cs_pq(1:2,j)' ; [-cs_pq(2,j),cs_pq(1,j)] ]...
            * bs(cs_pq((3:4),j),1);
    end
    % ------ Anfang Loesung xs berechnen ------
    ks = zeros(st.kd(2),1); % Koeffizientenvektor initialisieren
    % Rueckwaertssubstitution
    ks(st.kd(2)) = bs(st.kd(2))/ Hs(st.kd(2),st.kd(2));
    for j = st.kd(2)-1:-1:1
        ks(j,1) = (bs(j,1) - Hs(j,j+1:st.kd(2))*ks(j+1:st.kd(2),1)) / ...
                                                                   Hs(j,j);
    end
    % x fuer jeweilige Stufe ausrechnen
%     Giv{i,1}
%     Hb = [H(1:st.kd(2),1:st.kd(2)),b(1:st.kd(2))]
%     Hsbs = [Hs(1:st.kd(2),1:st.kd(2)),bs(1:st.kd(2)),ks]
    xs = Q( : , 1:st.kd(2)) * ks; % x = Q*l
    return
end
% Krylovdimensionen vor und nach der Stufe
kd = st.kd(i:i+1);
% Berechnung d.Indices fuer alle Berechnungsabschnitte (BA) in dieser Stufe
[indH,indH_all] = ind_pos( [st.b_exp(:,i)';st.kd(i:i+1)] , st.posQ);

% Indices fuer Zeilen in indH bzw. indH_all
zHi  = 1; % Spaltenindex H
zf0  = 2; % 1ter ~=0 Eintrag
zfGD = 3; % Index 1te Element GD
znGD = 4; % Anzahl GD
z1d  = 5; % +1 auf diag
zcQ  = 6; % Spaltenindices von Q
zgQ  = 7; % Q(:,i) erzeugen (ja/nein)=(1,0)
zl0  = 8; % letzter ~=0 Eintrag

% Spalten ab 2.Stufe Spalten in H vervollstï¿½ndigen u. 
% evt. zugehï¿½rige Spalten in Q erzeugen
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% disp('-----------------------------------------------------------------')
% disp('------------ 1. Spalten von H und Q vervollstï¿½ndigen ------------')
% disp('-----------------------------------------------------------------')
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
if i > 1
    iH = indH{1,1};
    for j = 1:size(iH,2)
        h_in = H(1:iH(zf0,j)-1,iH(zHi,j));
        posQ = iH(zcQ,j);
        [Q(:,posQ),dimre,H(iH(zf0,j):iH(zl0,j),iH(zHi,j))]...
           = orthov(Q,posQ,iH(zgQ,j)==1,Q(:,iH(zcQ,j)),h_in ,...   
              [iH(zf0,j),1],[iH(zl0,j) , iH(zl0,j)-iH(zgQ,j)],reTOL,dimre,1);
%            orthov(Q, posQ     ,genQ        ,v             ,h_in ,             
%              m            ,n                               ,TOL  ,dimre);
%              
    end
end

% Spalten in Q und H die vollstï¿½ndig erzeugt werden
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% disp('-----------------------------------------------------------------')
% disp('------------ 2. Spalten von H und Q normal berechnen ------------')
% disp('-----------------------------------------------------------------')
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

iH = indH{1,2};
sziH = size(iH,2);
if i == 1; Giv_i = cell(1,kd(2)); end % cellarray fuer GD in 1.Stufe

for j = 1:sziH
    
    % in diesem Teilabschnitt sind in den jeweiligen Spalten in H noch
    % keine Elemente vorhanden deshalb h_in = [];
    h_in = [];
    
    % Vektor der zu Q orthogonalisiert werden soll    
    [v,par{1:npar(2).Jacv}] = funcAv(rhs,t_m,u_m,norm_u_m,f_m,...
                                   Q(:,iH(zHi,j)),par{1:npar(1).Jacv});
    
    posQ = iH(zcQ,j);
    [Q(:,posQ),dimre,H(1:iH(zl0,j),iH(zHi,j))]...
         = orthov(Q,posQ,iH(zgQ,j)==1,v ,h_in ,[iH(zf0,j),1],[iH(zl0,j) , iH(zl0,j)-iH(zgQ,j)],reTOL,dimre,2);
    %      orthov(Q,posQ     ,genQ        ,v ,h_in ,m            ,n
    %      ,TOL);
    
    % nur in erster Stufe das Residuum ausrechnen
    if i == 1
        
        [res,H,Hs,bs,Giv_i] = residuum(H,Q,delta,Hs,bs,iH(zHi,j),iH(zl0,j),Giv_i);
%       [res,H,Hs,bs,Giv  ] = residuum(H,Q,delta,Hs,bs,   j     ,   k     , Giv)

        % Residuum fuer Nachortho.-Fehlerplot speichern
        res_out(iH(zHi,j),1) = norm(res);

        if j < sziH
            % Abrruchtest des Krylovraums der ersten Stufe --> [1] (2.4)
            st1TOL = h_m*sqrt(1/n*sum((res./(atol+rtol.*abs(u_m))).^2));

            if (st1TOL < theta && j > 3)
                % letztes Q unnormiert in richtige Spalte schreiben
%                 Q(:,iH(zcQ,j)) = Q(:,iH(zcQ,j)) * H(iH(zl0,j),iH(zHi,j));
                Q(:,iH(zcQ,j)+1) = Q(:,iH(zcQ,j)) * H(iH(zl0,j),iH(zHi,j));
                Qrej = Q(:,iH(zcQ,j)+1);

                Giv_i(:,j) = [];  % nicht benoetigte GD loeschen
                kd(2)      = j;   % Krydim der 1. Stufe wenn theta erreicht
                st.kd(2)   = j;
                break             % Abbruch da Toleranz theta erfuellt
            end

            % Berechnung der neuen Givensrotationen fuer die k-te Spalte, aber
            % nur wenn letzte Spalte fuer diese Stufe noch nicht erreicht
            [Hs(1:iH(zl0,j),iH(zHi,j)),Giv_i] = ...
                givrot(Hs(1:iH(zl0,j),iH(zHi,j)),iH(zHi,j),Giv_i);   
        else
            Qrej = Q(:,iH(zcQ,j));
        end    
    end
end

% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% disp('-----------------------------------------------------------------')
% disp('---------- 3. Matrix H zur Lsg des GS vervollstï¿½ndigen ----------')
% disp('-----------------------------------------------------------------')
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
iH = indH{1,3};
if i>1
    % vervollstï¿½ndigen von H mit unvollstï¿½ndigen Spalten in Q
    for j = 1:size(iH,2);
        
        h_in = [];
        [v,par{1:npar(2).Jacv}] = funcAv(rhs,t_m,u_m,norm_u_m,f_m,...
                                   Q(:,iH(zHi,j)),par{1:npar(1).Jacv});
        
        posQ = iH(zcQ,j);
        [Q(:,posQ),dimre,H(1:iH(zl0,j),iH(zHi,j))]...
        = orthov(Q,posQ,iH(zgQ,j)==1,v ,h_in ,[iH(zf0,j),1],[iH(zl0,j) , iH(zl0,j)-iH(zgQ,j)],reTOL,dimre,3);
        % orthov(Q,posQ     ,genQ        ,v               ,h_in ,m
        % ,n                              ,TOL);
        
    end
    
    transH = indH_all([zHi zf0 zl0],:); % ï¿½bertragungsindices von H nach Hs
    delHs  = indH_all([zHi znGD zfGD zl0],:); % Lï¿½schungsindices in Hs
    
    % H in Hsolve ï¿½bertragen
    for j = 1:size(transH,2)
        k = transH(1,j);                   % Spaltenindex
        rows = transH(2,j):transH(3,j);    % diesen Zeilen uebertragen
        Hs(rows,k) = -delta*H(rows,k);     % uebertragung
    end
        
    % Spalten wo 1 add. werden muss u. Vortrans. durchgefuehrt werden muss
    k = transH(1,indH_all(z1d,:)==1);
    
    % eine Eins zu den entsprechenden Diag.elementen in Hsolve add.                                    
    for j = k
        Hs(j,j) = Hs(j,j) + 1;
    end

    % vortransformieren mit alten GD
    cs_pq = [Giv{1:i-1,1}]; % alte GD fuer Hsolve aus Giv auslesen    
    for j = 1:size(cs_pq,2)
       Hs(cs_pq(3:4,j),k) = [ cs_pq(1:2,j)' ; [-cs_pq(2,j),cs_pq(1,j)] ]...
                                *Hs(cs_pq(3:4,j),k);
    end
    
    % neue Indices fï¿½r Givensdrehungen dieser Stufe Erzeugen
    szdelHs = size(delHs,2);
    Giv_i = cell(1,szdelHs); % Givensdrehungsmatrix fï¿½r diese Stufe
    for j = 1:szdelHs
        Giv_i{1,j} = [ zeros(2,delHs(2,j))
                       delHs(1,j)*ones(1,delHs(2,j))
                       delHs(3,j):delHs(4,j)             ];
    end

    % Berechnung der neuen Givensrotationen fï¿½r diese Stufe        
    [Hs(1:kd(2),delHs(1,:)),Giv_i] = ...
        givrot(Hs(1:kd(2),delHs(1,:)),1:szdelHs,Giv_i);
    
    Giv{i,1} = [Giv_i{1,:}]; % Givensdrehungen fï¿½r diese Stufe speichern
        
    % b fuer das zu loesende Gleichungssystem transformieren     
%     bs1(1:kd(2),1) = Q(:,1:kd(2))'*b;
%     
%     [i kd(1)]
%     [(1:kd(2))', bs1(1:kd(2),1), bs(1:kd(2),1), bs1(1:kd(2),1)-bs(1:kd(2),1)]
%     errbs = norm(bs1(1:kd(2),1)-bs(1:kd(2),1))
    
    cs_pq = [Giv{1:i,1}]; % Drehungen fuer bsolve aus Giv auslesen   
    for j = 1:size(cs_pq,2);
       bs(cs_pq(3:4,j),1) = [ cs_pq(1:2,j)' ; [-cs_pq(2,j),cs_pq(1,j)] ]...
                            * bs(cs_pq((3:4),j),1);
    end
else
    Giv{1,1} = [Giv_i{:}];
%     Giv{1,1}
end

% ------ Anfang Loesung xs berechnen ------
ks = zeros(kd(2),1); % Koeffizientenvektor initialisieren
% Rueckwaertssubstitution
ks(kd(2)) = bs(kd(2))/ Hs(kd(2),kd(2));
for j = kd(2)-1:-1:1
    ks(j,1) = (bs(j,1) - Hs(j,j+1:kd(2))*ks(j+1:kd(2),1)) / Hs(j,j);
end
% x fuer jeweilige Stufe ausrechnen
% if true
%     Giv{i,1}
%     Hb = [H(1:st.kd(i+1),1:st.kd(i+1)),b(1:st.kd(i+1))]
%     Hsbs = [Hs(1:st.kd(i+1),1:st.kd(i+1)),bs(1:st.kd(i+1)),ks]
% end
xs = Q( : , 1:kd(2)) * ks; % x = Q*l

% ****** Ende Lï¿½sung xs berechnen ******


% #########################################################################
% =========================================================================
% -------------------------------------------------------------------------
% ------------------------ zusaetzliche Funktionen ------------------------
% -------------------------------------------------------------------------
% =========================================================================
% #########################################################################

% -------------------------------------------------------------------------
% -- Orthogonalisierung von v zu Q bzw. erzeugen von Elementen in H u. Q --
% -------------------------------------------------------------------------
                                    
function [v,dimre,h,b_in] = orthov(Q,posQ,genQ,v,h_in,m,n,TOL,dimre,cn)
           
% Orthogonalisierung von v zu Q mit evt. Reorthogonalisierung in Abh. v.TOL
% Erzeugung d.Elemente v.H u.wenn genQ == true Erzeugung d.Spalte posQ in Q

h    = zeros(n(1),2); % h initialisieren

if m(1) > 1;  h(1:m(1)-1) = h_in; end % schon berechnete hij in h uebertragen

quot   = 0;           % quot so setzen das Schleife ausgefuehrt wird
reorth = true;
i      = 0;             % Reorthogonalisierungszaehler init.
b_in   = false;

while reorth    
    
    i = i + 1;                       % Reorthogonalisierungszaehler
    
    if n(2)>0
        for j = n(2) : -1 : m(i)          % v zu Q orthogonalisieren
            h(j,i) = Q(:,j)' * v;
            v      = v - h(j,i)*Q(:,j);
        end
    else 
        TOL = false;
    end
    
    if genQ  % norm(v) berechnen und quot fuer reorthogonalisieren berechnen
        h(n(1),i) = norm(v);
        
        if i < 2 && TOL
            
            % wenn TOL == true immer Nachortho egal wie gross quot ist
            if islogical(TOL)   
                reorth = true;
            else
                quot   = sum(abs(h(1:n(2),i)))/h(n(1),i);
                reorth = quot > TOL;
            end
            
            if cn == 0 % wenn es der Berechnungsabschnitt fuer b_i ist
                if quot > 1e+10            % b_i schon in K_{i-1} ?
                    b_in   = true;         % ja
                    reorth = false;        % keine Nachortho erforderlich
                    dimre{1,2} = dimre{1,2} + 1; % Anzahl der b_i in K                    
                    genQ       = false;
                end
            end    
        else
            reorth = false; 
            if i==2
                dimre{1,1}(1,posQ) = posQ; % Dim. der Reotho. speichern
            end
        end                
    else
        reorth = false; % Nachorthogonalisierung ausschalten
    end
end

% Erzeugen von q wenn dies mit genQ vorgesehen ist und es kein b_i ist was
% schon drin ist
if genQ; v = v / h(n(1),i); end

% h fuer Ausgabe anpassen u. uebergeben, wenn BA nicht fuer b_i
h = h(m(1,1):n(1,1),1); 

% ?????????????????????????????????????????????????????????????????????????
% h durch reorthogonalisieren verbessern
% if i == 2
%     h(n(1),1) = 0;
%     h         = sum(h,2);
% end

% -------------------------------------------------------------------------
% --                           Residuum berechnen                        --
% -------------------------------------------------------------------------

function [res,H,Hs,bs,Giv] = residuum(H,Q,delta,Hs,bs,j,k,Giv)
% nur in der ersten Stufe das residuum ausrechnen

% [res,H,Hs,bs,Giv_i] =
% residuum(H,Q,delta,Hs,b,bs,iH(zHi,j),iH(zl0,j),Giv_i);

% Erzeugung von Hs und bs
Hs(1:k,j) = -delta*H(1:k,j);
Hs(j,j)   = Hs(j,j) + 1;

% vortransfomieren der rechten Seite mit alten Givensdrehungen
% bei b muessen nicht alle alten Givensdrehungen erneut angewendet
% werden sondern nur die, die im letzten Schritt dazu gekommen sind
if j>1
    
    cs_pq = Giv{1,j-1};
    bs(cs_pq(3:4,1),1) = [cs_pq(1:2,1)';[-cs_pq(2,1),cs_pq(1,1)]]...
                            *bs(cs_pq((3:4),1),1);

%   vortransfomieren der k-ten Spalte mit alten Givensdrehungen
%   alle alten Givensdrehungen muessen erst auf die k-te Spalte
%   angewendet werden
    cs_pq = [Giv{1,1:j-1}];
    for i = 1:j-1
        Hs(cs_pq(3:4,i),j) = [cs_pq(1:2,i)';[-cs_pq(2,i),cs_pq(1,i)]]...
                             *Hs(cs_pq(3:4,i),j);
    end
end

% Berechnen der letzten Komponente des Koeffizientenvektors l fï¿½r
% Residuum
k_j = bs(j)/Hs(j,j);

% Berechnung des Residuums
if j < k;  
    res      = -delta* H(k,j) * Q(:,k) * k_j;
    Giv{1,j} = [0 0 j k]';
else
    res      = NaN(length(Q(:,j)),1);
end

%-------------------------------------------------------------------------- 
% Givensrotationen fuer eine Matrix H in verallgemeinerter Hessenbergform, 
% wobei nur neue Spalten von H in Form von A betrachten werden
%--------------------------------------------------------------------------

function [A,Giv] = givrot(A,k,Giv)
% A ... Matrix in der Elemente ab Zeile k+1 zu Null gemacht werden
% m ... Anzahl der Zeilen von A
% k ... gibt die Nummer der Spalte von H an bzw. welche Givensrotationenen 
%       aus Giv  benï¿½tigt werden 
% Giv ... cell array der fï¿½r die vorzunehmenden Givensdrehungen die Indices
%         der zu benutzenden Zeilen enthï¿½lt, das k-te Element von Giv
%         gehï¿½rt zur k-ten Spalte von H

%[m,n] = size(A);   % Auslesen der Grï¿½ï¿½e von A    | c   s |
G = zeros(2);      % Anfangsgivensmatrix     G = |       |
                   %                             |-s   c |
mA=size(A,1);
% % 1. Anwenden der alten Givensrotationen
% for i = 2:cs_pq(1,1)+1
%     A(cs_pq(3:4),i),:) = [cs_pq(1:2,i);[-cs_pq(2,i),cs_pq(1,i)]]...
%                          *A(cs_pq((3:4),i),:);
% end
%if length(k)==1;n=1;else n=length(k)-1;end
n=length(k);
% 2. Erzeugen und Anwenden der neuen Givensrotationen
for j = 1:n        % Schleife ï¿½ber die Spalten von A    
    cs_pq = Giv{k(j)};
    m = size(cs_pq,2);
    i=1;
    while i<=m
    %for i = 1:m    % Schleife ï¿½ber die Elemente der jeweiligen Spalte 
                   % von A die Null werden sollen  

        if cs_pq(4,i) <= mA
            b = A(cs_pq(4,i) ,j); % Element das Null wird
            a = A(cs_pq(3,i),j); % Element das zum auslï¿½schen von b benutzt wird
            % Erzeugen der Givensrotation, dabei Fallunterscheidung zum
            % vermeiden von Unter- bzw. ï¿½berlauf
            if abs(a) < abs(b)
                tau = a/b;
                G(1,2) =  1 / sqrt(1 + tau^2);   % s
                G(2,1) = -G(1,2);                % -s
                G(1,1) =  G(1,2) * tau;          % c
                G(2,2) =  G(1,1);                % c
            else
                tau = b/a;
                G(1,1) =  1 / sqrt(1 + tau^2);   % c
                G(2,2) =  G(1,1);                % c
                G(1,2) =  G(1,1) * tau;          % s           
                G(2,1) = -G(1,2);                % -s
            end
            % Givensrotation speichern
            cs_pq( 1:2 , i) = G(1,:)';
            % Givensrotation auf Matrix A anwenden
            A( cs_pq(3:4,i), j:end ) = G*A( cs_pq(3:4,i) , j:end );
            i=i+1;
        else
            cs_pq(:,i) = [];
            disp('in Giv Spalte gelï¿½scht!!!!!!!!!!!!')
            m=m-1;
        end
    end
    Giv{k(j)}=cs_pq;
end

% -------------------------------------------------------------------------
% ---------------------------- Indexberechnung ----------------------------
% -------------------------------------------------------------------------

function [indH,indH_all,indb] = ind_pos(kd,posQ)
% berechnet Indices fï¿½r , Spaltenberechnung in H, GD, H --> Hs und
% Speicherung in Q

% indH: Zeilen-cellarray mit drei Spalten
%  1.Spalte: Indices fï¿½r Auffï¿½llen der vervollstï¿½ndigten Spalten
%  2.Spalte: Indices fï¿½r Berechnung mit normalen Arnoldi
%  3.Spalte: Indices fï¿½r vervollstï¿½ndigen von H

%  In jeder Spalte befindet sich eine Matrix mit 7 Zeilen:
%   1.Zeile: Spaltenindices v.H f.d.jew.Berechnungsabschnitt in d.jew.Stufe
%   2.Zeile: 1ter ~=0-Eintrag in jew.Spalte u.Stufe f.Berechnung H u.H-->Hs
%   3.Zeile: Index des 1ten Elements in jeder Spalte fï¿½r Givens
%   4.Zeile: Anzahl der GD in jeder Spalte im jew.Berechnungsabschnitt (BA)
%   5.Zeile: Addition einer 1 auf diag von Hs im jew. BA
%   6.Zeile: Spaltenindice fï¿½r Speicherung in Q
%   7.Zeile: Q(:,i) erzeugen ? (ja/nein) = (1/0)
%   8.Zeile: letzter ~=0-Eintrag in jew. Spalte u.Berechnungsabschnitt (BA)

% indH_all: eine aus d.3 Matrizen in indH f.d.jew.BA Zusammengefï¿½gte Matrix

% indb: Spaltenindice des aktuellen b_i in Q

% =========================================================================

if nargin<1; kd = [2 3; 10 15]; end

% indb erzeugen
indb = kd(2,1) + 1;

% indH initialisieren
n = 3; m = 8;
indH = cell(m,n);

% Spaltenindex von H fï¿½r den jew. Berechnungsabschnitt in der jew. Stufe
indH{1,1} = (diff(kd(:,1)) + 1 : kd(2,1))';
indH{1,2} = (kd(2,1) + 1 : diff(kd(:,2)) + (kd(1,1)==0)*1)';
indH{1,3} = (kd(2,1) + length(indH{1,2}) + 1 : kd(2,2))' ;

% erster Nichnulleintrag in jew. Spalte und zugehï¿½riger Stufe fï¿½r H --> Hs
indH{2,1} = ((kd(2,1)+1)*ones(1,kd(1,1)))';
indH{2,2} = (ones(1,length(indH{1,2})))';
indH{2,3} = (ones(1,length(indH{1,3})))';

% erster Index in jeder Spalte fï¿½r Givens
indH{3,1} = indH{2,1};
indH{3,2} = indH{1,2}+1;
indH{3,3} = indH{1,3}+1;

% letzter Nichtnulleintrag in zugehoeriger Spalte und Berechnungsabschnitt
for i=1:2;    indH{m,i}                      = indH{1,i} + kd(1,2);
              indH{m,i}(indH{m,i} > kd(2,2)) = kd(2,2);
end
indH{m,3} = (kd(2,2)*ones(1,length(indH{1,3})))';

% Anzahl der GD in jeder Spalte im jew. BA
indH{4,1} = indH{m,1} - indH{3,1} + 1;
indH{4,2} = indH{m,2} - indH{3,2} + 1;
indH{4,3} = indH{m,3} - indH{3,3} + 1;

% Addition einer 1 auf diag von Hs im jew. BA
indH{5,1} = zeros(length(indH{1,1}),1);
indH{5,2} =  ones(length(indH{1,2}),1);
indH{5,3} =  ones(length(indH{1,3}),1);

% Spaltenindice fï¿½r Speicherung in Q
indH{6,1} = posQ(indH{1,1});
indH{6,2} = posQ(indH{1,2});
indH{6,3} = posQ(indH{1,3});

% Q(:,i) erzeugen ?
indH{7,1} = indH{6,1} <= indH{m,1};
indH{7,2} = indH{6,2} <= indH{m,2};
indH{7,3} = indH{6,3} <= indH{m,3};

% Alle cell-Eintraege einer Spalte in einem zusammenfassen
for i = 1:n
    indH{1,i} = [indH{:,i}]';
end

% ueberfluessige cells abschneiden
indH = indH(1,:);

% alle BAe zusammengefasst
indH_all = [indH{1,1:3}];

% -------------------------------------------------------------------------
% -------------------------- Spaltenindices von Q -------------------------
% -------------------------------------------------------------------------

function pQ = colQ(kd,nkd,szK,b_in)
% Berechnet die Spaltenindices von Q in denen die mit Hilfe von H erzeugten
% Vektoren gespeichert werden und schreibt sie auf pQ

cQ    = (1:szK+sum(b_in))'; % alle Spaltenindices von Q
pb    = kd(1:nkd) + 1;      % Spaltenindices der b_i in Q
indpQ = true(size(cQ));     % logischer Indicevektor der Spalten in Q
indpQ(pb(b_in)) = false;    % Spalten der b_i rausstreichen
pQ    = cQ(indpQ);          % Spaltenindices d.mit Hilfe v.H erzeugten Q_i

% -------------------------------------------------------------------------
% ------------- approximation of Av by a difference quotient --------------
% -------------------------------------------------------------------------
function [Av varargout] = AvDQ(rhs,t,y,normy,fy,v,varargin)
par               = varargin;
nout              = nargout-1;
delta             = 10^-7 * max(10^-5 , normy);
[fdv par{1:nout}] = rhs(t,y + delta*v,par{:});
Av                = (fdv - fy)/delta;
varargout         = par;

function [ft varargout] = FdtDQ(rhs,t_m,u_m,f_m,varargin)
par              = varargin;
nout             = nargout-1;
dt               = sqrt(eps * max(1e-5 , abs(t_m) ) );
[ft,par{1:nout}] = rhs(t_m + dt , u_m, par{:});
ft               = (ft - f_m)/dt;
varargout        = par;

%--------------------------------------------------------------------------
%----------------------------- choose method ------------------------------
%--------------------------------------------------------------------------
function [aij,b,be,cij,ci,di,gij,kd,b_par]=rowcoe(method,aut,kd_max)
         
%----------------------------------GRK4T-----------------------------------
% --> [2] S.11
if  strcmpi(method,'grk4t')
    
    gij = [  0.231             0                  0                 0
            -0.27062966775244  0.231              0                 0
             0.31125448329409  0.00852445628482   0.231             0
             0.28281683204353 -0.45795948328073  -0.11120833333333  0.231];
         
    aij = [  0                 0                  0                 0
             0.462             0                  0                 0
            -0.08156681683272  0.96177515016606   0                 0
            -0.08156681683272  0.96177515016606   0                 0]';
    
    % Gewichte des Verfahrens
    b   = [  0.21748737165273  0.48622903799012   0  0.29628359035715]';
    % Gewichte der eigebetteten Methode
    be  = [ -0.71708850449933  1.77617912176104  -0.05909061726171  0]';
    
    % Maximal Krylov dimension increments per stage 
    if aut;        kd = [4 3 1 3]; % autonomous [4 7 8 11]
        % welche rechte Seite in Krylovraum eingefï¿½gt werden soll und welche
        % nicht, weil sie evtl. schon im Krylovraum enthalten ist
        b_in = [true true true false];
    else          kd = [4 5 1 5]; % nonautonomous [4 9 10 15]
        % welche rechte Seite in Krylovraum eingefï¿½gt werden soll und welche
        % nicht, weil sie evtl. schon im Krylovraum enthalten ist
        b_in = [true true true true];
    end
    
    % muss f_i in der jeweiligen Stufe berechnet werden oder nicht, da es
    % evtl. identisch zu einem alten f_i ist
    cf = [true true true false]; % compute f
    
%----------------------------------LSTAB-----------------------------------
elseif  strcmpi(method,'lstab')
    
    gij = [  0.57282000000000  0                 0                0
            -2.34201389131923  0.57282000000000  0                0
            -0.02735980356646  0.21380314735851  0.57282000000000 0
            -0.25909062216449 -0.19059462272997 -0.22803686381559 0.57282000000000];
        
    aij = [  0                0                0    0   0
             1.14564000000000 0                0    0   0
             0.52092209544722 0.13429476836837 0    0   0
             0.52092209544722 0.13429476836837 0    0   0]';

    b  = [ 0.32453574762832 0.04908429214667 0.00000000000000 0.62637996022502 ]';
    be = [ 0.61994881642181 0.19268272217757 0.18736846140061 0                ]';
    
    % Maximal Krylov dimension increments per stage 
    if aut;       kd = [4 3 1 3]; % autonomous [4 7 8 11]
        % welche rechte Seite in Krylovraum eingefuegt werden soll und welche
        % nicht, weil sie evtl. schon im Krylovraum enthalten ist
        b_in = [true true true false];
    else          kd = [4 5 1 5]; % nonautonomous [4 9 10 15]
        % welche rechte Seite in Krylovraum eingefuegt werden soll und welche
        % nicht, weil sie evtl. schon im Krylovraum enthalten ist
        b_in = [true true true true];
    end
    
    % muss f_i in der jeweiligen Stufe berechnet werden oder nicht, da es
    % evtl. identisch zu einem/disc1/home/beck/ROWMAP/Bilder alten f_i ist
    cf = [true true true false]; % compute f
    
%----------------------------------RODAS-----------------------------------
elseif  strcmpi(method,'rodas')

    gij = [  0.25                  0                      0                     0                      0                  0
            -0.3543                0.25                   0                     0                      0                  0
            -0.1336025052681756   -0.1289749473182447e-1  0.25                  0                      0                  0
             0.1526849173006467e1 -0.5336562887504573    -0.127939288425601e1   0.25                   0                  0
             0.6981190951785019e1 -0.2092930097006117e1  -0.5870067663032753e1  0.73180680825385       0.25               0
            -0.2080189494180936e1  0.5957623556766833     0.1701617798267262e1 -0.8851451983588043e-1 -0.3786761399271284 0.25 ];

    aij = [  0                     0                     0                     0                  0          0
             0.386                 0                     0                     0                  0          0
             0.1460747075254179    0.639252924745821e-1  0                     0                  0          0
            -0.3308115036677301    0.7111510251682848    0.2496604784994454    0                  0          0
            -0.4552557186318031e1  0.1710181363241332e1  0.4014347332103172e1 -0.1719715090264738 0          0
             0.2428633765466988e1 -0.3827487337647846   -0.185572033092958e1   0.5598352992273763 0.25       0]';
       
    % Gewichte des Verfahrens
    b  = [ 0.3484442712860515  0.2130136219118987 -0.1541025326623185...
           0.4713207793914958 -0.1286761399271284  0.25                  ]';    
    % Gewichte der eigebetteten Methode
    be = [ 0.2428633765466988e1 -0.3827487337647846 -0.185572033092958e1...
           0.5598352992273763    0.25                0                   ]';

    % Maximal Krylov dimension increments per stage
    if aut                      
        kd = [4  5  4  1  1  1];    % autonomous [4  5  4  5  6  7] , [ 4 9 13 18 24 31]       
    else
        kd = [4  5  4  1  1  1];    % nonautonomous [4  5  4  5  6  7] , [ 4 9 13 18 24 31]       
    end

    % welche rechte Seite in Krylovraum eingefï¿½gt werden soll und welche
    % nicht, weil sie evtl. schon im Krylovraum enthalten ist
     b_in = true(1,6);
%      b_in = [true false(1,5)]
%      b_in = [true true false(1,4)]
%     b_in = [true(1,3) false(1,3)]

    % muss f_i in der jeweiligen Stufe berechnet werden oder nicht, da es
    % evtl. identisch zu einem alten f_i ist
    cf = true(1,6); % compute f
    
%---------------------------------RODASP-----------------------------------
elseif  strcmpi(method,'rodasp')
    
    aij = [  0                    0                     0                    0                   0    0
             0.75                 0                     0                    0                   0    0
             0.08612040081415219  0.1238795991858478    0                    0                   0    0
             0.7749345355073236   0.1492651549508680   -0.2941996904581916   0                   0    0
             5.308746682646142    1.330892140037269    -5.374137811655562   -0.2655010110278497  0    0
            -1.764437648774483    -0.4747565572063027   2.369691846915802    0.6195023590649829  0.25 0];

    betha_ij = [  0                      0                     0                    0                    0                  0
                  0                      0                     0                    0                    0                  0
                 -0.049392              -0.014112              0                    0                    0                  0
                 -0.4820494693877561    -0.10087955555555556   0.9267290249433117   0                    0                  0
                 -1.764437648774483     -0.4747565572063027    2.369691846915802    0.6195023590649829   0                  0
                 -0.080368370789113464  -0.056490613592447572  0.4882856300427991   0.5057162114816189  -0.1071428571428569 0];

    gij = betha_ij - aij;
    gij(logical(eye(6))) = 0.25;

    b = (aij(6,:) + gij(6,:))';
    be = [aij(5,1:5) + gij(5,1:5) , 0]';

    aij = aij';
    
    % Maximal Krylov dimension increments per stage
    if aut                          
        kd = [4  5  4  1  1  1]; % autonomous
    else                            
        kd = [4  5  4  1  1  1]; % nonautonomous
    end

    % welche rechte Seite in Krylovraum eingefï¿½gt werden soll und welche
    % nicht, weil sie evtl. schon im Krylovraum enthalten ist
    b_in = true(1,6);  
    
    % muss f_i in der jeweiligen Stufe berechnet werden oder nicht, da es
    % evtl. identisch zu einem alten f_i ist
    cf = true(1,6); % compute f

else
    error('Rowcoe:wrongMethod',['The method "',method,'", which is',...
          ' choosen in \n opts.method is not implemented.']);
end

%--------------------------------------------------------------------------
% fï¿½r alle Verfahren gemeinsam
di  = sum(gij,2);    % Parameter fï¿½r f_t bei nichtautonomen Problemen
cij = (gij/gij(1,1))';  % schon an Berechnung angepasste Gammas
ci  = sum(aij);    % Knoten
be  = b - be;        % Differenz der Gewichte, fï¿½r Fehlerschï¿½tzung

% defining the maximal first Krylov Dimension
if kd(1) <= kd_max; kd(1) = kd_max;
else% %         kd = [4 2*ones(1,5).*(2:6)+1];
%         kd = [4 20 20 20 20 20];
%         kd = [4 5 9 13 17 21];
%         kd = [4  5  4  3  3  3];

    warning('ROWMAP:first_kd_to_small','opts.MaxIter is to small, replaced by four.');
end

% Erweiterung des Krylovraums mit neuen rechten Seiten in der jew. Stufe
% b_exp ... b_expansion
nkd        = length(kd);                    % Anzahl der Stufen
b_exp      = zeros(2,nkd);                  % initialisieren
b_exp(2,:) = cumsum(b_in);                  % Anz. b_i in Q nach d.jew. Stufe
b_exp(1,:) = cumsum([0 , b_in(1:nkd-1)]);   % Anz. b_i in Q vor d.jew. Stufe

% cellarray von Parametern fuer rechte Seiten und berechnung der
% Funktionswerte
b_par = {b_in,b_exp,cf};




% -------------------------------------------------------------------------
% ---------------------- Fehler in Q mit Anzeige NO -----------------------
% -------------------------------------------------------------------------

function gQx = errQplot(Q,res,delta,Kerr,Ki,Wi,c_ij,dimre,kd,iH,...
                        gQx,nstep,colgQ,gQlines,par,H_map)
% plottet den Fehler in Q mit evtl. NO fï¿½r jeden Schritt

s   = length(kd);
par = par{1};
J   = par{1,end};

% wenn Spalte von H doppelt in iH vorkommt, Spalte aus kleinerer Stufe im 
% Krylovraum, also die mit kleinerem Spaltenindex in iH, loeschen da in
% Stage die gleiche Spalte in verschiedenen Berechnungsabschnitten
% betrachtet wird und wir benoetigen nur die Indices die fuer die letzte
% Berechnung benoetigt wurden

iH = [iH{:}];
i = 1;
while i < size(iH,2)
    if iH(1,i) >= iH(1,i+1);  iH(:,i) = [];
        if i>1; i=i-1;end
    else i = i+1;    
    end
end

% Berechnung von H aus Q und J
H = Q(:,1:kd(1,s))'*J*Q(:,1:kd(1,s));

% Elemente in H Null setzen die eigentlich Null sein sollen
for j = 1:kd(1,s)
    H(iH(8,j)+1:end,j) = 0;
end

H = H_map;
mat2file(H(1:kd(s),1:kd(s)),kd,'Hstrat1.txt');

% -------------------------------------------------------------------------
% ------------------ Fehler in Q fuer jede Dim berechnen ------------------
% -------------------------------------------------------------------------

Qvgl = Q(:,1:kd(1,s))'*Q(:,1:kd(1,s)); % Q^T*Q = I (normalerweise)
errQ = zeros(1,kd(1,s)); % Vektor fuer Orthogonalitaetsfehler in Q
ERR  = zeros(1,kd(1,s)); % Vektor für Residuum

K = zeros(size(Ki));

fprintf('\n Berechnung des Fehlers der Loesung    ')
for i = 1:s-1
    fprintf('  %2.0f',i)
    k0 = K(:,1:i-1) * c_ij(1:i-1,i);
    for j = kd(1,i)+1 : kd(1,i+1)
        errQ(1,j) = norm(eye(j)-Qvgl(1:j,1:j)); % Fehler in Q
        Hi = H(1:j,1:j);                        % Matrix H fuer diese Stufe
        % ErrH = norm(Hi-H_map(1:j,1:j));       % Fehler in H_map im Vgl. zu H = Q^T*J*Q
        % fprintf('\n ErrH = %10.5e',ErrH)
        Qi = Q(:,1:j);                          % Matrix Q fuer diese Stufe
        % Berechnung der Lsg. des Krylov-GS dieser Stufe
        ki = Qi*linsolve( eye(size(Hi)) - delta*Hi  , Qi' * Wi(:,i) );  
        k = ki - k0;        
        % Fehler des zu loesenden Gleichungssystems        
        ERR(1,j) = norm( ( speye(size(J)) - delta*J ) * ki - Wi(:,i));
    end
    K(:,i) = k; % Speichern der Steigungswerte
end

mst = struct('errQ'  ,'x-',...
             'resMAP','o-',...
             'errGS' ,'ks-',...
             'kMAP'  ,'p-');
   
% mst = struct('errQ'  ,'o-',...
%              'resMAP','+-',...
%              'errGS' ,'kx--',...
%              'kMAP'  ,'h-');

msz = 10;
lwd = 1.5;
         
% errQ fuer einzelne Stufen in unterschiedlichen Farben plotten
fprintf('\n Plotten der Ergebnisse                ')
for j = 1:s-1
    fprintf('  %2.0f',j)
    x = (kd(1,j)+1 : kd(1,j+1)); % Krylovdimensionen dieser Stufe
    x1 = gQx + x;
    
    % plotte Fehler in Q
    semilogy(x1,errQ(1,x),mst.errQ,'Color',colgQ(j,:),'MarkerSize',...
        msz,'Linewidth',lwd) 
    
    if j == 1;
        if nstep == 0;hold on;end
        % Fehler der Lsg. des GS aus Residuum des MAP in der 1.Stufe (j==1)
        % plotten
        semilogy(x1,res(1:length(x1))',mst.resMAP,'Color',colgQ(j,:),...
            'MarkerSize',msz,'Linewidth',lwd)     
    end
    % hier berechneter Fehler des GS plotten
    semilogy(x1,ERR(1,x),mst.errGS,'MarkerSize',msz,'Linewidth',lwd)
    
    % Fehler mit im Map berechneten Werten fuer K am Ende jeder Stufe
    Err = norm( ( speye(size(J)) - delta*J ) * Kerr(:,j) - Wi(:,j));   
    semilogy(gQx + x(end),Err,mst.kMAP,'Color',colgQ(j,:),...
        'MarkerFaceColor',colgQ(j,:),'MarkerSize',msz,'Linewidth',lwd)
end


% % legende hinzufï¿½gen
% if nstep == 0
%     leg = cell(s-1,1);
%     for j = 1:s-1 leg{j,1} = [num2str(j),'. Stufe'];end
%     legend(leg,'location','NWO')
% end

legend({'errQ','res_{MAP}','errGS-k_{new}','errGS-k_{map}'},...
        'Fontsize',16,'location','N')
xlabel('Krylovdimension','Fontsize',20)
ylabel('error','Fontsize',20)
title('Fehler bei der Loesung der lin. GSe im MAP','Fontsize',20)
set(gca,'Linewidth',2,'Fontsize',18,'FontWeight','bold')

% Nachorthogonalisierungen anzeigen
semilogy(gQx + dimre,10^-16.5,'mp','LineWidth',lwd,'Markersize',msz)
% wenn gewuenscht Trennlinien zeichnen zwischen Stufen
yl = ylim;
if gQlines{1,1}
    for j = 2:s
        x = gQx + kd([j j])+1;
        semilogy(x,yl,gQlines{2,1})
    end
end

% wenn gewuenscht Trennlinien zeichnen zwischen
% Integrationsschritten
if gQlines{1,2}
    x = gQx + kd([s s])+0.5;
    semilogy(x,yl,gQlines{2,2},'LineWidth',1)
end

% x-Achsen Zaehler, damit Fehler fortlaufend geplotted werden
gQx = gQx + kd(1,s);