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 © 2009-2022 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=0;
end
format long e;
% set operation
oper = operatormanager('so_2');
% Problem data

n1=500;
alpha=.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=1;

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 = 0;
%opts.adi.accumulateDeltaK = 0;
opts.adi.compute_sol_fac = 1;

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);
fprintf(1,'mess_lradi took %6.2f seconds \n',t_elapsed);

if istest
    if min(outB.res)>=1e-1
       error('MESS: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

disp('size outB.Z:');
disp(size(outB.Z));
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.002105e+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.188758e-01 
ADI step:   68 normalized residual: 1.172664e+00 
ADI step:   70 normalized residual: 1.129371e+00 
ADI step:   71 normalized residual: 8.961118e-01 
ADI step:   73 normalized residual: 8.713570e-01 
ADI step:   75 normalized residual: 7.378360e-01 
ADI step:   77 normalized residual: 5.279292e-01 
ADI step:   78 normalized residual: 8.063871e-01 
ADI step:   80 normalized residual: 6.666642e-01 
ADI step:   81 normalized residual: 7.099135e-01 
ADI step:   83 normalized residual: 6.526697e-01 
ADI step:   85 normalized residual: 4.615173e-01 
ADI step:   87 normalized residual: 4.450099e-01 
ADI step:   88 normalized residual: 5.715111e-01 
ADI step:   90 normalized residual: 5.039108e-01 
ADI step:   92 normalized residual: 3.200317e-01 
ADI step:   93 normalized residual: 4.143070e-01 
ADI step:   95 normalized residual: 4.217265e-01 
ADI step:   97 normalized residual: 3.333128e-01 
ADI step:   98 normalized residual: 3.777849e-01 
ADI step:  100 normalized residual: 3.617540e-01 
ADI step:  101 normalized residual: 3.912230e-01 
ADI step:  103 normalized residual: 3.370220e-01 
ADI step:  105 normalized residual: 3.317460e-01 
ADI step:  107 normalized residual: 3.355299e-01 
ADI step:  108 normalized residual: 3.173677e-01 
ADI step:  110 normalized residual: 2.591600e-01 
ADI step:  111 normalized residual: 2.713324e-01 
ADI step:  113 normalized residual: 2.680951e-01 
ADI step:  115 normalized residual: 2.582824e-01 
ADI step:  117 normalized residual: 2.575210e-01 
ADI step:  118 normalized residual: 1.824294e-01 
ADI step:  120 normalized residual: 2.049555e-01 
ADI step:  121 normalized residual: 2.397734e-01 
ADI step:  123 normalized residual: 2.421403e-01 
ADI step:  124 normalized residual: 2.247502e-01 
ADI step:  126 normalized residual: 2.264447e-01 
ADI step:  128 normalized residual: 2.230960e-01 
ADI step:  130 normalized residual: 2.271362e-01 
ADI step:  132 normalized residual: 2.325429e-01 
ADI step:  133 normalized residual: 2.222306e-01 
ADI step:  135 normalized residual: 1.468320e-01 
ADI step:  137 normalized residual: 1.418196e-01 
ADI step:  138 normalized residual: 1.719007e-01 
ADI step:  139 normalized residual: 1.616633e-01 
ADI step:  141 normalized residual: 1.078697e-01 
ADI step:  143 normalized residual: 1.049704e-01 
ADI step:  144 normalized residual: 1.892963e-01 
ADI step:  146 normalized residual: 1.895386e-01 
ADI step:  148 normalized residual: 1.778348e-01 
ADI step:  149 normalized residual: 2.060114e-01 
ADI step:  151 normalized residual: 2.060300e-01 
ADI step:  153 normalized residual: 1.763440e-01 
ADI step:  155 normalized residual: 1.804357e-01 
ADI step:  157 normalized residual: 1.739263e-01 
ADI step:  158 normalized residual: 1.921374e-01 
ADI step:  160 normalized residual: 1.461589e-01 
ADI step:  162 normalized residual: 1.433190e-01 
ADI step:  163 normalized residual: 1.485836e-01 
ADI step:  165 normalized residual: 1.227408e-01 
ADI step:  167 normalized residual: 1.224209e-01 
ADI step:  168 normalized residual: 1.295107e-01 
ADI step:  170 normalized residual: 1.159963e-01 
ADI step:  172 normalized residual: 1.100151e-01 
ADI step:  173 normalized residual: 1.154309e-01 
ADI step:  174 normalized residual: 9.360058e-02 
ADI step:  176 normalized residual: 9.177159e-02 
ADI step:  178 normalized residual: 8.333038e-02 
ADI step:  179 normalized residual: 7.422130e-02 
ADI step:  181 normalized residual: 7.365669e-02 
ADI step:  183 normalized residual: 7.292316e-02 
ADI step:  185 normalized residual: 7.307011e-02 
ADI step:  187 normalized residual: 7.096115e-02 
ADI step:  188 normalized residual: 8.845479e-02 
ADI step:  190 normalized residual: 9.420441e-02 
ADI step:  192 normalized residual: 9.470731e-02 
ADI step:  193 normalized residual: 9.967799e-02 
ADI step:  194 normalized residual: 9.501825e-02 
ADI step:  196 normalized residual: 9.489766e-02 
ADI step:  198 normalized residual: 9.454401e-02 
ADI step:  200 normalized residual: 9.479503e-02 
ADI step:  202 normalized residual: 9.439792e-02 
ADI step:  203 normalized residual: 1.034206e-01 
ADI step:  204 normalized residual: 1.074042e-01 
ADI step:  206 normalized residual: 1.052745e-01 
ADI step:  208 normalized residual: 1.052426e-01 
ADI step:  210 normalized residual: 1.051079e-01 
ADI step:  212 normalized residual: 1.042939e-01 
ADI step:  213 normalized residual: 1.027145e-01 
ADI step:  215 normalized residual: 1.020585e-01 
ADI step:  217 normalized residual: 1.017780e-01 
ADI step:  218 normalized residual: 1.010942e-01 
ADI step:  220 normalized residual: 1.002537e-01 
ADI step:  222 normalized residual: 1.020959e-01 
ADI step:  223 normalized residual: 9.459831e-02 
ADI step:  225 normalized residual: 9.295336e-02 
ADI step:  226 normalized residual: 7.709465e-02 
ADI step:  228 normalized residual: 7.670378e-02 
ADI step:  229 normalized residual: 6.854768e-02 
ADI step:  231 normalized residual: 6.711718e-02 
ADI step:  233 normalized residual: 6.611763e-02 
ADI step:  234 normalized residual: 5.228905e-02 
ADI step:  236 normalized residual: 5.226955e-02 
ADI step:  238 normalized residual: 5.362616e-02 
ADI step:  239 normalized residual: 9.468063e-02 
ADI step:  241 normalized residual: 9.402312e-02 
ADI step:  243 normalized residual: 9.365705e-02 
ADI step:  244 normalized residual: 6.628590e-02 
ADI step:  246 normalized residual: 6.501904e-02 
ADI step:  248 normalized residual: 6.482750e-02 
ADI step:  249 normalized residual: 7.485161e-02 
ADI step:  251 normalized residual: 7.483962e-02 
ADI step:  253 normalized residual: 7.504306e-02 
ADI step:  254 normalized residual: 6.256013e-02 
ADI step:  256 normalized residual: 6.236421e-02 
ADI step:  258 normalized residual: 6.001722e-02 
ADI step:  259 normalized residual: 5.742161e-02 
ADI step:  261 normalized residual: 5.713425e-02 
ADI step:  263 normalized residual: 5.971008e-02 
ADI step:  265 normalized residual: 6.031225e-02 
ADI step:  267 normalized residual: 6.027308e-02 
ADI step:  268 normalized residual: 6.064504e-02 
ADI step:  270 normalized residual: 6.397064e-02 
ADI step:  272 normalized residual: 6.372679e-02 
ADI step:  273 normalized residual: 6.674236e-02 
ADI step:  274 normalized residual: 3.606668e-02 
ADI step:  276 normalized residual: 3.602915e-02 
ADI step:  278 normalized residual: 3.561930e-02 
ADI step:  279 normalized residual: 4.706202e-02 
ADI step:  281 normalized residual: 4.623834e-02 
ADI step:  283 normalized residual: 4.618432e-02 
ADI step:  284 normalized residual: 3.994970e-02 
ADI step:  286 normalized residual: 3.984032e-02 
ADI step:  288 normalized residual: 3.966268e-02 
ADI step:  289 normalized residual: 3.880402e-02 
ADI step:  291 normalized residual: 3.878791e-02 
ADI step:  293 normalized residual: 3.884510e-02 
ADI step:  295 normalized residual: 3.844234e-02 
ADI step:  297 normalized residual: 3.868438e-02 
ADI step:  298 normalized residual: 3.985688e-02 
ADI step:  299 normalized residual: 5.061408e-02 
ADI step:  301 normalized residual: 5.055444e-02 
mess_lradi took   0.55 seconds 
size outB.Z:
        3002         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.sigma.fmin = 1e-4;
opts.sigma.fmax = 1e0;
opts.sigma.nsample = 200;
if istest
    opts.sigma.info = 1;
else
    opts.sigma.info = 2;
end
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
if istest
    if max(err) > 1000
        error('MESS: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