...
 
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 ## version 2.0
### Added
- New RADI iteration for AREs - New RADI iteration for AREs
- New splitting methods for autonomous DREs - New splitting methods for autonomous DREs
- New splitting and BDF methods for non-autonomous DREs - New splitting and BDF methods for non-autonomous DREs
- New operator manager only requires non-empty functions and replaces - New operator manager only requires non-empty functions and replaces
non-existent ones with a general `mess_do_nothing` function 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 - improved Riccati iteration
- updated minimum required/recommended Matlab and octave versions - updated minimum required/recommended Matlab and Octave versions
(see `DEPENDENCIES.md`) (see `DEPENDENCIES.md`)
- unified function interfaces for top level calls - unified function interfaces for top level calls
- unified handling of low rank updated operators. Now always A+UV' is - 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 - 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 * `opts.adi.shifts` has moved to `opts.shifts` such that also RADI
can use it independent of ADI can use it independent of ADI
* opts.norm now determines the norm for all methods rather than * opts.norm now determines the norm for all methods rather than
having to consistently specifiy the same norm in each substructure having to consistently specifiy the same norm in each substructure
* initial feedbacks for the Riccati solvers are now stored in the * initial feedbacks for the Riccati solvers are now stored in the
`opts` structure for the method rather than `eqn` `opts` structure for the method rather than `eqn`
- The projection shift routine uses the flag `opts.shifts.implicitVtAV`. - The projection shift routine uses the flag `opts.shifts.implicitVtAV`.
Default is `true`. If set to `false` A*V is computed explicitly. 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 - redesign of the demos
- turned scripts into actual demo functions - turned scripts into actual demo functions
- new demos for indefinite AREs and H-infinity control - 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 - CI testing
- demos serve as system tests - demos serve as system tests
- additional unit tests for the smaller building blocks and backend routines - additional unit tests for the smaller building blocks and backend routines
---
## version 1.0.1 ## version 1.0.1
- Minor consistency and bug fixes and improved integrity of metafiles. - Minor consistency and bug fixes and improved integrity of metafiles.
- updated documentation - updated documentation
- Removed replacements directory since its content was not needed for - 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 ## version 1.0
Compared to the predecessor lyapack a couple of things have changed. Compared to the predecessor lyapack a couple of things have changed.
- The user supplied functions are now managed by an operator manager - The user supplied functions are now managed by an operator manager
- The low rank ADI now has: - The low rank ADI now has:
- optimized treatment of E matrices in generalized equations - optimized treatment of E matrices in generalized equations
- more choices for shift selection, including completely automatic - 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 - improved stopping criteria based on low rank factors of the current residual
- automatic generation of real low rank factors also for complex shifts - automatic generation of real low rank factors also for complex shifts
- The Newton-Kleinman iteration features: - The Newton-Kleinman iteration features:
...@@ -56,4 +86,5 @@ Compared to the predecessor lyapack a couple of things have changed. ...@@ -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 - The SRM routine for balanced truncation is only available for
none-DAE systems. Still, DAE versions are included in the none-DAE systems. Still, DAE versions are included in the
corresponding DEMOS. 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 ...@@ -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. version and a sample BibTeX entry can be found below.
##DOI ##DOI
The DOI for version 2.0 is The DOI for version 2.0.1 is
[10.5281/zenodo.3368844](http://doi.org/10.5281/zenodo.3368844) [10.5281/zenodo.3606345](http://doi.org/10.5281/zenodo.3606345)
##BibTeX ##BibTeX
@Misc{SaaKB19-mmess-2.0, ```
@Misc{SaaKB19-mmess-2.0.1,
key = {MMESS}, key = {MMESS},
author = {Saak, J. and K\"{o}hler, M. and Benner, P.}, 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}, Solvers library},
howpublished = {DOI:10.5281/zenodo.3368844}, howpublished = {DOI:10.5281/zenodo.3606345},
month = aug, month = feb,
year = 2019, year = 2020,
note = {see also:\url{www.mpi-magdeburg.mpg.de/projects/mess}} note = {see also:\url{www.mpi-magdeburg.mpg.de/projects/mess}}
} }
```
name: Matrix Equations Sparse Solvers name: Matrix Equations Sparse Solvers
shortname: M.E.S.S. shortname: M.E.S.S.
version: 2.0 version: 2.0.1
release-date: 2019-08-23 release-date: 2020-02-??
id: 10.5281/zenodo.3368844 id: 10.5281/zenodo.3606345
id-type: doi id-type: doi
authors: MESS developer community authors: MESS developer community
copyright holders: Jens Saak, Martin Köhler, Peter Benner copyright holders: Jens Saak, Martin Köhler, Peter Benner
...@@ -15,6 +15,6 @@ repository: private development / public releases ...@@ -15,6 +15,6 @@ repository: private development / public releases
repository-type: git repository-type: git
languages: Matlab languages: Matlab
dependencies: GNU Octave >= 4.0, MATLAB >= 2014a dependencies: GNU Octave >= 4.0, MATLAB >= 2014a
systems: Linux, Windows systems: Linux, Windows, MacOS
website: https://gitlab.mpi-magdeburg.mpg.de/mess/mmess-releases 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 keywords: symmetric matrix equations, LR-ADI, Newton Kleinman, BDF methods, Rosenbrock methods, splitting methods, Riccati iteration, balanced trunation, IRKA
...@@ -7,7 +7,12 @@ ...@@ -7,7 +7,12 @@
- Dr. Jens Saak - Dr. Jens Saak
- Martin Köhler - 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 # Version 2.0
- Bjoern Baran (BDF methods for non-autonomous DREs, system tests) - Bjoern Baran (BDF methods for non-autonomous DREs, system tests)
- Patrick Kuerschner (RADI) - Patrick Kuerschner (RADI)
...@@ -17,7 +22,6 @@ ...@@ -17,7 +22,6 @@
- Steffen Werner (RADI, improved Operator Manager, - Steffen Werner (RADI, improved Operator Manager,
improved Riccati iteration) improved Riccati iteration)
---
# Version 1.0 & 1.0.1 # Version 1.0 & 1.0.1
## Student Assistants and Interns ## Student Assistants and Interns
- Bjoern Baran (LDL^T based Algorithms and Differential Equations) - Bjoern Baran (LDL^T based Algorithms and Differential Equations)
......
...@@ -31,7 +31,7 @@ function LQR_DAE1(istest) ...@@ -31,7 +31,7 @@ function LQR_DAE1(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
if nargin<1, istest=0; end if nargin<1, istest=0; end
...@@ -91,7 +91,7 @@ outnm = mess_lrnm(eqn, opts, oper); ...@@ -91,7 +91,7 @@ outnm = mess_lrnm(eqn, opts, oper);
toc; toc;
if istest if istest
if min(outnm.res)>=opts.nm.res_tol 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 end
else else
figure(1); figure(1);
...@@ -128,7 +128,7 @@ toc; ...@@ -128,7 +128,7 @@ toc;
if istest if istest
if min(outradi.res)>=opts.radi.res_tol 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 end
else else
figure(2); figure(2);
......
...@@ -28,7 +28,7 @@ function bt_mor_DAE1_tol(istest) ...@@ -28,7 +28,7 @@ function bt_mor_DAE1_tol(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
...@@ -81,7 +81,7 @@ toc; ...@@ -81,7 +81,7 @@ toc;
if istest if istest
if min(outB.res)>=opts.adi.res_tol if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure; figure;
...@@ -103,7 +103,7 @@ toc; ...@@ -103,7 +103,7 @@ toc;
if istest if istest
if min(outC.res)>=opts.adi.res_tol if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure; figure;
...@@ -161,11 +161,11 @@ end ...@@ -161,11 +161,11 @@ end
opts.sigma.fmin=-3; opts.sigma.fmin=-3;
opts.sigma.fmax=4; 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 istest
if max(err)>5e-3 if max(err)>5e-3
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure; 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) ...@@ -35,7 +35,7 @@ function LQR_DAE2(problem,lvl,re,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% Set operations %% Set operations
oper = operatormanager('dae_2'); oper = operatormanager('dae_2');
...@@ -79,6 +79,8 @@ switch lower(problem) ...@@ -79,6 +79,8 @@ switch lower(problem)
end end
eqn.C=mat.mat_v.C{lvl}; eqn.C=mat.mat_v.C{lvl};
eqn.st=mat.mat_mg.nv(lvl); eqn.st=mat.mat_mg.nv(lvl);
otherwise
error('input ''problem'' must be either ''NSE'' or ''Stokes''');
end end
%% %%
% First we run the Newton-ADI Method % First we run the Newton-ADI Method
...@@ -156,7 +158,7 @@ outradi = mess_lrradi(eqn, opts, oper); ...@@ -156,7 +158,7 @@ outradi = mess_lrradi(eqn, opts, oper);
toc; toc;
if not(istest) if not(istest)
figure(2); figure();
semilogy(outradi.res); semilogy(outradi.res);
title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM'); title('0= C^TC + A^TXM + M^TXA -M^TXBB^TXM');
xlabel('number of iterations'); xlabel('number of iterations');
...@@ -165,10 +167,10 @@ end ...@@ -165,10 +167,10 @@ end
%% compare %% compare
if istest if istest
if min(outnm.res)>=opts.nm.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 innacurate result'); end if min(outradi.res)>=opts.radi.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end
else else
figure(3); figure();
ls_nm=[outnm.adi.niter]; ls_nm=[outnm.adi.niter];
ls_radi=1:outradi.niter; ls_radi=1:outradi.niter;
......
...@@ -24,7 +24,7 @@ function bt_mor_DAE2(problem,lvl,re,istest) ...@@ -24,7 +24,7 @@ function bt_mor_DAE2(problem,lvl,re,istest)
% P. Benner, J. Saak, M. M. Uddin, Balancing based model reduction for % P. Benner, J. Saak, M. M. Uddin, Balancing based model reduction for
% structured index-2 unstable descriptor systems with application to flow % structured index-2 unstable descriptor systems with application to flow
% control, Numerical Algebra, Control and Optimization 6 (1) (2016) 1–20. % 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 % This program is free software; you can redistribute it and/or modify
...@@ -41,7 +41,7 @@ function bt_mor_DAE2(problem,lvl,re,istest) ...@@ -41,7 +41,7 @@ function bt_mor_DAE2(problem,lvl,re,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
% ADI tolerance and maximum iteration number % ADI tolerance and maximum iteration number
...@@ -91,6 +91,8 @@ switch lower(problem) ...@@ -91,6 +91,8 @@ switch lower(problem)
st=eqn.st; st=eqn.st;
eqn.haveE=1; eqn.haveE=1;
n=size(eqn.E_,1); n=size(eqn.E_,1);
otherwise
error('input ''problem'' must be either ''NSE'' or ''Stokes''');
end end
%% %%
eqn.type='N'; eqn.type='N';
...@@ -116,7 +118,7 @@ outB = mess_lradi(eqn,opts,oper); ...@@ -116,7 +118,7 @@ outB = mess_lradi(eqn,opts,oper);
toc; toc;
if not(istest) if not(istest)
figure(1); figure();
semilogy(outB.res); semilogy(outB.res);
title('AXM^T + MXA^T = -BB^T'); title('AXM^T + MXA^T = -BB^T');
xlabel('number of iterations'); xlabel('number of iterations');
...@@ -156,7 +158,7 @@ outC = mess_lradi(eqn, opts, oper); ...@@ -156,7 +158,7 @@ outC = mess_lradi(eqn, opts, oper);
toc; toc;
if not(istest) if not(istest)
figure(2); figure();
semilogy(outC.res); semilogy(outC.res);
title('A^TXM + M^TXA = -C^TC'); title('A^TXM + M^TXA = -C^TC');
xlabel('number of iterations'); xlabel('number of iterations');
...@@ -217,12 +219,12 @@ end ...@@ -217,12 +219,12 @@ end
opts.sigma.fmin=-3; opts.sigma.fmin=-3;
opts.sigma.fmax=4; opts.sigma.fmax=4;
err = mess_sigma_plot(eqn, opts, oper, ROM); out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err;
toc; toc;
%% %%
if istest 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 else
figure; figure;
semilogy(hsv); semilogy(hsv);
......
function closed_step(eqn,Ar,Br,Cr,problem,K,Kr,istest) 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 % 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) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
x0=zeros(size(eqn.A_,1),1); x0=zeros(size(eqn.A_,1),1);
xr0=zeros(size(Ar,1),1); xr0=zeros(size(Ar,1),1);
...@@ -39,13 +51,17 @@ relerr=abs(abserr./y); ...@@ -39,13 +51,17 @@ relerr=abs(abserr./y);
%% %%
if istest if istest
if max(abserr)>=1e-6 maerr = max(abserr);
error('MESS:TEST:accuracy','unexpectedly innacurate result'); if maerr>=1e-6
error('MESS:TEST:accuracy',['unexpectedly inaccurate result ' ...
'in closed loop simulation. Maximum ' ...
'absolute error %e > 1e-6'], maerr);
end end
else else
colors=['y','m','c','r','g','b','k']; colors=['y','m','c','r','g','b','k'];
figure(20); figure();
hold on; hold on;
for j=1:size(eqn.C,1) for j=1:size(eqn.C,1)
plot(T(range),y(j,range),colors(j)); plot(T(range),y(j,range),colors(j));
...@@ -63,10 +79,8 @@ else ...@@ -63,10 +79,8 @@ else
'out4','out4 red','out5','out5 red','Location','EastOutside'); 'out4','out4 red','out5','out5 red','Location','EastOutside');
end end
hold off; hold off;
figure(20);
%% %%
figure(21); figure();
for j=1:size(eqn.C,1) for j=1:size(eqn.C,1)
semilogy(T(range),abserr(j,range),colors(j)); semilogy(T(range),abserr(j,range),colors(j));
if j==1, hold on; end if j==1, hold on; end
...@@ -82,7 +96,7 @@ else ...@@ -82,7 +96,7 @@ else
end end
hold off; hold off;
figure(22); figure();
for j=1:size(eqn.C,1) for j=1:size(eqn.C,1)
semilogy(T(range),relerr(j,range),colors(j)); 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) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
[L,U,P,Q] = lu(E-tau*A); [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, ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
[L,U,P,Q] = lu(E-tau*A); [L,U,P,Q] = lu(E-tau*A);
......
function open_step(eqn,Ar,Br,Cr,problem,istest) 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 % 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 % 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) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
x0=zeros(size(eqn.A_,1),1); x0=zeros(size(eqn.A_,1),1);
xr0=zeros(size(Ar,1),1); xr0=zeros(size(Ar,1),1);
...@@ -38,8 +49,11 @@ abserr=abs(y-yr); ...@@ -38,8 +49,11 @@ abserr=abs(y-yr);
relerr=abs(abserr./y); relerr=abs(abserr./y);
%% %%
if istest if istest
if max(abserr)>=1e-6 maerr = max(abserr);
error('MESS:TEST:accuracy','unexpectedly innacurate result'); if maerr>=1e-6
error('MESS:TEST:accuracy',['unexpectedly inaccurate result ' ...
'in open loop simulation. Maximum ' ...
'absolute error %e > 1e-6'], maerr);
end end
else else
......
...@@ -21,7 +21,7 @@ function LQR_DAE2_SO(istest) ...@@ -21,7 +21,7 @@ function LQR_DAE2_SO(istest)
% ARGESIM-Reports, Vienna Univ. of Technology, ARGE Simulation News, % ARGESIM-Reports, Vienna Univ. of Technology, ARGE Simulation News,
% Vienna, Austria, 2009, pp. 1232–1243, iSBN/ISSN: 978-3-901608-35-3. % 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 % truncation for symmetric systems, IFAC Proceedings Volumes (7th
% Vienna International Conference on Mathematical Modelling) 45 (2) % Vienna International Conference on Mathematical Modelling) 45 (2)
% (2012) 758–762. https://doi.org/10.3182/20120215-3-AT-3016.00134. % (2012) 758–762. https://doi.org/10.3182/20120215-3-AT-3016.00134.
...@@ -47,7 +47,7 @@ function LQR_DAE2_SO(istest) ...@@ -47,7 +47,7 @@ function LQR_DAE2_SO(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
...@@ -121,7 +121,7 @@ toc; ...@@ -121,7 +121,7 @@ toc;
if istest if istest
if min(outnm.res)>=opts.nm.res_tol if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -161,7 +161,7 @@ toc; ...@@ -161,7 +161,7 @@ toc;
if istest if istest
if min(outradi.res)>=opts.radi.res_tol if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
......
...@@ -65,7 +65,7 @@ function BT_DAE3_SO(model, tol, max_ord, maxiter, istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% Input checks %% Input checks
...@@ -265,7 +265,7 @@ eqnu.haveE = 1; ...@@ -265,7 +265,7 @@ eqnu.haveE = 1;
operu = operatormanager('so_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; toc;
%% final accuracy test used in the continuous integration system or %% final accuracy test used in the continuous integration system or
...@@ -274,7 +274,7 @@ if istest ...@@ -274,7 +274,7 @@ if istest
% the errors are not always perfect in this example, but let's see % the errors are not always perfect in this example, but let's see
% wether they are "good enough"... % wether they are "good enough"...
if (max(err) > 50*tol) 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, ... 'for %s %g %d %d (%g)'], model, tol, ...
max_ord, maxiter,max(err)); max_ord, maxiter,max(err));
end end
......
...@@ -25,8 +25,7 @@ function LQR_DAE3_SO(model,istest) ...@@ -25,8 +25,7 @@ function LQR_DAE3_SO(model,istest)
% Reduction of Large-Scale Systems, volume 45 of Lecture Notes in % Reduction of Large-Scale Systems, volume 45 of Lecture Notes in
% Computational Science and Engineering, pages 83–115. Springer-Verlag, % Computational Science and Engineering, pages 83–115. Springer-Verlag,
% Berlin/Heidelberg, 2005. % Berlin/Heidelberg, 2005.
%
%
% This program is free software; you can redistribute it and/or modify % 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 % it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or % the Free Software Foundation; either version 2 of the License, or
...@@ -41,7 +40,7 @@ function LQR_DAE3_SO(model,istest) ...@@ -41,7 +40,7 @@ function LQR_DAE3_SO(model,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
narginchk(0,2); narginchk(0,2);
...@@ -111,7 +110,7 @@ toc; ...@@ -111,7 +110,7 @@ toc;
if istest if istest
if min(outnm.res)>=opts.nm.res_tol if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -154,7 +153,7 @@ toc; ...@@ -154,7 +153,7 @@ toc;
if istest if istest
if min(outradi.res)>=opts.radi.res_tol if min(outradi.res)>=opts.radi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); 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}, ... model={{'Stykel_large', 1e-4, 200, 500}, ...
{'Truhar_Veselic', 1e-4, 250, 300}, ... {'Truhar_Veselic', 1e-4, 250, 300}, ...
{'TV2', 1e-3, 300, 300},... {'TV2', 1e-3, 300, 300},...
......
...@@ -49,7 +49,7 @@ function [Er,Ar,Br,Cr] = IRKA_FDM(n0,r,istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
%% %%
if nargin<1 if nargin<1
......
...@@ -33,7 +33,7 @@ function LQR_FDM_unstable(n0, n_unst, istest) ...@@ -33,7 +33,7 @@ function LQR_FDM_unstable(n0, n_unst, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
narginchk(0,3); narginchk(0,3);
if nargin<1, n0 = 20; end if nargin<1, n0 = 20; end
...@@ -136,7 +136,7 @@ for type=['T','N'] ...@@ -136,7 +136,7 @@ for type=['T','N']
%% print output %% print output
if istest if istest
if min(outB.res)>=opts.nm.res_tol if min(outB.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -165,7 +165,7 @@ for type=['T','N'] ...@@ -165,7 +165,7 @@ for type=['T','N']
%% print output %% print output
if istest if istest
if min(outB2.res)>=opts.nm.res_tol if min(outB2.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
......
...@@ -55,7 +55,7 @@ function [Ar, Br, Cr] = bt_mor_FDM_tol(tol,n0,shifts,istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
narginchk(0,4); narginchk(0,4);
...@@ -120,7 +120,7 @@ toc; ...@@ -120,7 +120,7 @@ toc;
if istest if istest
if min(outB.res)>=opts.adi.res_tol if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -143,7 +143,7 @@ toc; ...@@ -143,7 +143,7 @@ toc;
if istest if istest
if min(outC.res)>=opts.adi.res_tol if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
...@@ -179,11 +179,11 @@ ROM.B = Br; ...@@ -179,11 +179,11 @@ ROM.B = Br;
ROM.C = Cr; ROM.C = Cr;
ROM.E = eye(size(ROM.A,1)); 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 istest
if max(err)>tol if max(err)>tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure; figure;
......
...@@ -123,7 +123,7 @@ toc; ...@@ -123,7 +123,7 @@ toc;
if istest if istest
if min(outC.res)>=opts.nm.res_tol if min(outC.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -145,7 +145,7 @@ toc; ...@@ -145,7 +145,7 @@ toc;
if istest if istest
if min(outB.res)>=opts.nm.res_tol if min(outB.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
...@@ -182,7 +182,7 @@ ROM.B = Br; ...@@ -182,7 +182,7 @@ ROM.B = Br;
ROM.C = Cr; ROM.C = Cr;
ROM.E = eye(size(ROM.A,1)); 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. % Report.
...@@ -190,6 +190,6 @@ err = mess_sigma_plot(eqn, opts, oper, ROM); ...@@ -190,6 +190,6 @@ err = mess_sigma_plot(eqn, opts, oper, ROM);
if istest if istest
if max(err)>=tol if max(err)>=tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
end end
...@@ -30,7 +30,7 @@ function out = LQR_LTV_smallscale_BDF(k, istest) ...@@ -30,7 +30,7 @@ function out = LQR_LTV_smallscale_BDF(k, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
if nargin < 1 if nargin < 1
k = 2; k = 2;
......
...@@ -43,7 +43,7 @@ function out = LQR_LTV_smallscale_splitting(method, istest) ...@@ -43,7 +43,7 @@ function out = LQR_LTV_smallscale_splitting(method, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % 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 ...@@ -3,7 +3,7 @@ the starting point for your own scripts. The directories serve the
following purposes. following purposes.
**models/** **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/** **FDM/ Rail/**
Both demonstrate the use of the *"default"* operators, i.e., Both demonstrate the use of the *"default"* operators, i.e.,
...@@ -32,4 +32,3 @@ following purposes. ...@@ -32,4 +32,3 @@ following purposes.
**RI/** **RI/**
Demonstrates the solver for Riccati equations with indefinite quadratic terms. Demonstrates the solver for Riccati equations with indefinite quadratic terms.
...@@ -25,7 +25,7 @@ function DEMO_RI_GE_DAE2(istest) ...@@ -25,7 +25,7 @@ function DEMO_RI_GE_DAE2(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
if nargin<1, istest=0; end if nargin<1, istest=0; end
...@@ -181,7 +181,7 @@ fprintf(1, 'RADI -> set tolerance vs. real residual: %e | %e\n', ... ...@@ -181,7 +181,7 @@ fprintf(1, 'RADI -> set tolerance vs. real residual: %e | %e\n', ...
if istest if istest
assert(relerrnm < opts.ri.res_tol, ... assert(relerrnm < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result'); 'MESS:TEST:accuracy','unexpectedly inaccurate result');
assert(relerrradi < opts.ri.res_tol, ... assert(relerrradi < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result'); 'MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
...@@ -25,7 +25,7 @@ function DEMO_RI_GE_T_N(istest) ...@@ -25,7 +25,7 @@ function DEMO_RI_GE_T_N(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
if nargin<1, istest=0; end if nargin<1, istest=0; end
...@@ -133,7 +133,7 @@ fprintf(1, 'Filter -> set tolerance vs. real residual: %e | %e\n', ... ...@@ -133,7 +133,7 @@ fprintf(1, 'Filter -> set tolerance vs. real residual: %e | %e\n', ...
if istest if istest
assert(relControl < opts.ri.res_tol, ... assert(relControl < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result'); 'MESS:TEST:accuracy','unexpectedly inaccurate result');
assert(relFilter < opts.ri.res_tol, ... assert(relFilter < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result'); 'MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
...@@ -25,7 +25,7 @@ function DEMO_RI_T_HYBRID(istest) ...@@ -25,7 +25,7 @@ function DEMO_RI_T_HYBRID(istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
if nargin<1, istest=0; end if nargin<1, istest=0; end
...@@ -129,5 +129,5 @@ fprintf(1, '\nset tolerance vs. real residual: %e | %e\n', ... ...@@ -129,5 +129,5 @@ fprintf(1, '\nset tolerance vs. real residual: %e | %e\n', ...
if istest if istest
assert(relerr < opts.ri.res_tol, ... assert(relerr < opts.ri.res_tol, ...
'MESS:TEST:accuracy','unexpectedly innacurate result'); 'MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
...@@ -58,7 +58,7 @@ function HINFR_rail(k, istest) ...@@ -58,7 +58,7 @@ function HINFR_rail(k, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
...@@ -117,7 +117,7 @@ toc; ...@@ -117,7 +117,7 @@ toc;
%% Residual behavior. %% Residual behavior.
if istest if istest
if min(out.res) >= opts.ri.res_tol if min(out.res) >= opts.ri.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
......
...@@ -55,7 +55,7 @@ function [Er,Ar,Br,Cr] = IRKA_rail(k,r,istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
%% %%
if nargin<1 if nargin<1
......
...@@ -73,7 +73,7 @@ function LQR_rail(k,shifts,inexact,Galerkin,istest) ...@@ -73,7 +73,7 @@ function LQR_rail(k,shifts,inexact,Galerkin,istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
...@@ -148,7 +148,7 @@ toc; ...@@ -148,7 +148,7 @@ toc;
if istest if istest
if min(outnm.res)>=opts.nm.res_tol if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -190,7 +190,7 @@ toc; ...@@ -190,7 +190,7 @@ toc;
if istest if istest
if min(outnm.res)>=opts.nm.res_tol if min(outnm.res)>=opts.nm.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
......
...@@ -54,7 +54,7 @@ function LQR_rail_BDF(k) ...@@ -54,7 +54,7 @@ function LQR_rail_BDF(k)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
narginchk(0,1); narginchk(0,1);
......
...@@ -54,7 +54,7 @@ function LQR_rail_Rosenbrock(k) ...@@ -54,7 +54,7 @@ function LQR_rail_Rosenbrock(k)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
narginchk(0,1); narginchk(0,1);
if nargin<1, k=2; end if nargin<1, k=2; end
......
...@@ -64,7 +64,7 @@ function out = LQR_rail_splitting(k, exp_action, method,istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
if nargin < 1 if nargin < 1
......
...@@ -78,7 +78,7 @@ function Lyapunov_rail_LDL_ADI(k,shifts,implicit,istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
%% %%
narginchk(0,4); narginchk(0,4);
...@@ -134,7 +134,7 @@ toc; ...@@ -134,7 +134,7 @@ toc;
if istest if istest
if min(out.res)>=opts.adi.res_tol if min(out.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -162,7 +162,7 @@ toc; ...@@ -162,7 +162,7 @@ toc;
if istest if istest
if min(out.res)>=opts.adi.res_tol if min(out.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
...@@ -176,9 +176,19 @@ disp(size(out1.Z)); ...@@ -176,9 +176,19 @@ disp(size(out1.Z));
%% Difference of Lyapunov solutions %% Difference of Lyapunov solutions
if k<3 if k<3
err = norm(out.Z * out.Z' - out1.Z * out1.D * out1.Z') / norm(out.Z * out.Z'); % This is mainly for consistency checking on our continuous
fprintf('Relative difference between solution with and without LDL^T: \t %g\n', err); % integration tests.
if err>1e-14 % NEVER FORM SUCH DYADIC PRODUCTS IN PRODUCTION CODE!!!
error('MESS:TEST:accuracy','unexpectedly innacurate result'); err = norm(out.Z * out.Z' - out1.Z * out1.D * out1.Z') /...
end 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 end
\ No newline at end of file
...@@ -67,7 +67,7 @@ function [Ar, Br, Cr] = bt_mor_rail_tol(k,tol,shifts,istest) ...@@ -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/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
narginchk(0,4); narginchk(0,4);
...@@ -117,7 +117,7 @@ toc; ...@@ -117,7 +117,7 @@ toc;
% residual norm plot % residual norm plot
if istest if istest
if min(outB.res)>=opts.adi.res_tol if min(outB.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -138,7 +138,7 @@ toc; ...@@ -138,7 +138,7 @@ toc;
% residual norm plot % residual norm plot
if istest if istest
if min(outC.res)>=opts.adi.res_tol if min(outC.res)>=opts.adi.res_tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(2); figure(2);
...@@ -176,11 +176,11 @@ end ...@@ -176,11 +176,11 @@ end
opts.sigma.fmin=-3; opts.sigma.fmin=-3;
opts.sigma.fmax=4; 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 istest
if max(err)>tol if max(err)>tol
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure; figure;
......
...@@ -51,7 +51,7 @@ function BT_TripleChain(version, istest) ...@@ -51,7 +51,7 @@ function BT_TripleChain(version, istest)
% along with this program; if not, see <http://www.gnu.org/licenses/>. % along with this program; if not, see <http://www.gnu.org/licenses/>.
% %
% Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others
% 2009-2019 % 2009-2020
% %
narginchk(0,2) narginchk(0,2)
...@@ -59,7 +59,7 @@ narginchk(0,2) ...@@ -59,7 +59,7 @@ narginchk(0,2)
if nargin==0, version = 'FO'; end if nargin==0, version = 'FO'; end
if nargin<2, istest=0; end if nargin<2, istest=0; end
format longe; format long e;
%% set operation %% set operation
oper = operatormanager('so_1'); oper = operatormanager('so_1');
...@@ -111,7 +111,7 @@ toc; ...@@ -111,7 +111,7 @@ toc;
if istest if istest
if min(outB.res)>=1e-1 if min(outB.res)>=1e-1
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end end
else else
figure(1); figure(1);
...@@ -134,7 +134,7 @@ toc; ...@@ -134,7 +134,7 @@ toc;
if istest if istest
if min(outC.res)>=1e-1 if min(outC.res)>=1e-1
error('MESS:TEST:accuracy','unexpectedly innacurate result'); error('MESS:TEST:accuracy','unexpectedly inaccurate result');
end