A collection of functions has been implemented as a MatrixEquation.jl Julia package to solve several classes of Lyapunov, Sylvester and Riccati matrix equations. The goal was to demonstrate that programs written in the Julia language can achieve high computational performance, which is comparable with the performance of efficient structure exploiting Fortran implementations, as those available in the Systems and Control Library SLICOT.

The available functions cover both standard and generalized continuous and discrete Lyapunov, Sylvester and Riccati equations for both real and complex data. The functions for the solution of Lyapunov and Sylvester equations rely on efficient structure exploiting solvers for which the input data are in Schur or generalized Schur forms. A comprehensive set of Lyapunov and Sylvester operators has been implemented, which allow the estimation of condition numbers of these operators. The implementations of Riccati equation solvers employ orthogonal Schur vectors based methods and their extensions to linear matrix pencil based reduction approaches. The calls of all functions with adjoint (in complex case) or transposed (in real case) arguments are fully supported by appropriate computational algorithms, thus the matrix copying operations are mostly avoided.

The current version of the package includes the following main functions:

Solution of Lyapunov equations

lyapc Solution of the continuous Lyapunov equations AX+XA'+C = 0 and A'X+XA+C = 0 and of their generalized versions AXE'+EXA'+C = 0 and A'XE+E'XA+C = 0.

lyapd Solution of discrete Lyapunov equations AXA' - X +C = 0 and AXA' - X +C = 0 and of their generalized versions AXA' - EXE'+C = 0 and A'XA - E'XE+C = 0.

plyapc Computation of the upper triangular factor U of the positive semi-definite solution X = UU' of the continuous Lyapunov equation AX+XA'+BB' = 0 or of the positive semi-definite solution X = U'U of the continuous Lyapunov equation A'X+XA+B'B = 0 and of their generalized versions AXE' + EXA' + BB' = 0 and A'XE + E'XA + B'B = 0.

plyapd Computation of the upper triangular factor U of the positive semi-definite solution X = UU' of the discrete Lyapunov equation AXA' - X+BB' = 0 or of the positive semi-definite solution X = U'U of the discrete Lyapunov equation A'XA - X+B'B = 0 and of their generalized versions AXA' - EXE' + BB' = 0 and A'XA - E'XE + B'B = 0.

plyaps Efficient computation of the above solutions for A in a Schur form or the pair (A,E) in a generalized Schur form.

Solution of algebraic Riccati equations

arec Solution of the continuous Riccati equation A'X+XA-XRX+Q = 0 and of the control oriented continuous Riccati equation A'X+XA-(XB+S)R^-1(B'X+S')+Q = 0.

garec Solution of the generalized continuous Riccati equation A'XE+E'XA-(A'XB+S)R^-1(B'XA+S')+Q = 0.

ared Solution of the discrete Riccati equation A'XA - X - (A'XB+S)(R+B'XB)^-1(B'XA+S') + Q = 0.

gared Solution of the generalized discrete Riccati equation A'XA - E'XE - (A'XB+S)(R+B'XB)^-1(B'XA+S') + Q = 0.

Solution of Sylvester equations and systems

sylvc Solution of the (continuous) Sylvester equation AX+XB = C.

sylvd Solution of the (discrete) Sylvester equation AXB+X = C.

gsylv Solution of the generalized Sylvester equation AXB+CXD = E.

sylvsys Solution of the Sylvester system of matrix equations AX+YB = C, DX+YE = F.

dsylvsys Solution of the dual Sylvester system of matrix equations AX+DY = C, XB+YE = F.

Norm, condition and separation estimation of linear operators

opnorm1 Computation of the 1-norm of a linear operator.

opnorm1est Estimation of 1-norm of a linear operator.

oprcondest Estimation of the reciprocal 1-norm condition number of a linear operator.

opsepest Estimation of the separation of a linear operator.