**Purpose**

To compute the residual matrix R for a continuous-time or discrete-time Riccati equation and/or the "closed-loop system" matrix op(C), using the formulas R = op(A)'*X + X*op(A) +/- X*G*X + Q, C = op(A) +/- G*X, or R = op(A)'*X*op(E) + op(E)'*X*op(A) +/- op(E)'*X*G*X*op(E) + Q, C = op(A) +/- G*X*op(E), or R = op(A)'*X*op(E) + op(E)'*X*op(A) +/- H*K + Q, C = op(A) +/- B*K, in the continuous-time case, or the formulas R = op(A)'*X*op(A) - X +/- op(A)'*X*G*X*op(A) + Q, C = op(A) +/- G*X*op(A), or R = op(A)'*X*op(A) - op(E)'*X*op(E) +/- op(A)'*X*G*X*op(A) + Q, C = op(A) +/- G*X*op(A), or R = op(A)'*X*op(A) - op(E)'*X*op(E) +/- H*K + Q, C = op(A) +/- B*K, in the discrete-time case, where X, G, and Q are symmetric matrices, A, E, H, K, B are general matrices, and op(W) is one of op(W) = W or op(W) = W'. _-1 Instead of the symmetric N-by-N matrix G, G = B*R *B', the N-by-M _-1 matrix D, D = B*L , such that G = D*D', may be given on entry. _ _ _ _ The matrix R, with R = L'*L, is a weighting matrix of the optimal _ _ problem, if DICO = 'C', or it is R = B'*X*B + Rd, if DICO = 'D', _ _ _ with Rd a similar weighting matrix; L is a Cholesky factor of R, _ _ if R is positive definite. If R is not positive definite, which may happen in the discrete-time case, a UdU' or LdL' factorization is used to compute the matrices H and K. If M = 0, the residual matrix of a (generalized) Lyapunov or Stein equation is computed. To this end, set JOBG = 'D' and JOB = 'R' (since op(C) = A in this case). Optionally, the quadratic term in the formulas for R is specified as H*K, where H = L + op(E)'*X*B, if DICO = 'C', or H = L + op(A)'*X*B, if DICO = 'D', and _-1 K = R *H', with L an N-by-M matrix. This is useful, e.g., for DICO = 'D', _ _ when L <> 0 and/or Rd is singular, hence R might be numerically indefinite; it might be indefinite in the first iterations of Newton's algorithm. Depending on JOB, part or all of the matrices H, K, and B should be given in such a case. _ If R is positive definite, the quadratic term can be specified as F*F', and the second term in the formulas for C is D*F', where _-1 F = H*L . The matrices F and/or D should be given. This option is not useful when L = 0, unless F and D are available. If DICO = 'C', the computational problem with L <> 0 is equivalent with one with L = 0 after replacing _-1 _-1 A := A - B*R* L', Q := Q - L*R* L'. _ _ These formulas, with R replaced by Rd, can also be used in the _ discrete-time case, if Rd is nonsingular and well-conditioned with respect to inversion. 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.

SUBROUTINE SG02CW( DICO, JOB, JOBE, FLAG, JOBG, UPLO, TRANS, N, M, $ A, LDA, E, LDE, G, LDG, X, LDX, F, LDF, K, LDK, $ XE, LDXE, R, LDR, C, LDC, NORMS, DWORK, LDWORK, $ INFO ) C .. Scalar Arguments .. CHARACTER DICO, FLAG, JOB, JOBE, JOBG, TRANS, UPLO INTEGER INFO, LDA, LDC, LDE, LDF, LDG, LDK, LDR, LDWORK, $ LDX, LDXE, M, N C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), C(LDC,*), DWORK(*), E(LDE,*), $ F(LDF,*), G(LDG,*), K(LDK,*), NORMS(*), $ R(LDR,*), X(LDX,*), XE(LDXE,*)

**Mode Parameters**

DICO CHARACTER*1 Specifies the type of the Riccati equation, as follows: = 'C': continuous-time algebraic Riccati equation; = 'D': discrete-time algebraic Riccati equation. JOB CHARACTER*1 Specifies which results must be computed, as follows: = 'A': Both (all) matrices R and C must be computed; = 'R': The matrix R only must be computed; = 'C': The matrix C only must be computed; = 'N': The matrices R and C 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. FLAG CHARACTER*1 Specifies which sign is used, as follows: = 'P': The plus sign is used; = 'M': The minus sign is used. JOBG CHARACTER*1 Specifies how the quadratic term in the formulas for R is defined, as follows: = 'G': The matrix G is given; = 'D': The matrix D is given; = 'F': The matrix F is given; = 'H': The matrices H and K are given. UPLO CHARACTER*1 Specifies which triangles of the symmetric matrices X, G (if JOBG = 'G'), and Q (if JOB <> 'C') 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'.

N (input) INTEGER The order of the matrices A, E, Q, X, C and R. N >= 0. M (input) INTEGER If JOBG <> 'G', the number of columns of the matrices D, F, and/or B, H, and K'. M >= 0. If JOBG = 'G', the value of M is meaningless. A (input) DOUBLE PRECISION array, dimension (LDA,N) The leading N-by-N part of this array must contain the matrix A. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N). E (input) DOUBLE PRECISION array, dimension (LDE,*) If JOBE = 'G' and (JOB <> 'C' or (DICO = 'C' and (JOBG = 'G' or JOBG = 'D'))), the leading N-by-N part of this array must contain the matrix E. If JOBE = 'I' or (JOB = 'C' and (DICO = 'D' or JOBG = 'F' or JOBG = 'H')), this array is not referenced. LDE INTEGER The leading dimension of array E. LDE >= MAX(1,N), if JOBE = 'G' and (JOB <> 'C' or (DICO = 'C' and (JOBG = 'G' or JOBG = 'D'))); LDE >= 1, if JOBE = 'I' or (JOB = 'C' and (DICO = 'D' or JOBG = 'F' or JOBG = 'H')). G (input/works.) DOUBLE PRECISION array, dimension (LDG,*) If JOBG = 'G', 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 G. The other strictly triangular part is not referenced. If DICO = 'D', (JOB = 'R' or JOB = 'B'), and JOBG = 'G', the diagonal elements of this array are modified internally, but are restored on exit. If JOBG = 'D' or (JOBG = 'F' and JOB <> 'R' and JOB <> 'B'), the leading N-by-M part of this array must contain the matrix D, so that G = D*D'. If JOBG = 'H' and JOB <> 'R' and JOB <> 'B', the leading N-by-M part of this array must contain the matrix B. If (JOBG = 'F' or JOBG = 'H') and JOB = 'R' or JOB = 'B', this array is not referenced. LDG INTEGER The leading dimension of array G. LDG >= MAX(1,N), if JOBG = 'G' or JOBG = 'D' or (JOB <> 'R' and JOB <> 'B'); LDG >= 1, if (JOBG = 'F' or JOBG = 'H') and (JOB = 'R' or JOB = 'B'). X (input/works.) DOUBLE PRECISION array, dimension (LDX,N) The leading N-by-N part of this array must contain the symmetric matrix X, and it is unchanged on exit. If DICO = 'D', JOBE = 'G' and JOB <> 'C', the diagonal elements of this array are modified internally, but they are restored on exit. The full matrix X should be input if DICO = 'C', JOBE = 'I', and the conditions in the lines of the table below are satisfied JOBG JOB LDWORK ---------------------------------------------- 'F','H' 'A','R' LDWORK < N*N 'G' 'A','R','N' LDWORK < 2*N*N 'G' 'C' LDWORK < N*N 'G' 'B' LDWORK < 3*N*N 'D' 'R' (M<=N, LDWORK < N*N) or (M> N, LDWORK < 3*N*N) 'D' 'A' (M<=N, LDWORK < N*N) or (LDWORK >= N*N and LDWORK < 2*N*N) ---------------------------------------------- For all the other cases, including when the optimal length of the workspace array DWORK is used, only the relevant upper or lower triangular part (depending on UPLO) of this array must be input, and the other strictly triangular part is not referenced. LDX INTEGER The leading dimension of array X. LDX >= MAX(1,N). F (input) DOUBLE PRECISION array, dimension (LDF,*) If JOBG = 'F', the leading N-by-M part of this array must contain the matrix F. If JOBG = 'H', the leading N-by-M part of this array must contain the matrix H. If JOBG = 'G' or JOBG = 'D', this array is not referenced. LDF INTEGER The leading dimension of array F. LDF >= MAX(1,N), if JOBG = 'F' or JOBG = 'H'; LDF >= 1, if JOBG = 'G' or JOBG = 'D'. K (input) DOUBLE PRECISION array, dimension (LDK,*) If JOBG = 'H', the leading M-by-N part of this array must contain the matrix K. If JOBG <> 'H', this array is not referenced. LDK INTEGER The leading dimension of array K. LDK >= MAX(1,M), if JOBG = 'H'; LDK >= 1, if JOBG <> 'H'. XE (input) DOUBLE PRECISION array, dimension (LDXE,*) If (JOBG = 'F' or JOBG = 'H'), JOB <> 'C', DICO = 'C', and JOBE = 'G', the leading N-by-N part of this array must contain the matrix product X*E, if TRANS = 'N', or E*X, if TRANS = 'T' or 'C'. If (JOBG = 'F' or JOBG = 'H'), JOB <> 'C', and DICO = 'D', the leading N-by-N part of this array must contain the matrix product X*A, if TRANS = 'N', or A*X, if TRANS = 'T' or 'C'. These matrix products are needed for computing F or H. If JOBG = 'G' or JOBG = 'D' or JOB = 'C' or (DICO = 'C' and JOBE = 'I') this array is not referenced. LDXE INTEGER The leading dimension of array XE. LDXE >= MAX(1,N), if (JOBG = 'F' or JOBG = 'H'), JOB <> 'C', and either DICO = 'C' and JOBE = 'G', or DICO = 'D'; LDXE >= 1, if JOBG = 'G' or JOBG = 'D' or JOB = 'C' or (DICO = 'C' and JOBE = 'I'). R (input/output) DOUBLE PRECISION array, dimension (LDR,*) On entry, if JOB <> 'C', 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, if JOB <> 'C' and INFO = 0, 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. If JOB = 'C', this array is not referenced. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,N), if JOB <> 'C'; LDR >= 1, if JOB = 'C'. C (output) DOUBLE PRECISION array, dimension (LDC,*) If JOB <> 'R' and JOB <> 'B' and INFO = 0, the leading N-by-N part of this array contains the matrix op(C). If JOB = 'R' or JOB = 'B', this array is not referenced. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,N), if JOB <> 'R' and JOB <> 'B'; LDC >= 1, if JOB = 'R' or JOB = 'B'. NORMS (output) DOUBLE PRECISION array, dimension (LN) If JOB = 'N' or JOB = 'B', LN = 2 or 3, 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'; NORMS(2) contains the Frobenius norm of the matrix product X*G*X, if JOBE = 'I', or of the matrix product V = op(E)'*X*G*X*op(E), if JOBE = 'G' (for JOBG = 'G' or JOBG = 'D'), or of V = F*F', if JOBG = 'F', or of V = H*K, if JOBG = 'H'. If DICO = 'D', NORMS(1) contains the Frobenius norm of the matrix op(A)'*X*op(A); NORMS(2) contains the Frobenius norm of the matrix product V = op(A)'*X*G*X*op(A), if JOBG = 'G' or JOBG = 'D', or of V = F*F', if JOBG = 'F', or of V = H*K, if JOBG = 'H'; if JOBE = 'G', NORMS(3) contains the Frobenius norm of the matrix product op(E)'*X*op(E). If JOB <> 'N' and JOB <> 'B', this array is not referenced.

DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = -30, 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 JOBG JOB v ----------------------------------------------- 'C' 'F','H' 'A','C','R' 0 'C' 'F','H' 'N' a*N*N 'C' 'F','H' 'B' N*N 'C' 'G' 'A','C' a*N*N 'C' 'G' 'N','R' (a+1)*N*N 'C' 'G' 'B' (a+2)*N*N 'C' 'D' 'A' N*MIN(M,(a+1)*N) 'C' 'D' 'C' N*MIN(N,M) 'C' 'D' 'N' N*(N+MIN(a*N,M)) 'C' 'D' 'B' N*(N+MIN(N+a*N,M)) 'C' 'D' 'R' N*MIN(a*N+M,(a+2)*N) ----------------------------------------------- 'D' 'F','H' 'A','C' 0 'D' 'F','H' 'N','R' a*N*N 'D' 'F','H' 'B' (a+1)*N*N 'D' 'G' 'A','C' N*N 'D' 'G' 'N','R' 2*N*N 'D' 'G' 'B' 3*N*N 'D' 'D' 'A','N' N*MIN(MAX(N,M),2*N) 'D' 'D' 'B' N*(N+MAX(N,M)) 'D' 'D' 'C' N*MIN(N,M) 'D' 'D' 'R' N*MIN(3*N,N+M) ----------------------------------------------- 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. This evaluation assumes that only the specified triangle of the array X is always used, and the other strict triangle is not referenced. 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. This evaluation assumes that full matrix is given in the array X, when needed (see the table at the description of the array X).

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.

