function BT_sym_TripleChain(variant, istest)
%
% Computes a reduced order model (ROM) for the triple chain example of
% Truhar and Veselic [1] via Balanced truncation, e.g. [2], exploiting
% symmetry of the control system, i.e. equality of the system Gramians.
% In comparison to BT_TripleChain this demo exploits, that the first order
% form of the system has symmetric E and A and B=C', such that the two
% Lyapunov equations coincide.
%
% Usage:   BT_sym_TripleChain(version)
%
% Input:
%
% variant  Decides the Balanced Truncation version to use.
%          Possible values:
%          'FO' for reduction of the first order form to first order form
%          'VV' velocity-velocity balancing of the second order form to
%               second order form.
%          'PP' position-position balancing of the second order form to
%               second order form.
%          'PV' position-velocity balancing of the second order form to
%               second order form.
%          'VP' velocity-position balancing of the second order form to
%               second order form.%
%
% 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] N. Truhar and K. Veselic, An efficient method for estimating the
%     optimal dampers’ viscosity for linear vibrating systems using
%     Lyapunov equation, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 18–39.
%     https://doi.org/10.1137/070683052
%
% [2] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Vol.
%     6 of Adv. Des. Control, SIAM Publications, Philadelphia, PA, 2005.
%     https://doi.org/10.1137/1.9780898718713

%
% 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, 2);

if nargin == 0
    variant = 'FO';
end
if nargin < 2
    istest = false;
end
format long e;
% set operation
opts = struct();
[oper, opts] = operatormanager(opts, 'so_2');
% Problem data

n1      = 500;
alpha   = 0.002;
Beta    = alpha;
v       = 5;

[eqn.M_, eqn.E_, eqn.K_] = ...
    triplechain_MSD(n1, alpha, Beta, v);

s  = size(eqn.K_, 1);

B = ones(s, 1);
O = zeros(s, 1);
Cp = ones(1, s);
Cv = O';
eqn.B = [B; O];
eqn.C = [Cp, Cv];

eqn.haveE = true;

ADI tolerances and maximum iteration number

opts.adi.maxiter = 300;
opts.adi.res_tol = 1e-10;
opts.norm        = 'fro';

% opts.adi.rel_diff_tol    = 1e-16;
opts.adi.rel_diff_tol    = 0;
opts.adi.info            = 1;
% opts.adi.accumulateK = false;
% opts.adi.accumulateDeltaK = false;
opts.adi.compute_sol_fac = true;

eqn.type = 'N';

Heuristic shift parameters via projection

opts.shifts.num_desired  = 5;

opts.shifts.info         = 0;
opts.shifts.method       = 'projection';

Compute Gramian Factor (one is enough, since the Gramians are equal)

t_mess_lradi = tic;
outB = mess_lradi(eqn, opts, oper);
t_elapsed = toc(t_mess_lradi);
mess_fprintf(opts, 'mess_lradi took %6.2f seconds \n', t_elapsed);

if istest
    if min(outB.res) >= 1e-1
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outB.res, 'LineWidth', 3);
    title('0 = A X ^T + E X A^T - BB^T');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end

