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
%             (1-4, i.e. 1357-79841Dofs)
%             (optinal, defaults to 1)
%
% 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 largescale 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).
%
% [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 Lect. Notes Comput. Sci. Eng., 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 program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
%               2009-2020
%
narginchk(0,4);
if nargin<1, k=1; end
if nargin<2, shifts='wachspress'; end
if nargin<3, implicit=0; end
if nargin<4, istest=0; end

set operation

oper = operatormanager('default');
% Problem data
eqn = getrail(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.shifts.p=mess_para(eqn,opts,oper);
opts.norm = 'fro';
fprintf('########################\n');
fprintf('# ADI with ZZ^T:        \n');
fprintf('########################\n');
tic;
out = mess_lradi(eqn, opts, oper);
toc;

if istest
    if min(out.res)>=opts.adi.res_tol
       error('MESS:TEST:accuracy','unexpectedly inaccurate result');
   end
else
    figure(1);
    semilogy(out.res);
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size out.Z:');
disp(size(out.Z));
########################
# ADI with ZZ^T:        
########################
ADI step:    1 normalized residual: 9.295015e-01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 8.668178e-01 relative change in Z: 6.761705e-01
ADI step:    3 normalized residual: 8.051179e-01 relative change in Z: 5.486908e-01
ADI step:    4 normalized residual: 7.409268e-01 relative change in Z: 4.816916e-01
ADI step:    5 normalized residual: 6.727096e-01 relative change in Z: 4.393321e-01
ADI step:    6 normalized residual: 6.002287e-01 relative change in Z: 4.078640e-01
ADI step:    7 normalized residual: 5.241929e-01 relative change in Z: 3.813269e-01
ADI step:    8 normalized residual: 4.461882e-01 relative change in Z: 3.570955e-01
ADI step:    9 normalized residual: 3.688459e-01 relative change in Z: 3.338816e-01
ADI step:   10 normalized residual: 2.958357e-01 relative change in Z: 3.108154e-01
ADI step:   11 normalized residual: 2.311456e-01 relative change in Z: 2.874086e-01
ADI step:   12 normalized residual: 1.776061e-01 relative change in Z: 2.639929e-01
ADI step:   13 normalized residual: 1.355442e-01 relative change in Z: 2.419540e-01
ADI step:   14 normalized residual: 1.028509e-01 relative change in Z: 2.230617e-01
ADI step:   15 normalized residual: 7.660017e-02 relative change in Z: 2.079796e-01
ADI step:   16 normalized residual: 5.480322e-02 relative change in Z: 1.954637e-01
ADI step:   17 normalized residual: 3.697245e-02 relative change in Z: 1.833006e-01
ADI step:   18 normalized residual: 2.336999e-02 relative change in Z: 1.699510e-01
ADI step:   19 normalized residual: 1.399895e-02 relative change in Z: 1.554995e-01
ADI step:   20 normalized residual: 8.249446e-03 relative change in Z: 1.414048e-01
ADI step:   21 normalized residual: 5.025404e-03 relative change in Z: 1.287738e-01
ADI step:   22 normalized residual: 3.181876e-03 relative change in Z: 1.165575e-01
ADI step:   23 normalized residual: 2.001619e-03 relative change in Z: 1.024183e-01
ADI step:   24 normalized residual: 1.205763e-03 relative change in Z: 8.534499e-02
ADI step:   25 normalized residual: 7.172900e-04 relative change in Z: 6.705462e-02
ADI step:   26 normalized residual: 4.694384e-04 relative change in Z: 5.153971e-02
ADI step:   27 normalized residual: 3.527571e-04 relative change in Z: 4.274872e-02
ADI step:   28 normalized residual: 2.854815e-04 relative change in Z: 4.066369e-02
ADI step:   29 normalized residual: 2.370529e-04 relative change in Z: 4.160086e-02
ADI step:   30 normalized residual: 1.962232e-04 relative change in Z: 4.324807e-02
ADI step:   31 normalized residual: 1.568807e-04 relative change in Z: 4.518725e-02
ADI step:   32 normalized residual: 1.170897e-04 relative change in Z: 4.679393e-02
ADI step:   33 normalized residual: 7.879545e-05 relative change in Z: 4.699500e-02
ADI step:   34 normalized residual: 4.589439e-05 relative change in Z: 4.478056e-02
ADI step:   35 normalized residual: 2.187953e-05 relative change in Z: 3.950392e-02
ADI step:   36 normalized residual: 7.902182e-06 relative change in Z: 3.136163e-02
ADI step:   37 normalized residual: 1.961263e-06 relative change in Z: 2.177405e-02
ADI step:   38 normalized residual: 3.071648e-07 relative change in Z: 1.294172e-02
ADI step:   39 normalized residual: 3.094205e-08 relative change in Z: 6.510069e-03
ADI step:   40 normalized residual: 1.953366e-09 relative change in Z: 2.662204e-03
ADI step:   41 normalized residual: 4.715330e-11 relative change in Z: 7.698935e-04
ADI step:   42 normalized residual: 4.709828e-13 relative change in Z: 1.230127e-04
Elapsed time is 0.076513 seconds.
size out.Z:
        1357         252

Set LDL fields

opts.LDL_T = 1;
eqn.S = diag([4,4,9,9,16,16]);
eqn.G = eqn.C' * diag([4,4,9,9,16,16].^(-0.5));
fprintf('########################\n');
fprintf('# ADI with LDL^T:       \n');
fprintf('########################\n');
opts.shifts.p=mess_para(eqn,opts,oper);
tic;
out1 = mess_lradi(eqn, opts, oper);
toc;

if istest
    if min(out.res)>=opts.adi.res_tol
       error('MESS:TEST:accuracy','unexpectedly inaccurate result');
   end
else
    figure(2);
    semilogy(out1.res);
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size out1.Z:');
disp(size(out1.Z));
########################
# ADI with LDL^T:       
########################
ADI step:    1 normalized residual: 9.295015e-01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 8.668178e-01 relative change in Z: 6.845616e-01
ADI step:    3 normalized residual: 8.051179e-01 relative change in Z: 5.622064e-01
ADI step:    4 normalized residual: 7.409268e-01 relative change in Z: 4.996237e-01
ADI step:    5 normalized residual: 6.727096e-01 relative change in Z: 4.615544e-01
ADI step:    6 normalized residual: 6.002287e-01 relative change in Z: 4.342783e-01
ADI step:    7 normalized residual: 5.241929e-01 relative change in Z: 4.115193e-01
ADI step:    8 normalized residual: 4.461882e-01 relative change in Z: 3.900572e-01
ADI step:    9 normalized residual: 3.688459e-01 relative change in Z: 3.679366e-01
ADI step:   10 normalized residual: 2.958357e-01 relative change in Z: 3.438585e-01
ADI step:   11 normalized residual: 2.311456e-01 relative change in Z: 3.173376e-01
ADI step:   12 normalized residual: 1.776061e-01 relative change in Z: 2.891777e-01
ADI step:   13 normalized residual: 1.355442e-01 relative change in Z: 2.616553e-01
ADI step:   14 normalized residual: 1.028509e-01 relative change in Z: 2.376189e-01
ADI step:   15 normalized residual: 7.660017e-02 relative change in Z: 2.184736e-01
ADI step:   16 normalized residual: 5.480322e-02 relative change in Z: 2.028716e-01
ADI step:   17 normalized residual: 3.697245e-02 relative change in Z: 1.877678e-01
ADI step:   18 normalized residual: 2.336999e-02 relative change in Z: 1.707601e-01
ADI step:   19 normalized residual: 1.399895e-02 relative change in Z: 1.516709e-01
ADI step:   20 normalized residual: 8.249446e-03 relative change in Z: 1.325671e-01
ADI step:   21 normalized residual: 5.025404e-03 relative change in Z: 1.159418e-01
ADI step:   22 normalized residual: 3.181876e-03 relative change in Z: 1.021666e-01
ADI step:   23 normalized residual: 2.001619e-03 relative change in Z: 8.940539e-02
ADI step:   24 normalized residual: 1.205763e-03 relative change in Z: 7.627503e-02
ADI step:   25 normalized residual: 7.172900e-04 relative change in Z: 6.385690e-02
ADI step:   26 normalized residual: 4.694384e-04 relative change in Z: 5.519119e-02
ADI step:   27 normalized residual: 3.527571e-04 relative change in Z: 5.247123e-02
ADI step:   28 normalized residual: 2.854815e-04 relative change in Z: 5.461069e-02
ADI step:   29 normalized residual: 2.370529e-04 relative change in Z: 5.889037e-02
ADI step:   30 normalized residual: 1.962232e-04 relative change in Z: 6.357819e-02
ADI step:   31 normalized residual: 1.568807e-04 relative change in Z: 6.783239e-02
ADI step:   32 normalized residual: 1.170897e-04 relative change in Z: 7.069701e-02
ADI step:   33 normalized residual: 7.879545e-05 relative change in Z: 7.095376e-02
ADI step:   34 normalized residual: 4.589439e-05 relative change in Z: 6.735728e-02
ADI step:   35 normalized residual: 2.187953e-05 relative change in Z: 5.906471e-02
ADI step:   36 normalized residual: 7.902182e-06 relative change in Z: 4.646421e-02
ADI step:   37 normalized residual: 1.961263e-06 relative change in Z: 3.180240e-02
ADI step:   38 normalized residual: 3.071648e-07 relative change in Z: 1.851713e-02
ADI step:   39 normalized residual: 3.094205e-08 relative change in Z: 9.120345e-03
ADI step:   40 normalized residual: 1.953366e-09 relative change in Z: 3.697835e-03
ADI step:   41 normalized residual: 4.715330e-11 relative change in Z: 1.072806e-03
ADI step:   42 normalized residual: 4.709828e-13 relative change in Z: 1.717290e-04
Elapsed time is 0.074660 seconds.
size out1.Z:
        1357         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');
    fprintf(['Relative difference between solution with and without ' ...
             'LDL^T: \t %g\n'], err);
    if err>1e-12
        if implicit
            shifts = [shifts '(implicit)'];
        end
        error('MESS: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: 	 4.11452e-16