subspace_pspa: A subspace method for computing the pseudospectral abscissa.


This package contains the MATLAB code used in the paper “Subspace methods for computing the pseudospectral abscissa and the stability radius”.


D. Kressner, EPF Lausanne
B. Vandereycken, EPF Lausanne


The code should work on any reasonably recent version of MATLAB. It was tested with R2011b, R2013b (Mac OS X) and R2010b (Linux).
Although parts of the code is MEX, binary MEX files are provided so there is hopefully no need to compile them. If needed, you can try to compile them using install_mex.


Download and unpack the package.
Run the startup file to adjust your MATLAB path.
Try out example_dense, example_sparse, difficult_example_dense, difficult_example_sparse. You should get something like the following output.

>> startup
>> example_dense
Iteration  1. Current = -7.816030706939740e-02.
              Absolute error = 1.34e-04.  Relative error = 1.18e-06.
Iteration  2. Current = -7.802629355658958e-02. Relative change = 1.34e-04.
              Absolute error = 3.49e-09.  Relative error = 3.09e-11.
Iteration  3. Current = -7.802629006654559e-02. Relative change = 3.49e-09.
              Absolute error = 2.08e-15.  Relative error = 1.84e-17.
Iteration  4. Current = -7.802629006654559e-02. Relative change = 0.00e+00.
-- Tolerance on relative change satisfied at iteration  4. Current = -7.802629006654559e-02.
License and use

The code is GPLv3 licensed. In addition, it is research code and not intended for production use. If you publish a paper using this code, a reference to the above paper would be appreciated.
The code contains parts of PROPACK, EigTool, the criss-cross algorithm, PSAPSR, and eigsplus all of which may have their own licenses.


All comments, bugs and requests are most welcome at .

Version history

Obligatory marketing picture 1

The rectangular pseudospectra Lambda_varepsilon(AV_k, V_K) of the subspaces V_k generated by subspace_pspa for iterations k = {1, 2, 3, 6, 8, 10} (from left to right, top to bottom) applied to the Grcar matrix for varepsilon = 10^{-2}. In addition, the thin lines show the square pseudospectra Lambda_varepsilon(A).

Rectangular pseudospectra

Obligatory marketing picture 2

Iterates of subspace_pspa using options.tol_svd=-1 (blue discs, always eigenvalue computations) and options.tol_svd=Inf (red circles, no eigenvalue computations) for the Grcar (left figure) and the Orr-Sommerfeld (right figure) matrices of size n = 100. The default options.tol_svd=0.1 only performs 2 to 3 eigenvalue computations (not shown).

Rectangular pseudospectra