Commit 6857de14 authored by Tim Mitchell's avatar Tim Mitchell
Browse files

ROSTAPACK v2.2

parent 607913e7
......@@ -2,6 +2,31 @@
ROSTAPACK: RObust STAbility PACKage
## Version: 2.2 --- 2019-09-06
### Added
- specValSet now uses Hessenberg factorizations by default to make evaluations
of the norm of transfer function much faster. This results in specValSet
being noticeably faster as well.
### Changed
- specValSet: E must be an invertible matrix
- specValSet: better tolerance for getting imaginary eigenvalues when using eig
- specValSetOptions: opts.use_permutations has been replaced by opts.solve_type
- specValSetOptions: opts.ignore_infinite has been removed
### Documentation
- Various documentation typos fixed
### Maintenance
- Some minor code cleanup and refactoring of private subroutines
## Version: 2.1 --- 2019-04-24
### Fixed
......@@ -18,7 +43,7 @@ ROSTAPACK: RObust STAbility PACKage
### Added
- New specValSet routine for exact computation of the spectral value set
abscissa and radius (or the pseudospectal abscissa and radius)
abscissa and radius (or the pseudospectral abscissa and radius)
## Version: 1.0.1 --- 2019-02-13
......
ROSTAPACK: RObust STAbility PACKage
Version 2.1
Version 2.2
Licensed under the AGPLv3, Copyright (C) 2014-2019 Tim Mitchell
ROSTAPACK (pronounced rost-a-pack, with rost rhyming with cost, frost, lost)
......@@ -13,7 +13,7 @@ These measures include:
- the spectral value set abscissa and radius
- the pseudospectral abscissa and radius.
Version 2.1 includes both exact and approximation methods for the spectral
Version 2.2 includes both exact and approximation methods for the spectral
value set (or pseudospectral) abscissa and radius measures. It also includes
scalable routines for approximating the H-infinity norm, complex stability
radius and real stability radius. Corresponding fast exact methods for
......
No preview for this file type
name: ROSTAPACK: RObust STAbility PACKage
shortname: ROSTAPACK
version: 2.1
release-date: 2019-04-24
version: 2.2
release-date: 2019-09-06
authors: Tim Mitchell
topic: Numerical Linear Algebra
type: Software Package / Toolbox
......
function [eta,loc,info] = specValSet(A,varargin)
% specValSet:
% Given a linear dynamical system, in state-space matrix form,
%
% E xdot = Ax + Bu
% = Cx + Du,
%
% specValSet computes the epsilon spectral value set abscissa or
% radius. When B=C=E=I and D=0, with n=m=p, then specValSet computes
% the epsilon pseudospectral abscissa or radius.
% the epsilon pseudospectral abscissa or radius. E should be
% invertible.
%
% Computationally, specValSet requires O(n^3) work and O(n^2) memory
% per iteration. Every iteration, a dense matrix pencil of order 2n
......@@ -35,15 +40,17 @@ function [eta,loc,info] = specValSet(A,varargin)
% indeed sparse and are provided in sparse format, there can be
% additional efficiency benefits in the horizontal/radial
% searches. In the general case, the fast searches require
% applying (zE-A)^{-1} for complex scalars z, which is done via
% LUs and then doing backsolves. When B=C=I, D=0, and n=m=p, an
% alterate more efficient formulation is used by default, where
% only zE-A is required (and not its inverse), thus avoiding the
% LUs and backsolves. However, as this alterate form involves
% computing the smallest singular value, numerical accuracy can
% be adversely affected for some very poorly scaled matrices. In
% this case, one can set opts.force_two_norm to true to force the
% the computation involving (zE-A)^{-1}.
% applying (zE-A)^{-1} for complex scalars z, which by default is
% done via Hessenberg factorizations. If A and E are large but
% then this can optionally be done with LU factorizations. When
% B=C=I, D=0, and n=m=p, an alterate more efficient formulation
% is used by default, where only zE-A is required (and not its
% inverse), thus avoiding any factorizations and backsolves.
% However, as this alterate form involves computing the smallest
% singular value, numerical accuracy can be adversely affected
% for some very poorly scaled matrices. In this case, one can
% set opts.force_two_norm to true to force the the computation
% involving (zE-A)^{-1}.
%
% epsilon [required: real value in [0,1/norm(D))]
% Specifies the perturbation level of the spectral value set.
......@@ -247,10 +254,13 @@ function [eta,loc,info] = specValSet(A,varargin)
[is_real,A,B,C,D,E] = isRealSystem(A,B,C,D,E);
% Compute spectrum to find initial starting eigenvalue.
[d0,~,d_co,d_u] = getEigenvalues(A,B,C,D,E,spec_opts);
[eta,loc] = cplx2EtaLoc(d0);
[d0,~,d_co,d_u,inf_co,inf_u] = getEigenvalues(A,B,C,D,E,spec_opts);
if inf_co || inf_u
error('E must be invertible');
end
[eta,loc] = cplx2EtaLoc(d0);
if is_real
[d0,loc] = ensureInUpperHalf(d0,loc);
[d0,loc] = ensureInUpperHalf(d0,loc);
end
if epsilon == 0 || isinf(eta) || isempty(d0)
......@@ -419,6 +429,9 @@ function [eta,loc,info] = specValSet(A,varargin)
for j = 1:total
l = mps(j);
[e,l] = getExtremalPoint(A,B,C,D,E,epsilon,eta,l,sub_opts);
if is_real
l = abs(l);
end
if e > e_max
e_max = e;
l_max = l;
......@@ -429,6 +442,8 @@ function [eta,loc,info] = specValSet(A,varargin)
% of the given angle, the radial search will not be drawn
% from the midpoint of this arc cross section but from its
% opposite point, i.e. e*exp(1i*(l + pi)).
% However, if the problem has real symmetry, the directions
% will be drawn in the upper half plane.
plotSearch(eta,e,l)
end
end
......@@ -707,7 +722,7 @@ function [eta,loc,info] = specValSet(A,varargin)
% radius case, with a modification to ensure the perturbation is at
% least large enough to make a relative change beyond machine
% precision.
fn = @(e) ntfFn(e,l);
fn = @(e) ntfRootFn(e,l);
[f,df] = fn(e);
delta = newtonStep(f,df);
% Make a tolerance a bit above machine precision (1e-15 for 64 bit)
......@@ -733,9 +748,13 @@ function [eta,loc,info] = specValSet(A,varargin)
end
function opts = getSubOpts()
opts.discrete_time = discrete;
opts.ham_symp_tol = sqrt(eps)*max(norm(toFull(A)),epsilon);
opts.is_symmetric = is_real;
opts.discrete_time = discrete;
nEA = norm(toFull(A));
if ~isempty(E)
nEA = cond(toFull(E))*nEA;
end
opts.ham_symp_tol = sqrt(eps)*max(nEA,epsilon);
opts.is_symmetric = is_real;
end
function plotIntervals(intervals,freqs)
......@@ -764,7 +783,7 @@ end
function o = getSpectrumOptions(opts)
o.discrete_time = opts.discrete_time;
o.ignore_infinite = opts.ignore_infinite;
o.ignore_infinite = false;
o.ignore_unperturbable = opts.ignore_unperturbable;
o.find_unperturbable = opts.ignore_unperturbable || opts.plotting;
end
......@@ -776,7 +795,7 @@ function [d,y] = ensureInUpperHalf(d,y)
if y < 0
y = abs(y);
d = conj(d);
d = conj(d);
end
end
......
......@@ -28,13 +28,6 @@ function opts = specValSetOptions(user_opts,varargin)
% abscissa. If you wish to instead compute the spectral value
% set radius, set .discrete_time to true.
%
% .ignore_infinite [logical | {true}]
% By default, specValSet computes the epsilon spectral value set
% abscissa|radius, ignoring any infinite eigenvalues of (A,E). If
% one wishes the epsilon spectral value set abscissa|radius to be
% computed as infinity whenever any eigenvalue of (A,E) is infinite,
% set this to true.
%
% .ignore_unperturbable [logical | {true}]
% By default, specValSet computes the epsilon spectral value set
% abscissa|radius by initializing at a globably rightmost|outermost
......@@ -135,10 +128,12 @@ function opts = specValSetOptions(user_opts,varargin)
% not be changed but some problems with relatively large m,p, using
% one of these two first-order methods could be faster.
%
% .use_permutations [logical | {false}]
% One can optionally have the LUs done with permutations by setting
% this to true. Whether or not doing so is more efficient will
% depend on the particular problem. A and E should be sparse.
% .solve_type [any value in [0,1,2] | {0}]
% Sets how inverses of zE-A should be applied:
% 0: via upper triangular Hessenberg factorization
% 1: via LU
% 2: via LU with permutations
% 0 is generally recommended for small-scale systems.
%
% .force_two_norm [logical | {false}]
% Only relevant when opts.fast_search is true. If B=C=I, D=0, and
......@@ -231,7 +226,6 @@ function opts = specValSetOptions(user_opts,varargin)
validator.setLogical('discrete_time');
% Starting options
validator.setLogical('ignore_infinite');
validator.setLogical('ignore_unperturbable');
if validator.isSpecified('warm_start')
validator.setDimensioned('warm_start',1,1);
......@@ -243,12 +237,11 @@ function opts = specValSetOptions(user_opts,varargin)
validator.setRealInIntervalOO('tol',0,inf);
% Performance options
%validator.setRealInIntervalCC('auto_tune',0,1);
validator.setLogical('fast_search');
validator.setLogical('vertical_search_first');
validator.setIntegerInRange('bracket_order',1,2);
validator.setRealInIntervalCC('root_order',1,2);
validator.setLogical('use_permutations');
validator.setIntegerInRange('solve_type',0,2);
validator.setLogical('force_two_norm');
validator.setIntegerInRange('random_directions',1,inf);
validator.setRealInIntervalCC('safeguard_width',0,1);
......@@ -266,7 +259,6 @@ end
function default_opts = getDefaults()
default_opts = struct( ...
'discrete_time', false, ...
'ignore_infinite', true, ...
'ignore_unperturbable', true, ...
'warm_start', [], ...
'maxit', 100, ...
......@@ -275,7 +267,7 @@ function default_opts = getDefaults()
'vertical_search_first', false, ...
'bracket_order', 2, ...
'root_order', 2, ...
'use_permutations', false, ...
'solve_type', 0, ...
'force_two_norm', false, ...
'random_directions', 3, ...
'safeguard_width', 0.75, ...
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment