Purpose
To compute the Cholesky factor U of the matrix X, H X = op(U) * op(U), which is the solution of either the generalized c-stable continuous-time Lyapunov equation H H op(A) * X * op(E) + op(E) * X * op(A) 2 H = - SCALE * op(B) * op(B), (1) or the generalized d-stable discrete-time Lyapunov equation H H op(A) * X * op(A) - op(E) * X * op(E) 2 H = - SCALE * op(B) * op(B), (2) without first finding X and without the need to form the matrix op(B)**H * op(B). op(K) is either K or K**H 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 SG03BZ( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, $ LDQ, Z, LDZ, B, LDB, SCALE, ALPHA, BETA, DWORK, $ ZWORK, LZWORK, INFO ) C .. Scalar Arguments .. DOUBLE PRECISION SCALE INTEGER INFO, LDA, LDB, LDE, LDQ, LDZ, LZWORK, M, N CHARACTER DICO, FACT, TRANS C .. Array Arguments .. COMPLEX*16 A(LDA,*), ALPHA(*), B(LDB,*), BETA(*), E(LDE,*), $ Q(LDQ,*), Z(LDZ,*), ZWORK(*) DOUBLE PRECISION DWORK(*)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 (complex) 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 conjugate transposed equation is to be solved or not: = 'N': op(A) = A, op(E) = E; = 'C': op(A) = A**H, op(E) = E**H.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. If M = 0, A and E are unchanged on exit, and Q, Z, ALPHA and BETA are not set. A (input/output) COMPLEX*16 array, dimension (LDA,N) On entry, if FACT = 'F', then the leading N-by-N upper triangular 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 triangular matrix. The elements below the upper triangular part of the array A are used as workspace. If FACT = 'N', then the leading N-by-N part of this array must contain the matrix A. On exit, if FACT = 'N', the leading N-by-N upper triangular part of this array contains the generalized Schur factor A_s of the matrix A. (A_s is an upper triangular matrix.) If FACT = 'F', the leading N-by-N upper triangular part of this array is unchanged. LDA INTEGER The leading dimension of the array A. LDA >= MAX(1,N). E (input/output) COMPLEX*16 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). E_s must be an upper triangular matrix. The elements below the upper triangular part of the array E are used as workspace. If FACT = 'N', then the leading N-by-N part of this array must contain the coefficient matrix E of the equation. On exit, if FACT = 'N', the leading N-by-N upper triangular part of this array contains the generalized Schur factor E_s of the matrix E. (E_s is an upper triangular matrix.) If FACT = 'F', the leading N-by-N upper triangular part of this array is unchanged. LDE INTEGER The leading dimension of the array E. LDE >= MAX(1,N). Q (input/output) COMPLEX*16 array, dimension (LDQ,N) On entry, if FACT = 'F', then the leading N-by-N part of this array must contain the unitary matrix Q from the generalized Schur factorization (see definitions (3) and (4) in section METHOD), or an identity matrix (if the original equation has upper triangular matrices A and E). If FACT = 'N', Q need not be set on entry. On exit, if FACT = 'N', the leading N-by-N part of this array contains the unitary matrix Q from the generalized Schur factorization. If FACT = 'F', this array is unchanged. LDQ INTEGER The leading dimension of the array Q. LDQ >= MAX(1,N). Z (input/output) COMPLEX*16 array, dimension (LDZ,N) On entry, if FACT = 'F', then the leading N-by-N part of this array must contain the unitary matrix Z from the generalized Schur factorization (see definitions (3) and (4) in section METHOD), or an identity matrix (if the original equation has upper triangular matrices A and E). If FACT = 'N', Z need not be set on entry. On exit, if FACT = 'N', the leading N-by-N part of this array contains the unitary matrix Z from the generalized Schur factorization. If FACT = 'F', this array is unchanged. LDZ INTEGER The leading dimension of the array Z. LDZ >= MAX(1,N). B (input/output) COMPLEX*16 array, dimension (LDB,N1) On entry, if TRANS = 'C', 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, if INFO = 0, 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)**H * 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 = 'C', 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. ALPHA (output) COMPLEX*16 arrays, dimension (N) BETA If INFO = 0, 5, 6, or 7, then ALPHA(j)/BETA(j), j = 1, ... , N, are the eigenvalues of the matrix pencil A - lambda * E (the diagonals of the complex Schur form). All BETA(j) are non-negative real numbers. ALPHA will be always less than and usually comparable with norm(A) in magnitude, and BETA always less than and usually comparable with norm(B).Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK), where LDWORK = 0, if MIN(M,N) = 0 or FACT = 'F' and N <= 1; else, LDWORK = N-1, if FACT = 'F' and DICO = 'C'; LDWORK = MAX(N-1,10), if FACT = 'F' and DICO = 'D'; LDWORK = 8*N, if FACT = 'N'. ZWORK COMPLEX*16 array, dimension (LZWORK) On exit, if INFO = 0, ZWORK(1) returns the optimal value of LZWORK. On exit, if INFO = -21, ZWORK(1) returns the minimum value of LZWORK. LZWORK INTEGER The dimension of the array ZWORK. LZWORK >= MAX(1,3*N-3,2*N). For good performance, LZWORK should be larger. If LZWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the ZWORK array, returns this value as the first entry of the ZWORK array, and no error message related to LZWORK is issued by XERBLA.Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 4: FACT = 'N' and the pencil A - lambda * E cannot be reduced to generalized Schur form: LAPACK routine ZGGES 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 ZSTEIN utilized to factorize M3 failed to converge in the discrete-time case (see section METHOD for SLICOT Library routine SG03BS). 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 complex generalized Schur form A_s - lambda * E_s by means of unitary transformations (QZ-algorithm): A_s = Q**H * A * Z (upper triangular), (3) E_s = Q**H * 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 unitary 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 = 'C' and M < N, B_s is the upper triangular matrix from the RQ-factorization ( Q_B O ) (B_s O ) * ( ) = ( Q**H * B O ). ( O I ) 4. If TRANS = 'C' and M >= N, B_s is the upper triangular matrix from the (rectangular) RQ-factorization ( B_s O ) * Q_B = Q**H * 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 H H op(A_s) op(U_s) op(U_s) op(E_s) H H + op(E_s) op(U_s) op(U_s) op(A_s) H = - op(B_s) op(B_s) (5) or to the reduced discrete-time equation H H op(A_s) op(U_s) op(U_s) op(A_s) H H - op(E_s) op(U_s) op(U_s) op(E_s) H = - 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. We compute U11, U12**H, and U22 in three steps. Step I: From (5) and (7) we get the 1-by-1 equation H H H H A11 * U11 * U11 * E11 + E11 * U11 * U11 * A11 H = - 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 5). Furthermore, the auxiliary scalars M1 and M2 defined as follows M1 = A11 / E11 , M2 = B11 / E11 / U11 , are computed in a numerically reliable way. Step II: The generalized Sylvester equation H H H H A22 * U12 + E22 * U12 * M1 = H H H H H - B12 * M2 - A12 * U11 - E12 * U11 * M1 is solved for U12**H, as a linear system of order N-1. Step III: It can be shown that H H H H A22 * U22 * U22 * E22 + E22 * U22 * U22 * A22 = H H - B22 * B22 - y * y (8) holds, where y is defined as H H H H H 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**H ) where Q_B_tilde is a unitary matrix of order N, then H H H - B22 * B22 - y * y = - B22_tilde * B22_tilde. Replacing the right hand side in (8) by the term - B22_tilde**H * B22_tilde leads to a reduced generalized Lyapunov equation like (5), but of dimension N-1. 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**H (if TRANS = 'N'), or U * Q_U = Z * U_s (if TRANS = 'C'), where Q_U is an N-by-N unitary matrix. Again, the unitary 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)/3Further 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.Example
Program Text
NoneProgram Data
NoneProgram Results
None