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)