[mZ, nZ] = size(outB.Z);
mess_fprintf(opts, 'size outB.Z: %d x %d\n\n', mZ, nZ);
ADI step:    1 normalized residual: 1.270108e+00   
ADI step:    2 normalized residual: 3.427418e+00   
ADI step:    4 normalized residual: 4.477922e+02   
ADI step:    6 normalized residual: 3.943500e+02   
ADI step:    8 normalized residual: 6.085237e+02   
ADI step:    9 normalized residual: 7.293151e+03   
ADI step:   11 normalized residual: 7.332737e+03   
ADI step:   13 normalized residual: 7.483298e+03   
ADI step:   14 normalized residual: 5.877023e+02   
ADI step:   16 normalized residual: 7.272517e+02   
ADI step:   18 normalized residual: 8.594421e+02   
ADI step:   20 normalized residual: 2.043734e+02   
ADI step:   21 normalized residual: 8.008062e+02   
ADI step:   23 normalized residual: 7.946889e+02   
ADI step:   25 normalized residual: 6.376582e+01   
ADI step:   27 normalized residual: 6.430936e+01   
ADI step:   28 normalized residual: 9.394251e+01   
ADI step:   29 normalized residual: 7.753989e+01   
ADI step:   31 normalized residual: 8.002104e+00   
ADI step:   33 normalized residual: 7.641858e+00   
ADI step:   35 normalized residual: 6.705939e+00   
ADI step:   36 normalized residual: 5.442093e+00   
ADI step:   38 normalized residual: 5.539283e+00   
ADI step:   40 normalized residual: 4.767961e+00   
ADI step:   42 normalized residual: 4.926510e+00   
ADI step:   43 normalized residual: 5.954219e+00   
ADI step:   44 normalized residual: 4.564030e+00   
ADI step:   46 normalized residual: 5.300509e+00   
ADI step:   48 normalized residual: 2.859722e+00   
ADI step:   50 normalized residual: 1.514706e+00   
ADI step:   51 normalized residual: 3.250111e+00   
ADI step:   53 normalized residual: 2.174333e+00   
ADI step:   55 normalized residual: 1.813608e+00   
ADI step:   57 normalized residual: 1.797216e+00   
ADI step:   58 normalized residual: 1.571941e+00   
ADI step:   60 normalized residual: 1.090202e+00   
ADI step:   62 normalized residual: 1.067706e+00   
ADI step:   63 normalized residual: 1.210877e+00   
ADI step:   65 normalized residual: 1.156980e+00   
ADI step:   66 normalized residual: 7.188757e-01   
ADI step:   68 normalized residual: 1.172664e+00   
ADI step:   70 normalized residual: 1.129371e+00   
ADI step:   71 normalized residual: 8.961106e-01   
ADI step:   73 normalized residual: 8.713581e-01   
ADI step:   75 normalized residual: 7.378372e-01   
ADI step:   77 normalized residual: 5.279300e-01   
ADI step:   78 normalized residual: 8.063868e-01   
ADI step:   80 normalized residual: 6.666618e-01   
ADI step:   81 normalized residual: 7.099150e-01   
ADI step:   83 normalized residual: 6.526703e-01   
ADI step:   85 normalized residual: 4.615184e-01   
ADI step:   87 normalized residual: 4.450107e-01   
ADI step:   88 normalized residual: 5.715073e-01   
ADI step:   90 normalized residual: 5.039067e-01   
ADI step:   92 normalized residual: 3.200372e-01   
ADI step:   93 normalized residual: 4.142961e-01   
ADI step:   95 normalized residual: 4.217176e-01   
ADI step:   97 normalized residual: 3.332969e-01   
ADI step:   98 normalized residual: 3.777805e-01   
ADI step:  100 normalized residual: 3.617469e-01   
ADI step:  101 normalized residual: 3.912357e-01   
ADI step:  103 normalized residual: 3.370402e-01   
ADI step:  105 normalized residual: 3.317560e-01   
ADI step:  107 normalized residual: 3.355225e-01   
ADI step:  108 normalized residual: 3.173618e-01   
ADI step:  110 normalized residual: 2.591786e-01   
ADI step:  111 normalized residual: 2.713283e-01   
ADI step:  113 normalized residual: 2.680890e-01   
ADI step:  115 normalized residual: 2.582752e-01   
ADI step:  117 normalized residual: 2.575143e-01   
ADI step:  118 normalized residual: 1.824167e-01   
ADI step:  120 normalized residual: 2.048510e-01   
ADI step:  121 normalized residual: 2.396396e-01   
ADI step:  123 normalized residual: 2.420217e-01   
ADI step:  124 normalized residual: 2.250442e-01   
ADI step:  126 normalized residual: 2.267417e-01   
ADI step:  128 normalized residual: 2.233053e-01   
ADI step:  130 normalized residual: 2.276675e-01   
ADI step:  132 normalized residual: 2.321377e-01   
ADI step:  133 normalized residual: 2.234471e-01   
ADI step:  135 normalized residual: 1.473972e-01   
ADI step:  137 normalized residual: 1.421866e-01   
ADI step:  138 normalized residual: 1.709875e-01   
ADI step:  139 normalized residual: 1.586470e-01   
ADI step:  141 normalized residual: 9.900877e-02   
ADI step:  143 normalized residual: 9.614724e-02   
ADI step:  144 normalized residual: 1.857306e-01   
ADI step:  146 normalized residual: 1.861647e-01   
ADI step:  148 normalized residual: 1.863501e-01   
ADI step:  149 normalized residual: 1.995936e-01   
ADI step:  151 normalized residual: 1.992495e-01   
ADI step:  153 normalized residual: 1.787107e-01   
ADI step:  155 normalized residual: 1.668872e-01   
ADI step:  157 normalized residual: 1.534941e-01   
ADI step:  158 normalized residual: 1.706631e-01   
ADI step:  160 normalized residual: 1.645401e-01   
ADI step:  162 normalized residual: 1.596682e-01   
ADI step:  163 normalized residual: 1.581303e-01   
ADI step:  165 normalized residual: 1.586407e-01   
ADI step:  167 normalized residual: 1.645811e-01   
ADI step:  168 normalized residual: 1.726095e-01   
ADI step:  170 normalized residual: 1.726114e-01   
ADI step:  172 normalized residual: 1.692694e-01   
ADI step:  173 normalized residual: 1.691214e-01   
ADI step:  175 normalized residual: 1.703841e-01   
ADI step:  177 normalized residual: 1.661153e-01   
ADI step:  178 normalized residual: 1.627408e-01   
ADI step:  180 normalized residual: 1.603131e-01   
ADI step:  182 normalized residual: 1.609220e-01   
ADI step:  183 normalized residual: 1.654816e-01   
ADI step:  185 normalized residual: 1.652809e-01   
ADI step:  187 normalized residual: 1.635957e-01   
ADI step:  188 normalized residual: 1.638068e-01   
ADI step:  190 normalized residual: 1.645263e-01   
ADI step:  192 normalized residual: 1.648147e-01   
ADI step:  193 normalized residual: 1.640379e-01   
ADI step:  195 normalized residual: 1.610816e-01   
ADI step:  197 normalized residual: 1.279600e-01   
ADI step:  198 normalized residual: 1.263336e-01   
ADI step:  199 normalized residual: 1.137194e-01   
ADI step:  201 normalized residual: 1.132135e-01   
ADI step:  203 normalized residual: 1.045018e-01   
ADI step:  204 normalized residual: 1.395207e-01   
ADI step:  206 normalized residual: 1.353857e-01   
ADI step:  208 normalized residual: 1.356013e-01   
ADI step:  209 normalized residual: 7.226735e-02   
ADI step:  211 normalized residual: 7.212180e-02   
ADI step:  213 normalized residual: 7.147190e-02   
ADI step:  214 normalized residual: 1.019081e-01   
ADI step:  216 normalized residual: 1.019346e-01   
ADI step:  218 normalized residual: 1.015503e-01   
ADI step:  219 normalized residual: 9.629021e-02   
ADI step:  221 normalized residual: 9.773470e-02   
ADI step:  223 normalized residual: 9.797324e-02   
ADI step:  224 normalized residual: 7.189089e-02   
ADI step:  226 normalized residual: 7.215250e-02   
ADI step:  228 normalized residual: 7.139122e-02   
ADI step:  229 normalized residual: 7.682273e-02   
ADI step:  231 normalized residual: 7.550024e-02   
ADI step:  233 normalized residual: 7.565256e-02   
ADI step:  234 normalized residual: 8.643170e-02   
ADI step:  236 normalized residual: 8.579344e-02   
ADI step:  238 normalized residual: 8.526905e-02   
ADI step:  239 normalized residual: 7.558314e-02   
ADI step:  241 normalized residual: 7.507395e-02   
ADI step:  243 normalized residual: 7.458020e-02   
ADI step:  244 normalized residual: 9.377813e-02   
ADI step:  246 normalized residual: 9.378112e-02   
ADI step:  248 normalized residual: 9.372578e-02   
ADI step:  249 normalized residual: 9.792557e-02   
ADI step:  251 normalized residual: 9.748247e-02   
ADI step:  253 normalized residual: 9.725190e-02   
ADI step:  254 normalized residual: 8.359994e-02   
ADI step:  256 normalized residual: 8.365475e-02   
ADI step:  258 normalized residual: 8.316054e-02   
ADI step:  259 normalized residual: 7.309216e-02   
ADI step:  261 normalized residual: 7.298501e-02   
ADI step:  263 normalized residual: 7.186703e-02   
ADI step:  265 normalized residual: 7.128906e-02   
ADI step:  266 normalized residual: 6.935223e-02   
ADI step:  268 normalized residual: 6.973855e-02   
ADI step:  270 normalized residual: 6.975180e-02   
ADI step:  272 normalized residual: 6.951310e-02   
ADI step:  273 normalized residual: 7.076928e-02   
ADI step:  274 normalized residual: 7.909230e-02   
ADI step:  276 normalized residual: 7.930569e-02   
ADI step:  278 normalized residual: 7.928774e-02   
ADI step:  279 normalized residual: 4.677406e-02   
ADI step:  281 normalized residual: 4.644850e-02   
ADI step:  283 normalized residual: 4.551248e-02   
ADI step:  284 normalized residual: 6.268944e-02   
ADI step:  286 normalized residual: 6.265200e-02   
ADI step:  288 normalized residual: 6.149662e-02   
ADI step:  289 normalized residual: 6.020292e-02   
ADI step:  291 normalized residual: 5.981246e-02   
ADI step:  293 normalized residual: 5.946440e-02   
ADI step:  294 normalized residual: 4.795946e-02   
ADI step:  296 normalized residual: 4.789738e-02   
ADI step:  298 normalized residual: 4.765180e-02   
ADI step:  299 normalized residual: 5.558819e-02   
ADI step:  301 normalized residual: 5.557317e-02   
Warning: LR-ADI reached 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/TripleChain/BT_sym_TripleChain.m',107)">BT_sym_TripleChain (line 107)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mat-eqn-solvers/mess_lradi.m',865)">mess_lradi (line 865)</a>
mess_lradi took   0.58 seconds   
size outB.Z: 3002 x 301  
  
