Purpose
To compute the residual matrix R for a continuous-time or discrete-time "reduced" Lyapunov equation, using the formulas R = op(A)'*X + X*op(A) + Q, or R = op(A)'*X*op(E) + op(E)'*X*op(A) + Q, in the continuous-time case, or the formulas R = op(A)'*X*op(A) - X + Q, or R = op(A)'*X*op(A) - op(E)'*X*op(E) + Q, in the discrete-time case, where X and Q are symmetric matrices, A is in upper real Schur form, E is upper triangular, and op(W) is op(W) = W or op(W) = W'. Optionally, the Frobenius norms of the product terms defining the denominator of the relative residual are also computed. The norms of Q and X are not computed.Specification
SUBROUTINE SG02CV( DICO, JOB, JOBE, UPLO, TRANS, N, A, LDA, E, $ LDE, X, LDX, R, LDR, NORMS, DWORK, LDWORK, $ INFO ) C .. Scalar Arguments .. CHARACTER DICO, JOB, JOBE, TRANS, UPLO INTEGER INFO, LDA, LDE, LDR, LDWORK, LDX, N C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), DWORK(*), E(LDE,*), NORMS(*), $ R(LDR,*), X(LDX,*)Arguments
Mode Parameters
DICO CHARACTER*1 Specifies the type of the Lyapunov equation, as follows: = 'C': continuous-time Lyapunov equation; = 'D': discrete-time Lyapunov equation. JOB CHARACTER*1 Specifies which results must be computed, as follows: = 'R': The matrix R only must be computed; = 'N': The matrix R and the norms must be computed; = 'B': The matrix R and the norms must be computed. JOBE CHARACTER*1 Specifies whether E is a general or an identity matrix, as follows: = 'G': The matrix E is general and is given; = 'I': The matrix E is assumed identity and is not given. UPLO CHARACTER*1 Specifies which triangles of the symmetric matrices X and Q are given, as follows: = 'U': The upper triangular part is given; = 'L': The lower triangular part is given. TRANS CHARACTER*1 Specifies the form of op(W) to be used in the formulas above, as follows: = 'N': op(W) = W; = 'T': op(W) = W'; = 'C': op(W) = W'.Input/Output Parameters
N (input) INTEGER The order of the matrices A, E, Q, X, and R. N >= 0. A (input) DOUBLE PRECISION array, dimension (LDA,N) The leading N-by-N upper part of this array must contain the upper real Schur matrix A. If TRANS = 'N' and (DICO = 'D' or (JOB = 'R' and JOBE = 'G')), the entries 3, 4,..., N of the first column are modified internally, but are restored on exit. Otherwise, the part of this array below the first subdiagonal is not referenced. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N). E (input) DOUBLE PRECISION array, dimension (LDE,*) If JOBE = 'G', the leading N-by-N upper triangular part of this array must contain the upper triangular matrix E. The strictly lower triangular part of this array is not referenced. If JOBE = 'I', this array is not referenced. LDE INTEGER The leading dimension of array E. LDE >= MAX(1,N), if JOBE = 'G'; LDE >= 1, if JOBE = 'I'. X (input/works.) DOUBLE PRECISION array, dimension (LDX,N) On entry, if UPLO = 'U', the leading N-by-N upper triangular part of this array must contain the upper triangular part of the symmetric matrix X and the strictly lower triangular part of the array is not referenced. On entry, if UPLO = 'L', the leading N-by-N lower triangular part of this array must contain the lower triangular part of the symmetric matrix X and the strictly upper triangular part of the array is not referenced. If DICO = 'D' or (JOB = 'R' and JOBE = 'G'), the diagonal elements of this array are modified internally, but they are restored on exit. LDX INTEGER The leading dimension of array X. LDX >= MAX(1,N). R (input/output) DOUBLE PRECISION array, dimension (LDR,*) On entry, the leading N-by-N upper or lower triangular part (depending on UPLO) of this array must contain the upper or lower triangular part, respectively, of the matrix Q. The other strictly triangular part is not referenced. On exit, the leading N-by-N upper or lower triangular part (depending on UPLO) of this array contains the upper or lower triangular part, respectively, of the matrix R. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,N). NORMS (output) DOUBLE PRECISION array, dimension (LN) If JOB = 'N' or JOB = 'B', LN = 1 or 2, if (DICO = 'C' or JOBE = 'I'), or (DICO = 'D' and JOBE = 'G'), respectively. If DICO = 'C', NORMS(1) contains the Frobenius norm of the matrix op(A)'*X (or of X*op(A)), if JOBE = 'I', or of the matrix op(A)'*X*op(E) (or of op(E)'*X*op(A)), if JOBE = 'G'. If DICO = 'D', NORMS(1) contains the Frobenius norm of the matrix op(A)'*X*op(A); if JOBE = 'G', NORMS(2) contains the Frobenius norm of the matrix op(E)'*X*op(E). If JOB <> 'N', this array is not referenced.Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = -17 or if LDWORK = -2 on input, then DWORK(1) returns the minimum value of LDWORK. On exit, if INFO = 0, or if LDWORK = -1 on input, then DWORK(1) returns the optimal value of LDWORK. LDWORK The length of the array DWORK. LDWORK >= MAX(v,1), with v specified in the following table, where a = 1, if JOBE = 'G'; a = 0, if JOBE = 'I'. DICO JOB v ---------------------------- 'C' 'R' a*N*N 'C' 'N','B' N*N ---------------------------- 'D' 'R' N*N 'D' 'N','B' 2*N*N ---------------------------- If LDWORK = -1, an optimal 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 is issued by XERBLA. If LDWORK = -2, a minimal workspace query is assumed; the routine only calculates the minimal size of the DWORK array, returns this value as the first entry of the DWORK array, and no error message is issued by XERBLA.Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.Method
The matrix expressions are efficiently evaluated, using symmetry. If JOB = 'N' or JOB = 'B', then: If DICO = 'C', the matrices op(op(A)'*X*op(E)) or op(X*op(A)), are efficiently computed. If DICO = 'D', the matrices op(A)'*X*op(A) and op(E)'*X*op(E), if JOBE = 'G', are efficiently computed. The results are used to evaluate R and the norms. If JOB = 'R', then the needed parts of the intermediate results are obtained and used to evaluate R.Numerical Aspects
The calculations are backward stable. The algorithm requires approximately a*N^3 operations, where ^ denotes the power operator, and a = 1, if DICO = 'C' and JOB <> 'R' and JOBE = 'G'; a = 1/2, otherwise. An "operation" includes a multiplication, an addition, and some address calculations.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None