function [Ar, Br, Cr] = bt_mor_FDM_tol(tol,n0,shifts,istest)
% bt_mor_FDM_tol computes a reduced order model via the standard Lyapunov
% balanced truncation (see e.g. [1]) for a finite difference discretized
% convection diffusion model on the unit square described in [2].
%
% Usage:
%    [Ar, Br, Cr] = bt_mor_FDM_tol(tol,max_ord,n0,test)
%
% Inputs
%
% tol         truncation tolerance for the Hankel singular values
%             (optional; defalts to 1e-6)
%
% n0          n0^2 gives the dimension of the original model, i.e. n0 is
%             the number of degrees of freedom, i.e. grid points, per
%             spatial direction
%             (optional; defaults to 50)
%
% shifts      shift selection used in ADI;  possible choices:
%               'heur'       :   Penzl heuristic shifts
%               'projection' :   projection shifts using the last columns
%                                of the solution factor
%             (optional, defaults to 'heur')
%
% istest      flag to determine whether this demo runs as a CI test or
%             interactive demo
%             (optional, defaults to 0, i.e. interactive demo)
%
% Outputs
%
% Ar, Br, Cr  the reduced orde system matrices.
%
% References
% [1] 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.
%
% [2] T. Penzl, Lyapack Users Guide, Tech. Rep. SFB393/00-33,
%     Sonderforschungsbereich 393 Numerische Simulation auf massiv
%     parallelen Rechnern, TU Chemnitz, 09107 Chemnitz, Germany,
%     available from http://www.tu-chemnitz.de/sfb393/sfb00pr.html. (2000).


% 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);
% BT tolerance and maximum order for the ROM
if nargin<1, tol=1e-6; end
if nargin<2, n0 = 50; end
if nargin<3, shifts='heur'; end
if nargin<4, istest = 0; end

% ADI tolerance and maximum iteration number
opts.adi.maxiter = 100;
opts.adi.res_tol = 1e-9;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 1;
opts.norm = 'fro';

operations

oper = operatormanager('default');

% Problem data
eqn.A_ = fdm_2d_matrix(n0,'10*x','100*y','0');
eqn.B = fdm_2d_vector(n0,'.1<x<=.3');
eqn.C = fdm_2d_vector(n0,'.7<x<=.9');eqn.C=eqn.C';
eqn.haveE = 0;

% problem dimension
n=oper.size(eqn, opts);
%Heuristic Parameters via basic Arnoldi or Projection shifts
opts.shifts.method = shifts;
switch opts.shifts.method
    case 'heur'
        opts.shifts.num_desired=25;
        opts.shifts.num_Ritz=50;
        opts.shifts.num_hRitz=25;
        opts.shifts.b0=ones(n,1);
    case 'projection'
        opts.shifts.method = 'projection';
        opts.shifts.num_desired = 10;
    otherwise
        warning('MESS:Test:invalid_parameter_fallback',...
            'invalid shift selection, falling back to "heur"');
        opts.shifts.num_desired=25;
        opts.shifts.num_Ritz=50;
        opts.shifts.num_hRitz=25;
        opts.shifts.b0=ones(n,1);
        opts.shifts.method = 'heur';
end

opts.shifts.p=mess_para(eqn,opts,oper);

