Commit 08a7e1c4 authored by Jens Saak's avatar Jens Saak

updates for version 2.0 from dev repository

parent 7083c9a7
## 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));