The huge amount of theoretical results available in systems and control theory has led to a variety of methods and techniques used both in industry and academia. Although based on sound theoretical results, these methods often fail when applied to real-life problems, which frequently tend to be ill-posed or highly dimensional. Such failures are often due to the lack of numerical robustness when implemented in finite-precision computations.
SLICOT Library (Subroutine Library in Control Theory) has been designed and developed with the main aim of achieving the robustness and efficiency needed in solving possibly ill-conditioned, ill-scaled, and/or large-scale real-life problems. A special emphasis has been put on structure-preserving algorithms. The main advantage of such algorithms is that the structural properties of a problem are preserved during computations. This allows the computed result to be interpreted as the exact solution of the original problem with perturbed input data, which may otherwise not be the case. This not only increases the reliability of the returned results, but often also improves their accuracy.

The following results illustrate the performance (in terms of efficiency, reliability, and accuracy) which can be obtained using SLICOT components, in comparison with equivalent computations performed by some MATLAB 5.3 functions included in the MATLAB nucleus or in the MATLAB Control Toolbox. Results obtained on newer MATLAB versions are published in related NICONET reports and papers authored by the NICONET team. The results show that some SLICOT computations are several times faster than MATLAB computations, at comparable or better accuracy; moreover, less memory is required by SLICOT routines, because the problem structure is fully exploited whenever possible. The calculations have been performed on a SUN Ultra 2 Creator 2200 workstation with 128 M RAM under SunOS 5.5, by calling from MATLAB the SLICOT gateways produced by the NAGWare Gateway Generator of the Numerical Algorithms Group Ltd. (NAG). SLICOT routines have been compiled with f77 compiler using options -O4 -native -u. SUN performance libraries have not been used for obtaining the results. The data have been randomly generated from a uniform distribution in the interval (0,1). The timing values reported are the execution times in seconds exclusively spent in the corresponding SLICOT or MATLAB code portions.

Fast Fourier transforms

Tables 1 and 2 display comparative results for the SLICOT Fast Fourier transform routines DG01MD and DG01ND (for complex and real sequences, respectively), and the corresponding MATLAB function fft. Random sequences X of length n were generated. Besides better efficiency, the accuracy of SLICOT routines shows an improvement which can be two orders of magnitude. The accuracy was measured by computing the distance between the original sequence and the inverse Fourier transform of the transformed sequence, norm(X - ifft( fft( X ) ))/norm(X), which should theoretically be zero, and similarly for the SLICOT routines. Note that DG01ND is 2-4 times faster than MATLAB.

Table 1. Comparison between DG01MD and MATLAB results
 TimeRelative error





0.00 7.54e-16 8.14e-15
4096 0.02 0.03 1.36e-15 2.49e-14
16384 0.08 0.13 2.32e-15 1.10e-13
65536 0.80 1.22 5.04e-15 4.77e-13
Table 2. Comparison between DG01ND and MATLAB results
 TimeRelative error
1024 0.00 0.01 6.44e-16 6.56e-15
4096 0.01 0.02 1.44e-15 2.00e-14
16384 0.03 0.11 2.46e-15 8.21e-14
65536 0.23 0.88 4.22e-15 4.26e-13

Matrix exponential

Tables 3 and 4 display comparative results for the SLICOT subroutines MB05MD and MB05OD, for the computation of the matrix exponential, and the equivalent MATLAB functions expm3 and expm1, respectively. The relative errors reported are those corresponding to the differences between the results obtained from SLICOT and MATLAB. SLICOT codes are faster than MATLAB calculations, at a comparable accuracy. An increase of over 30 times in speed has been achieved for n = 256 for upper triangular matrices. This remarkable speeding up was possible because the balancing strategy implemented in MB05OD can reduce the matrix norm for triangular matrices and save matrix multiplications, while expm1 cannot.

Table 3. Comparison between MB05MD and MATLAB expm3 results
 TimeRelative error
