Contents

function Lyapunov_rail_LDL_ADI(k, shifts, implicit, istest)
% Computes the solution of the generalized Lyapunov equation via the
% low-rank ADI iteration in both ZZ^T [1,2] and LDL^T[3] formulation for the
% selective cooling of Steel profiles application described in [4,5,6].
%
% Usage:      Lyapunov_rail_LDL_ADI(k,shifts,istest)
%
% Inputs:
%
% k           refinement level of the model to use
%             (0 - 5, i.e. 109 - 79841 Dofs)
%             (optional, defaults to 2, i.e. 1357 Dofs)
%
% shifts      ADI shift selection strategy. Possible values:
%              'heur'        Penzl's heuristic shifts
%              'Wachspress'  Wachspress shifts, optimally solving the dense
%                            shift selection problem.
%              'projection'  projection shifts [7]
%             (optional, defaults to 'heur')
%
% implicit    projection shifts with implicit computation of projected A
%             matrix
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%
% References:
%
% [1] P. Benner, J.-R. Li, T. Penzl, Numerical solution of large-scale
%     Lyapunov equations, Riccati equations, and linear-quadratic optimal
%     control problems, Numer. Lin. Alg. Appl. 15 (9) (2008) 755–777.
%     https://doi.org/10.1002/nla.622
%
% [2] P. Kürschner, Efficient low-rank solution of large-scale matrix
%     equations, Dissertation, Otto-von-Guericke-Universität, Magdeburg,
%     Germany, shaker Verlag, ISBN 978-3-8440-4385-3 (Apr. 2016).
%     URL http://hdl.handle.net/11858/00-001M-0000-0029-CE18-2
%
% [3] N. Lang, H. Mena, J. Saak, On the benefits of the LDLT factorization
%     for large-scale differential matrix equation solvers, Linear Algebra
%     Appl. 480 (2015) 44–71.
%     https://doi.org/10.1016/j.laa.2015.04.006
%
% [4] J. Saak, Effiziente numerische Lösung eines
%     Optimalsteuerungsproblems für die Abkühlung von Stahlprofilen,
%     Diplomarbeit, Fachbereich 3/Mathematik und Informatik, Universität
%     Bremen, D-28334 Bremen (Sep. 2003).
%     https://doi.org/10.5281/zenodo.1187040
%
% [5] P. Benner, J. Saak, A semi-discretized heat transfer model for
%     optimal cooling of steel profiles, in: P. Benner, V. Mehrmann, D.
%     Sorensen (Eds.), Dimension Reduction of Large-Scale Systems, Vol. 45
%     of Lecture Notes in Computational Science and Engineering,
%     Springer-Verlag, Berlin/Heidelberg, Germany, 2005, pp. 353–356.
%     https://doi.org/10.1007/3-540-27909-1_19
%
% [6] J. Saak, Efficient numerical solution of large scale algebraic matrix
%     equations in PDE control and model order reduction, Dissertation,
%     Technische Universität Chemnitz, Chemnitz, Germany (Jul. 2009).
%     URL http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642
%
% [7] P. Benner, P. Kürschner, J. Saak, Self-generating and efficient shift
%     parameters in ADI methods for large Lyapunov and Sylvester equations,
%     Electron. Trans. Numer. Anal. 43 (2014) 142–162.
%     URL http://etna.mcs.kent.edu/volumes/2011-2020/vol43/abstract.php?vol=43&pages=142-162

%
% 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)
%
narginchk(0, 4);
if nargin < 1
    k = 2;
end
if nargin < 2
    shifts = 'wachspress';
end
if nargin < 3
    implicit = 0;
end
if nargin < 4
    istest = false;
end

set operation

opts = struct();
[oper, opts] = operatormanager(opts, 'default');
% Problem data
eqn = mess_get_linear_rail(k);
% ADI tolerances and maximum iteration number
opts.adi.maxiter      = 100;
opts.adi.res_tol      = 1e-12;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info         = 1;

