function KSM_Rail(k, type, transformed, eqtype, SISO, space, istest)
% KSM_RAIL solves Lyapunov or Riccati equations for the revised
% Oberwolfach steel profile benchmark problem.
%
% INPUTS:
% k            refinement level for the rail data. See
%              mess_get_linear_rail for allowed values.
%              (optional, default: 1)
%
% type         selects Lyapunov ('LE') or Riccati ('CARE') equation
%
% transformed  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 < 7
    istest = false;
end

if nargin < 6
    space = 'EK';
end

if nargin < 1
    k = 1;
end
% fetch rail model of appropriate size
eqn = mess_get_linear_rail(k);
eqn.haveE = true;

if nargin < 5
    SISO = false;
end
% For the SISO test pick one column of B and one row of C
if SISO
    eqn.B = sum(eqn.B, 2);
    eqn.C = sum(eqn.C, 1);
end
if nargin < 3
    transformed = false;
end
if transformed
    % Note that this is for demonstration purposes only and should never be
    % done explicitly, otherwise.
    eqn.A_ = eqn.E_ \ eqn.A_;
    eqn.B = eqn.E_ \ eqn.B;
    eqn.haveE = false;
    eqn.E_ = speye(size(eqn.A_, 1));
end

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.

opts = struct();
[oper, opts] = operatormanager(opts, 'state_space_transformed_default');
if nargin < 4
    eqtype = 'N';
end
eqn.type = eqtype;
opts.KSM.maxiter = 100;
opts.KSM.res_tol = 1e-6;
opts.KSM.comp_res = 1;
opts.KSM.space = space;
switch space
    case 'RK'
        opts.KSM.type_shifts = 'real';
end

if nargin < 2
    type = 'CARE';
elseif strcmp(type, 'LE') && strcmp(type, 'CARE')
    mess_err(opts, 'illegal_input', ...
             'type must be either ''LE'' or ''CARE''');
end
opts.KSM.type_eqn = type;
opts.KSM.CARE_shifts = 'Ritz';
opts.KSM.symmetric = true;
opts.KSM.trunc_tol = 1e-17;
opts.KSM.trunc_info = 1;
opts.KSM.explicit_proj = false;

opts.KSM.info = 1;

opts.norm = 'fro';

mess_fprintf(opts, '\n');
mess_fprintf(opts, 'Computing solution factors\n');

t_KSM = tic;
[out, eqn, opts, ~] = mess_KSM(eqn, opts, oper);
toc(t_KSM);
Computing solution factors
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   1 is:  6.093519e-14 (absolute)  6.036974e-01 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   2 is:  2.496149e-14 (absolute)  2.472986e-01 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   3 is:  9.403720e-15 (absolute)  9.316458e-02 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   4 is:  3.403934e-15 (absolute)  3.372348e-02 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   5 is:  1.144504e-15 (absolute)  1.133883e-02 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   6 is:  3.820469e-16 (absolute)  3.785017e-03 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   7 is:  1.153273e-16 (absolute)  1.142571e-03 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   8 is:  3.221288e-17 (absolute)  3.191396e-04 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration   9 is:  9.134052e-18 (absolute)  9.049293e-05 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration  10 is:  1.831665e-18 (absolute)  1.814668e-05 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration  11 is:  3.428914e-19 (absolute)  3.397096e-06 (normalized)  
Warning: dense Newton method stopped by maximum iteration number. Results may be
inaccurate! 
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/package/package.m',13)">package (line 13)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/_release/publish_demos.m',18)">publish_demos (line 18)</a>
↳ In <a href="matlab:opentoline('/matlab/R2020b/toolbox/matlab/codetools/publish.p',0)">publish (line 0)</a>
↳ In <a href="matlab:opentoline('',21)">evalmxdom (line 21)</a>
↳ In <a href="matlab:opentoline('',109)">instrumentAndRun (line 109)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/DEMOS/Rail/KSM_Rail.m',120)">KSM_Rail (line 120)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_KSM.m',322)">mess_KSM (line 322)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_solve_projected_eqn.m',347)">mess_solve_projected_eqn (line 347)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/helpers/mess_dense_nm.m',159)">mess_dense_nm (line 159)</a>
KSM residual at iteration  12 is:  6.388142e-20 (absolute)  6.328863e-07 (normalized)  
Elapsed time is 8.324799 seconds.