16 0.01 0.01 5.84e-15
32 0.03 0.04 1.08e-14
64 0.18 0.33 1.99e-14
128 1.21 2.60 1.82e-14
256 9.29 34.04 1.43e-13
Table 4. Comparison between MB05OD and MATLAB expm1 results
 TimeRelative error
16 0.00 0.00 4.84e-17
32 0.00 0.02 7.67e-17
64 0.02 0.11 1.60e-16
128 0.08 0.86 2.14e-16
256 0.61 19.01 1.99e-16

A = triu(rand(n,n),n/2)

Controllability staircase form

Tables 5 and 6 demonstrate the behaviour of the SLICOT routine AB01ND for computing the controllability staircase form of a linear system in state-space form given by a matrix pair (A,B), with A n-by-n, B n-by-m, and MATLAB function ctrbf provided by the Control Toolbox. In Table 6, m is taken as 1. It should be mentioned that ctrbf uses singular value decompositions for rank determination while AB01ND employs rank-revealing QR factorizations. The relative residuals have been computed using the following formulas
max(norm(z'*A*z - a)/norm(A), norm(z'*B - b)/norm(B));
max(norm(A - z*a*z')/norm(A), norm(B - z*b)/norm(B));
for SLICOT and MATLAB, respectively, where A and B are the given system matrices, a and b are the computed matrices in controllability staircase form, and z is the computed orthogonal transformation matrix. The SLICOT routine is up to ten times faster than ctrbf in the multi-input case and up to 60 times faster in the single-input case.

Table 5. Comparison between AB01ND and MATLAB results
 TimeRelative residuals
n m



16 2



4.75e-16 6.38e-16
32 4 0.01 0.04 3.92e-16 1.67e-15
64 8 0.04 0.22 6.05e-16 8.50e-16
128 16 0.32 1.44 6.46e-16 2.13e-15
256 32 2.56 25.10 1.35e-15 3.82e-15
Table 6. Comparison between AB01ND and MATLAB results, m=1
 TimeRelative residuals
16 0.00 0.03 4.94e-16 4.04e-16
32 0.00 0.17 7.01e-16 9.30e-16
64 0.04 1.29 5.83e-16 1.23e-15
128 0.32 19.12 7.41e-16 3.58e-15

Lyapunov and Riccati equations

Tables 7 and 8 give performance results for SLICOT's general Lyapunov solver SB03MD, which solves both continuous- and discrete-time equations, and MATLAB functions lyap and dlyap. The SLICOT routine can solve equations defined by
op(A)' X + X op(A) = sigma C,(1)
op(A)' X op(A) - X = sigma C,(2)
where op(A) = A or A' (the transpose of A), and sigma is a scale factor, set less than or equal to 1 to avoid overflow in computing X. Moreover, an estimate of the condition number is returned. Specifically, Table 7 reports results for continuous-time problems with op(A) = A, generated using the MATLAB statements
A = rand(n,n); C = rand(n,n); C = C + C';
so the solutions were not known. Therefore, only relative residuals are given in Table 7.

Table 7. Comparison between SB03MD and MATLAB results (continuous-time case)
 TimeRelative residuals
16 0.01 0.05 3.69e-15 3.59e-15
32 0.03 0.12 3.54e-15 8.96e-15
64 0.21 0.70 3.56e-15 1.96e-14
128 1.50 5.38 1.73e-14 2.86e-14

Similar results have been obtained for the transposed case, or for discrete-time equations. SLICOT residuals were usually less then MATLAB residuals in our tests. The speed-up has been even larger when compared with the previous MATLAB 4.2 versions of lyap and dlyap.
Table 8 reports results for discrete-time problems with op(A) = A'. The problems were generated using the MATLAB statements
A = rand(n,n); X0 = rand(n,n) + n*eye(n); X0 = X0 + X0';
C = A*X0*A' - X0; C = (C + C')/2;
so the relative errors could be computed.

Table 8. Comparison between SB03MD and MATLAB results (discrete-time case)
 TimeRelative errors
16 0.01 0.03 4.01e-13 2.22e-13
32 0.03 0.11 1.63e-13 1.14e-12
64 0.22 0.70 3.48e-12 9.04e-12
128 1.73 5.49 1.51e-11 5.82e-11

Table 9 gives performance results for the SLICOT Riccati solver SB02OD, which solves both continuous- and discrete-time equations, and the MATLAB functions care and dare, all based on the generalized Schur vector approaches. The SLICOT solver can also compute the anti-stabilizing solutions. The problems solved were generated using the MATLAB statements
A = rand(n,n); B = rand(n,m); G = B*inv(R)*B'; G = (G + G')/2;
(with m = n), and the weighting matrices R and Q have been computed to have given random (non-negative) eigenvalues, using a gateway for the test routine DLATMS from LAPACK. Table 9 reports the timings and relative residuals for computing stabilizing solutions for continuous-time problems.