switch upper(variant)
    case 'FO'

Compute first order ROM

        opts.srm.max_ord = 150;
        opts.srm.tol = eps;
        opts.srm.info = 1;

        [TL, TR] = mess_square_root_method(eqn, opts, oper, outB.Z, outB.Z);

        ROM.E = eye(size(TL, 2));
        ROM.A = TL' * oper.mul_A(eqn, opts, 'N', TR, 'N');
        ROM.B = TL' * eqn.B;
        ROM.C = eqn.C * TR;
        ROM.D = [];
reduced system order: 150  (max possible/allowed: 301/150)  
  
    case 'VV'
        U = outB.Z(1:s, :);
        V = outB.Z(1:s, :);
    case 'PP'
        U = outB.Z(s + 1:end, :);
        V = outB.Z(s + 1:end, :);
    case 'PV'
        U = outB.Z(s + 1:end, :);
        V = outB.Z(1:s, :);
    case 'VP'
        U = outB.Z(1:s, :);
        V = outB.Z(s + 1:end, :);
end
if not(strcmp(variant, 'FO'))
    max_ord = 75;
    tol = eps;
    inform = 1;

    [TL, TR] = square_root_method_SO(eqn.M_, max_ord, tol, inform, U, V);

    ROM.M = eye(size(TL, 2));
    ROM.E = TL' * (eqn.E_ * TR);
    ROM.K = TL' * (eqn.K_ * TR);
    ROM.B = TL' * B;
    ROM.Cv = Cv * TR;
    ROM.Cp = Cp * TR;
end

plot results

opts.tf_plot.fmin = 1e-4;
opts.tf_plot.fmax = 1e0;
opts.tf_plot.nsample = 200;
if istest
    opts.tf_plot.info = 1;
else
    opts.tf_plot.info = 2;
end

opts.tf_plot.type = 'sigma';

out = mess_tf_plot(eqn, opts, oper, ROM);
err = out.err;
if istest
    if max(err) > 1000
        mess_err(opts, 'TEST:accuracy', 'unexpectedly inaccurate result %g', max(err));
    end
end
Computing TFMs of original and reduced order systems and MOR errors
 Step  20 / 200 Step  40 / 200 Step  60 / 200 Step  80 / 200 Step 100 / 200 Step 120 / 200 Step 140 / 200 Step 160 / 200 Step 180 / 200 Step 200 / 200