disp('(Initial) ADI shifts');
disp(opts.shifts.p);
(Initial) ADI shifts
     -2.065258190191611e+04 + 0.000000000000000e+00i
     -1.937292407469257e+04 + 0.000000000000000e+00i
     -1.749376770662731e+04 + 0.000000000000000e+00i
     -1.330511046698838e+04 + 0.000000000000000e+00i
     -9.263930276803778e+03 + 0.000000000000000e+00i
     -6.671163618267795e+03 + 0.000000000000000e+00i
     -1.698501218018387e+03 + 0.000000000000000e+00i
     -4.649900486215604e+02 + 0.000000000000000e+00i
     -2.603557887074633e+02 + 0.000000000000000e+00i
     -1.918393347707617e+02 + 0.000000000000000e+00i
     -1.112802038431261e+02 + 0.000000000000000e+00i
     -4.453715272674594e+03 + 9.448801986471070e+02i
     -4.453715272674594e+03 - 9.448801986471070e+02i
     -3.511454699403188e+03 + 1.017535965450999e+03i
     -3.511454699403188e+03 - 1.017535965450999e+03i
     -1.941507072350244e+03 + 7.426412079916519e+02i
     -1.941507072350244e+03 - 7.426412079916519e+02i
     -9.796882754126804e+02 + 9.470372612793864e+02i
     -9.796882754126804e+02 - 9.470372612793864e+02i
     -8.124319370443204e+02 + 4.557291003561135e+02i
     -8.124319370443204e+02 - 4.557291003561135e+02i
     -3.712015209114301e+02 + 2.393803494728078e+02i
     -3.712015209114301e+02 - 2.393803494728078e+02i
     -1.392548614286558e+02 + 2.618919038961343e+01i
     -1.392548614286558e+02 - 2.618919038961343e+01i

%observability
eqn.type='N';
tic;
outB  = mess_lradi(eqn, opts, oper);
toc;