Table 9. Comparison between SB02OD and MATLAB results (continuous-time case, m = n)
 TimeRelative residuals
16 0.07 0.37 1.83e-14 3.81e-14
32 0.30 0.86 4.99e-14 1.13e-13
64 2.80 4.01 7.21e-12 2.10e-11
128 24.00 79.96 2.44e-11 6.72e-11

Table 10 gives performance results for the SLICOT Riccati solver SB02MD, which solves both continuous- and discrete-time equations, using the Schur vectors approach. For discrete-time problems, SB02MD includes an option to use the inverse of the symplectic 2n-by-2n matrix, which is always faster than the standard approach.

Table 10. Comparison between SB02MD and MATLAB results (continuous-time case, m = n/2)
 TimeRelative errors
16 0.03 0.22 3.03e-14 8.90e-14
32 0.14 1.48 1.67e-13 7.31e-13
64 1.02 3.89 3.10e-13 3.85e-13
128 7.78 35.86 4.66e-11 9.97e-11

Discrete-time problems can be solved even faster by SB02MD, but the accuracy depends on the conditioning of the matrix A.

Structure-Preserving Algorithms

Hankel singular values

The use of structure-preserving algorithms not only increases the reliability of the returned results, but often also improves their accuracy, as shown below for an important computational problem.
The controllability and observability Gramians Pc, Po, respectively, of a stable state-space realization (A,B,C) of a continuous-time linear time-invariant system, where A is n-by-n, B is n-by-m, and C is p-by-n, are given by the solutions of the stable Lyapunov equations
A Pc + Pc A' = - B B' (3)
A' Po + Po A = - C' C (4)
which are non-negative definite. The nonnegative square roots of the eigenvalues of the product PcPo, called the Hankel singular values, play a fundamental role in finding balanced realizations and in model reduction. To guarantee the symmetry and definiteness of the computed Gramians, the special problem structure should be taken into account. Besides this, exploiting the structure usually results in a reduction of necessary computational operations and memory requirements. The SLICOT Library routine SB03OD solves for X = op(U)' op(U) either continuous-time or discrete-time Lyapunov equations,
op(A)' X + X op(A) = -sigma*sigma op(M)' op(M) (5)
op(A)' X op(A) - X = -sigma*sigma op(M)' op(M) (6)
where A is square and stable (in the continuous- or discrete-time sense, respectively), M is rectangular, U is upper triangular, and sigma is a scale factor, set less than or equal to 1 to avoid overflow in computing X. The routine uses the same real Schur form of A for all different computations, and optionally A can be given in such a form on input.
The data have been generated as
A = rand(n,n) - n*eye(n); B = rand(n,m); C = rand(l,n);
Table 11 gives comparative results using the gateway for the SLICOT routine SB03OD and MATLAB function lyap, for solving the two Lyapunov equations (3) and (4). It is apparent that the speed-up compared with MATLAB is about 5. Besides improved efficiency, the accuracy of the SLICOT routine is sometimes better. More important, the Hankel singular values can be obtained as the singular values of the triangular matrix V U, where U and V are the Cholesky factors computed by SB03OD for the solutions of the equations (5), for op(K) = K' and op(K) = K, respectively. On the contrary, MATLAB calculations need to be based on sqrt(eig(X*Y)), where X = Pc and Y = Po are the solutions of the Lyapunov equations (3) and (4), respectively.

