Purpose
To solve for X = op(U)'*op(U) either the stable non-negative
definite continuous-time Lyapunov equation
2
op(A)'*X + X*op(A) = -scale *op(B)'*op(B) (1)
or the convergent non-negative definite discrete-time Lyapunov
equation
2
op(A)'*X*op(A) - X = -scale *op(B)'*op(B) (2)
where op(K) = K or K' (i.e., the transpose of the matrix K), A is
an N-by-N matrix, op(B) is an M-by-N matrix, U is an upper
triangular matrix containing the Cholesky factor of the solution
matrix X, X = op(U)'*op(U), and scale is an output scale factor,
set less than or equal to 1 to avoid overflow in X. If matrix B
has full rank then the solution matrix X will be positive-definite
and hence the Cholesky factor U will be nonsingular, but if B is
rank deficient then X may be only positive semi-definite and U
will be singular.
In the case of equation (1) the matrix A must be stable (that
is, all the eigenvalues of A must have negative real parts),
and for equation (2) the matrix A must be convergent (that is,
all the eigenvalues of A must lie inside the unit circle).
Specification
SUBROUTINE SB03OD( DICO, FACT, TRANS, N, M, A, LDA, Q, LDQ, B,
$ LDB, SCALE, WR, WI, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, FACT, TRANS
INTEGER INFO, LDA, LDB, LDQ, LDWORK, M, N
DOUBLE PRECISION SCALE
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*), Q(LDQ,*), WI(*),
$ WR(*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of Lyapunov equation to be solved as
follows:
= 'C': Equation (1), continuous-time case;
= 'D': Equation (2), discrete-time case.
FACT CHARACTER*1
Specifies whether or not the real Schur factorization
of the matrix A is supplied on entry, as follows:
= 'F': On entry, A and Q contain the factors from the
real Schur factorization of the matrix A;
= 'N': The Schur factorization of A will be computed
and the factors will be stored in A and Q.
TRANS CHARACTER*1
Specifies the form of op(K) to be used, as follows:
= 'N': op(K) = K (No transpose);
= 'T': op(K) = K**T (Transpose).
Input/Output Parameters
N (input) INTEGER
The order of the matrix A and the number of columns in
matrix op(B). N >= 0.
M (input) INTEGER
The number of rows in matrix op(B). M >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the matrix A. If FACT = 'F', then A contains
an upper quasi-triangular matrix S in Schur canonical
form; the elements below the upper Hessenberg part of the
array A are not referenced.
On exit, the leading N-by-N upper Hessenberg part of this
array contains the upper quasi-triangular matrix S in
Schur canonical form from the Shur factorization of A.
The contents of array A is not modified if FACT = 'F'.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
Q (input or output) DOUBLE PRECISION array, dimension
(LDQ,N)
On entry, if FACT = 'F', then the leading N-by-N part of
this array must contain the orthogonal matrix Q of the
Schur factorization of A.
Otherwise, Q need not be set on entry.
On exit, the leading N-by-N part of this array contains
the orthogonal matrix Q of the Schur factorization of A.
The contents of array Q is not modified if FACT = 'F'.
LDQ INTEGER
The leading dimension of array Q. LDQ >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
if TRANS = 'N', and dimension (LDB,max(M,N)), if
TRANS = 'T'.
On entry, if TRANS = 'N', the leading M-by-N part of this
array must contain the coefficient matrix B of the
equation.
On entry, if TRANS = 'T', the leading N-by-M part of this
array must contain the coefficient matrix B of the
equation.
On exit, the leading N-by-N part of this array contains
the upper triangular Cholesky factor U of the solution
matrix X of the problem, X = op(U)'*op(U).
If M = 0 and N > 0, then U is set to zero.
LDB INTEGER
The leading dimension of array B.
LDB >= MAX(1,N,M), if TRANS = 'N';
LDB >= MAX(1,N), if TRANS = 'T'.
SCALE (output) DOUBLE PRECISION
The scale factor, scale, set less than or equal to 1 to
prevent the solution overflowing.
WR (output) DOUBLE PRECISION array, dimension (N)
WI (output) DOUBLE PRECISION array, dimension (N)
If FACT = 'N', and INFO >= 0 and INFO <= 2, WR and WI
contain the real and imaginary parts, respectively, of
the eigenvalues of A.
If FACT = 'F', WR and WI are not referenced.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, or INFO = 1, DWORK(1) returns the
optimal value of LDWORK.
LDWORK INTEGER
The length of the array DWORK.
If M > 0, LDWORK >= MAX(1,4*N + MIN(M,N));
If M = 0, LDWORK >= 1.
For optimum performance LDWORK should sometimes be larger.
If LDWORK = -1, then a workspace query is assumed; the
routine only calculates the optimal size of the DWORK
array, returns this value as the first entry of the DWORK
array, and no error message related to LDWORK is issued by
XERBLA.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if the Lyapunov equation is (nearly) singular
(warning indicator);
if DICO = 'C' this means that while the matrix A
(or the factor S) has computed eigenvalues with
negative real parts, it is only just stable in the
sense that small perturbations in A can make one or
more of the eigenvalues have a non-negative real
part;
if DICO = 'D' this means that while the matrix A
(or the factor S) has computed eigenvalues inside
the unit circle, it is nevertheless only just
convergent, in the sense that small perturbations
in A can make one or more of the eigenvalues lie
outside the unit circle;
perturbed values were used to solve the equation;
= 2: if FACT = 'N' and DICO = 'C', but the matrix A is
not stable (that is, one or more of the eigenvalues
of A has a non-negative real part), or DICO = 'D',
but the matrix A is not convergent (that is, one or
more of the eigenvalues of A lies outside the unit
circle); however, A will still have been factored
and the eigenvalues of A returned in WR and WI.
= 3: if FACT = 'F' and DICO = 'C', but the Schur factor S
supplied in the array A is not stable (that is, one
or more of the eigenvalues of S has a non-negative
real part), or DICO = 'D', but the Schur factor S
supplied in the array A is not convergent (that is,
one or more of the eigenvalues of S lies outside the
unit circle);
= 4: if FACT = 'F' and the Schur factor S supplied in
the array A has two or more consecutive non-zero
elements on the first sub-diagonal, so that there is
a block larger than 2-by-2 on the diagonal;
= 5: if FACT = 'F' and the Schur factor S supplied in
the array A has a 2-by-2 diagonal block with real
eigenvalues instead of a complex conjugate pair;
= 6: if FACT = 'N' and the LAPACK Library routine DGEES
has failed to converge. This failure is not likely
to occur. The matrix B will be unaltered but A will
be destroyed.
Method
The method used by the routine is based on the Bartels and Stewart
method [1], except that it finds the upper triangular matrix U
directly without first finding X and without the need to form the
normal matrix op(B)'*op(B).
The Schur factorization of a square matrix A is given by
A = QSQ',
where Q is orthogonal and S is an N-by-N block upper triangular
matrix with 1-by-1 and 2-by-2 blocks on its diagonal (which
correspond to the eigenvalues of A). If A has already been
factored prior to calling the routine however, then the factors
Q and S may be supplied and the initial factorization omitted.
If TRANS = 'N', the matrix B is factored as (QR factorization)
_ _ _ _ _
B = P ( R ), M >= N, B = P ( R Z ), M < N,
( 0 )
_ _
where P is an M-by-M orthogonal matrix and R is a square upper
_ _ _ _ _
triangular matrix. Then, the matrix B = RQ, or B = ( R Z )Q (if
M < N) is factored as
_ _
B = P ( R ), M >= N, B = P ( R Z ), M < N.
If TRANS = 'T', the matrix B is factored as (RQ factorization)
_
_ _ ( Z ) _
B = ( 0 R ) P, M >= N, B = ( _ ) P, M < N,
( R )
_ _
where P is an M-by-M orthogonal matrix and R is a square upper
_ _ _ _ _
triangular matrix. Then, the matrix B = Q'R, or B = Q'( Z' R' )'
(if M < N) is factored as
_ _
B = ( R ) P, M >= N, B = ( Z ) P, M < N.
( R )
These factorizations are utilised to either transform the
continuous-time Lyapunov equation to the canonical form
2
op(S)'*op(V)'*op(V) + op(V)'*op(V)*op(S) = -scale *op(F)'*op(F),
or the discrete-time Lyapunov equation to the canonical form
2
op(S)'*op(V)'*op(V)*op(S) - op(V)'*op(V) = -scale *op(F)'*op(F),
where V and F are upper triangular, and
F = R, M >= N, F = ( R Z ), M < N, if TRANS = 'N';
( 0 0 )
F = R, M >= N, F = ( 0 Z ), M < N, if TRANS = 'T'.
( 0 R )
The transformed equation is then solved for V, from which U is
obtained via the QR factorization of V*Q', if TRANS = 'N', or
via the RQ factorization of Q*V, if TRANS = 'T'.
References
[1] Bartels, R.H. and Stewart, G.W.
Solution of the matrix equation A'X + XB = C.
Comm. A.C.M., 15, pp. 820-826, 1972.
[2] Hammarling, S.J.
Numerical solution of the stable, non-negative definite
Lyapunov equation.
IMA J. Num. Anal., 2, pp. 303-325, 1982.
Numerical Aspects
3 The algorithm requires 0(N ) operations and is backward stable.Further Comments
The Lyapunov equation may be very ill-conditioned. In particular,
if A is only just stable (or convergent) then the Lyapunov
equation will be ill-conditioned. A symptom of ill-conditioning
is "large" elements in U relative to those of A and B, or a
"small" value for scale. A condition estimate can be computed
using SLICOT Library routine SB03MD.
SB03OD routine can be also used for solving "unstable" Lyapunov
equations, i.e., when matrix A has all eigenvalues with positive
real parts, if DICO = 'C', or with moduli greater than one,
if DICO = 'D'. Specifically, one may solve for X = op(U)'*op(U)
either the continuous-time Lyapunov equation
2
op(A)'*X + X*op(A) = scale *op(B)'*op(B), (3)
or the discrete-time Lyapunov equation
2
op(A)'*X*op(A) - X = scale *op(B)'*op(B), (4)
provided, for equation (3), the given matrix A is replaced by -A,
or, for equation (4), the given matrices A and B are replaced by
inv(A) and B*inv(A), if TRANS = 'N' (or inv(A)*B, if TRANS = 'T'),
respectively. Although the inversion generally can rise numerical
problems, in case of equation (4) it is expected that the matrix A
is enough well-conditioned, having only eigenvalues with moduli
greater than 1. However, if A is ill-conditioned, it could be
preferable to use the more general SLICOT Lyapunov solver SB03MD.
Example
Program Text
* SB03OD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D0 )
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX
PARAMETER ( NMAX = 20, MMAX = 20 )
INTEGER LDA, LDB, LDQ, LDX, LDWORK
PARAMETER ( LDA = NMAX, LDB = MAX( MMAX,NMAX ),
$ LDQ = NMAX, LDX = NMAX )
PARAMETER ( LDWORK = 4*NMAX+MIN(MMAX,NMAX) )
* .. Local Scalars ..
DOUBLE PRECISION SCALE, TEMP
INTEGER I, INFO, J, K, M, N
CHARACTER*1 DICO, FACT, TRANS
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,LDB), DWORK(LDWORK),
$ Q(LDQ,NMAX), WR(NMAX), WI(NMAX), X(LDX,NMAX)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL SB03OD
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M, DICO, FACT, TRANS
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
IF ( LSAME( FACT, 'F' ) ) READ ( NIN, FMT = * )
$ ( ( Q(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) M
ELSE
IF ( LSAME( TRANS, 'N' ) ) THEN
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,M )
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
END IF
* Find the Cholesky factor U.
CALL SB03OD( DICO, FACT, TRANS, N, M, A, LDA, Q, LDQ, B,
$ LDB, SCALE, WR, WI, DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 20 J = 1, N
WRITE ( NOUT, FMT = 99996 ) ( B(I,J), I = 1,J )
20 CONTINUE
* Form the solution matrix X = op(U)'*op(U).
IF ( LSAME( TRANS, 'N' ) ) THEN
DO 80 I = 1, N
DO 60 J = I, N
TEMP = ZERO
DO 40 K = 1, I
TEMP = TEMP + B(K,I)*B(K,J)
40 CONTINUE
X(I,J) = TEMP
X(J,I) = TEMP
60 CONTINUE
80 CONTINUE
ELSE
DO 140 I = 1, N
DO 120 J = I, N
TEMP = ZERO
DO 100 K = J, N
TEMP = TEMP + B(I,K)*B(J,K)
100 CONTINUE
X(I,J) = TEMP
X(J,I) = TEMP
120 CONTINUE
140 CONTINUE
END IF
WRITE ( NOUT, FMT = 99995 )
DO 160 J = 1, N
WRITE ( NOUT, FMT = 99996 ) ( X(I,J), I = 1,N )
160 CONTINUE
WRITE ( NOUT, FMT = 99992 ) SCALE
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB03OD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB03OD = ',I2)
99997 FORMAT (' The transpose of the Cholesky factor U is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The solution matrix X = op(U)''*op(U) is ')
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' M is out of range.',/' M = ',I5)
99992 FORMAT (/' Scaling factor = ',F8.4)
END
Program Data
SB03OD EXAMPLE PROGRAM DATA 4 5 C N N -1.0 37.0 -12.0 -12.0 -1.0 -10.0 0.0 4.0 2.0 -4.0 7.0 -6.0 2.0 2.0 7.0 -9.0 1.0 2.5 1.0 3.5 0.0 1.0 0.0 1.0 -1.0 -2.5 -1.0 -1.5 1.0 2.5 4.0 -5.5 -1.0 -2.5 -4.0 3.5Program Results
SB03OD EXAMPLE PROGRAM RESULTS The transpose of the Cholesky factor U is 1.0000 3.0000 1.0000 2.0000 -1.0000 1.0000 -1.0000 1.0000 -2.0000 1.0000 The solution matrix X = op(U)'*op(U) is 1.0000 3.0000 2.0000 -1.0000 3.0000 10.0000 5.0000 -2.0000 2.0000 5.0000 6.0000 -5.0000 -1.0000 -2.0000 -5.0000 7.0000 Scaling factor = 1.0000