...
 
Commits (2)
## version 2.0.1
### Changed
- many function headers and help texts got improved/completed
### Fixed
- DAE_1 usfs failed for certain systems with non-symmetric A.
- LTV BDF could break in certain situations and was not following the
general naming scheme for some variables.
- mess_res2_norms would break when more than 4 output arguments were requested
## version 2.0
### Added
- 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
- CI testing
- demos serve as system tests
- additional unit tests for the smaller building blocks and backend routines
### Changed
- improved Riccati iteration
- updated minimum required/recommended Matlab and octave versions
- 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)
used. (Note the sign of the update and the transposition in V)
- major updates in the MOR routines
- some resturcturing in the opts structure.
- some restructuring 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`.
`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*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
### Fixed
- several consistency updates and bug fixes
- general code cleaning and pretty printing
## version 1.0.1
### Changed
- updated documentation
- Removed replacements directory since its content was not needed for
Matlab after release 2010b and Octave after 4.0.
### Fixed
- Minor consistency and bug fixes and improved integrity of metafiles.
- CI testing
- demos serve as system tests
- 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
- updated documentation
- Removed replacements directory since its content was not needed for
Matlab after release 2010b and Octave after 4.0.
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
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:
......@@ -56,4 +86,5 @@ Compared to the predecessor lyapack a couple of things have changed.
- 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
- A tangential IRKA implementation for non-DAE systems was added
......@@ -10,17 +10,19 @@ 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 2.0 is
[10.5281/zenodo.3368844](http://doi.org/10.5281/zenodo.3368844)
The DOI for version 2.0.1 is
[10.5281/zenodo.3606345](http://doi.org/10.5281/zenodo.3606345)
##BibTeX
@Misc{SaaKB19-mmess-2.0,
```
@Misc{SaaKB19-mmess-2.0.1,
key = {MMESS},
author = {Saak, J. and K\"{o}hler, M. and Benner, P.},
title = {{M-M.E.S.S.}-2.0 -- The Matrix Equations Sparse
title = {{M-M.E.S.S.}-2.0.1 -- The Matrix Equations Sparse
Solvers library},
howpublished = {DOI:10.5281/zenodo.3368844},
month = aug,
year = 2019,
howpublished = {DOI:10.5281/zenodo.3606345},
month = feb,
year = 2020,
note = {see also:\url{www.mpi-magdeburg.mpg.de/projects/mess}}
}
```
name: Matrix Equations Sparse Solvers
shortname: M.E.S.S.
version: 2.0
release-date: 2019-08-23
id: 10.5281/zenodo.3368844
version: 2.0.1
release-date: 2020-02-??
id: 10.5281/zenodo.3606345
id-type: doi
authors: MESS developer community
copyright holders: Jens Saak, Martin Köhler, Peter Benner
......@@ -15,6 +15,6 @@ repository: private development / public releases
repository-type: git
languages: Matlab
dependencies: GNU Octave >= 4.0, MATLAB >= 2014a
systems: Linux, Windows
systems: Linux, Windows, MacOS
website: https://gitlab.mpi-magdeburg.mpg.de/mess/mmess-releases
keywords: symmetric matrix equations, LR-ADI, Newton Kleinman, BDF methods, Rosenbrock methods, splitting methods, Riccati iteration, balanced trunation, IRKA
......@@ -7,7 +7,12 @@
- Dr. Jens Saak
- Martin Köhler
---
# Version 2.0.1
- Bjoern Baran (DRE method fixes)
- Christian Himpe (code review and documentation fixes)
- Jens Saak (improved MOR functions, partial release automation)
- Steffen Werner (fixed DAE_1 usfs)
# Version 2.0
- Bjoern Baran (BDF methods for non-autonomous DREs, system tests)
- Patrick Kuerschner (RADI)
......@@ -17,7 +22,6 @@
- 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)
......
......@@ -31,7 +31,7 @@ function LQR_DAE1(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
if nargin<1, istest=0; end
......@@ -91,7 +91,7 @@ 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');
error('MESS:TEST:accuracy','unexpectedly inaccurate result in LRNM');
end
else
figure(1);
......@@ -128,7 +128,7 @@ toc;
if istest
if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result in RADI');
error('MESS:TEST:accuracy','unexpectedly inaccurate result in RADI');
end
else
figure(2);
......
......@@ -28,7 +28,7 @@ function bt_mor_DAE1_tol(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
......@@ -81,7 +81,7 @@ toc;
if istest
if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure;
......@@ -103,7 +103,7 @@ toc;
if istest
if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure;
......@@ -161,11 +161,11 @@ end
opts.sigma.fmin=-3;
opts.sigma.fmax=4;
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
if istest
if max(err)>5e-3
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure;
......
function IRKA_mor_Stokes(istest)
% Computes a standard ROM by implicitly running IRKA for the equivalent
% projected system on the hidden manifold.
%
% Inputs:
% istest flag to determine whether this demo runs as a CI test or
% interactive demo
% (optional, defaults to 0, i.e. interactive demo)
%
% 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
%
% IRKA tolerance and maximum iteration number
opts.irka.maxiter = 150;
opts.irka.r = 30;
opts.irka.flipeig = 0;
opts.irka.h2_tol = 1e-6;
opts.irka.init = 'logspace';
if nargin < 1, istest = 0; end
if istest
opts.irka.info = 1;
else
opts.irka.info = 1;
end
oper = operatormanager('dae_2');
%% Problem data
nin = 5;
nout = 5;
nx = 10;
ny = 10;
[eqn.E_, eqn.A_, eqn.Borig, eqn.Corig] = ...
stokes_ind2(nin, nout, nx, ny);
n = size(eqn.E_, 1);
eqn.haveE = 1;
st = full(sum(diag(eqn.E_)));
eqn.st = st;
eqn.B = eqn.Borig(1:st, :);
eqn.C = eqn.Corig(:, 1:st);
%% Compute reduced system matrices
tic;
[ROM.E, ROM.A, ROM.B, ROM.C, S, b, c, V, W] = ...
mess_tangential_irka(eqn, opts, oper);
toc;
%%
tic;
%% Evaluate the ROM quality
% while the Gramians are computed exploiting the DAE structure, due to the
% construction of the function handles we can not do so for the transfer
% function. Therfore we need to extend the matrices B and C and call the
% 'default' usfs for unstructured computation:
eqn.B = eqn.Borig;
eqn.C = eqn.Corig;
oper = operatormanager('default');
if istest
opts.sigma.info = 0;
else
opts.sigma.info = 2;
end
opts.sigma.fmin = -3;
opts.sigma.fmax = 4;
out = mess_sigma_plot(eqn, opts, oper, ROM);
toc;
%%
maerr = max(abs(out.err));
mrerr = max(abs(out.relerr));
if istest && (maerr>1e-6 || mrerr>1e-4)
error('MESS:TEST:accuracy',['unexpectedly inaccurate result.\n' ...
'max. abs err: %e (allowed 1e-6)\n' ...
'max rel err: %e (allowed 1e-4)'], maerr, mrerr);
end
%%
problem = 'Stokes';
fprintf(['\nComputing open loop step response of original and ', ...
'reduced-order systems and time domain MOR errors\n']);
open_step(eqn, ROM.A, ROM.B, ROM.C, problem, istest);
%%
fprintf('\nComputing ROM based feedback\n');
if exist('care', 'file')
[~, ~, Kr] = care(ROM.A, ROM.B, ROM.C'*ROM.C, eye(size(ROM.B, 2)));
else
Y = care_nwt_fac([], ROM.A, ROM.B, ROM.C, 1e-12, 50);
Kr = (Y * ROM.B)' * Y;
end
K = [Kr * W' * eqn.E_(1:st, 1:st), zeros(size(Kr, 1), n-st)];
%%
fprintf(['\nComputing closed loop step response of original and ', ...
'reduced-order systems and time domain MOR errors\n']);
closed_step(eqn, ROM.A, ROM.B, ROM.C, problem, K, Kr, istest);
......@@ -35,7 +35,7 @@ function LQR_DAE2(problem,lvl,re,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%% Set operations
oper = operatormanager('dae_2');
......@@ -79,6 +79,8 @@ switch lower(problem)
end
eqn.C=mat.mat_v.C{lvl};
eqn.st=mat.mat_mg.nv(lvl);
otherwise
error('input ''problem'' must be either ''NSE'' or ''Stokes''');
end
%%
% First we run the Newton-ADI Method
......@@ -156,7 +158,7 @@ outradi = mess_lrradi(eqn, opts, oper);
toc;
if not(istest)
figure(2);
figure();
semilogy(outradi.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations');
......@@ -165,10 +167,10 @@ 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
if min(outnm.res)>=opts.nm.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
if min(outradi.res)>=opts.radi.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
else
figure(3);
figure();
ls_nm=[outnm.adi.niter];
ls_radi=1:outradi.niter;
......
......@@ -24,7 +24,7 @@ function bt_mor_DAE2(problem,lvl,re,istest)
% P. Benner, J. Saak, M. M. Uddin, Balancing based model reduction for
% structured index-2 unstable descriptor systems with application to flow
% control, Numerical Algebra, Control and Optimization 6 (1) (2016) 1–20.
% https://doir.org/10.3934/naco.2016.6.1.
% https://doi.org/10.3934/naco.2016.6.1.
% This program is free software; you can redistribute it and/or modify
......@@ -41,7 +41,7 @@ function bt_mor_DAE2(problem,lvl,re,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
% ADI tolerance and maximum iteration number
......@@ -91,6 +91,8 @@ switch lower(problem)
st=eqn.st;
eqn.haveE=1;
n=size(eqn.E_,1);
otherwise
error('input ''problem'' must be either ''NSE'' or ''Stokes''');
end
%%
eqn.type='N';
......@@ -116,7 +118,7 @@ outB = mess_lradi(eqn,opts,oper);
toc;
if not(istest)
figure(1);
figure();
semilogy(outB.res);
title('AXM^T + MXA^T = -BB^T');
xlabel('number of iterations');
......@@ -156,7 +158,7 @@ outC = mess_lradi(eqn, opts, oper);
toc;
if not(istest)
figure(2);
figure();
semilogy(outC.res);
title('A^TXM + M^TXA = -C^TC');
xlabel('number of iterations');
......@@ -217,12 +219,12 @@ end
opts.sigma.fmin=-3;
opts.sigma.fmax=4;
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
toc;
%%
if istest
if max(err)>=opts.srm.tol, error('MESS:TEST:accuracy','unexpectedly innacurate result'); end
if max(err)>=opts.srm.tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
else
figure;
semilogy(hsv);
......
function closed_step(eqn,Ar,Br,Cr,problem,K,Kr,istest)
% Simple validation of the DAE2 MESS closed loop example via a basic
% step response computation
%
% ope_step(eqn,Ar,Br,Cr,problem, istest)
%
% INPUTS:
% eqn The original system equation structure
% Ar,Br,Cr The reduced order system matrices
% K,Kr full and reduced order feedback matrices
% problem 'NSE' or 'Stokes' switching between the Stokes demo or the
% linearized Navier-Stokes-Equation
% istest flag to determine whether this demo runs as a CI test or
% interactive demo
%
% This program is free software; you can redistribute it and/or modify
......@@ -16,7 +28,7 @@ function closed_step(eqn,Ar,Br,Cr,problem,K,Kr,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
x0=zeros(size(eqn.A_,1),1);
xr0=zeros(size(Ar,1),1);
......@@ -39,13 +51,17 @@ relerr=abs(abserr./y);
%%
if istest
if max(abserr)>=1e-6
error('MESS:TEST:accuracy','unexpectedly innacurate result');
maerr = max(abserr);
if maerr>=1e-6
error('MESS:TEST:accuracy',['unexpectedly inaccurate result ' ...
'in closed loop simulation. Maximum ' ...
'absolute error %e > 1e-6'], maerr);
end
else
colors=['y','m','c','r','g','b','k'];
figure(20);
figure();
hold on;
for j=1:size(eqn.C,1)
plot(T(range),y(j,range),colors(j));
......@@ -63,10 +79,8 @@ else
'out4','out4 red','out5','out5 red','Location','EastOutside');
end
hold off;
figure(20);
%%
figure(21);
figure();
for j=1:size(eqn.C,1)
semilogy(T(range),abserr(j,range),colors(j));
if j==1, hold on; end
......@@ -82,7 +96,7 @@ else
end
hold off;
figure(22);
figure();
for j=1:size(eqn.C,1)
semilogy(T(range),relerr(j,range),colors(j));
......
......@@ -33,7 +33,7 @@ function [y,yr] = impeuler(E,A,B,C,Er,Ar,Br,Cr,tau,tmin,tmax,x0,xr0,alpha)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
[L,U,P,Q] = lu(E-tau*A);
......
......@@ -34,7 +34,7 @@ function [y,yr] = impeuler_closed(E,A,B,C,Er,Ar,Br,Cr,K,Kr,tau,tmin,tmax,x0,xr0,
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
[L,U,P,Q] = lu(E-tau*A);
......
function open_step(eqn,Ar,Br,Cr,problem,istest)
% Simple validation of the DAE2 MESS open loop example via a basic
% step response computation
%
% ope_step(eqn,Ar,Br,Cr,problem, istest)
%
% INPUTS:
% eqn The original system equation structure
% Ar,Br,Cr The reduced order system matrices
% problem 'NSE' or 'Stokes' switching between the Stokes demo or the
% linearized Navier-Stokes-Equation
% istest flag to determine whether this demo runs as a CI test or
% interactive demo
%
% 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
......@@ -16,7 +27,7 @@ function open_step(eqn,Ar,Br,Cr,problem,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
x0=zeros(size(eqn.A_,1),1);
xr0=zeros(size(Ar,1),1);
......@@ -38,8 +49,11 @@ abserr=abs(y-yr);
relerr=abs(abserr./y);
%%
if istest
if max(abserr)>=1e-6
error('MESS:TEST:accuracy','unexpectedly innacurate result');
maerr = max(abserr);
if maerr>=1e-6
error('MESS:TEST:accuracy',['unexpectedly inaccurate result ' ...
'in open loop simulation. Maximum ' ...
'absolute error %e > 1e-6'], maerr);
end
else
......
......@@ -21,7 +21,7 @@ function LQR_DAE2_SO(istest)
% ARGESIM-Reports, Vienna Univ. of Technology, ARGE Simulation News,
% Vienna, Austria, 2009, pp. 1232–1243, iSBN/ISSN: 978-3-901608-35-3.
%
% [3] P. Benner, P. Krschner, J. Saak, Improved second-order balanced
% [3] P. Benner, P. Kürschner, J. Saak, Improved second-order balanced
% truncation for symmetric systems, IFAC Proceedings Volumes (7th
% Vienna International Conference on Mathematical Modelling) 45 (2)
% (2012) 758–762. https://doi.org/10.3182/20120215-3-AT-3016.00134.
......@@ -47,7 +47,7 @@ function LQR_DAE2_SO(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
......@@ -121,7 +121,7 @@ toc;
if istest
if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -161,7 +161,7 @@ toc;
if istest
if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......
......@@ -65,7 +65,7 @@ function BT_DAE3_SO(model, tol, max_ord, maxiter, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%% Input checks
......@@ -265,7 +265,7 @@ eqnu.haveE = 1;
operu = operatormanager('so_1');
err = mess_sigma_plot(eqnu, opts, operu, ROM);
out = mess_sigma_plot(eqnu, opts, operu, ROM); err = out.err;
toc;
%% final accuracy test used in the continuous integration system or
......@@ -274,7 +274,7 @@ if istest
% the errors are not always perfect in this example, but let's see
% wether they are "good enough"...
if (max(err) > 50*tol)
error('MESS:TEST:accuracy', ['unexpectedly innacurate result ' ...
error('MESS:TEST:accuracy', ['unexpectedly inaccurate result ' ...
'for %s %g %d %d (%g)'], model, tol, ...
max_ord, maxiter,max(err));
end
......
......@@ -25,8 +25,7 @@ function LQR_DAE3_SO(model,istest)
% Reduction of Large-Scale Systems, volume 45 of Lecture Notes in
% Computational Science and Engineering, pages 83–115. Springer-Verlag,
% Berlin/Heidelberg, 2005.
%
%
% 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
......@@ -41,7 +40,7 @@ function LQR_DAE3_SO(model,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
narginchk(0,2);
......@@ -111,7 +110,7 @@ toc;
if istest
if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -154,7 +153,7 @@ toc;
if istest
if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......
% Driver script for the results presented in:
%
% J. Saak, M. Voigt, Model reduction of constrained mechanical systems
% in M-M.E.S.S., IFAC-PapersOnLine 9th Vienna International Conference
% on Mathematical Modelling MATHMOD 2018, Vienna, Austria, 21–23
% February 2018 51 (2) (2018) 661–666.
% https://doi.org/10.1016/j.ifacol.2018.03.112.
%
% 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
%
model={{'Stykel_large', 1e-4, 200, 500}, ...
{'Truhar_Veselic', 1e-4, 250, 300}, ...
{'TV2', 1e-3, 300, 300},...
......
......@@ -49,7 +49,7 @@ function [Er,Ar,Br,Cr] = IRKA_FDM(n0,r,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%%
if nargin<1
......
......@@ -33,7 +33,7 @@ function LQR_FDM_unstable(n0, n_unst, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
narginchk(0,3);
if nargin<1, n0 = 20; end
......@@ -136,7 +136,7 @@ for type=['T','N']
%% print output
if istest
if min(outB.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -165,7 +165,7 @@ for type=['T','N']
%% print output
if istest
if min(outB2.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......
......@@ -55,7 +55,7 @@ function [Ar, Br, Cr] = bt_mor_FDM_tol(tol,n0,shifts,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
narginchk(0,4);
......@@ -120,7 +120,7 @@ toc;
if istest
if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -143,7 +143,7 @@ toc;
if istest
if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......@@ -179,11 +179,11 @@ ROM.B = Br;
ROM.C = Cr;
ROM.E = eye(size(ROM.A,1));
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
if istest
if max(err)>tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure;
......
......@@ -123,7 +123,7 @@ toc;
if istest
if min(outC.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -145,7 +145,7 @@ toc;
if istest
if min(outB.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......@@ -182,7 +182,7 @@ ROM.B = Br;
ROM.C = Cr;
ROM.E = eye(size(ROM.A,1));
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
%%
% Report.
......@@ -190,6 +190,6 @@ err = mess_sigma_plot(eqn, opts, oper, ROM);
if istest
if max(err)>=tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
end
......@@ -30,7 +30,7 @@ function out = LQR_LTV_smallscale_BDF(k, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
if nargin < 1
k = 2;
......
......@@ -43,7 +43,7 @@ function out = LQR_LTV_smallscale_splitting(method, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
......
......@@ -3,7 +3,7 @@ the starting point for your own scripts. The directories serve the
following purposes.
**models/**
Contains the actual benchmark model data or functios for their generation
Contains the actual benchmark model data or functions for their generation
**FDM/ Rail/**
Both demonstrate the use of the *"default"* operators, i.e.,
......@@ -32,4 +32,3 @@ following purposes.
**RI/**
Demonstrates the solver for Riccati equations with indefinite quadratic terms.
......@@ -25,7 +25,7 @@ function DEMO_RI_GE_DAE2(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
if nargin<1, istest=0; end
......@@ -181,7 +181,7 @@ fprintf(1, 'RADI -> set tolerance vs. real residual: %e | %e\n', ...
if istest
assert(relerrnm < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result');
'MESS:TEST:accuracy','unexpectedly inaccurate result');
assert(relerrradi < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result');
'MESS:TEST:accuracy','unexpectedly inaccurate result');
end
......@@ -25,7 +25,7 @@ function DEMO_RI_GE_T_N(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
if nargin<1, istest=0; end
......@@ -133,7 +133,7 @@ fprintf(1, 'Filter -> set tolerance vs. real residual: %e | %e\n', ...
if istest
assert(relControl < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result');
'MESS:TEST:accuracy','unexpectedly inaccurate result');
assert(relFilter < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result');
'MESS:TEST:accuracy','unexpectedly inaccurate result');
end
......@@ -25,7 +25,7 @@ function DEMO_RI_T_HYBRID(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
if nargin<1, istest=0; end
......@@ -129,5 +129,5 @@ fprintf(1, '\nset tolerance vs. real residual: %e | %e\n', ...
if istest
assert(relerr < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result');
'MESS:TEST:accuracy','unexpectedly inaccurate result');
end
......@@ -58,7 +58,7 @@ function HINFR_rail(k, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
......@@ -117,7 +117,7 @@ toc;
%% Residual behavior.
if istest
if min(out.res) >= opts.ri.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......
......@@ -55,7 +55,7 @@ function [Er,Ar,Br,Cr] = IRKA_rail(k,r,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%%
if nargin<1
......
......@@ -73,7 +73,7 @@ function LQR_rail(k,shifts,inexact,Galerkin,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
......@@ -148,7 +148,7 @@ toc;
if istest
if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -190,7 +190,7 @@ toc;
if istest
if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......
......@@ -54,7 +54,7 @@ function LQR_rail_BDF(k)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
narginchk(0,1);
......
......@@ -54,7 +54,7 @@ function LQR_rail_Rosenbrock(k)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
narginchk(0,1);
if nargin<1, k=2; end
......
......@@ -64,7 +64,7 @@ function out = LQR_rail_splitting(k, exp_action, method,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
if nargin < 1
......
......@@ -78,7 +78,7 @@ function Lyapunov_rail_LDL_ADI(k,shifts,implicit,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
narginchk(0,4);
......@@ -134,7 +134,7 @@ toc;
if istest
if min(out.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -162,7 +162,7 @@ toc;
if istest
if min(out.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......@@ -176,9 +176,19 @@ disp(size(out1.Z));
%% Difference of Lyapunov solutions
if k<3
err = norm(out.Z * out.Z' - out1.Z * out1.D * out1.Z') / norm(out.Z * out.Z');
fprintf('Relative difference between solution with and without LDL^T: \t %g\n', err);
if err>1e-14
error('MESS:TEST:accuracy','unexpectedly innacurate result');
end
% This is mainly for consistency checking on our continuous
% integration tests.
% NEVER FORM SUCH DYADIC PRODUCTS IN PRODUCTION CODE!!!
err = norm(out.Z * out.Z' - out1.Z * out1.D * out1.Z') /...
norm(out.Z * out.Z');
fprintf(['Relative difference between solution with and without ' ...
'LDL^T: \t %g\n'], err);
if err>1e-12
if implicit
shifts = [shifts '(implicit)'];
end
error('MESS:TEST:accuracy',...
['unexpectedly inaccurate result relative difference',...
' %e > 1e-12 in case %s'],err,shifts);
end
end
\ No newline at end of file
......@@ -67,7 +67,7 @@ function [Ar, Br, Cr] = bt_mor_rail_tol(k,tol,shifts,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
narginchk(0,4);
......@@ -117,7 +117,7 @@ toc;
% residual norm plot
if istest
if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -138,7 +138,7 @@ toc;
% residual norm plot
if istest
if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......@@ -176,11 +176,11 @@ end
opts.sigma.fmin=-3;
opts.sigma.fmax=4;
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
if istest
if max(err)>tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure;
......
......@@ -51,7 +51,7 @@ function BT_TripleChain(version, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
narginchk(0,2)
......@@ -59,7 +59,7 @@ narginchk(0,2)
if nargin==0, version = 'FO'; end
if nargin<2, istest=0; end
format longe;
format long e;
%% set operation
oper = operatormanager('so_1');
......@@ -111,7 +111,7 @@ toc;
if istest
if min(outB.res)>=1e-1
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -134,7 +134,7 @@ toc;
if istest
if min(outC.res)>=1e-1
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......@@ -201,10 +201,10 @@ if istest
else
opts.sigma.info = 2;
end
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
if istest
if max(err) > 1000
error('MESS:TEST:accuracy','unexpectedly innacurate result %g',max(err));
error('MESS:TEST:accuracy','unexpectedly inaccurate result %g',max(err));
end
end
......@@ -52,7 +52,7 @@ function BT_sym_TripleChain(version,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
......@@ -64,7 +64,7 @@ end
if nargin<2
istest=0;
end
format longe;
format long e;
% set operation
oper = operatormanager('so_2');
% Problem data
......@@ -116,7 +116,7 @@ toc;
if istest
if min(outB.res)>=1e-1
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -183,10 +183,10 @@ if istest
else
opts.sigma.info = 2;
end
err = mess_sigma_plot(eqn, opts, oper, ROM);
out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
if istest
if max(err) > 1000
error('MESS:TEST:accuracy','unexpectedly innacurate result %g', max(err));
error('MESS:TEST:accuracy','unexpectedly inaccurate result %g', max(err));
end
end
......@@ -62,7 +62,7 @@ function LQR_TripleChain(n1, usfs, shifts, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%%
narginchk(0,4);
......@@ -154,7 +154,7 @@ toc;
if istest
if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(1);
......@@ -192,7 +192,7 @@ toc;
if istest
if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result');
error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end
else
figure(2);
......@@ -210,7 +210,7 @@ if istest
nrmNM=norm(outnm.K,'fro');
if nrm/nrmNM >= 1e-9
error('MESS:TEST:accuracy',...
'unexpectedly innacurate result: ||K_NM - K_RADI||_F / ||K_NM||_F=%g',nrm/nrmNM);
'unexpectedly inaccurate result: ||K_NM - K_RADI||_F / ||K_NM||_F=%g',nrm/nrmNM);
end
else
figure(3);
......
......@@ -34,7 +34,7 @@ function eqn = getrail(k)
% along with this program; if not, see <http://www.gnu.org/licenses/>.
%
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019
% 2009-2020
%
%% set path
......
% Simple wrapper script to generate the version of the triple chain
% oscillator model from
% [1] N.Truhar and K.Veselic
% An efficient method for estimating the optimal dampers' viscosity for
% linear vibrating systems using Lyapunov equations
% SIAM J. Matrix Anal. Appl. vol.31 no.1 pp 18-39
%
% with the parametrization used in
%
% [2] J. Saak, Efficient numerical solution of large scale algebraic matrix equations
% in PDE control and model order reduction, Dissertation,
% Technische Universität Chemnitz, Chemnitz, Germany (Jul. 2009).
% URL http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642
%
% 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
%
%%
n1=150;
alpha=0.01;
......@@ -5,7 +36,7 @@ beta=alpha;
v=5e0;
%%
[M,D,K]=triplechain_MSD(n1,alpha,beta,v);
[M,E,K]=triplechain_MSD(n1,alpha,beta,v);
B=ones(3*n1+1,1);
Cp=B';
Cv=zeros(size(Cp));
......
function [M,D,K]=triplechain_MSD(n1,alpha,beta,v)
% function [M,D,K]=triplechain_MSD(n1,alpha,beta,v)
function [M,E,K]=triplechain_MSD(n1,alpha,beta,v)
% function [M,E,K]=triplechain_MSD(n1,alpha,beta,v)
%
% Generates mass spring damper system of three coupled mass sprong damper
% chains as in [1,example 2] with proportional damping. The resulting
......@@ -7,7 +7,7 @@ function [M,D,K]=triplechain_MSD(n1,alpha,beta,v)
%
% Output:
%
% M,D,K mass , damping and stiffness matrices of the system
% M,E,K mass , damping and stiffness matrices of the system
%
% Input:
%
......@@ -21,6 +21,26 @@ function [M,D,K]=triplechain_MSD(n1,alpha,beta,v)
% linear vibrating systems using Lyapunov equations
% SIAM J. Matrix Anal. Appl. vol.31 no.1 pp 18-39
%
%
% 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
%
m1=1;
m2=2;
m3=3;
......@@ -59,9 +79,10 @@ K(3*n1+1, 3*n1+1)=k1+k2+k3+k0;
%fractional damping is used in [1] but needs full matrices due to the sqrtm
%SM=spdiags([sqrt(m1)*ones(1,n1) sqrt(m2)*ones(1,n1) sqrt(m3)*ones(1,n1)...
% sqrt(m0)]',0,3*n1+1,3*n1+1);
%D=.02*(2*SM*sqrtm((SM\K)/SM)*SM);
%E=.02*(2*SM*sqrtm((SM\K)/SM)*SM);
D=alpha*M+beta*K;
D(1,1)=D(1,1)+v;
D(n1,n1)=D(n1,n1)+v;
D(2*n1+1,2*n1+1)=D(2*n1+1,2*n1+1)+v;
\ No newline at end of file
% here we want internal damping based on sparse matrices via Rayleigh damping:
E=alpha*M+beta*K;
E(1,1)=E(1,1)+v;
E(n1,n1)=E(n1,n1)+v;
E(2*n1+1,2*n1+1)=E(2*n1+1,2*n1+1)+v;
\ No newline at end of file
Dependencies of M-M.E.S.S. 2.0
============================================================================
==============================
Matlab R2014a and above, or Octave 4.0 and above.
Matlab R2014a and above, or Octave 4.0 and above.
Some functions can benefit from the Control Systems Toolbox in either
of the above, but fallbacks exist in case it is not available.
Some functions can benefit from the Control Systems Toolbox in Matlab or the
Control package in Octave, but fallbacks exist in case those are not available.
Matlab R2017a has improved handling negative definite matrices with "backslash".
We recommend using this version or later ones for optimal performance.
Matlab R2017a has improved handling of negative definite matrices with in the
"backslash" operator. We recommend using this version or later ones for optimal
performance.
Note that Matlab R2017a R2017b contain a bug in "backslash" that can cause
extraordinarily slow computations with certain block structured
matrices. Use spparms('usema57', 0); to fix this, or upgrade to at
least R2017b update 5.
Note that Matlab R2017a and R2017b also contain a bug in "backslash" that can
cause extraordinarily slow computations with certain block structured matrices.
Use
```
spparms('usema57', 0);
```
to fix this, or upgrade to at least R2017b update 5.
......@@ -15,8 +15,23 @@ Riccati equations.
### ZIP-archive based downloads
Simply unpack the content of the archive to some folder, change to