...
 
Commits (2)
## version 2.0
- New RADI iteration for AREs
- New splitting methods for autonomous DREs
- New splitting and BDF methods for non-autonomous DREs
- New operator manager only requires non-empty functions and replaces
non-existent ones with a general `mess_do_nothing` function
- improved Riccati iteration
- updated minimum required/recommended Matlab and octave versions
(see `DEPENDENCIES.md`)
- unified function interfaces for top level calls
- unified handling of low rank updated operators. Now always A+UV' is
used. (Note the sign of the update and the transposition in V)
- major updates in the MOR routines
- some resturcturing in the opts structure.
* `opts.adi.shifts` has moved to `opts.shifts` such that also RADI
can use it independent of ADI
* opts.norm now determines the norm for all methods rather than
having to consistently specifiy the same norm in each substructure
* initial feedbacks for the Riccati solvers are now stored in the
`opts` structure for the method rather than `eqn`
- The projection shift routine uses the flag `opts.shifts.implicitVtAV`.
Default is `true`. If set to `false` $A\cdot{}V$ is computed explicitly.
- several consistency updates and bug fixes
- general code cleaning and pretty printing
- redesign of the demos
- turned scripts into actual demo functions
- new demos for indefinite AREs and H-infinity control
- CI testing
- demos serve as system tests
- additional unit tests for the smaller building blocks and backend routines
---
## version 1.0.1
- Minor consistency and bug fixes and improved integrity of metafiles.
- updated documentation
- Removed replacements directory since its content was not needed for Matlab
after release 2010b and Octave after 4.0.
- Removed replacements directory since its content was not needed for
Matlab after release 2010b and Octave after 4.0.
---
## version 1.0
Compared to the predecessor lyapack a couple of things have changed.
- The user supplied functions are now managed by an operator manager
- The low rank ADI now has:
- optimized treatment of E matrices in generalized equations
- more choices for shift selection, including completely automatic generation of shifts
- more choices for shift selection, including completely automatic
generation of shifts
- improved stopping criteria based on low rank factors of the current residual
- automatic generation of real low rank factors also for complex shifts
- The Newton-Kleinman iteration features:
......@@ -19,5 +53,7 @@ Compared to the predecessor lyapack a couple of things have changed.
- Examples have been extended
- The Riccati iteration for H-infinity Riccati equations was added
- DSPMR has not yet been ported to the new infrastructure
- The SRM routine for balanced truncation is only available for none-DAE systems. Still, DAE versions are included in the corresponding DEMOS.
- The SRM routine for balanced truncation is only available for
none-DAE systems. Still, DAE versions are included in the
corresponding DEMOS.
- A tangential IRKA implementation for none-DAE systems was added
M.M.E.S.S. - The Matrix Equation Sparse Solver Library for MATLAB and Octave
M-M.E.S.S. - The Matrix Equation Sparse Solver Library for MATLAB and Octave
============================================================================
##General
......@@ -6,20 +6,21 @@ please refer to the M-M.E.S.S. toolbox by at least its acronym and the
version number. The core Copyright holders are Jens Saak, Martin
Köhler and Peter Benner. They should be given as the authors in a
citation. The project is frequently updated on Zenodo and every stable
release will receive a DOI for proper Citation there. A sample BibTeX
entry can be found below.
release will receive a DOI for proper citation there. The DOI for this
version and a sample BibTeX entry can be found below.
##DOI
The DOI for version 1.0.1 is
[10.5281/zenodo.50575](dx.doi.org/10.5281/zenodo.50575).
The DOI for version 2.0 is
[10.5281/zenodo.3368844](http://doi.org/10.5281/zenodo.3368844)
##BibTeX
@Misc{SaaKB16-mmess-1.0.1,
@Misc{SaaKB19-mmess-2.0,
key = {MMESS},
author = {Saak, J. and K\"{o}hler, M. and Benner, P.},
title = {M-M.E.S.S.-1.0.1 -- The Matrix Equations Sparse Solvers library},
howpublished = {DOI:10.5281/zenodo.50575},
month = apr,
year = 2016,
title = {{M-M.E.S.S.}-2.0 -- The Matrix Equations Sparse
Solvers library},
howpublished = {DOI:10.5281/zenodo.3368844},
month = aug,
year = 2019,
note = {see also:\url{www.mpi-magdeburg.mpg.de/projects/mess}}
}
name: Matrix Equations Sparse Solvers
shortname: M.E.S.S.
version: 1.0.1
release-date: 2016-04-27
doi: 10.5281/zenodo.50575
authors: Jens Saak, Martin Köhler, Peter Benner
orcids: 0000-0001-5567-963, none, 0000-0003-3362-4103
version: 2.0
release-date: 2019-08-27
id: 10.5281/zenodo.3368844
id-type: doi
authors: MESS developer community
copyright holders: Jens Saak, Martin Köhler, Peter Benner
orcids: 0000-0001-5567-963, 0000-0003-2338-9904, 0000-0003-3362-4103
topic: large sparse matrix equations
type: Toolbox
license: GPL v2
......@@ -12,7 +14,7 @@ license-type: Open
repository: private development / public releases
repository-type: git
languages: Matlab
dependencies: GNU Octave >= 4.0, MATLAB >= 2011b
dependencies: GNU Octave >= 4.0, MATLAB >= 2014a
systems: Linux, Windows
website: https://gitlab.mpi-magdeburg.mpg.de/mess/mmess-releases
keywords: symmetric matrix equations, LR-ADI, Newton Kleinman, BDF, Rosenbrock methods
keywords: symmetric matrix equations, LR-ADI, Newton Kleinman, BDF methods, Rosenbrock methods, splitting methods, Riccati iteration, balanced trunation, IRKA
# Core Team
## Scientific Advisor
- Prof. Dr. Peter Benner
## Core Team
## Core Developers
- Dr. Jens Saak
- Martin Köhler
---
# Version 2.0
- Bjoern Baran (BDF methods for non-autonomous DREs, system tests)
- Patrick Kuerschner (RADI)
- Jens Saak (improved MOR functions, test framework, unit and
system tests, code and toolbox restructuring)
- Tony Stillfjord (splitting schemes for Differential Riccati equations)
- Steffen Werner (RADI, improved Operator Manager,
improved Riccati iteration)
---
# Version 1.0 & 1.0.1
## Student Assistants and Interns
- Bjoern Baran (LDL^T based Algorithms and Differential Equations)
- Maximilian Behr (Operator Manager, DAE function handles)
......@@ -15,9 +29,9 @@
- Patrick Kürschner (experimental prototype codes for:
adaptive shifts,
residual factor based algorithms,
non-symmetric equations)
non-symmetric equations, RADI)
- Norman Lang (experimental prototypes codes for:
LDL^T based algorithms,
Differential Lyapunov and Riccati equations)
- Heiko Weichelt (experimental prototype codes for:
inexact Newton with linesearch)
\ No newline at end of file
inexact Newton with linesearch)
function LQR_DAE1(istest)
% Computes a Riccati feedback control for the proper
% index-1 System BIPS98_606 from https://sites.google.com/site/rommes/software
% following the ideas introduced in [1] for Lyapunov equations using he
% Newton-ADI iteration.
%
% Input:
% istest decides whether the function runs as an interactive demo or a
% continuous integration test. (optional; defaults to 0, i.e.
% interactive demo)
%
% References:
%[1] F. Freitas, J. Rommes, N. Martins, Gramian-based reduction method
% applied to large sparse power system descriptor models, IEEE Trans.
% Power Syst. 23 (3) (2008) 1258–1270. doi:10.1109/TPWRS.2008.926693.
%
%
......@@ -15,65 +31,127 @@
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016
% 2009-2019
%
clear , close all
%%
if nargin<1, istest=0; end
%%
% set operation
oper = operatormanager('dae_1');
%% Problem data
load('bips98_606.mat')
fname = sprintf('%s/../models/BIPS/bips98_606.mat',...
fileparts(mfilename('fullpath')));
Bips = load(fname);
% from https://sites.google.com/site/rommes/software
p = find(diag(E));
np = find(diag(E) == 0);
p = find(diag(Bips.E));
np = find(diag(Bips.E) == 0);
pp = [p;np];
eqn.A_ = A(pp, pp);
eqn.E_ = E(pp, pp);
eqn.B = b(pp, :);
eqn.C = c( : , pp);
eqn.A_ = Bips.A(pp, pp);
eqn.E_ = Bips.E(pp, pp);
eqn.B = Bips.b(pp, :);
eqn.C = 0.01*Bips.c( : , pp);
eqn.st = length(p);
eqn.haveE = 1;
clear Bips;
%% Turn off close to singular warnings
% (this model is really badly conditioned)
orig_warnstate = warning('OFF','MATLAB:nearlySingularMatrix');
%%
opts.norm = 'fro';
% ADI tolerances and maximum iteration number
opts.adi.maxiter = 300;
opts.adi.restol = 1e-12;
opts.adi.rctol = 1e-16;
opts.adi.info = 1;
opts.adi.res_tol = 1e-12;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 0;
opts.adi.projection.freq=0;
eqn.type='N';
%%
% n=size(eqn.A_, 1);
opts.adi.shifts.l0=25;
opts.adi.shifts.kp=50;
opts.adi.shifts.km=25;
opts.adi.shifts.method = 'projection';
opts.adi.shifts.l0= 9;
% opts.adi.shifts.b0=ones(n,1);
opts.shifts.num_desired=25;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.method = 'projection';
opts.shifts.num_desired= 9;
%%
% Newton tolerances and maximum iteration number
opts.nm.maxiter = 30;
opts.nm.restol = 1e-10;
opts.nm.rctol = 1e-16;
opts.nm.res_tol = 1e-10;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 1;
opts.nm.linesearch = 1;
opts.nm.accumulateRes = 1;
opts.nm.norm = 'fro';
%%
tic
[ZB,out]=mess_lrnm(eqn,opts,oper);
toc
figure(1)
out.nm.res
semilogy(out.nm.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations');
ylabel('normalized residual norm');
pause(1)
disp('size ZB:')
size(ZB)
tic;
outnm = mess_lrnm(eqn, opts, oper);
toc;
if istest
if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result in LRNM');
end
else
figure(1);
semilogy(outnm.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations');
ylabel('normalized residual norm');
pause(1);
end
disp('size outnm.Z:');
disp(size(outnm.Z));
%% Lets try RADI
opts.norm = 2;
% RADI-MESS settings
opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1);
opts.shifts.method = 'gen-ham-opti';
opts.shifts.naive_update_mode = false;
% .. Suggest false (smart update is faster; convergence is the same).
opts.radi.compute_sol_fac = 1;
opts.radi.get_ZZt = 1;
opts.radi.compute_res = 0;
opts.radi.maxiter = 500;
opts.radi.res_tol = opts.nm.res_tol;
opts.radi.rel_diff_tol = 0;
opts.radi.info = 1;
tic;
outradi = mess_lrradi( eqn, opts, oper );
toc;
if istest
if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result in RADI');
end
else
figure(2);
semilogy(outradi.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations');
ylabel('normalized residual norm');
end
disp('size outradi.Z:');
disp(size(outradi.Z));
%% compare
if not(istest)
figure(3);
ls_nm=[outnm.adi.niter];
ls_radi=1:outradi.niter;
semilogy(cumsum(ls_nm),outnm.res,'k--',ls_radi,outradi.res,'b-');
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of solves with A+p*M');
ylabel('normalized residual norm');
legend('LR-NM','RADI');
end
%% reset warning state
warning(orig_warnstate');
function bt_mor_DAE1_tol(istest)
% Computes a reduced order model via Balanced Truncation for the proper
% index-1 System BIPS98_606 from https://sites.google.com/site/rommes/software
% following the method suggested in [1]
%
% Input:
% istest decides whether the function runs as an interactive demo or a
% continuous integration test. (optional; defaults to 0, i.e.
% interactive demo)
%
% References:
%[1] F. Freitas, J. Rommes, N. Martins, Gramian-based reduction method
% applied to large sparse power system descriptor models, IEEE Trans.
% Power Syst. 23 (3) (2008) 1258–1270. doi:10.1109/TPWRS.2008.926693.
%
% 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
......@@ -12,165 +28,151 @@
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016
% 2009-2019
%
clear , close all
%%
if nargin<1, istest=0; end
%%
% set operation
% set operation manager for the Gramian computations
oper = operatormanager('dae_1');
%% Problem data
load('bips98_606.mat')
%% Read problem data
fname = sprintf('%s/../models/BIPS/bips98_606.mat',...
fileparts(mfilename('fullpath')));
Bips = load(fname);
% from https://sites.google.com/site/rommes/software
% Note that we here use the alpha shift suggested by Rommes and coauthors,
% but do not perform the shift back.
p = find(diag(E));
np = find(diag(E) == 0);
% Note that we here use the alpha shift suggested by Rommes and coauthors
% [1], but do not perform the shift back.
p = find(diag(Bips.E));
np = find(diag(Bips.E) == 0);
pp = [p;np];
eqn.A_ = A(pp, pp)-0.05*E(pp,pp);
eqn.E_ = E(pp, pp);
eqn.B = b(pp, :);
eqn.C = c( : , pp);
eqn.A_ = Bips.A(pp, pp)-0.05*Bips.E(pp,pp);
eqn.E_ = Bips.E(pp, pp);
eqn.B = Bips.b(pp, :);
eqn.C = Bips.c( : , pp);
eqn.st = length(p);
eqn.haveE = 1;
%%
% BT tolerance and maximum order for the ROM
tol=1e-3;
max_ord=250;
%% Turn off close to singular warnings
% (this model is really badly conditioned)
orig_warnstate = warning('OFF','MATLAB:nearlySingularMatrix');
%%
% ADI tolerances and maximum iteration number
opts.adi.maxiter = 300;
opts.adi.restol = 1e-12;
opts.adi.rctol = 1e-16;
opts.adi.res_tol = 1e-10;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 1;
opts.adi.projection.freq=0;
opts.adi.norm = 'fro';
opts.norm = 'fro';
%%
opts.adi.shifts.kp = 50;
opts.adi.shifts.km = 35;
opts.adi.shifts.method = 'projection';
opts.adi.shifts.l0= 20;
opts.shifts.method = 'projection';
opts.shifts.num_desired= 20;
opts.adi.shifts.p=mess_para(eqn,opts,oper);
opts.shifts.p=mess_para(eqn,opts,oper);
disp(opts.adi.shifts.p);
disp(opts.shifts.p);
%%
eqn.type='N';
tic
[ZB,out]=mess_lradi(eqn,opts,oper);
toc
figure
disp('normalize residuals')
semilogy(out.res);
title('0= BB^T + AXM^T + MXA^T');
xlabel('number of iterations');
ylabel('normalized residual norm');
tic;
outB = mess_lradi(eqn, opts, oper);
toc;
if istest
if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
end
else
figure;
semilogy(outB.res);
title('0= BB^T + AXM^T + MXA^T');
xlabel('number of iterations');
ylabel('normalized residual norm');
end
disp('size ZB:')
size(ZB)
disp('size outB.Z:');
disp(size(outB.Z));
%%
eqn.type='T';
tic
[ZC,out]=mess_lradi(eqn,opts,oper);
toc
figure
disp('normalize residuals')
semilogy(out.res);
title('0= C^TC + A^TXE + E^TXA');
xlabel('number of iterations');
ylabel('normalized residual norm');
pause(1)
disp('size ZC:')
size(ZC)
%%
[U0,S0,V0] = svd(ZC'*(eqn.E_(1:eqn.st,1:eqn.st)*ZB),0);
s0=diag(S0);
ks=length(s0);
k=ks;
while (sum(s0(k-1:ks))<tol/2)&&(k>2)
k=k-1;
tic;
outC = mess_lradi(eqn, opts, oper);
toc;
if istest
if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
end
else
figure;
semilogy(outC.res);
title('0= C^TC + A^TXE + E^TXA');
xlabel('number of iterations');
ylabel('normalized residual norm');
pause(1);
end
k0=k
r= min([max_ord k0]);
fprintf(1,'reduced system order: %d\n\n',r);
sigma_0=diag(S0);
sigma_r=diag(S0(1:r,1:r));
disp('size outC.Z:');
disp(size(outC.Z));
VB = ZB*V0(:,1:r);
VC = ZC*U0(:,1:r);
%% Compute reduced system matrices
% Perform Square Root Method (SRM)
SB = VB*diag(ones(r,1)./sqrt(sigma_r));
SC = VC*diag(ones(r,1)./sqrt(sigma_r));
b1 = SC'*(eqn.A_(1:eqn.st,1:eqn.st))*SB;
b2 = SC'*(eqn.A_(1:eqn.st,eqn.st+1:end));
a1 = eqn.A_(eqn.st+1:end,1:eqn.st)*SB;
% BT tolerance and maximum order for the ROM
opts.srm.tol=1e-3;
opts.srm.max_ord=250;
% SRM verbosity
if istest
opts.srm.info=1;
else
opts.srm.info=2;
end
Ar = b1 - b2*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\a1);%+a1*eye(r,r);
Br = SC'*eqn.B(1:eqn.st,:) - b2*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\eqn.B(eqn.st+1:end,:));
Cr = eqn.C(:,1:eqn.st)*SB - eqn.C(:,eqn.st+1:end)*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\a1);
Dr = -eqn.C(:,eqn.st+1:end)*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\eqn.B(eqn.st+1:end,:));
%%
nsample = 200;
w = logspace(-3, 4, nsample);
tr1=zeros(1,nsample); tr2=tr1; err=tr1; relerr=tr1;
fprintf(['Computing TFMs of original and reduced order systems and ' ...
'MOR errors\n'])
Ir=eye(r);
for k=1:nsample
if mod(k,10)==0, fprintf('\r Step %3d / %3d',k,nsample), end
g1 = eqn.C / (1i*w(k)*eqn.E_ - eqn.A_) * eqn.B;
g2 = Cr / (1i*w(k)*Ir - Ar) * Br + Dr;
err(k) = max(svds(g1-g2));
tr1(k) = max(svds(g1));
tr2(k) = max(svds(g2));
relerr(k)=err(k)/tr1(k);
%The actual SRM
[TL,TR,hsv] = mess_square_root_method(eqn,opts,oper,outB.Z,outC.Z);
% compute ROM matrices
B1 = TL'*(eqn.A_(1:eqn.st,1:eqn.st))*TR;
B2 = TL'*(eqn.A_(1:eqn.st,eqn.st+1:end));
A1 = eqn.A_(eqn.st+1:end,1:eqn.st)*TR;
ROM.A = B1 - B2*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\A1);
ROM.B = TL'*eqn.B(1:eqn.st,:) - ...
B2*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\eqn.B(eqn.st+1:end,:));
ROM.C = eqn.C(:,1:eqn.st)*TR - ...
eqn.C(:,eqn.st+1:end)*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\A1);
ROM.D = -eqn.C(:,eqn.st+1:end)*(eqn.A_(eqn.st+1:end,eqn.st+1:end)\...
eqn.B(eqn.st+1:end,:));
ROM.E = eye(size(ROM.A));
%% Evaluate the ROM quality
% while the Gramians are computed on the hidden manifold, we need to do the
% frequency domain computations without (implicitly) using the Schur
% complement (due to the construction of the function handles)
oper = operatormanager('default');
if istest
opts.sigma.info=0;
else
opts.sigma.info=2;
end
clear I
fprintf('\n\n');
figure
subplot(2,1,1);
loglog(w, err);
title('absolute model reduction error')
xlabel('\omega')
ylabel('\sigma_{max}(G(j\omega) - G_r(j\omega))')
subplot(2,1,2);
loglog(w, relerr);
title('relative model reduction error')
xlabel('\omega')
ylabel(['\sigma_{max}(G(j\omega) - G_r(j\omega)) / \' ...
'sigma_{max}(G(j\omega))'])
figure
loglog(w, tr1)
hold on
loglog(w, tr2, 'r--')
legend({'original system','reduced system'})
xlabel('\omega')
ylabel('\sigma_{max}(G(j\omega))')
title('Transfer functions of original and reduced systems')
hold off
figure
semilogy(diag(S0));
title('Computed Hankel singular values');
xlabel('index');
ylabel('magnitude');
opts.sigma.fmin=-3;
opts.sigma.fmax=4;
err = mess_sigma_plot(eqn, opts, oper, ROM);
if istest
if max(err)>5e-3
error('MESS:TEST:accuracy','unexpectedly innacurate result');
end
else
figure;
semilogy(hsv);
title('Computed Hankel singular values');
xlabel('index');
ylabel('magnitude');
end
%% reset warning state
warning(orig_warnstate);
\ No newline at end of file
function LQR_DAE2(problem,lvl,re,istest)
% Computes a stabilizing feedback by implicitly solving the generalized
% Riccati equation for the equivalent projected system on the hidden manifold.
%
% Inputs:
% problem either 'Stokes' or 'NSE' to choose the Stokes demo or the
% linearized Navier-Stokes-Equation. (required)
%
% lvl discretization level 1 through 5
% (optional, only used in 'NSE' case, default: 1)
%
% re Reynolds number 300, 400, or 500
% (optional, only used in 'NSE' case, default: 500)
%
% istest flag to determine whether this demo runs as a CI test or
% interactive demo
% (optional, defaults to 0, i.e. interactive demo)
%
% Note that the 'NSE' option requires additional data available in a
% separate 270MB archive and at least the 5th discretization level needs a
% considerable amount of main memory installed in your machine.
%
% This program is free software; you can redistribute it and/or modify
......@@ -15,21 +35,19 @@
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016
% 2009-2019
%
clear all, close all
%%
% set operation
oper = operatormanager('dae_2');
%%
% Problem data
problem='stokes'
% problem='NSE'
% lvl=1
% re=300
if nargin<2, lvl=1; end
if nargin<3, re=500; end
if nargin<4, istest=0; end
switch problem
switch lower(problem)
case 'stokes'
nin = 5;
nout = 5;
......@@ -41,49 +59,51 @@ switch problem
eqn.st=st;
eqn.B=eqn.B(1:st,:);
eqn.C=eqn.C(:,1:st);
case 'NSE'
case 'nse'
try
load(sprintf('mat_nse_re_%d',re))
load(sprintf('%s/../models/NSE/mat_nse_re_%d',...
fileparts(mfilename('fullpath')),re),'mat');
catch
error('The files mat_nse_re_300.mat, mat_nse_re_400.mat and mat_nse_re_500.mat ar available for dowload in a separate archive (270MB each). Please fetch them from the MESS download mage and unpack them into the DEMOS/models/NSE folder.')
error(['The files mat_nse_re_300.mat, mat_nse_re_400.mat and ', ...
'mat_nse_re_500.mat are available for dowload in a ', ...
'separate archive (270MB each). Please fetch them from the ', ...
'MESS download page and unpack them into the ', ...
'DEMOS/models/NSE folder.']);
end
eqn.A_=mat.mat_v.fullA{lvl};
eqn.E_=mat.mat_v.E{lvl};
eqn.haveE=1;
eqn.B=mat.mat_v.B{lvl};
if re>200
eqn.K0=mat.mat_v.Feed_0{lvl};
eqn.haveUV=1;
opts.nm.K0=mat.mat_v.Feed_0{lvl}';
opts.radi.K0 = opts.nm.K0;
end
eqn.C=mat.mat_v.C{lvl};
eqn.st=mat.mat_mg.nv(lvl);
st=eqn.st;
end
k=2;
%%
% First we run the Newton-ADI Method
opts.norm = 2;
% ADI tolerances and maximum iteration number
opts.adi.maxiter = 300;
opts.adi.restol = 1e-12;
opts.adi.rctol = 1e-16;
opts.adi.res_tol = 1e-12;
opts.adi.rel_diff_tol = 1e-16;
opts.adi.info = 1;
opts.adi.projection.freq=0;
opts.adi.LDL_T=0;
eqn.type='T';
%%
n=size(eqn.A_, 1);
opts.adi.shifts.l0=25;
opts.adi.shifts.kp=50;
opts.adi.shifts.km=25;
opts.adi.shifts.method = 'projection';
opts.adi.shifts.l0= 9;
opts.adi.shifts.b0=ones(n,1);
opts.shifts.num_desired=5;%*nout;
opts.shifts.num_Ritz=50;
opts.shifts.num_hRitz=25;
opts.shifts.method = 'projection';
opts.shifts.b0=ones(n,1);
%%
% Newton tolerances and maximum iteration number
opts.nm.maxiter = 20;
opts.nm.restol = 1e-10;
opts.nm.rctol = 1e-16;
opts.nm.res_tol = 1e-10;
opts.nm.rel_diff_tol = 1e-16;
opts.nm.info = 1;
opts.nm.projection.freq=0;
opts.nm.projection.ortho=1;
......@@ -93,19 +113,69 @@ opts.nm.linesearch = 1;
opts.nm.inexact = 'superlinear';
opts.nm.tau = 0.1;
opts.nm.accumulateRes = 1;
opts.nm.norm = 'fro';
%%
tic
[ZB,out]=mess_lrnm(eqn,opts,oper);
toc
figure(1)
out.nm.res
semilogy(out.nm.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations');
ylabel('normalized residual norm');
pause(1)
%% use low-rank Newton-Kleinman-ADI
tic;
outnm = mess_lrnm(eqn, opts, oper);
toc;
if not(istest)
figure(1);
disp(outnm.res);
semilogy(outnm.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of newton iterations');
ylabel('normalized residual norm');
pause(1);
end
disp('size outnm.Z:');
disp(size(outnm.Z));
%% Lets try the RADI method and compare
opts.norm = 2;
% RADI-MESS settings
opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1);
opts.shifts.num_desired = opts.shifts.num_desired;
% choose either of the three shift methods, here
opts.shifts.method = 'gen-ham-opti';
% opts.shifts.method = 'heur';
% opts.shifts.method = 'projection';
opts.shifts.naive_update_mode = false; % .. Suggest false (smart update is faster; convergence is the same).
opts.shifts.info = 0;
opts.radi.compute_sol_fac = 1; % Turned on for numerical stability reasons.
opts.radi.get_ZZt = 0;
opts.radi.maxiter = opts.adi.maxiter;
opts.radi.res_tol = opts.nm.res_tol;
opts.radi.rel_diff_tol = 0;
opts.radi.info = 1;
disp('size ZB:')
size(ZB)
tic;
outradi = mess_lrradi(eqn, opts, oper);
toc;
if not(istest)
figure(2);
semilogy(outradi.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations');
ylabel('normalized residual norm');
end
%% compare
if istest
if min(outnm.res)>=opts.nm.res_tol, error('MESS:TEST:accuracy','unexpectedly innacurate result'); end
if min(outradi.res)>=opts.radi.res_tol, error('MESS:TEST:accuracy','unexpectedly innacurate result'); end
else
figure(3);
ls_nm=[outnm.adi.niter];
ls_radi=1:outradi.niter;
semilogy(cumsum(ls_nm),outnm.res,'k--',ls_radi,outradi.res,'b-');
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of solves with A+p*M');
ylabel('normalized residual norm');
legend('LR-NM','RADI');
end
This diff is collapsed.
function closed_step(eqn,Ar,Br,Cr,problem,K,Kr,istest)
%
%
......@@ -15,7 +16,7 @@
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016
% 2009-2019
%
x0=zeros(size(eqn.A_,1),1);
xr0=zeros(size(Ar,1),1);
......@@ -28,64 +29,73 @@ T=tmin:tau:tmax;
ntau=floor((tmax-tmin)/tau);
range=1:floor(ntau/500):ntau;
%%
tic
[y,yr] = impeuler_closed(eqn.E_,eqn.A_,eqn.B,eqn.C,eye(size(Ar)),Ar,Br,Cr,K,Kr,tau,tmin,tmax,x0,xr0,alpha);
toc
tic;
[y,yr] = impeuler_closed(eqn.E_,eqn.A_,eqn.B,eqn.C,eye(size(Ar)),Ar,Br,...
Cr,K,Kr,tau,tmin,tmax,x0,xr0,alpha);
toc;
%%
abserr=abs(y-yr);
relerr=abs(abserr./y);
%%
colors=['y','m','c','r','g','b','k'];
figure(20)
hold on
for j=1:size(eqn.C,1)
plot(T(range),y(j,range),colors(j));
plot(T(range),yr(j,range),strcat(colors(j),'--'));
end
xlabel('time')
ylabel('magnitude of outputs')
title('step response')
if strcmp(problem,'NSE')
legend('out1','out1 red','out2','out2 red','out3','out3 red','out4','out4 red','out5','out5 red','out6','out6 red','out7','out7 red','Location','EastOutside')
else
legend('out1','out1 red','out2','out2 red','out3','out3 red','out4','out4 red','out5','out5 red','Location','EastOutside')
end
hold off
figure(20)
%%
figure(21)
for j=1:size(eqn.C,1)
semilogy(T(range),abserr(j,range),colors(j));
if j==1, hold on; end
end
xlabel('time')
ylabel('magnitude')
title('absolute error')
if strcmp(problem,'NSE')
legend('out1','out2','out3','out4','out5','out6','out7','Location','EastOutside')
else
legend('out1','out2','out3','out4','out5','Location','EastOutside')
end
hold off
figure(21)
figure(22)
for j=1:size(eqn.C,1)
semilogy(T(range),relerr(j,range),colors(j));
if j==1, hold on; end
end
xlabel('time')
ylabel('magnitude')
title('relative error')
if strcmp(problem,'NSE')
legend('out1','out2','out3','out4','out5','out6','out7','Location','EastOutside')
if istest
if max(abserr)>=1e-6
error('MESS:TEST:accuracy','unexpectedly innacurate result');
end
else
legend('out1','out2','out3','out4','out5','Location','EastOutside')
end
hold off
figure(22)
colors=['y','m','c','r','g','b','k'];
figure(20);
hold on;
for j=1:size(eqn.C,1)
plot(T(range),y(j,range),colors(j));
plot(T(range),yr(j,range),strcat(colors(j),'--'));
end
xlabel('time');
ylabel('magnitude of outputs');
title('step response');
if strcmp(problem,'NSE')
legend('out1','out1 red','out2','out2 red','out3','out3 red',...
'out4','out4 red','out5','out5 red','out6','out6 red',...
'out7','out7 red','Location','EastOutside');
else
legend('out1','out1 red','out2','out2 red','out3','out3 red',...
'out4','out4 red','out5','out5 red','Location','EastOutside');
end
hold off;
figure(20);
%%
figure(21);
for j=1:size(eqn.C,1)
semilogy(T(range),abserr(j,range),colors(j));
if j==1, hold on; end
end
xlabel('time');
ylabel('magnitude');
title('absolute error');
if strcmp(problem,'NSE')
legend('out1','out2','out3','out4','out5','out6','out7',...
'Location','EastOutside');
else
legend('out1','out2','out3','out4','out5','Location','EastOutside');
end
hold off;
figure(22);
for j=1:size(eqn.C,1)
semilogy(T(range),relerr(j,range),colors(j));
if j==1, hold on; end
end
xlabel('time');
ylabel('magnitude');
title('relative error');
if strcmp(problem,'NSE')
legend('out1','out2','out3','out4','out5','out6','out7',...
'Location','EastOutside');
else
legend('out1','out2','out3','out4','out5','Location','EastOutside');
end
hold off;
end
\ No newline at end of file
%
%
% 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, 2010, 2011, 2012, 2013, 2014, 2015, 2016
%
function [y,yr] = impeuler(E,A,B,C,Er,Ar,Br,Cr,tau,tmin,tmax,x0,xr0,alpha)
% Simple implicit Euler implementation for validation of the DAE2
% MESS open loop example via a basic step response computation
......@@ -36,7 +17,25 @@ function [y,yr] = impeuler(E,A,B,C,Er,Ar,Br,Cr,tau,tmin,tmax,x0,xr0,alpha)
% y,yr outputs of the full and reduced systems in [tmin,tmax]
%
% Author: Jens Saak
%
% 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-2019
%
[L,U,P,Q] = lu(E-tau*A);
[Lr,Ur,Pr] = lu(Er-tau*Ar);
......@@ -44,7 +43,7 @@ function [y,yr] = impeuler(E,A,B,C,Er,Ar,Br,Cr,tau,tmin,tmax,x0,xr0,alpha)
y=zeros(size(C,1),ntau);
yr=y;
for i=1:ntau
if ~mod(i,ceil(ntau/10)),
if not(mod(i,ceil(ntau/10)))
fprintf('\r Implicit Euler step %d / %d',i,ntau);
end
if i<0.1*ntau
......
%
%
% 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, 2010, 2011, 2012, 2013, 2014, 2015, 2016
%
function [y,yr] = impeuler_closed(E,A,B,C,Er,Ar,Br,Cr,K,Kr,tau,tmin,tmax,x0,xr0,alpha)
% Simple implicit Euler implementation for validation of the DAE2
% MESS closed loop example via a basic step response computation
......@@ -37,7 +18,25 @@ function [y,yr] = impeuler_closed(E,A,B,C,Er,Ar,Br,Cr,K,Kr,tau,tmin,tmax,x0,xr0,
% y,yr outputs of the full and reduced systems in [tmin,tmax]
%
% Author: Jens Saak
%
% 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-2019
%
[L,U,P,Q] = lu(E-tau*A);
[Lr,Ur,Pr] = lu(Er-tau*(Ar-Br*Kr));
......@@ -47,7 +46,7 @@ function [y,yr] = impeuler_closed(E,A,B,C,Er,Ar,Br,Cr,K,Kr,tau,tmin,tmax,x0,xr0,
AinvtB=Q*(U\(L\(P*(tau*B))));
tK=(eye(size(K,1))+K*AinvtB)\K;
for i=1:ntau
if ~mod(i,ceil(ntau/10)), fprintf('\r Implicit Euler step %d / %d',i,ntau); end
if not(mod(i,ceil(ntau/10))), fprintf('\r Implicit Euler step %d / %d',i,ntau); end
if i<0.1*ntau
Ex=E*x0;
AinvEx=Q*(U\(L\(P*Ex)));
......
function open_step(eqn,Ar,Br,Cr,problem,istest)
%
%
......@@ -15,7 +16,7 @@
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016
% 2009-2019
%
x0=zeros(size(eqn.A_,1),1);
xr0=zeros(size(Ar,1),1);
......@@ -28,64 +29,74 @@ T=tmin:tau:tmax;
ntau=floor((tmax-tmin)/tau);
range=1:floor(ntau/500):ntau;
%%
tic
[y,yr] = impeuler(eqn.E_,eqn.A_,eqn.B,eqn.C,eye(size(Ar)),Ar,Br,Cr,tau,tmin,tmax,x0,xr0,alpha);
toc
tic;
[y,yr] = impeuler(eqn.E_,eqn.A_,eqn.B,eqn.C,eye(size(Ar)),Ar,Br,Cr,tau,...
tmin,tmax,x0,xr0,alpha);
toc;
%%
abserr=abs(y-yr);
relerr=abs(abserr./y);
%%
colors=['y','m','c','r','g','b','k'];
figure(10)
hold on
for j=1:size(eqn.C,1)
plot(T(range),y(j,range),colors(j));
plot(T(range),yr(j,range),strcat(colors(j),'--'));
end
xlabel('time')
ylabel('magnitude of outputs')
title('step response')
if strcmp(problem,'NSE')
legend('out1','out1 red','out2','out2 red','out3','out3 red','out4','out4 red','out5','out5 red','Location','EastOutside')
%%
if istest
if max(abserr)>=1e-6
error('MESS:TEST:accuracy','unexpectedly innacurate result');
end
else
legend('out1','out1 red','out2','out2 red','out3','out3 red','out4','out4 red','out5','out5 red','Location','EastOutside')
end
hold off
figure(10)
%%
figure(11)
for j=1:size(eqn.C,1)
semilogy(T(range),abserr(j,range),colors(j));
if j==1, hold on; end
end
xlabel('time')