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; defaults 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 order 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 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, 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

controllability

eqn.type = 'N';
t_mess_lradi = tic;
outB  = mess_lradi(eqn, opts, oper);
t_elapsed1 = toc(t_mess_lradi);
fprintf(1,'mess_lradi took %6.2f seconds \n' ,t_elapsed1);

if istest
    if min(outB.res)>=opts.adi.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outB.res,'LineWidth',3);
    title('A X + 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: 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
mess_lradi took   0.32 seconds 
size outB.Z:
        2500          25

observability

eqn.type = 'T';
t_mess_lradi = tic;
outC = mess_lradi(eqn, opts, oper);
t_elapsed2 = toc(t_mess_lradi);
fprintf(1, 'mess_lradi took %6.2f seconds \n', t_elapsed2);

if istest
    if min(outC.res) >= opts.adi.res_tol
        error('MESS:TEST:accuracy', 'unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(outC.res, 'LineWidth', 3);
    title('A^T X + X A = -C^T C');
    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
mess_lradi took   0.65 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, 'LineWidth', 3);
    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 =

   1.0e+03 *

  Columns 1 through 7

   -0.0883   -0.0010   -0.0945   -0.0390   -0.0027   -0.0201   -0.0079
    0.0010   -0.0001   -0.2738   -0.0053   -0.0003   -0.0025   -0.0010
   -0.0945    0.2738   -0.6449   -0.4721   -0.0354   -0.2609   -0.1039
   -0.0390    0.0053   -0.4721   -1.1982   -0.2092   -1.1249   -0.4791
    0.0027   -0.0003    0.0354    0.2092   -0.0249   -0.6521   -0.1604
   -0.0201    0.0025   -0.2609   -1.1249    0.6521   -3.3650   -2.2758
   -0.0079    0.0010   -0.1039   -0.4791    0.1604   -2.2758   -3.1028
   -0.0007    0.0001   -0.0087   -0.0410    0.0121   -0.2386   -0.8482
   -0.0019    0.0002   -0.0253   -0.1184    0.0364   -0.6412   -1.4211
    0.0001   -0.0000    0.0012    0.0054   -0.0017    0.0298    0.0718
   -0.0010    0.0001   -0.0129   -0.0604    0.0185   -0.3301   -0.7647
    0.0001   -0.0000    0.0008    0.0036   -0.0011    0.0199    0.0465

  Columns 8 through 12

    0.0007   -0.0019   -0.0001   -0.0010   -0.0001
    0.0001   -0.0002   -0.0000   -0.0001   -0.0000
    0.0087   -0.0253   -0.0012   -0.0129   -0.0008
    0.0410   -0.1184   -0.0054   -0.0604   -0.0036
    0.0121   -0.0364   -0.0017   -0.0185   -0.0011
    0.2386   -0.6412   -0.0298   -0.3301   -0.0198
    0.8481   -1.4211   -0.0718   -0.7647   -0.0464
   -0.0562    0.3902    0.0136    0.1698    0.0099
   -0.3902   -2.9133   -0.6835   -2.6582   -0.1892
    0.0136    0.6835   -0.0100   -0.2759   -0.0122
   -0.1698   -2.6582    0.2759   -6.4306   -1.5055
    0.0099    0.1896   -0.0123    1.5069   -0.0476