this shouldn't be computed for large problems

if k > 4
    return
end

X = out.Z * out.D * out.Z';
mess_fprintf(opts, '\n');  % \n at the beginning of a string breaks
% spellchecking

mess_fprintf(opts, ...
             ['The final internal and real normalized residual norms ', ...
              'after %d iterations are:\n'], out.niter);

if strcmp(eqn.type, 'T')
    if  strcmp(opts.KSM.type_eqn, 'CARE')
        mess_fprintf(opts, '\n');
        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');
        mess_fprintf(opts, ['CARE Frobenius norm: %10g (real) \t %10g ', ...
                            '(internal)\n'], ...
                     nrm, out.res(end));

        gnrm = norm(eqn.A_' * X * eqn.E_ + ...
                    eqn.E_' * X * eqn.A_ - ...
                    eqn.E_' * (X * eqn.B) * (eqn.B' * X) * eqn.E_ + ...
                    eqn.C' * eqn.C, 'fro') / ...
               norm(eqn.C' * eqn.C, 'fro');
        mess_fprintf(opts, ['gCARE Frobenius norm: %10g (real) \t %10g ', ...
                            '(internal)\n'], ...
                     gnrm, out.res(end));
    else
        mess_fprintf(opts, '\n');
        nrm = norm(eqn.A_' * X + X * eqn.A_ + eqn.C' * eqn.C, 'fro') / ...
              norm(eqn.C' * eqn.C, 'fro');
        mess_fprintf(opts, ['LE Frobenius norm: %10g (real) \t %10g ', ...
                            '(internal)\n'], ...
                     nrm, out.res(end));

        gnrm = norm(eqn.A_' * X * eqn.E_ + ...
                    eqn.E_' * X * eqn.A_ + ...
                    eqn.C' * eqn.C, 'fro') / ...
               norm(eqn.C' * eqn.C, 'fro');
        mess_fprintf(opts, ['gLE Frobenius norm: %10g (real) \t %10g ', ...
                            '(internal)\n'], ...
                     gnrm, out.res(end));
    end
else
    if  strcmp(opts.KSM.type_eqn, 'CARE')
        mess_fprintf(opts, '\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');
        mess_fprintf(opts, ['CARE Frobenius norm: %10g (real) \t %10g ', ...
                            '(internal)\n'], ...
                     nrm, out.res(end));

        gnrm = norm(eqn.A_ * X * eqn.E_' + ...
                    eqn.E_ * X * eqn.A_' - ...
                    eqn.E_ * (X * eqn.C') * (eqn.C * X) * eqn.E_' + ...
                    eqn.B * eqn.B', 'fro') / ...
               norm(eqn.B * eqn.B', 'fro');
        mess_fprintf(opts, ['gCARE Frobenius norm: %10g (real) \t %10g ', ...
                            '(internal)\n'], ...
                     gnrm, out.res(end));
    else
        mess_fprintf(opts, '\n');
        nrm = norm(eqn.A_ * X + X * eqn.A_' + eqn.B * eqn.B', 'fro') / ...
              norm(eqn.B * eqn.B', 'fro');
        mess_fprintf(opts, ['LE Frobenius norm: %10g (real) \t %10g ' ...
                            '(internal)\n'], ...
                     nrm, out.res(end));

        gnrm = norm(eqn.A_ * X * eqn.E_' + ...
                    eqn.E_ * X * eqn.A_' + ...
                    eqn.B * eqn.B', 'fro') / ...
               norm(eqn.B * eqn.B', 'fro');
        mess_fprintf(opts, ...
                     'gLE Frobenius norm: %10g (real) \t %10g (internal)\n', ...
                     gnrm, out.res(end));
    end
end

if istest
    if not(transformed)
        if (abs(gnrm - out.res(end)) / gnrm) > 1e-4
            mess_err(opts, 'failure', 'test failed');
        end
    else
        if (abs(nrm - out.res(end)) / nrm) > 1e-4
            mess_err(opts, 'failure', 'test failed');
        end
    end
end
The final internal and real normalized residual norms after 12 iterations are:  

CARE Frobenius norm:    4062.14 (real) 	 6.32886e-07 (internal)  
gCARE Frobenius norm: 6.32886e-07 (real) 	 6.32886e-07 (internal)