if istest
    if min(outB.res)>=opts.adi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outB.res);
    title('AX + XA^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: 8.428057e-01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 7.378016e-01 relative change in Z: 6.917507e-01
ADI step:    3 normalized residual: 6.508479e-01 relative change in Z: 5.637997e-01
ADI step:    4 normalized residual: 5.614350e-01 relative change in Z: 5.169589e-01
ADI step:    5 normalized residual: 4.639724e-01 relative change in Z: 4.948088e-01
ADI step:    6 normalized residual: 3.674147e-01 relative change in Z: 4.645217e-01
ADI step:    7 normalized residual: 1.830519e-01 relative change in Z: 5.907596e-01
ADI step:    8 normalized residual: 4.914772e-02 relative change in Z: 5.625659e-01
ADI step:    9 normalized residual: 1.439251e-02 relative change in Z: 3.708279e-01
ADI step:   10 normalized residual: 2.640914e-03 relative change in Z: 2.116408e-01
ADI step:   11 normalized residual: 1.583523e-04 relative change in Z: 7.793154e-02
ADI step:   13 normalized residual: 1.066545e-04 relative change in Z: 8.368038e-03
ADI step:   15 normalized residual: 7.634587e-05 relative change in Z: 7.751726e-03
ADI step:   17 normalized residual: 4.229441e-05 relative change in Z: 8.133399e-03
ADI step:   19 normalized residual: 1.675468e-05 relative change in Z: 6.178842e-03
ADI step:   21 normalized residual: 2.303907e-06 relative change in Z: 4.177396e-03
ADI step:   23 normalized residual: 1.919842e-07 relative change in Z: 1.872833e-03
ADI step:   25 normalized residual: 1.238977e-10 relative change in Z: 6.576621e-04
Elapsed time is 0.110372 seconds.
size outB.Z:
        2500          25

%controllability
eqn.type = 'T';
tic;
outC = mess_lradi(eqn, opts, oper);
toc;

if istest
    if min(outC.res)>=opts.adi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(outC.res);
    title('A^TX + XA = -C^TC');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
disp('size outC.Z:');
disp(size(outC.Z));
ADI step:    1 normalized residual: 9.739434e-01 relative change in Z: 1.000000e+00
ADI step:    2 normalized residual: 9.600200e-01 relative change in Z: 7.149146e-01
ADI step:    3 normalized residual: 9.507972e-01 relative change in Z: 5.989262e-01
ADI step:    4 normalized residual: 9.434197e-01 relative change in Z: 5.644110e-01
ADI step:    5 normalized residual: 9.378728e-01 relative change in Z: 5.588782e-01
ADI step:    6 normalized residual: 9.355761e-01 relative change in Z: 5.491134e-01
ADI step:    7 normalized residual: 9.533528e-01 relative change in Z: 7.354621e-01
ADI step:    8 normalized residual: 1.080148e+00 relative change in Z: 8.167288e-01
ADI step:    9 normalized residual: 1.086516e+00 relative change in Z: 7.398403e-01
ADI step:   10 normalized residual: 4.755716e-01 relative change in Z: 5.745268e-01
ADI step:   11 normalized residual: 5.254634e-02 relative change in Z: 2.652358e-01
ADI step:   13 normalized residual: 3.985457e-02 relative change in Z: 4.242675e-02
ADI step:   15 normalized residual: 2.784426e-02 relative change in Z: 4.005181e-02
ADI step:   17 normalized residual: 1.387804e-02 relative change in Z: 4.044519e-02
ADI step:   19 normalized residual: 4.202789e-03 relative change in Z: 2.932766e-02
ADI step:   21 normalized residual: 2.110828e-04 relative change in Z: 1.612450e-02
ADI step:   23 normalized residual: 7.179830e-07 relative change in Z: 3.008242e-03
ADI step:   25 normalized residual: 2.305240e-07 relative change in Z: 1.640317e-04
ADI step:   26 normalized residual: 1.970927e-07 relative change in Z: 3.081957e-05
ADI step:   27 normalized residual: 1.682114e-07 relative change in Z: 2.941062e-05
ADI step:   28 normalized residual: 1.425346e-07 relative change in Z: 2.854049e-05
ADI step:   29 normalized residual: 1.163449e-07 relative change in Z: 2.983788e-05
ADI step:   30 normalized residual: 8.919103e-08 relative change in Z: 3.178993e-05
ADI step:   31 normalized residual: 6.415080e-08 relative change in Z: 3.227039e-05
ADI step:   32 normalized residual: 2.271836e-08 relative change in Z: 4.511281e-05
ADI step:   33 normalized residual: 2.096728e-09 relative change in Z: 3.126206e-05
ADI step:   34 normalized residual: 8.476203e-10 relative change in Z: 7.747729e-06
Elapsed time is 0.159134 seconds.
size outC.Z:
        2500          34

opts.srm.tol=tol;
opts.srm.info=1;
[TL, TR, hsv, eqn, opts, ~] = mess_square_root_method(eqn,opts,oper,...
    outB.Z,outC.Z);

Ar = TL'*(eqn.A_*TR);
Br = TL'*eqn.B;
Cr = eqn.C*TR;

opts.sigma.nsample = 200;  % 200 frequency samples
opts.sigma.fmin = -3;      % min. frequency 1e-3
opts.sigma.fmax = 4;       % max. frequency 1e4
if istest
    opts.sigma.info=1;     % no output
else
    opts.sigma.info = 2;   % show messages and plots
end
ROM.A = Ar;
ROM.B = Br;
ROM.C = Cr;
ROM.E = eye(size(ROM.A,1));

out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;

if istest
    if max(err)>tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure;
    semilogy(hsv);
    title('Computed Hankel singular values');
    xlabel('index');
    ylabel('magnitude');
end
reduced system order: 12  (max possible/allowed: 25/2500)

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


ans =

  Columns 1 through 3

    -8.829057200166768e+01    -1.019807420789014e+00    -9.445425120873438e+01
     1.019809700641548e+00    -6.032814822854384e-02    -2.738130320613079e+02
    -9.445425124221941e+01     2.738130476572438e+02    -6.449423333932912e+02
    -3.895979695274298e+01     5.269097785916864e+00    -4.720956862483179e+02
     2.657259610089089e+00    -3.229743304497948e-01     3.544604944344591e+01
    -2.007136244229511e+01     2.505374132701651e+00    -2.609151425920568e+02
    -7.941125513916729e+00     9.846878060831060e-01    -1.038824212154909e+02
    -6.650291192945058e-01     8.230925023975408e-02    -8.715254979717750e+00
    -1.933456615594407e+00     2.394417042909453e-01    -2.532332230991804e+01
     8.840248681283010e-02    -1.094672911079475e-02     1.157997721452221e+00
    -9.859754967001821e-01     1.220951580802673e-01    -1.291468089680958e+01
     5.933504024874833e-02    -7.346286808055215e-03     7.772053996853828e-01

  Columns 4 through 6

    -3.895979674387746e+01    -2.657256449973046e+00    -2.007136239791246e+01
    -5.269087839516895e+00    -3.229719701910339e-01    -2.505368460790542e+00
    -4.720956816023990e+02    -3.544601011733786e+01    -2.609151414699180e+02
    -1.198154522101892e+03    -2.091589392939177e+02    -1.124891937041254e+03
     2.091591199188712e+02    -2.488425796556966e+01    -6.521257009031148e+02
    -1.124891937831908e+03     6.521260437433831e+02    -3.364971173014024e+03
    -4.790746202810799e+02     1.603658658652016e+02    -2.275787585528893e+03
    -4.100898102733039e+01     1.212810319522633e+01    -2.385831669273557e+02
    -1.183748028283387e+02     3.640648114391952e+01    -6.412041706538694e+02
     5.421138574529287e+00    -1.652442001380978e+00     2.982625078550474e+01
    -6.041914951522367e+01     1.849139038149388e+01    -3.300697867036573e+02
     3.636575821763832e+00    -1.112034912791116e+00     1.990042621069726e+01

  Columns 7 through 9

    -7.941124327467590e+00     6.650193716657503e-01    -1.933454879063501e+00
    -9.846840714258981e-01     8.230685197982039e-02    -2.394416150829136e-01
    -1.038824070795128e+02     8.715128575640250e+00    -2.532329868077609e+01
    -4.790745158538591e+02     4.100836444178677e+01    -1.183747120669450e+02
    -1.603657970798759e+02     1.212800776525438e+01    -3.640633077964564e+01
    -2.275786982086108e+03     2.385800241693510e+02    -6.412036399833823e+02
    -3.102754701000881e+03     8.481464517018636e+02    -1.421087701775793e+03
    -8.481508882672374e+02    -5.616941834441856e+01     3.902170010744658e+02
    -1.421090573132168e+03    -3.902234082395105e+02    -2.913330645316168e+03
     7.184392791618765e+01     1.358213962932752e+01     6.834788023398814e+02
    -7.647495486709025e+02    -1.698252188863975e+02    -2.658153174937496e+03
     4.653399280207658e+01     9.935020544111936e+00     1.896467105169808e+02

  Columns 10 through 12

    -8.839465760239280e-02    -9.859739496303747e-01    -5.917363441388186e-02
    -1.094541017677039e-02    -1.220951463707010e-01    -7.327480813125683e-03
    -1.157895568422158e+00    -1.291466004288380e+01    -7.750897632465272e-01
    -5.420651529498241e+00    -6.041906542795733e+01    -3.626712140989104e+00
    -1.652327004379968e+00    -1.849128882615374e+01    -1.108880341142835e+00
    -2.982361743621728e+01    -3.300693060137281e+02    -1.984645671967601e+01
    -7.183851286384869e+01    -7.647474632338972e+02    -4.640645751367969e+01
     1.358163344690188e+01     1.698209970079201e+02     9.903309425863881e+00
    -6.834610280793736e+02    -2.658152870223816e+03    -1.892032426369110e+02
    -9.978790743996990e+00    -2.758794823644630e+02    -1.220734908508392e+01
     2.759007134037782e+02    -6.430581526719097e+03    -1.505522035506696e+03
    -1.226347412651648e+01     1.506870839135303e+03    -4.757197654714582e+01