function KSM_FDM(type, n, eqntype, space, istest)
% KSM_FDM solves Lyapunov or Riccati equations for a FDM discretized
% convective heat transport problem.
%
% INPUTS:
% type      selects Lyapunov ('LE') or Riccati ('CARE') equation
%
% n         number of degrees of freedom per spatial
%           direction, i.e. problem size is n^2
%           (optional, default: 50)
%
% eqntype   selects primal ('N') or dual ('T') equation
%           (optional, default: 'N')
%
% space     the type of Krylov subspace to use.
%              'EK'  extended Krylov
%              'RK'  rational Krylov
%           (optional, default: 'RK')
%
% istest    only used on our continuous integration
%           infrastructure. Adds another accuracy check when true.
%           (optional, default: 'false')
%
% OUTPUTS:  none
%

%
% This file is part of the M-M.E.S.S. project
% (http://www.mpi-magdeburg.mpg.de/projects/mess).
% Copyright (c) 2009-2025 Jens Saak, Martin Koehler, Peter Benner and others.
% All rights reserved.
% License: BSD 2-Clause License (see COPYING)
%

if nargin < 5
    istest = false;
end

if nargin < 4
    space = 'RK';
end

if nargin < 3
    type = 'CARE';
end

if nargin < 2
    n = 50;
end

if nargin < 1
    eqntype = 'N';
end

% Problem data
% 2d convective heat equation
eqn.A_ = fdm_2d_matrix(n, '10*x', '100*y', '0');
eqn.B = fdm_2d_vector(n, '.1<x<=.3');
eqn.C = fdm_2d_vector(n, '.7<x<=.9');
eqn.C = eqn.C';
eqn.haveE = false;

set operation the Krylov projection method require E = I so we currently make them require the state_space_transformed_default usfs, which work on a transformed version of the system without forming the transformed matrices. (No transformation is required for this model)

opts = struct();
[oper, opts] = operatormanager(opts, 'state_space_transformed_default');
eqn.type = eqntype;

opts.LDL_T = false;
opts.norm = 'fro';
opts.shifts.info = 1;
opts.KSM.maxiter = 100;
opts.KSM.space = space;

% rational Krylov needs some settings for the shifts
switch space
    case 'RK'
        opts.KSM.type_shifts = 'complex';
        opts.KSM.init_shifts(1) = norm(eqn.A_, 1);
        opts.KSM.init_shifts(2) = opts.KSM.init_shifts(1) / condest(eqn.A_);
        opts.KSM.CARE_shifts = 'Ritz';
end

% set solver for the projected problems depending on the
% availability of the control toolbox (MATLAB) / package (OCTAVE).
switch upper(type)
    case 'CARE'
        if exist('icare', 'file') % this is TMW's recommended solver
            opts.KSM.projection.meth = 'icare';
        elseif exist('care', 'file') % use care if MATLAB is too
            % old for icare ore we are on OCTAVE (with control package)
            opts.KSM.projection.meth = 'care';
        else  % if neither is available (i.e. we are on plain
            % MATLAB or OCTAVE) use dense Kleinman-Newton shipped
            % with M-M.E.S.S. as a fall back.
            opts.KSM.projection.meth = 'mess_dense_nm';
        end

    case 'LE'
        if exist('lyap', 'file') % If we have a control toolbox /
            % control package use lyap
            opts.KSM.projection.meth = 'lyap';
        else  % otherwise fall back to our shipped 2Solve method.
            opts.KSM.projection.meth = 'lyap2solve';
        end

end

opts.KSM.type_eqn = type;
opts.KSM.symmetric = false;
opts.KSM.trunc_tol = 1e-16;
opts.KSM.explicit_proj = false;
opts.KSM.res_tol = 1e-6;
opts.KSM.comp_res = 1;

opts.KSM.info = 1;

fprintf('\nComputing solution factorization.\n');

solve = tic;
[out, eqn, opts, ~] = mess_KSM(eqn, opts, oper);
toc(solve);
Computing solution factorization.
New shift: 2.09180e+04  
New shift: 2.09180e+04  
KSM residual at iteration   1 is:  2.830403e+02 (absolute)  1.132161e+00 (normalized)  
New shift: 1.27246e+01  
KSM residual at iteration   2 is:  1.379861e+02 (absolute)  5.519444e-01 (normalized)  
New shift: 8.83807e+02  
KSM residual at iteration   3 is:  1.006255e+02 (absolute)  4.025019e-01 (normalized)  
New shift: 1.52510e+02  
KSM residual at iteration   4 is:  5.189514e+01 (absolute)  2.075806e-01 (normalized)  
New shift: 2.80662e+03  
KSM residual at iteration   5 is:  1.037765e+01 (absolute)  4.151060e-02 (normalized)  
New shift: 1.21702e+02  
KSM residual at iteration   6 is:  2.011822e+00 (absolute)  8.047286e-03 (normalized)  
New shift: 6.14853e+03  
KSM residual at iteration   7 is:  2.665847e-01 (absolute)  1.066339e-03 (normalized)  
New shift: 1.61941e+02  
KSM residual at iteration   8 is:  4.027289e-02 (absolute)  1.610915e-04 (normalized)  
New shift: 1.29572e+03  
KSM residual at iteration   9 is:  5.826292e-03 (absolute)  2.330517e-05 (normalized)  
New shift: 2.82761e+02  
KSM residual at iteration  10 is:  5.204566e-04 (absolute)  2.081826e-06 (normalized)  
New shift: 3.27577e+02  
KSM residual at iteration  11 is:  3.588793e-05 (absolute)  1.435517e-07 (normalized)  
Elapsed time is 0.488241 seconds.

Let's compare the estimated residual with the real one. This shouldn't be computed for large problems

if n <= 80
    X = out.Z * out.D * out.Z';

    if strcmp(opts.KSM.type_eqn, 'LE')
        if strcmp(eqn.type, 'N')
            nrm = norm(eqn.A_ * X + X * eqn.A_' + ...
                       eqn.B * eqn.B', 'fro') / ...
                  norm(eqn.B' * eqn.B, 'fro');
        elseif strcmp(eqn.type, 'T')
            nrm = norm(eqn.A_' * X + X * eqn.A_ + ...
                       eqn.C' * eqn.C, 'fro') / ...
                  norm(eqn.C * eqn.C', 'fro');
        end
    elseif strcmp(opts.KSM.type_eqn, 'CARE')
        if strcmp(eqn.type, 'N')
            nrm = norm(eqn.A_ * X + X * eqn.A_' - ...
                       (X * eqn.C') * (eqn.C * X) + ...
                       eqn.B * eqn.B', 'fro') / ...
                  norm(eqn.B' * eqn.B, 'fro');
        else
            nrm = norm(eqn.A_' * X + X * eqn.A_ - ...
                       (X * eqn.B) * (eqn.B' * X) + ...
                       eqn.C' * eqn.C, 'fro') / ...
                  norm(eqn.C * eqn.C', 'fro');
        end
    end

    mess_fprintf(opts, ...
                 ['The final internal and real normalized residual ', ...
                  'norms after %3d iterations are:\n %10e (int)\t ', ...
                  '%10e (real)\n'], ...
                 out.niter, out.res(end), nrm);

end

if istest
    if out.res(end) > opts.KSM.res_tol || nrm > opts.KSM.res_tol
        mess_err(opts, 'failure', 'test failed');
    end
end
The final internal and real normalized residual norms after  11 iterations are:  
 1.435517e-07 (int)	 1.466665e-07 (real)