Table 11. Comparison between SB03OD and MATLAB results (m = p = n/2)
 TimeRelative residuals in X
16 0.01 0.06 3.17e-14 2.97e-14
32 0.04 0.24 8.41e-14 7.23e-14
64 0.24 1.30 1.48e-13 1.82e-13
128 1.68 9.44 4.13e-13 4.21e-13

Similar results have been obtained for the relative residuals in Y.
For one sample of the largest example, two computed eigenvalues of the matrix X Y were negative (close to the machine accuracy), and four were complex conjugate. In other words, in extreme cases, MATLAB could give physically meaningless negative or complex Hankel "singular values". This, of course, does not happen for well-conditioned problems.

Eigenvalues of Hamiltonian matrices

Table 12 displays timing results for the SLICOT Hamiltonian eigenvalue solver based on the square-reduced form (MB04ZD and MB03SD) and MATLAB function eig. There is no MATLAB function exploiting the Hamiltonian structure.

Table 12. Comparison between MB04ZD + MB03SD and MATLAB results
16 0.00 0.01
32 0.04 0.09
64 0.26 0.60
128 2.10 4.54
256 17.68 94.84

For random matrices, MATLAB often gave eigenvalues very close to the imaginary axis, and which do not not appear in pairs lambda, and -lambda, as the structure ensures. An example is given below. Let G and Q be symmetric matrices, and consider the MATLAB statements
A = [ ...
8.156785069508218e-01 9.583801041682488e-01 3.592539466930961e-01
3.786356003299222e-01 3.034683538610545e-01 7.876353002957823e-01
4.773517592106414e-01 8.360131178618799e-01 7.242888156356292e-01 ];

Q = [ ...
8.890570332721912e-02 9.990740343443876e-01 8.732331409873815e-01
9.990740343443876e-01 6.033984620815241e-01 5.861126862288001e-01
8.732331409873815e-01 5.861126862288001e-01 7.904409030391265e-01 ];

G = [ ...
5.951660134153334e-01 9.187384627530358e-01 6.419067121426090e-01
9.187384627530358e-01 2.943561392390958e-01 4.831226444677942e-01
6.419067121426090e-01 4.831226444677942e-01 7.635187432100927e-01 ];

H = [ A G; Q -A'];


ans =
6.938893903907228e-17 + 2.317956004861914e-01i
6.938893903907228e-17 - 2.317956004861914e-01i

Condition estimates for Lyapunov and Riccati equations

Many linear systems analysis and design problems (for instance, stability analysis, transfer-function matrix factorization, linear-quadratic or robust optimization problems), involve the solution of either continuous-time or discrete-time matrix Lyapunov and Riccati equations. In such applications, it is very useful to have easily computable estimates of the problem conditioning and error bounds for the solutions. Such estimators have been developed and are included in the SLICOT Library.
Table 13 displays estimates for the "separation" of the matrices A' and A, sepd(A',A), which is needed to compute the condition number for discrete-time matrix Lyapunov and Riccati equations. Similarly, the separation of the matrices A' and -A, sep(A',-A), is used for continuous-time equations. For A an n-by-n matrix with real elements, sep(.,.) and sepd(.,.) are defined in terms of the smallest singular value of a large matrix, by

