function BT_sym_TripleChain(version,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:
%
% version  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.
%
% [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 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,2);

if nargin==0
    version = '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)

tic;
outB = mess_lradi(eqn, opts, oper);
toc;

if istest
    if min(outB.res)>=1e-1
       error('MESS:TEST:accuracy','unexpectedly inaccurate result');
   end
else
    figure(1);
    semilogy(outB.res);
    title('0= AXM^T + MXA^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.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.813607e+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.961108e-01 
ADI step:   73 normalized residual: 8.713580e-01 
ADI step:   75 normalized residual: 7.378370e-01 
ADI step:   77 normalized residual: 5.279300e-01 
ADI step:   78 normalized residual: 8.063868e-01 
ADI step:   80 normalized residual: 6.666620e-01 
ADI step:   81 normalized residual: 7.099150e-01 
ADI step:   83 normalized residual: 6.526702e-01 
ADI step:   85 normalized residual: 4.615182e-01 
ADI step:   87 normalized residual: 4.450105e-01 
ADI step:   88 normalized residual: 5.715083e-01 
ADI step:   90 normalized residual: 5.039076e-01 
ADI step:   92 normalized residual: 3.200358e-01 
ADI step:   93 normalized residual: 4.142990e-01 
ADI step:   95 normalized residual: 4.217200e-01 
ADI step:   97 normalized residual: 3.333007e-01 
ADI step:   98 normalized residual: 3.777807e-01 
ADI step:  100 normalized residual: 3.617455e-01 
ADI step:  101 normalized residual: 3.912320e-01 
ADI step:  103 normalized residual: 3.370354e-01 
ADI step:  105 normalized residual: 3.317509e-01 
ADI step:  107 normalized residual: 3.355219e-01 
ADI step:  108 normalized residual: 3.173621e-01 
ADI step:  110 normalized residual: 2.591779e-01 
ADI step:  111 normalized residual: 2.713245e-01 
ADI step:  113 normalized residual: 2.680856e-01 
ADI step:  115 normalized residual: 2.582790e-01 
ADI step:  117 normalized residual: 2.575180e-01 
ADI step:  118 normalized residual: 1.824067e-01 
ADI step:  120 normalized residual: 2.048517e-01 
ADI step:  121 normalized residual: 2.396445e-01 
ADI step:  123 normalized residual: 2.420309e-01 
ADI step:  124 normalized residual: 2.250282e-01 
ADI step:  126 normalized residual: 2.267304e-01 
ADI step:  128 normalized residual: 2.233017e-01 
ADI step:  130 normalized residual: 2.275887e-01 
ADI step:  132 normalized residual: 2.323182e-01 
ADI step:  133 normalized residual: 2.233032e-01 
ADI step:  135 normalized residual: 1.474431e-01 
ADI step:  137 normalized residual: 1.422854e-01 
ADI step:  138 normalized residual: 1.711563e-01 
ADI step:  139 normalized residual: 1.590640e-01 
ADI step:  141 normalized residual: 1.009243e-01 
ADI step:  143 normalized residual: 9.804670e-02 
ADI step:  144 normalized residual: 1.874438e-01 
ADI step:  146 normalized residual: 1.877888e-01 
ADI step:  148 normalized residual: 1.842147e-01 
ADI step:  149 normalized residual: 2.021664e-01 
ADI step:  151 normalized residual: 2.020127e-01 
ADI step:  153 normalized residual: 1.791933e-01 
ADI step:  155 normalized residual: 1.766039e-01 
ADI step:  157 normalized residual: 1.592584e-01 
ADI step:  158 normalized residual: 1.749497e-01 
ADI step:  160 normalized residual: 1.720773e-01 
ADI step:  162 normalized residual: 1.675864e-01 
ADI step:  163 normalized residual: 1.559620e-01 
ADI step:  165 normalized residual: 1.566123e-01 
ADI step:  167 normalized residual: 1.591720e-01 
ADI step:  168 normalized residual: 1.588365e-01 
ADI step:  169 normalized residual: 1.641616e-01 
ADI step:  171 normalized residual: 1.639970e-01 
ADI step:  173 normalized residual: 1.451533e-01 
ADI step:  175 normalized residual: 1.383603e-01 
ADI step:  177 normalized residual: 1.234213e-01 
ADI step:  178 normalized residual: 1.053449e-01 
ADI step:  179 normalized residual: 1.262741e-01 
ADI step:  181 normalized residual: 1.191211e-01 
ADI step:  183 normalized residual: 1.181099e-01 
ADI step:  184 normalized residual: 6.863730e-02 
ADI step:  186 normalized residual: 6.824738e-02 
ADI step:  188 normalized residual: 6.207974e-02 
ADI step:  189 normalized residual: 9.478628e-02 
ADI step:  191 normalized residual: 9.501622e-02 
ADI step:  193 normalized residual: 9.264676e-02 
ADI step:  194 normalized residual: 9.925469e-02 
ADI step:  196 normalized residual: 9.861238e-02 
ADI step:  198 normalized residual: 9.929500e-02 
ADI step:  199 normalized residual: 1.043346e-01 
ADI step:  201 normalized residual: 1.040538e-01 
ADI step:  203 normalized residual: 1.035106e-01 
ADI step:  204 normalized residual: 7.429947e-02 
ADI step:  206 normalized residual: 7.597314e-02 
ADI step:  208 normalized residual: 7.518157e-02 
ADI step:  209 normalized residual: 9.231182e-02 
ADI step:  211 normalized residual: 8.773386e-02 
ADI step:  213 normalized residual: 8.725790e-02 
ADI step:  214 normalized residual: 6.214829e-02 
ADI step:  216 normalized residual: 6.135600e-02 
ADI step:  218 normalized residual: 5.987776e-02 
ADI step:  219 normalized residual: 1.004795e-01 
ADI step:  221 normalized residual: 1.003918e-01 
ADI step:  223 normalized residual: 9.785594e-02 
ADI step:  224 normalized residual: 9.320641e-02 
ADI step:  226 normalized residual: 9.283261e-02 
ADI step:  228 normalized residual: 8.572871e-02 
ADI step:  229 normalized residual: 9.553589e-02 
ADI step:  231 normalized residual: 9.657063e-02 
ADI step:  233 normalized residual: 9.249383e-02 
ADI step:  234 normalized residual: 7.123043e-02 
ADI step:  236 normalized residual: 6.958544e-02 
ADI step:  238 normalized residual: 7.414030e-02 
ADI step:  239 normalized residual: 5.894823e-02 
ADI step:  241 normalized residual: 6.038939e-02 
ADI step:  243 normalized residual: 5.916657e-02 
ADI step:  244 normalized residual: 8.411748e-02 
ADI step:  246 normalized residual: 8.425901e-02 
ADI step:  248 normalized residual: 8.243929e-02 
ADI step:  249 normalized residual: 7.637549e-02 
ADI step:  251 normalized residual: 7.632137e-02 
ADI step:  253 normalized residual: 7.705582e-02 
ADI step:  255 normalized residual: 7.636337e-02 
ADI step:  256 normalized residual: 7.639694e-02 
ADI step:  258 normalized residual: 7.923463e-02 
ADI step:  259 normalized residual: 4.060082e-02 
ADI step:  261 normalized residual: 4.059216e-02 
ADI step:  263 normalized residual: 3.871716e-02 
ADI step:  264 normalized residual: 4.788324e-02 
ADI step:  266 normalized residual: 4.702063e-02 
ADI step:  268 normalized residual: 4.721126e-02 
ADI step:  269 normalized residual: 5.602991e-02 
ADI step:  271 normalized residual: 5.538506e-02 
ADI step:  273 normalized residual: 5.433885e-02 
ADI step:  274 normalized residual: 3.250868e-02 
ADI step:  276 normalized residual: 3.236078e-02 
ADI step:  278 normalized residual: 3.151680e-02 
ADI step:  279 normalized residual: 4.724736e-02 
ADI step:  281 normalized residual: 4.753528e-02 
ADI step:  283 normalized residual: 4.491330e-02 
ADI step:  284 normalized residual: 4.464692e-02 
ADI step:  286 normalized residual: 4.490678e-02 
ADI step:  287 normalized residual: 5.450309e-02 
ADI step:  288 normalized residual: 5.716043e-02 
ADI step:  289 normalized residual: 7.086736e-02 
ADI step:  291 normalized residual: 7.097008e-02 
ADI step:  293 normalized residual: 6.753440e-02 
ADI step:  294 normalized residual: 4.793382e-02 
ADI step:  296 normalized residual: 4.775572e-02 
ADI step:  298 normalized residual: 4.622585e-02 
ADI step:  299 normalized residual: 5.323988e-02 
ADI step:  301 normalized residual: 5.316499e-02 
Elapsed time is 0.364077 seconds.
size outB.Z:
        3002         301

switch upper(version)
    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(version,'FO'))
    max_ord = 75;
    tol = eps;
    info = 1;

    [TL,TR] = square_root_method_SO(eqn.M_, max_ord, tol, info, 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