Purpose
To compute the Cholesky factor U of the matrix X,
T
X = op(U) * op(U),
which is the solution of either the generalized
c-stable continuous-time Lyapunov equation
T T
op(A) * X * op(E) + op(E) * X * op(A)
2 T
= - SCALE * op(B) * op(B), (1)
or the generalized d-stable discrete-time Lyapunov equation
T T
op(A) * X * op(A) - op(E) * X * op(E)
2 T
= - SCALE * op(B) * op(B), (2)
without first finding X and without the need to form the matrix
op(B)**T * op(B).
op(K) is either K or K**T for K = A, B, E, U. A and E are N-by-N
matrices, op(B) is an M-by-N matrix. The resulting matrix U is an
N-by-N upper triangular matrix with non-negative entries on its
main diagonal. SCALE is an output scale factor set to avoid
overflow in U.
In the continuous-time case (1) the pencil A - lambda * E must be
c-stable (that is, all eigenvalues must have negative real parts).
In the discrete-time case (2) the pencil A - lambda * E must be
d-stable (that is, the moduli of all eigenvalues must be smaller
than one).
Specification
SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q,
$ LDQ, Z, LDZ, B, LDB, SCALE, ALPHAR, ALPHAI,
$ BETA, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
DOUBLE PRECISION SCALE
INTEGER INFO, LDA, LDB, LDE, LDQ, LDWORK, LDZ, M, N
CHARACTER DICO, FACT, TRANS
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
$ BETA(*), DWORK(*), E(LDE,*), Q(LDQ,*), Z(LDZ,*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies which type of the equation is considered:
= 'C': Continuous-time equation (1);
= 'D': Discrete-time equation (2).
FACT CHARACTER*1
Specifies whether the generalized real Schur
factorization of the pencil A - lambda * E is supplied
on entry or not:
= 'N': Factorization is not supplied;
= 'F': Factorization is supplied.
TRANS CHARACTER*1
Specifies whether the transposed equation is to be solved
or not:
= 'N': op(A) = A, op(E) = E;
= 'T': op(A) = A**T, op(E) = E**T.
Input/Output Parameters
N (input) INTEGER
The order of the matrix A. N >= 0.
M (input) INTEGER
The number of rows in the matrix op(B). M >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, if FACT = 'F', then the leading N-by-N upper
Hessenberg part of this array must contain the
generalized Schur factor A_s of the matrix A (see
definition (3) in section METHOD). A_s must be an upper
quasitriangular matrix. The elements below the upper
Hessenberg part of the array A are not referenced.
If FACT = 'N', then the leading N-by-N part of this
array must contain the matrix A.
On exit, the leading N-by-N part of this array contains
the generalized Schur factor A_s of the matrix A. (A_s is
an upper quasitriangular matrix.)
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,N).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, if FACT = 'F', then the leading N-by-N upper
triangular part of this array must contain the
generalized Schur factor E_s of the matrix E (see
definition (4) in section METHOD). The elements below the
upper triangular part of the array E are not referenced.
If FACT = 'N', then the leading N-by-N part of this
array must contain the coefficient matrix E of the
equation.
On exit, the leading N-by-N part of this array contains
the generalized Schur factor E_s of the matrix E. (E_s is
an upper triangular matrix.)
LDE INTEGER
The leading dimension of the array E. LDE >= MAX(1,N).
Q (input/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 from
the generalized Schur factorization (see definitions (3)
and (4) in section METHOD).
If FACT = 'N', Q need not be set on entry.
On exit, the leading N-by-N part of this array contains
the orthogonal matrix Q from the generalized Schur
factorization.
LDQ INTEGER
The leading dimension of the array Q. LDQ >= MAX(1,N).
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
On entry, if FACT = 'F', then the leading N-by-N part of
this array must contain the orthogonal matrix Z from
the generalized Schur factorization (see definitions (3)
and (4) in section METHOD).
If FACT = 'N', Z need not be set on entry.
On exit, the leading N-by-N part of this array contains
the orthogonal matrix Z from the generalized Schur
factorization.
LDZ INTEGER
The leading dimension of the array Z. LDZ >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,N1)
On entry, if TRANS = 'T', the leading N-by-M part of this
array must contain the matrix B and N1 >= MAX(M,N).
If TRANS = 'N', the leading M-by-N part of this array
must contain the matrix B and N1 >= N.
On exit, the leading N-by-N part of this array contains
the Cholesky factor U of the solution matrix X of the
problem, X = op(U)**T * op(U).
If M = 0 and N > 0, then U is set to zero.
LDB INTEGER
The leading dimension of the array B.
If TRANS = 'T', LDB >= MAX(1,N).
If TRANS = 'N', LDB >= MAX(1,M,N).
SCALE (output) DOUBLE PRECISION
The scale factor set to avoid overflow in U.
0 < SCALE <= 1.
ALPHAR (output) DOUBLE PRECISION array, dimension (N)
ALPHAI (output) DOUBLE PRECISION array, dimension (N)
BETA (output) DOUBLE PRECISION array, dimension (N)
If INFO = 0, 3, 5, 6, or 7, then
(ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, are the
eigenvalues of the matrix pencil A - lambda * E.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The dimension of the array DWORK.
LDWORK >= MAX(1,4*N,6*N-6), if FACT = 'N';
LDWORK >= MAX(1,2*N,6*N-6), if FACT = 'F'.
For good performance, LDWORK should 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: the pencil A - lambda * E is (nearly) singular;
perturbed values were used to solve the equation
(but the reduced (quasi)triangular matrices A and E
are unchanged);
= 2: FACT = 'F' and the matrix contained in the upper
Hessenberg part of the array A is not in upper
quasitriangular form;
= 3: FACT = 'F' and there is a 2-by-2 block on the main
diagonal of the pencil A_s - lambda * E_s whose
eigenvalues are not conjugate complex;
= 4: FACT = 'N' and the pencil A - lambda * E cannot be
reduced to generalized Schur form: LAPACK routine
DGEGS (or DGGES) has failed to converge;
= 5: DICO = 'C' and the pencil A - lambda * E is not
c-stable;
= 6: DICO = 'D' and the pencil A - lambda * E is not
d-stable;
= 7: the LAPACK routine DSYEVX utilized to factorize M3
failed to converge in the discrete-time case (see
section METHOD for SLICOT Library routine SG03BU).
This error is unlikely to occur.
Method
An extension [2] of Hammarling's method [1] to generalized
Lyapunov equations is utilized to solve (1) or (2).
First the pencil A - lambda * E is reduced to real generalized
Schur form A_s - lambda * E_s by means of orthogonal
transformations (QZ-algorithm):
A_s = Q**T * A * Z (upper quasitriangular) (3)
E_s = Q**T * E * Z (upper triangular). (4)
If the pencil A - lambda * E has already been factorized prior to
calling the routine however, then the factors A_s, E_s, Q and Z
may be supplied and the initial factorization omitted.
Depending on the parameters TRANS and M the N-by-N upper
triangular matrix B_s is defined as follows. In any case Q_B is
an M-by-M orthogonal matrix, which need not be accumulated.
1. If TRANS = 'N' and M < N, B_s is the upper triangular matrix
from the QR-factorization
( Q_B O ) ( B * Z )
( ) * B_s = ( ),
( O I ) ( O )
where the O's are zero matrices of proper size and I is the
identity matrix of order N-M.
2. If TRANS = 'N' and M >= N, B_s is the upper triangular matrix
from the (rectangular) QR-factorization
( B_s )
Q_B * ( ) = B * Z,
( O )
where O is the (M-N)-by-N zero matrix.
3. If TRANS = 'T' and M < N, B_s is the upper triangular matrix
from the RQ-factorization
( Q_B O )
(B_s O ) * ( ) = ( Q**T * B O ).
( O I )
4. If TRANS = 'T' and M >= N, B_s is the upper triangular matrix
from the (rectangular) RQ-factorization
( B_s O ) * Q_B = Q**T * B,
where O is the N-by-(M-N) zero matrix.
Assuming SCALE = 1, the transformation of A, E and B described
above leads to the reduced continuous-time equation
T T
op(A_s) op(U_s) op(U_s) op(E_s)
T T
+ op(E_s) op(U_s) op(U_s) op(A_s)
T
= - op(B_s) op(B_s) (5)
or to the reduced discrete-time equation
T T
op(A_s) op(U_s) op(U_s) op(A_s)
T T
- op(E_s) op(U_s) op(U_s) op(E_s)
T
= - op(B_s) op(B_s). (6)
For brevity we restrict ourself to equation (5) and the case
TRANS = 'N'. The other three cases can be treated in a similar
fashion.
We use the following partitioning for the matrices A_s, E_s, B_s
and U_s
( A11 A12 ) ( E11 E12 )
A_s = ( ), E_s = ( ),
( 0 A22 ) ( 0 E22 )
( B11 B12 ) ( U11 U12 )
B_s = ( ), U_s = ( ). (7)
( 0 B22 ) ( 0 U22 )
The size of the (1,1)-blocks is 1-by-1 (iff A_s(2,1) = 0.0) or
2-by-2.
We compute U11 and U12**T in three steps.
Step I:
From (5) and (7) we get the 1-by-1 or 2-by-2 equation
T T T T
A11 * U11 * U11 * E11 + E11 * U11 * U11 * A11
T
= - B11 * B11.
For brevity, details are omitted here. See [2]. The technique
for computing U11 is similar to those applied to standard
Lyapunov equations in Hammarling's algorithm ([1], section 6).
Furthermore, the auxiliary matrices M1 and M2 defined as
follows
-1 -1
M1 = U11 * A11 * E11 * U11
-1 -1
M2 = B11 * E11 * U11
are computed in a numerically reliable way.
Step II:
The generalized Sylvester equation
T T T T
A22 * U12 + E22 * U12 * M1 =
T T T T T
- B12 * M2 - A12 * U11 - E12 * U11 * M1
is solved for U12**T.
Step III:
It can be shown that
T T T T
A22 * U22 * U22 * E22 + E22 * U22 * U22 * A22 =
T T
- B22 * B22 - y * y (8)
holds, where y is defined as
T T T T T T
y = B12 - ( E12 * U11 + E22 * U12 ) * M2 .
If B22_tilde is the square triangular matrix arising from the
(rectangular) QR-factorization
( B22_tilde ) ( B22 )
Q_B_tilde * ( ) = ( ),
( O ) ( y**T )
where Q_B_tilde is an orthogonal matrix of order N, then
T T T
- B22 * B22 - y * y = - B22_tilde * B22_tilde.
Replacing the right hand side in (8) by the term
- B22_tilde**T * B22_tilde leads to a reduced generalized
Lyapunov equation of lower dimension compared to (5).
The recursive application of the steps I to III yields the
solution U_s of the equation (5).
It remains to compute the solution matrix U of the original
problem (1) or (2) from the matrix U_s. To this end we transform
the solution back (with respect to the transformation that led
from (1) to (5) (from (2) to (6)) and apply the QR-factorization
(RQ-factorization). The upper triangular solution matrix U is
obtained by
Q_U * U = U_s * Q**T (if TRANS = 'N')
or
U * Q_U = Z * U_s (if TRANS = 'T')
where Q_U is an N-by-N orthogonal matrix. Again, the orthogonal
matrix Q_U need not be accumulated.
References
[1] Hammarling, S.J.
Numerical solution of the stable, non-negative definite
Lyapunov equation.
IMA J. Num. Anal., 2, pp. 303-323, 1982.
[2] Penzl, T.
Numerical solution of generalized Lyapunov equations.
Advances in Comp. Math., vol. 8, pp. 33-48, 1998.
Numerical Aspects
The number of flops required by the routine is given by the
following table. Note that we count a single floating point
arithmetic operation as one flop.
| FACT = 'F' FACT = 'N'
---------+--------------------------------------------------
M <= N | (13*N**3+6*M*N**2 (211*N**3+6*M*N**2
| +6*M**2*N-2*M**3)/3 +6*M**2*N-2*M**3)/3
|
M > N | (11*N**3+12*M*N**2)/3 (209*N**3+12*M*N**2)/3
Further Comments
The Lyapunov equation may be very ill-conditioned. In particular, if DICO = 'D' and the pencil A - lambda * E has a pair of almost reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost degenerate pair of eigenvalues, then the Lyapunov equation will be ill-conditioned. Perturbed values were used to solve the equation. A condition estimate can be obtained from the routine SG03AD. When setting the error indicator INFO, the routine does not test for near instability in the equation but only for exact instability.Example
Program Text
* SG03BD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX
PARAMETER ( NMAX = 20 )
INTEGER LDA, LDB, LDE, LDQ, LDZ
PARAMETER ( LDA = NMAX, LDB = NMAX, LDE = NMAX,
$ LDQ = NMAX, LDZ = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( 1, 4*NMAX, 6*NMAX-6 ) )
* .. Local Scalars ..
CHARACTER*1 DICO, FACT, TRANS
DOUBLE PRECISION SCALE
INTEGER I, INFO, J, N, M
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX),
$ B(LDB,NMAX), BETA(NMAX), DWORK(LDWORK),
$ E(LDE,NMAX), Q(LDQ,NMAX), Z(LDZ,NMAX)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL SG03BD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. 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 = 99995 ) N
ELSE IF ( M.LT.0 .OR. M.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) M
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
IF ( LSAME( FACT, 'F' ) ) THEN
READ ( NIN, FMT = * ) ( ( Q(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( Z(I,J), J = 1,N ), I = 1,N )
END IF
IF ( LSAME( FACT, 'T' ) ) THEN
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,M )
END IF
* Find the Cholesky factor U of the solution matrix.
CALL SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, LDQ,
$ Z, LDZ, B, LDB, SCALE, ALPHAR, ALPHAI, BETA,
$ DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 ) SCALE
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,N )
20 CONTINUE
END IF
END IF
STOP
*
99999 FORMAT (' SG03BD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SG03BD = ',I2)
99997 FORMAT (' SCALE = ',F8.4,//' The Cholesky factor U of the solution
$ matrix is')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' N is out of range.',/' N = ',I5)
99994 FORMAT (/' M is out of range.',/' M = ',I5)
END
Program Data
SG03BD EXAMPLE PROGRAM DATA 3 1 C N N -1.0 3.0 -4.0 0.0 5.0 -2.0 -4.0 4.0 1.0 2.0 1.0 3.0 2.0 0.0 1.0 4.0 5.0 1.0 2.0 -1.0 7.0Program Results
SG03BD EXAMPLE PROGRAM RESULTS SCALE = 1.0000 The Cholesky factor U of the solution matrix is 1.6003 -0.4418 -0.1523 0.0000 0.6795 -0.2499 0.0000 0.0000 0.2041