function IRKA_TripleChain(n1, usfs, istest)
% IRKA_TripleChain computes first order IRKA based ROM for the triple chain
%
% Usage:      IRKA_TripleChain(n1, oper, istest)
%
% Inputs:
%
% n1          length of a single chain in the model
%             (optional; defaults to 1000)
%
% usfs        the set of user supplied functions to use, or the indication to
%             call IRKA with matrices.
%             possible values: 'so_1', 'so_2', 'default', 'matrices'
%             (optional; defaults to 'so_1')
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%

%
% 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, 3);
if nargin < 1
    n1 = 1000;
end
if nargin < 2
    usfs = 'so_1';
end
if nargin < 3
    istest = 0;
end

set usfs, unless we want to use the matrices

opts = struct();
if strcmp(usfs, 'matrices')
    [oper, opts] = operatormanager(opts, 'so_2');
else
    [oper, opts] = operatormanager(opts, usfs);
end

% Problem data
alpha = .002;
Beta  = .002;
v     = 5;

matrices = false;
switch usfs
    case 'so_1'
        [eqn.M_, eqn.E_, eqn.K_] = ...
            triplechain_MSD(n1, alpha, Beta, v);
        s  = size(eqn.K_, 1);
        eqn.B = [zeros(s, 1); ones(s, 1)];
        eqn.C = [ones(1, s) zeros(1, s)];
    case 'so_2'
        [eqn.M_, eqn.E_, eqn.K_] = ...
            triplechain_MSD(n1, alpha, Beta, v);
        s  = size(eqn.K_, 1);
        eqn.B = [ones(s, 1); zeros(s, 1)];
        eqn.C = eqn.B';
    case 'default'
        [M_, E_, K_] = triplechain_MSD(n1, alpha, Beta, v);
        s = size(K_, 1);
        eqn.A_ = [sparse(s, s), -K_; -K_, -E_];
        eqn.E_ = [-K_, sparse(s, s); sparse(s, s), M_];
        eqn.B = [zeros(s, 1); ones(s, 1)];
        eqn.C = [ones(1, s) zeros(1, s)];
        clear M_ E_ K_ s;
    case 'matrices'
        % setup matrices for the IRKA call
        [M, E, K] = triplechain_MSD(n1, alpha, Beta, v);
        s = size(K, 1);
        B = ones(s, 1);
        Cp = B';
        Cv = zeros(1, s);
        % We will need eqn for tf_plot below
        eqn = struct('M_', M, 'E_', E, 'K_', K, ...
                     'B', [B; zeros(s, 1)], ...
                     'C', [Cp Cv]);
        matrices = true;
end
eqn.haveE = 1;
opts.irka = struct('r', 30, ...
                   'maxiter', 100, ...
                   'shift_tol', 1e-3, ...
                   'h2_tol', 1e-8, ...
                   'num_prev_shifts', 5, ...
                   'num_prev_roms', 5, ...
                   'flipeig', false, ...
                   'init', 'subspace');
if matrices
    [Er, Ar, Br, Cr, Dr, outinfo] = ...
        mess_tangential_irka(M, E, K, B, Cp, Cv, [], opts);
else
    [Er, Ar, Br, Cr, Dr, outinfo] = ...
        mess_tangential_irka(eqn, opts, oper);
end
ROM.E = Er;
ROM.A = Ar;
ROM.B = Br;
ROM.C = Cr;
ROM.D = Dr;

opts.tf_plot.nsample = 400;  % 400 frequency samples
opts.tf_plot.fmin = -4;      % min. frequency 1e-3
opts.tf_plot.fmax = 0;       % max. frequency 1e4
if istest
    opts.tf_plot.info = 0;
else
    opts.tf_plot.info = 3;
end
opts.tf_plot.type = 'sigma';

[out, ~, opts, ~] = mess_tf_plot(eqn, opts, oper, ROM);

if istest
    if isequal(outinfo.term_flag, 'maxiter')
        mess_err(opts, 'TEST:accuracy', ...
                 'terminated with maximum number of iterations');
    end
    if any(out.relerr > 1)
        mess_err(opts, 'TEST:accuracy', ...
                 'unexpectedly inaccurate results');

    end
end
Warning: IRKA step 1 : 15 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   1, rel. chg. shifts = 2.729673e+23 , rel. H2-norm chg. ROM = 1.000000e+00  
Warning: IRKA step 2 : 11 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   2, rel. chg. shifts = 7.228632e+03 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 3 : 11 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   3, rel. chg. shifts = 1.000000e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 4 : 5 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   4, rel. chg. shifts = 1.000000e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 5 : 4 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   5, rel. chg. shifts = 1.000233e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 6 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   6, rel. chg. shifts = 1.000323e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 7 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   7, rel. chg. shifts = 1.000197e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 8 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   8, rel. chg. shifts = 1.034142e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 9 : 4 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step   9, rel. chg. shifts = 6.143890e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 10 : 3 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  10, rel. chg. shifts = 2.112838e+00 , rel. H2-norm chg. ROM = Inf  
IRKA step  11, rel. chg. shifts = 3.429764e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 12 : 4 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  12, rel. chg. shifts = 8.136740e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 13 : 4 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  13, rel. chg. shifts = 6.144289e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 14 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  14, rel. chg. shifts = 4.452914e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 15 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  15, rel. chg. shifts = 5.367510e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 16 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  16, rel. chg. shifts = 1.631884e+00 , rel. H2-norm chg. ROM = Inf  
Warning: IRKA step 17 : 2 non-stable reduced eigenvalues detected.
 
↳ 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/IRKA_TripleChain.m',106)">IRKA_TripleChain (line 106)</a>
↳ In <a href="matlab:opentoline('/builds/mess/mmess/mor/mess_tangential_irka.m',369)">mess_tangential_irka (line 369)</a>
IRKA step  17, rel. chg. shifts = 2.278493e+00 , rel. H2-norm chg. ROM = Inf  
IRKA step  18, rel. chg. shifts = 1.016258e+01 , rel. H2-norm chg. ROM = 1.226916e-03  
IRKA step  19, rel. chg. shifts = 8.854805e+00 , rel. H2-norm chg. ROM = 6.931831e-04  
IRKA step  20, rel. chg. shifts = 4.613083e+00 , rel. H2-norm chg. ROM = 3.863685e-05  
IRKA step  21, rel. chg. shifts = 5.687914e+00 , rel. H2-norm chg. ROM = 6.879632e-06  
IRKA step  22, rel. chg. shifts = 4.613082e+00 , rel. H2-norm chg. ROM = 1.412944e-06  
IRKA step  23, rel. chg. shifts = 1.542359e-01 , rel. H2-norm chg. ROM = 6.332403e-07  
IRKA step  24, rel. chg. shifts = 4.070484e-04 , rel. H2-norm chg. ROM = 4.232296e-07  
IRKA terminated due to relative change of shifts criterion.

Computing TFMs of original and reduced order systems and MOR errors
 Step  40 / 400 Step  80 / 400 Step 120 / 400 Step 160 / 400 Step 200 / 400 Step 240 / 400 Step 280 / 400 Step 320 / 400 Step 360 / 400 Step 400 / 400

end