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.
Specification
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,*)
Arguments
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'.
Input/Output Parameters
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.
Workspace
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).
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, 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).Numerical Aspects
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.
Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None