function [Ar, Br, Cr] = lqgbt_mor_FDM(tol,max_ord,n0,istest)
% LQGBT_MOR_FDM computes a reduced order model via the lienar quadratic
% Gaussian balanced truncation [1] for a finite difference discretized
% convection diffusion model on the unit square described in [2].
%
% Usage:
%    [Ar, Br, Cr] = lqgbt_mor_FDM(tol,max_ord,n0,test)
%
% Inputs
%
% tol         truncation tolerance for the LQG characteristic values
%
% max_ord     maximum allowed order for the reuced order model
%
% n0          n0^2 gives the dimension of the original model, i.e. n0 is
%             the number of degrees of freedom per spatial direction
%
% 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] D. Mustafa, K. Glover, Controller design by H∞ -balanced truncation,
%     IEEE Trans. Autom. Control 36 (6) (1991) 668–682.
%     https://doi.org/10.1109/9.86941.
%
% [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-2017
%
narginchk(0,4);

LQGBT Tolerance and maximum order of the ROM.

if nargin<1
    tol     = 1e-6;
end
if nargin<2
    max_ord = 250;
end
if nargin<3
    n0 = 60; % n0 = number of grid points in either space direction;
             % n = n0^2 is the problem dimension!
             % (Change n0 to generate problems of different size.)
end
if nargin<4
    istest=0;
end
% 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<=.8'), fdm_2d_vector(n0,'.8<x<=.9')];
eqn.C  = eqn.C';

Set operations.

oper = operatormanager('default');

ADI options.

opts.adi.maxiter          = 200;
opts.adi.res_tol           = 1.0e-10;
opts.adi.rel_diff_tol            = 0;
opts.adi.info             = 0;
opts.adi.accumulateK      = 1;
opts.adi.accumulateDeltaK = 0;
opts.adi.compute_sol_fac  = 1;

n = oper.size(eqn, opts);
opts.shifts.num_desired = 25;
opts.shifts.num_Ritz = 50;
opts.shifts.num_hRitz = 25;

opts.shifts.info   = 0;
opts.shifts.method = 'heur';
opts.shifts.b0     = ones(n,1);
opts.norm          = 'fro';
% Newton options.
opts.nm.maxiter       = 25;
opts.nm.res_tol        = 1e-10;
opts.nm.rel_diff_tol         = 1e-16;
opts.nm.info          = 1;
opts.nm.accumulateRes = 1;
opts.nm.linesearch    = 1;
opts.nm.inexact       = 'quadratic';
opts.nm.tau           = 0.1;

Solve the filter Riccati equation. A*X + X*A' - X*C'*C*X + B*B' = 0

tic;
eqn.type = 'N';
outC = mess_lrnm(eqn, opts, oper);
toc;

if istest
    if min(outC.res)>=opts.nm.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure(1);
    semilogy(outC.res);
    title('AX + XA^T - XC^TCX + BB^T = 0');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
NM step:    1  normalized residual: 8.637261e-01
               relative change in K: 1.000000e+00
               number of ADI steps: 10 


NM step:    2  normalized residual: 1.247428e-01
               relative change in K: 7.069622e-01
               number of ADI steps: 9 


NM step:    3  normalized residual: 1.105634e-02
               relative change in K: 2.214896e-01
               number of ADI steps: 13 


NM step:    4  normalized residual: 2.578772e-04
               relative change in K: 3.403586e-02
               number of ADI steps: 20 


NM step:    5  normalized residual: 2.177720e-07
               relative change in K: 9.882326e-04
               number of ADI steps: 26 


NM step:    6  normalized residual: 3.531935e-11
               relative change in K: 6.630748e-07
               number of ADI steps: 24 

Elapsed time is 2.584722 seconds.

Solve the regulator Riccati equation. A'*X + X*A - X*B*B'*X + C'*C = 0

tic;
eqn.type = 'T';
outB = mess_lrnm(eqn, opts, oper);
toc;

if istest
    if min(outB.res)>=opts.nm.res_tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
else
    figure(2);
    semilogy(outB.res);
    title('A^TX + XA - XBB^TX + C^TC = 0');
    xlabel('number of iterations');
    ylabel('normalized residual norm');
    pause(1);
end
		Using line search (res: 3.850120e+00)
		lambda: 4.243950e-01

NM step:    1  normalized residual: 3.249656e-01
               relative change in K: 1.000000e+00
               number of ADI steps: 8 


NM step:    2  normalized residual: 2.089833e-02
               relative change in K: 1.547483e-01
               number of ADI steps: 17 


NM step:    3  normalized residual: 8.639340e-05
               relative change in K: 1.010884e-02
               number of ADI steps: 23 


NM step:    4  normalized residual: 5.402815e-10
               relative change in K: 2.528390e-05
               number of ADI steps: 63 


NM step:    5  normalized residual: 3.669164e-11
               relative change in K: 9.263361e-11
               number of ADI steps: 40 

Elapsed time is 2.459279 seconds.

% Model reduction by square root method.

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

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

opts.sigma.nsample = 200;  % 200 frequency samples
opts.sigma.fmin = -2;      % min. frequency 1e-3
opts.sigma.fmax = 6;       % 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;
reduced system order: 14  (max possible/allowed: 21/250)

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

Report.

if istest
    if max(err)>=tol
        error('MESS:TEST:accuracy','unexpectedly inaccurate result');
    end
end
ans =

  Columns 1 through 3

    -1.379787850218477e+02    -2.098991930101151e+02     1.945158039664558e+01
    -2.098993425497630e+02    -6.657215444470442e+02     1.614371684602637e+02
    -1.946667317851269e+01    -1.615077870129067e+02    -1.295849155971231e+01
    -1.024027644293116e+02    -5.391996682124566e+02    -2.024614301379341e+02
    -1.744205843518507e+01    -9.779198798747869e+01    -2.606563434539390e+01
    -7.023755627850534e+01    -3.970330138923267e+02    -1.013966002010898e+02
    -9.150460288643462e+00    -5.209958465729899e+01    -1.286698240536472e+01
    -2.268014554594061e+01    -1.295551486684216e+02    -3.153584685792218e+01
     8.322966600542122e+00     4.759460210521294e+01     1.152776228085584e+01
     9.719557585030199e+00     5.562584502265517e+01     1.342564012423486e+01
     6.939465733257562e-01     3.972197438524026e+00     9.579607739711724e-01
     4.262413356774420e+00     2.439795081839154e+01     5.884498658659648e+00
     6.033862060430781e-01     3.453861992598300e+00     8.329305474459288e-01
     7.483603779847949e-01     4.283731030368058e+00     1.033040700245889e+00

  Columns 4 through 6

    -1.023861517535651e+02    -1.257265753365973e+01    -7.020968286808396e+01
    -5.390823440111986e+02    -6.913639688064477e+01    -3.968588556025126e+02
     2.024322765223696e+02     1.934226789051771e+01     1.013891097414623e+02
    -1.336635014831780e+03    -1.896681405832564e+02    -1.588626390628468e+03
    -3.875847108236453e+02    -1.128331790966220e+02    -3.076351123487116e+02
    -1.588739115888475e+03    -7.769856898243755e+02    -4.294402418865620e+03
    -2.244621207331687e+02    -9.202614780827804e+01    -1.214954771348800e+03
    -5.724844404607552e+02    -2.808235840601703e+02    -2.378021421495675e+03
     2.126703768177391e+02     1.033116404120387e+02     9.463701348687551e+02
     2.505551286842463e+02     1.248710273103652e+02     1.169204825284935e+03
     1.792397684706062e+01     8.924611522241307e+00     8.464215229317914e+01
     1.100708029188558e+02     5.495820945795410e+01     5.190283233035052e+02
     1.558633105562935e+01     7.786827973264802e+00     7.362940427766884e+01
     1.933216470099979e+01     9.656576562709146e+00     9.135211730811348e+01

  Columns 7 through 9

    -5.762907629269254e+00    -2.255820725763411e+01     6.899733044648654e+00
    -3.258465573456419e+01    -1.288411597600884e+02     3.946654083682273e+01
     7.930361340571264e+00     3.147434073436109e+01    -9.302300294728083e+00
    -1.239472326186882e+02    -5.702713128348192e+02     1.731975507251304e+02
    -9.072347262702132e+01    -1.696823799899645e+02     1.016987789518224e+02
     1.431333555576709e+02    -2.368544902722163e+03     7.363474880154822e+02
    -1.073835474461003e+02    -1.711611529017456e+02     1.626087911267472e+02
    -4.903154184481200e+02    -2.699668828733954e+03     5.687170287191307e+02
     1.724395322287991e+02     1.654522865596895e+03    -6.328592216076318e+02
     2.249870929749682e+02     2.121595681156901e+03    -1.319948145407470e+03
     1.606838746635216e+01     1.654458517805755e+02    -9.903282862547407e+01
     9.970848975561816e+01     9.983165364770227e+02    -6.327006238557901e+02
     1.415611734941333e+01     1.430680124160949e+02    -9.166796428898252e+01
     1.754733163546086e+01     1.778768773584253e+02    -1.134959682950343e+02

  Columns 10 through 12

     9.717602262181458e+00     9.592151167264153e-02     4.254669380276249e+00
     5.561342171295966e+01     5.527709909417748e-01     2.435518294052400e+01
    -1.342274925606812e+01    -9.403170850854337e-02    -5.854466430986755e+00
     2.505544701611049e+02     2.134927449869502e+00     1.097349851471159e+02
     8.807496306129146e+01     6.985852993597978e+00     4.189789960824827e+01
     1.169165331058345e+03     8.476629068149602e+00     5.169903579454201e+02
     1.351594033378031e+02     1.368531267303926e+01     6.718592439511617e+01
     2.112695595640415e+03    -4.524755620151915e+00     9.844808283609590e+02
    -9.998790056277833e+02    -5.978103395380602e+01    -5.400936052782504e+02
    -5.459930174060410e+03     8.362096707968609e+02    -4.039196379420825e+03
    -1.353328209723596e+03    -4.054835313754420e+01     3.537779672140710e+01
    -4.058055607424087e+03    -5.076752522830531e+02    -5.804564925906087e+03
    -6.541249424131895e+02    -7.508203931425328e+01    -1.482065562089189e+03
    -8.281040083041950e+02    -8.641576252209734e+01    -1.876005251054115e+03

  Columns 13 through 14

     3.877705080129896e-01     6.857053679354591e-01
     2.216788100672909e+00     3.926920029770919e+00
    -5.608939550648020e-01    -9.288816778976313e-01
     1.022678003802499e+01     1.757265460390456e+01
    -5.110050536178274e-01     9.067141593922733e+00
     4.885020466737890e+01     8.263724241086231e+01
    -2.700274831327105e+00     1.559093172725845e+01
     1.024288020593655e+02     1.546652355795167e+02
    -8.738817096754818e+00    -1.116285562998081e+02
    -4.078788012136729e+02    -7.494790145425161e+02
     5.361441726831008e+01    -4.365927575017179e+01
    -5.012711118639284e+02    -1.738258548102236e+03
    -3.612341476083874e+02     7.709498948108468e+00
    -8.983095857443517e+02    -1.917034827209667e+03