CHANGELOG.md

+## 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
+

CITATION.md

@@ -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}}
 }
+```

CODE

 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

CONTRIBUTORS.md

@@ -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)

DEMOS/DAE1/LQR_DAE1.m

@@ -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);

DEMOS/DAE1/bt_mor_DAE1_tol.m

@@ -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;

DEMOS/DAE2/IRKA_mor_Stokes.m

+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);

DEMOS/DAE2/LQR_DAE2.m

@@ -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;
     

DEMOS/DAE2/bt_mor_DAE2.m

@@ -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);

DEMOS/DAE2/private/closed_step.m

 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));

DEMOS/DAE2/private/impeuler.m

@@ -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);

DEMOS/DAE2/private/impeuler_closed.m

@@ -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);

DEMOS/DAE2/private/open_step.m

 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
 

DEMOS/DAE2_SO/LQR_DAE2_SO.m

@@ -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. Kürschner, 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);

DEMOS/DAE3_SO/BT_DAE3_SO.m

@@ -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

DEMOS/DAE3_SO/LQR_DAE3_SO.m

@@ -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);

DEMOS/DAE3_SO/runme_BT_DAE3_SO.m

+%  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},...

DEMOS/FDM/IRKA_FDM.m

@@ -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

DEMOS/FDM/LQR_FDM_unstable.m

@@ -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);

DEMOS/FDM/bt_mor_FDM_tol.m

@@ -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;

DEMOS/FDM/lqgbt_mor_FDM.m

@@ -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

DEMOS/LTV/LQR_LTV_smallscale_BDF.m

@@ -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;

DEMOS/LTV/LQR_LTV_smallscale_splitting.m

@@ -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
 %
 
 

DEMOS/README.md

@@ -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.
-  

DEMOS/RI/DEMO_RI_GE_DAE2.m

@@ -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

DEMOS/RI/DEMO_RI_GE_T_N.m

@@ -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

DEMOS/RI/DEMO_RI_T_HYBRID.m

@@ -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

DEMOS/Rail/HINFR_rail.m

@@ -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);

DEMOS/Rail/IRKA_rail.m

@@ -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 

DEMOS/Rail/LQR_rail.m

@@ -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);

DEMOS/Rail/LQR_rail_BDF.m

@@ -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);

DEMOS/Rail/LQR_rail_Rosenbrock.m

@@ -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

DEMOS/Rail/LQR_rail_splitting.m

@@ -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

DEMOS/Rail/Lyapunov_rail_LDL_ADI.m

@@ -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

DEMOS/Rail/bt_mor_rail_tol.m

@@ -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;

DEMOS/TripleChain/BT_TripleChain.m

@@ -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    
 

DEMOS/TripleChain/BT_sym_TripleChain.m

@@ -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    
 

DEMOS/TripleChain/LQR_TripleChain.m

@@ -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);

DEMOS/models/Data_Rail/getrail.m

@@ -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

DEMOS/models/TripleChain/example_from_Saak09.m

+% 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));

DEMOS/models/TripleChain/triplechain_MSD.m

-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