The matrix expressions are efficiently evaluated, using symmetry, common matrix subexpressions, and proper order of matrix multiplications. If JOB = 'N' or JOB = 'B', then: If DICO = 'C', the matrices op(op(A)'*X*op(E)) or op(X*op(A)), and V = op(E)'*X*G*X*op(E) or V = F*F' or V = H*K, are efficiently computed. If DICO = 'D', the matrices op(A)'*X*op(A), V = op(A)'*X*G*X*op(A) or V = F*F' or V = H*K, and op(E)'*X*op(E), if JOBE = 'G', are efficiently computed. The results are used to evaluate R, op(C) (if JOB = 'N'), and the norms. If JOB <> 'N', then the needed parts of the intermediate results are obtained and used to evaluate R and/or op(C).

The calculations are backward stable. 3 2 The algorithm requires approximately (a+b)N + cN M operations, where, a = 0, if JOBE = 'I', a = 1, if JOBE = 'G' and (DICO = 'C' or (DICO = 'D' and JOB = 'C')), a = 1.5, if JOBE = 'G' and DICO = 'D', and b and c are implicitly defined below. Specifically, the effort is approximately as follows (using ^ to denote the power operator) For DICO = 'C': JOBG JOB C R A, N, B 'G' (a+1)*N^3 (a+2)*N^3 (a+2)*N^3,(a+2.5)*N^3 'D' (a+2)*N^2*M (a+1)*N^3+1.5*N^2*M (a+1)*N^3+2.5*N^2*M 'F','H' N^2*M N^3+0.5*N^2*M N^3+1.5*N^2*M For DICO = 'D': JOBG JOB C R A, N, B 'G' 2*N^3 (a+3 )*N^3 (a+3 )*N^3 'D' 3*N^2*M (a+1.5)*N^3+1.5*N^2*M (a+1.5)*N^3+2.5*N^2*M 'F','H' N^2*M (a+0.5)*N^3+0.5*N^2*M (a+0.5)*N^3+1.5*N^2*M For JOBG <> 'G' and JOB = 'B', the effort reduces by N^2*M in both tables. An "operation" includes a multiplication, an addition, and some address calculations.

None

**Program Text**

None

None

None