eqn.type = 'T';

Heuristic shift parameters via basic Arnoldi

n = oper.size(eqn, opts);
opts.shifts.num_Ritz    = 50;
opts.shifts.num_hRitz   = 25;
opts.shifts.num_desired = 6;

opts.shifts.b0 = ones(n, 1);
switch lower(shifts)

    case 'heur'
        opts.shifts.method       = 'heur';

    case 'wachspress'
        opts.shifts.method       = 'wachspress';
        opts.shifts.wachspress   = 'T';

    case 'projection'
        opts.shifts.method       = 'projection';
        opts.shifts.implicitVtAV = implicit;
end
opts.norm = 'fro';
mess_fprintf(opts, '########################\n');
mess_fprintf(opts, '# ADI with ZZ^T:        \n');
mess_fprintf(opts, '########################\n');
t_mess_lradi = tic;
out = mess_lradi(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lradi);
mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n', t_elapsed1);

if istest
    if min(out.res) >= opts.adi.res_tol
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(out.res, 'LineWidth', 3);
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
[mZ, nZ] = size(out.Z);
mess_fprintf(opts, 'size out.Z: %d x %d\n\n', mZ, nZ);
########################
# ADI with ZZ^T:        
########################
ADI step:    1 normalized residual: 8.784598e-01 relative change in Z: 1.000000e+00  
ADI step:    2 normalized residual: 7.834838e-01 relative change in Z: 6.486596e-01  
ADI step:    3 normalized residual: 6.992446e-01 relative change in Z: 5.096197e-01  
ADI step:    4 normalized residual: 6.191398e-01 relative change in Z: 4.349779e-01  
ADI step:    5 normalized residual: 5.404000e-01 relative change in Z: 3.866261e-01  
ADI step:    6 normalized residual: 4.622520e-01 relative change in Z: 3.506752e-01  
ADI step:    7 normalized residual: 3.854888e-01 relative change in Z: 3.212854e-01  
ADI step:    8 normalized residual: 3.124516e-01 relative change in Z: 2.954727e-01  
ADI step:    9 normalized residual: 2.465196e-01 relative change in Z: 2.714156e-01  
ADI step:   10 normalized residual: 1.907576e-01 relative change in Z: 2.484019e-01  
ADI step:   11 normalized residual: 1.463406e-01 relative change in Z: 2.270182e-01  
ADI step:   12 normalized residual: 1.119302e-01 relative change in Z: 2.086859e-01  
ADI step:   13 normalized residual: 8.467511e-02 relative change in Z: 1.943956e-01  
ADI step:   14 normalized residual: 6.211017e-02 relative change in Z: 1.834507e-01  
ADI step:   15 normalized residual: 4.326477e-02 relative change in Z: 1.737083e-01  
ADI step:   16 normalized residual: 2.826568e-02 relative change in Z: 1.631056e-01  
ADI step:   17 normalized residual: 1.733668e-02 relative change in Z: 1.509104e-01  
ADI step:   18 normalized residual: 1.020709e-02 relative change in Z: 1.380123e-01  
ADI step:   19 normalized residual: 6.049633e-03 relative change in Z: 1.259393e-01  
ADI step:   20 normalized residual: 3.740716e-03 relative change in Z: 1.148802e-01  
ADI step:   21 normalized residual: 2.361170e-03 relative change in Z: 1.030612e-01  
ADI step:   22 normalized residual: 1.449919e-03 relative change in Z: 8.867504e-02  
ADI step:   23 normalized residual: 8.556644e-04 relative change in Z: 7.196930e-02  
ADI step:   24 normalized residual: 5.234750e-04 relative change in Z: 5.568756e-02  
ADI step:   25 normalized residual: 3.689986e-04 relative change in Z: 4.388619e-02  
ADI step:   26 normalized residual: 2.940349e-04 relative change in Z: 3.889762e-02  
ADI step:   27 normalized residual: 2.483663e-04 relative change in Z: 3.874220e-02  
ADI step:   28 normalized residual: 2.153700e-04 relative change in Z: 4.035949e-02  
ADI step:   29 normalized residual: 1.868929e-04 relative change in Z: 4.288721e-02  
ADI step:   30 normalized residual: 1.578141e-04 relative change in Z: 4.616489e-02  
ADI step:   31 normalized residual: 1.263773e-04 relative change in Z: 4.944642e-02  
ADI step:   32 normalized residual: 9.371539e-05 relative change in Z: 5.174894e-02  
ADI step:   33 normalized residual: 6.255849e-05 relative change in Z: 5.206152e-02  
ADI step:   34 normalized residual: 3.615248e-05 relative change in Z: 4.946079e-02  
ADI step:   35 normalized residual: 1.717309e-05 relative change in Z: 4.356664e-02  
ADI step:   36 normalized residual: 6.311398e-06 relative change in Z: 3.503761e-02  
ADI step:   37 normalized residual: 1.701167e-06 relative change in Z: 2.550287e-02  
ADI step:   38 normalized residual: 3.284777e-07 relative change in Z: 1.665063e-02  
ADI step:   39 normalized residual: 4.358497e-08 relative change in Z: 9.364006e-03  
ADI step:   40 normalized residual: 3.157861e-09 relative change in Z: 4.097720e-03  
ADI step:   41 normalized residual: 7.886436e-11 relative change in Z: 1.197046e-03  
ADI step:   42 normalized residual: 5.875199e-13 relative change in Z: 1.901405e-04  
mess_lradi took   0.37 seconds   
size out.Z: 1357 x 252  
  