sep(A', -A) = sigma_min (P), P = kron(I_n, A') + kron(A', I_n),
sepd(A', A) = sigma_min (P), P = kron(A', A') - I_n*n,

where kron(X,Y) denotes the Kronecker product of the matrices A and B, kron(X,Y) = (x_{ij} Y). The evaluation of the norms above by computing the minimal singular value sigma_min is not practical for medium-size or large-scale problems, due to the excessive computational effort and memory requirements. SLICOT routines exploit the problem structure.
The results in Table 13 have been computed by two approaches:
(a) solution of original Lyapunov equations;
(b) solution of reduced Lyapunov equations,
as well as by MATLAB function svd. The best results were obtained for approach (b), which can be two-four times faster than approach (a), and possibly thousands times faster than MATLAB svd.

Table 13. Comparison between SLICOT and MATLAB results for sepd
n (a) (b) MATLAB (a) (b) MATLAB
10 0.01 0.00 0.17 2.50e-2 2.81e-2 7.69e-2
20 0.02 0.01 6.48 6.00e-4 7.36e-4 1.83e-3
30 0.04 0.02 71.19 1.37e-3 3.47e-3 6.42e-3
40 0.08 0.03 411.71 3.08e-4 4.46e-4 6.31e-4


SLICOT computations are several times faster than MATLAB computations, at comparable or better accuracy. Moreover, the problem structure is better exploited by SLICOT routines. Even better efficiency is attainable by using machine-tailored BLAS libraries, which, together with LAPACK subroutines provide the basic algorithmic level of SLICOT, and enable to exploit the potential of modern high-performance computer architectures.

This email address is being protected from spambots. You need JavaScript enabled to view it., April 7, 1999

Control Software and Benchmarking

In the analysis of numerical methods and their implementation as numerical software it is extremely important to be able to test the correctness of the implementation as well as the performance of the method. This validation is one of the major steps in the construction of a software library, in particular if this library is used in practical applications.

In order to carry out such a test it is important to have tools that yield an evaluation of the performance of the method as well as the implementation with respect to correctness, accuracy and speed. Similar tools are needed to be able to compare different numerical methods, to test their robustness, and also to analyse the behaviour of the methods in extreme situation, i.e., on problems where the limit of the possible accuracy is reached.

In many application areas therefore benchmark collections have been created that can partially serve for this purpose. Such collections are heavily used. In our opinion, in order to have a fair evaluation and a comparison of methods and software, there should be a standardized set of examples on which newly developed methods and their implementation can be tested. It was one of the goals of WGS to create such testing and validation environments for the area of numerical methods in control, in particular it was planned to accompany the SLICOT library with benchmark collections for each of the major problem areas. In order to make such collections useful, it is important that they cover a wide range of problems and also problems difficult to solve in finite arithmetic are included. Such problems in particular drive the methods and their implementation to a limit. These are ideal test cases, since errors and failures usually occur only in extreme cases and these are often not covered by standard software validation procedures.

The SLICOT library currently contains such benchmark collections for:
- standard and generalized continuous-time and discrete-time systems models (see reports 2 and 3)
- continuous-time and discrete-time standard and generalized Lyapunov matrix equations (see reports 4 and 5)
- continuous-time and discrete-time algebraic Riccati matrix equations (see reports 6 and 7)
- identification (see report 8)
- model reduction of (high order) linear time invariant dynamical systems containing some useful "real world" examples reflecting current problems in applications can be found here (see also report 9)

Everybody is invited to submit such benchmark examples for the SLICOT benchmark collection. For submissions please contact This email address is being protected from spambots. You need JavaScript enabled to view it.

Available Reports

  1. Volker Mehrmann and Thilo Penzl: Benchmark collections in SLICOT SLICOT Working Note 1998-5, June 1998.
  2. Daniel Kressner, Volker Mehrmann and Thilo Penzl:  DTDSX - A collection of benchmark examples for state-space realizations of time-invariant discrete-time systems; SLICOT Working Note 1998-10: November 1998, revised June 1999.
  3. Daniel Kressner, Volker Mehrmann and Thilo Penzl: CTDSX - A collection of benchmark examples for state-space realizations of time-invariant continuous-time systems; SLICOT Working Note 1998-9: November 1998.
  4. Daniel Kressner, Volker Mehrmann and Thilo Penzl: DTLEX - A collection of benchmark examples for discrete-time Lyapunuv equations; SLICOT Working Note 1999-7: June 1999.
  5. Daniel Kressner, Volker Mehrmann and Thilo Penzl: CTLEX - A collection of benchmark examples for continuous-time Lyapunuv equations; SLICOT Working Note 1999-6: June 1999.
  6. Jörn Abels and Peter Benner: DAREX - A collection of benchmark examples for discrete-time algebraic Riccati equations (version 2.0); SLICOT Working Note 1999-16: December 1999.
  7. Jörn Abels and Peter Benner: CAREX - A collection of benchmark examples for continuous-time algebraic Riccati equations (version 2.0); SLICOT Working Note 1999-14: December 1999.
  8. Ad van den Boom, Ton Backx and Yucai Zhu: Benchmarks for identification; NICONET Report 1999-19: July 2000.
  9. Younès Chahlaoui and Paul Van Dooren: A collection of benchmark examples for model reduction of linear time invariant dynamical systems; SLICOT Working Note 2002-2: February 2002.

Some Matrix Collections Suitable for Benchmarking:

Ad van den Boom, March 4, 2002
Updated: Vasile Sima, August 31, 2005; June 15, 2006; June 8, 2009

The production of high quality, portable, and easy-to-use scientific software must heavily rely on rigorous implementation and documentation standards. To ensure such high quality demands for the SLICOT library, the Working Group on Software (WGS) produced in 1990 the

SLICOT Implementation and Documentation Standards,

to which all software submissions for the first two releases of the SLICOT library had to conform. With the new direction in the development of SLICOT, i.e. its transition to a copyrighted library entirely based on LAPACK and BLAS routines, also a change in the documentation and implementation standards was necessary. For this reason a new standard has been developed by the WGS.

The main changes in the new standard consist of introducing new calling sequences for subroutines to conform to LAPACK calling conventions, adopting of new data structures as well as enlarging the Library Index with new chapters and subchapters reflecting recent developments in the computational aspects of systems theory.

To get a compressed postscript version of a report describing the new standard click here. A new version of the Contributors Kit is also available.

Andras Varga,  March 21, 1997
Updated: Vasile Sima,  June 13, 2006

The SLICOT library is organised by chapters. Each chapter can be identified by a single letter. The following chapters are included:

The documentation of SLICOT is on-line available. It contains complete descriptions of each available subroutine. To see a sample documentation of a routine (AB01MD.f), click here.

Vasile Sima, April 28, 2005; updated: June 15, 2006; November 4, 2013

Please find below a list of Source codes for SLICOT software that can be downloaded via Obtaining SLICOT. SLICOT Version 4.5 Archives are available under GNU License policy and are directly and freely downloadable for everybody.
The Release Notes can be dowloaded below from this page.

SLICOT Source Code Archives

Main SLICOT Archives

  • slicot.tar.gz
    Contains the SLICOT Library source files and related example and documentation files for Unix platforms.
    Contains the SLICOT Library source files and related example and documentation files for Windows platforms.

SLICOT Version 4.5 Archives

  • slicot45.tar.gz
    Contains the SLICOT Library Version 4.5 source files and related example and documentation files for Unix platforms.
    Contains the SLICOT Library Version 4.5 source files and related example and documentation files for Windows platforms.

Parallel SLICOT Archives for Model Reduction

  • pab09a.tar.gz
    Contains source codes for model reduction of high order systems by direct methods.
  • plicmr.tar.gz
    Contains source codes for model reduction of high order systems by iterative methods.

Interface to Nonlinear Solvers Archives

  • ODESolver.tgz
    Contains interface codes for the ordinary differential equations solvers.
  • DAESolver.tgz
    Contains interface codes for the differential algebraic equations solvers.
  • EtcSolver.tgz
    Contains interface codes for:
    - the differential Riccati equations solver DRESOL;
    - the nonlinear algebraic equations solver KINSOL;
    - the optimization problems solver FSQP.

BLAS Level-3 Implementations of Common Solvers for (Quasi-)Triangular Generalized Lyapunov Equations

  • SLWN2014_1_code.tar.gz
    Contains Fortran 90 BLAS Level-3 codes for solving (quasi-)triangular generalized Lyapunov equations on Linux Ubuntu 12.10.

Release Notes

Vasile Sima, April 28, 2005, updated: October 18, 2009