Set LDL fields

opts.LDL_T = true;
eqn.T = diag([4, 4, 9, 9, 16, 16]);
eqn.W = eqn.C' * diag([4, 4, 9, 9, 16, 16].^(-0.5));
mess_fprintf(opts, '########################\n');
mess_fprintf(opts, '# ADI with LDL^T:       \n');
mess_fprintf(opts, '########################\n');
t_mess_lradi = tic;
out1 = mess_lradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lradi);
mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n', t_elapsed2);

if istest
    if min(out.res) >= opts.adi.res_tol
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(out1.res, 'LineWidth', 3);
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
[mZ, nZ] = size(out1.Z);
mess_fprintf(opts, 'size out1.Z: %d x %d\n\n', mZ, nZ);
########################
# ADI with LDL^T:       
########################
ADI step:    1 normalized residual: 8.784598e-01 relative change in Z: 1.000000e+00  
ADI step:    2 normalized residual: 7.834838e-01 relative change in Z: 6.607722e-01  
ADI step:    3 normalized residual: 6.992446e-01 relative change in Z: 5.292964e-01  
ADI step:    4 normalized residual: 6.191398e-01 relative change in Z: 4.611570e-01  
ADI step:    5 normalized residual: 5.404000e-01 relative change in Z: 4.187750e-01  
ADI step:    6 normalized residual: 4.622520e-01 relative change in Z: 3.877834e-01  
ADI step:    7 normalized residual: 3.854888e-01 relative change in Z: 3.614172e-01  
ADI step:    8 normalized residual: 3.124516e-01 relative change in Z: 3.358785e-01  
ADI step:    9 normalized residual: 2.465196e-01 relative change in Z: 3.090984e-01  
ADI step:   10 normalized residual: 1.907576e-01 relative change in Z: 2.808051e-01  
ADI step:   11 normalized residual: 1.463406e-01 relative change in Z: 2.526762e-01  
ADI step:   12 normalized residual: 1.119302e-01 relative change in Z: 2.277203e-01  
ADI step:   13 normalized residual: 8.467511e-02 relative change in Z: 2.083123e-01  
ADI step:   14 normalized residual: 6.211017e-02 relative change in Z: 1.939634e-01  
ADI step:   15 normalized residual: 4.326477e-02 relative change in Z: 1.815841e-01  
ADI step:   16 normalized residual: 2.826568e-02 relative change in Z: 1.679855e-01  
ADI step:   17 normalized residual: 1.733668e-02 relative change in Z: 1.518485e-01  
ADI step:   18 normalized residual: 1.020709e-02 relative change in Z: 1.342506e-01  
ADI step:   19 normalized residual: 6.049633e-03 relative change in Z: 1.176705e-01  
ADI step:   20 normalized residual: 3.740716e-03 relative change in Z: 1.036853e-01  
ADI step:   21 normalized residual: 2.361170e-03 relative change in Z: 9.150195e-02  
ADI step:   22 normalized residual: 1.449919e-03 relative change in Z: 7.935436e-02  
ADI step:   23 normalized residual: 8.556644e-04 relative change in Z: 6.700389e-02  
ADI step:   24 normalized residual: 5.234750e-04 relative change in Z: 5.656666e-02  
ADI step:   25 normalized residual: 3.689986e-04 relative change in Z: 5.095069e-02  
ADI step:   26 normalized residual: 2.940349e-04 relative change in Z: 5.092780e-02  
ADI step:   27 normalized residual: 2.483663e-04 relative change in Z: 5.452279e-02  
ADI step:   28 normalized residual: 2.153700e-04 relative change in Z: 5.966108e-02  
ADI step:   29 normalized residual: 1.868929e-04 relative change in Z: 6.548959e-02  
ADI step:   30 normalized residual: 1.578141e-04 relative change in Z: 7.148644e-02  
ADI step:   31 normalized residual: 1.263773e-04 relative change in Z: 7.680915e-02  
ADI step:   32 normalized residual: 9.371539e-05 relative change in Z: 8.032059e-02  
ADI step:   33 normalized residual: 6.255849e-05 relative change in Z: 8.063136e-02  
ADI step:   34 normalized residual: 3.615248e-05 relative change in Z: 7.639703e-02  
ADI step:   35 normalized residual: 1.717309e-05 relative change in Z: 6.709429e-02  
ADI step:   36 normalized residual: 6.311398e-06 relative change in Z: 5.379174e-02  
ADI step:   37 normalized residual: 1.701167e-06 relative change in Z: 3.903856e-02  
ADI step:   38 normalized residual: 3.284777e-07 relative change in Z: 2.544198e-02  
ADI step:   39 normalized residual: 4.358497e-08 relative change in Z: 1.430879e-02  
ADI step:   40 normalized residual: 3.157861e-09 relative change in Z: 6.268140e-03  
ADI step:   41 normalized residual: 7.886436e-11 relative change in Z: 1.832423e-03  
ADI step:   42 normalized residual: 5.875199e-13 relative change in Z: 2.910450e-04  
mess_lradi took   0.25 seconds   
size out1.Z: 1357 x 252  
  

Difference of Lyapunov solutions

if k < 3
    % This is mainly for consistency checking on our continuous
    % integration tests.
    % NEVER FORM SUCH DYADIC PRODUCTS IN PRODUCTION CODE!!!
    err = norm(out.Z * out.Z' - out1.Z * out1.D * out1.Z') / ...
          norm(out.Z * out.Z');
    mess_fprintf(opts, ...
                 ['Relative difference between solution with and without ' ...
                  'LDL^T: \t %g\n'], ...
                 err);
    if err > 1e-11
        if implicit
            shifts = [shifts '(implicit)'];
        end
        mess_err(opts, 'TEST:accuracy', ...
                 ['unexpectedly inaccurate result relative difference', ...
                  ' %e > 1e-12 in case %s'], ...
                 err, shifts);
    end
end
Relative difference between solution with and without LDL^T: 	 5.27035e-16