**Purpose**

To solve for X either the continuous-time algebraic Riccati equation -1 Q + A'X + XA - (L+XB)R (L+XB)' = 0 (1) or the discrete-time algebraic Riccati equation -1 X = A'XA - (L+A'XB)(R + B'XB) (L+A'XB)' + Q (2) where A, B, Q, R, and L are N-by-N, N-by-M, N-by-N, M-by-M and N-by-M matrices, respectively, such that Q = C'C, R = D'D and L = C'D; X is an N-by-N symmetric matrix. The routine also returns the computed values of the closed-loop spectrum of the system, i.e., the stable eigenvalues lambda(1), ..., lambda(N) of the corresponding Hamiltonian or symplectic pencil, in the continuous-time case or discrete-time case, respectively. -1 Optionally, matrix G = BR B' may be given instead of B and R. Other options include the case with Q and/or R given in a factored form, Q = C'C, R = D'D, and with L a zero matrix. The routine uses the method of deflating subspaces, based on reordering the eigenvalues in a generalized Schur matrix pair. A standard eigenproblem is solved in the continuous-time case if G is given.

SUBROUTINE SB02OD( DICO, JOBB, FACT, UPLO, JOBL, SORT, N, M, P, A, $ LDA, B, LDB, Q, LDQ, R, LDR, L, LDL, RCOND, X, $ LDX, ALFAR, ALFAI, BETA, S, LDS, T, LDT, U, $ LDU, TOL, IWORK, DWORK, LDWORK, BWORK, INFO ) C .. Scalar Arguments .. CHARACTER DICO, FACT, JOBB, JOBL, SORT, UPLO INTEGER INFO, LDA, LDB, LDL, LDQ, LDR, LDS, LDT, LDU, $ LDWORK, LDX, M, N, P DOUBLE PRECISION RCOND, TOL C .. Array Arguments .. LOGICAL BWORK(*) INTEGER IWORK(*) DOUBLE PRECISION A(LDA,*), ALFAI(*), ALFAR(*), B(LDB,*), BETA(*), $ DWORK(*), L(LDL,*), Q(LDQ,*), R(LDR,*), $ S(LDS,*), T(LDT,*), U(LDU,*), X(LDX,*)

**Mode Parameters**

DICO CHARACTER*1 Specifies the type of Riccati equation to be solved as follows: = 'C': Equation (1), continuous-time case; = 'D': Equation (2), discrete-time case. JOBB CHARACTER*1 Specifies whether or not the matrix G is given, instead of the matrices B and R, as follows: = 'B': B and R are given; = 'G': G is given. FACT CHARACTER*1 Specifies whether or not the matrices Q and/or R (if JOBB = 'B') are factored, as follows: = 'N': Not factored, Q and R are given; = 'C': C is given, and Q = C'C; = 'D': D is given, and R = D'D; = 'B': Both factors C and D are given, Q = C'C, R = D'D. UPLO CHARACTER*1 If JOBB = 'G', or FACT = 'N', specifies which triangle of the matrices G and Q (if FACT = 'N'), or Q and R (if JOBB = 'B'), is stored, as follows: = 'U': Upper triangle is stored; = 'L': Lower triangle is stored. JOBL CHARACTER*1 Specifies whether or not the matrix L is zero, as follows: = 'Z': L is zero; = 'N': L is nonzero. JOBL is not used if JOBB = 'G' and JOBL = 'Z' is assumed. SLICOT Library routine SB02MT should be called just before SB02OD, for obtaining the results when JOBB = 'G' and JOBL = 'N'. SORT CHARACTER*1 Specifies which eigenvalues should be obtained in the top of the generalized Schur form, as follows: = 'S': Stable eigenvalues come first; = 'U': Unstable eigenvalues come first.

N (input) INTEGER The actual state dimension, i.e. the order of the matrices A, Q, and X, and the number of rows of the matrices B and L. N >= 0. M (input) INTEGER The number of system inputs. If JOBB = 'B', M is the order of the matrix R, and the number of columns of the matrix B. M >= 0. M is not used if JOBB = 'G'. P (input) INTEGER The number of system outputs. If FACT = 'C' or 'D' or 'B', P is the number of rows of the matrices C and/or D. P >= 0. Otherwise, P is not used. A (input) DOUBLE PRECISION array, dimension (LDA,N) The leading N-by-N part of this array must contain the state matrix A of the system. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N). B (input) DOUBLE PRECISION array, dimension (LDB,*) If JOBB = 'B', the leading N-by-M part of this array must contain the input matrix B of the system. If JOBB = 'G', the leading N-by-N upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the upper triangular part or lower triangular part, respectively, of the matrix -1 G = BR B'. The stricly lower triangular part (if UPLO = 'U') or stricly upper triangular part (if UPLO = 'L') is not referenced. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). Q (input) DOUBLE PRECISION array, dimension (LDQ,N) If FACT = 'N' or 'D', the leading N-by-N upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the upper triangular part or lower triangular part, respectively, of the symmetric state weighting matrix Q. The stricly lower triangular part (if UPLO = 'U') or stricly upper triangular part (if UPLO = 'L') is not referenced. If JOBB = 'B', the triangular part of this array defined by UPLO is modified internally, but is restored on exit. If FACT = 'C' or 'B', the leading P-by-N part of this array must contain the output matrix C of the system. If JOBB = 'B', this part is modified internally, but is restored on exit. LDQ INTEGER The leading dimension of array Q. LDQ >= MAX(1,N) if FACT = 'N' or 'D', LDQ >= MAX(1,P) if FACT = 'C' or 'B'. R (input) DOUBLE PRECISION array, dimension (LDR,M) If FACT = 'N' or 'C', the leading M-by-M upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the upper triangular part or lower triangular part, respectively, of the symmetric input weighting matrix R. The stricly lower triangular part (if UPLO = 'U') or stricly upper triangular part (if UPLO = 'L') is not referenced. The triangular part of this array defined by UPLO is modified internally, but is restored on exit. If FACT = 'D' or 'B', the leading P-by-M part of this array must contain the direct transmission matrix D of the system. This part is modified internally, but is restored on exit. If JOBB = 'G', this array is not referenced. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,M) if JOBB = 'B' and FACT = 'N' or 'C'; LDR >= MAX(1,P) if JOBB = 'B' and FACT = 'D' or 'B'; LDR >= 1 if JOBB = 'G'. L (input) DOUBLE PRECISION array, dimension (LDL,M) If JOBL = 'N' (and JOBB = 'B'), the leading N-by-M part of this array must contain the cross weighting matrix L. This part is modified internally, but is restored on exit. If JOBL = 'Z' or JOBB = 'G', this array is not referenced. LDL INTEGER The leading dimension of array L. LDL >= MAX(1,N) if JOBL = 'N' and JOBB = 'B'; LDL >= 1 if JOBL = 'Z' or JOBB = 'G'. RCOND (output) DOUBLE PRECISION An estimate of the reciprocal of the condition number (in the 1-norm) of the N-th order system of algebraic equations from which the solution matrix X is obtained. X (output) DOUBLE PRECISION array, dimension (LDX,N) The leading N-by-N part of this array contains the solution matrix X of the problem. LDX INTEGER The leading dimension of array X. LDX >= MAX(1,N). ALFAR (output) DOUBLE PRECISION array, dimension (2*N) ALFAI (output) DOUBLE PRECISION array, dimension (2*N) BETA (output) DOUBLE PRECISION array, dimension (2*N) The generalized eigenvalues of the 2N-by-2N matrix pair, ordered as specified by SORT (if INFO = 0). For instance, if SORT = 'S', the leading N elements of these arrays contain the closed-loop spectrum of the system matrix A - BF, where F is the optimal feedback matrix computed based on the solution matrix X. Specifically, lambda(k) = [ALFAR(k)+j*ALFAI(k)]/BETA(k) for k = 1,2,...,N. If DICO = 'C' and JOBB = 'G', the elements of BETA are set to 1. S (output) DOUBLE PRECISION array, dimension (LDS,*) The leading 2N-by-2N part of this array contains the ordered real Schur form S of the first matrix in the reduced matrix pencil associated to the optimal problem, or of the corresponding Hamiltonian matrix, if DICO = 'C' and JOBB = 'G'. That is, (S S ) ( 11 12) S = ( ), (0 S ) ( 22) where S , S and S are N-by-N matrices. 11 12 22 Array S must have 2*N+M columns if JOBB = 'B', and 2*N columns, otherwise. LDS INTEGER The leading dimension of array S. LDS >= MAX(1,2*N+M) if JOBB = 'B', LDS >= MAX(1,2*N) if JOBB = 'G'. T (output) DOUBLE PRECISION array, dimension (LDT,2*N) If DICO = 'D' or JOBB = 'B', the leading 2N-by-2N part of this array contains the ordered upper triangular form T of the second matrix in the reduced matrix pencil associated to the optimal problem. That is, (T T ) ( 11 12) T = ( ), (0 T ) ( 22) where T , T and T are N-by-N matrices. 11 12 22 If DICO = 'C' and JOBB = 'G' this array is not referenced. LDT INTEGER The leading dimension of array T. LDT >= MAX(1,2*N+M) if JOBB = 'B', LDT >= MAX(1,2*N) if JOBB = 'G' and DICO = 'D', LDT >= 1 if JOBB = 'G' and DICO = 'C'. U (output) DOUBLE PRECISION array, dimension (LDU,2*N) The leading 2N-by-2N part of this array contains the right transformation matrix U which reduces the 2N-by-2N matrix pencil to the ordered generalized real Schur form (S,T), or the Hamiltonian matrix to the ordered real Schur form S, if DICO = 'C' and JOBB = 'G'. That is, (U U ) ( 11 12) U = ( ), (U U ) ( 21 22) where U , U , U and U are N-by-N matrices. 11 12 21 22 LDU INTEGER The leading dimension of array U. LDU >= MAX(1,2*N).

TOL DOUBLE PRECISION The tolerance to be used to test for near singularity of the original matrix pencil, specifically of the triangular factor obtained during the reduction process. If the user sets TOL > 0, then the given value of TOL is used as a lower bound for the reciprocal condition number of that matrix; a matrix whose estimated condition number is less than 1/TOL is considered to be nonsingular. If the user sets TOL <= 0, then a default tolerance, defined by TOLDEF = EPS, is used instead, where EPS is the machine precision (see LAPACK Library routine DLAMCH). This parameter is not referenced if JOBB = 'G'.

IWORK INTEGER array, dimension (LIWORK) LIWORK >= MAX(1,M,2*N) if JOBB = 'B', LIWORK >= MAX(1,2*N) if JOBB = 'G'. DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK. If JOBB = 'B' and N > 0, DWORK(2) returns the reciprocal of the condition number of the M-by-M lower triangular matrix obtained after compressing the matrix pencil of order 2N+M to obtain a pencil of order 2N. If INFO = 0 or INFO = 6, DWORK(3) returns the scaling factor used internally, which should multiply the submatrix Y2 to recover X from the first N columns of U (see METHOD). LDWORK INTEGER The length of the array DWORK. LDWORK >= MAX(3,6*N), if JOBB = 'G', DICO = 'C'; LDWORK >= MAX(7*(2*N+1)+16,16*N), if JOBB = 'G', DICO = 'D'; LDWORK >= MAX(7*(2*N+1)+16,16*N,2*N+M,3*M), if JOBB = 'B'. For optimum performance LDWORK should be larger. BWORK LOGICAL array, dimension (2*N)

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: if the computed extended matrix pencil is singular, possibly due to rounding errors; = 2: if the QZ (or QR) algorithm failed; = 3: if reordering of the (generalized) eigenvalues failed; = 4: if after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the (generalized) Schur form no longer satisfy the stability condition; this could also be caused due to scaling; = 5: if the computed dimension of the solution does not equal N; = 6: if a singular matrix was encountered during the computation of the solution matrix X.

The routine uses a variant of the method of deflating subspaces proposed by van Dooren [1]. See also [2], [3]. It is assumed that (A,B) is stabilizable and (C,A) is detectable. Under these assumptions the algebraic Riccati equation is known to have a unique non-negative definite solution. The first step in the method of deflating subspaces is to form the extended Hamiltonian matrices, dimension 2N + M given by discrete-time continuous-time |A 0 B| |I 0 0| |A 0 B| |I 0 0| |Q -I L| - z |0 -A' 0|, |Q A' L| - s |0 -I 0|. |L' 0 R| |0 -B' 0| |L' B' R| |0 0 0| Next, these pencils are compressed to a form (see [1]) lambda x A - B . f f This generalized eigenvalue problem is then solved using the QZ algorithm and the stable deflating subspace Ys is determined. If [Y1'|Y2']' is a basis for Ys, then the required solution is -1 X = Y2 x Y1 . A standard eigenvalue problem is solved using the QR algorithm in the continuous-time case when G is given (DICO = 'C', JOBB = 'G').

[1] Van Dooren, P. A Generalized Eigenvalue Approach for Solving Riccati Equations. SIAM J. Sci. Stat. Comp., 2, pp. 121-135, 1981. [2] Mehrmann, V. The Autonomous Linear Quadratic Control Problem. Theory and Numerical Solution. Lect. Notes in Control and Information Sciences, vol. 163, Springer-Verlag, Berlin, 1991. [3] Sima, V. Algorithms for Linear-Quadratic Optimization. Pure and Applied Mathematics: A Series of Monographs and Textbooks, vol. 200, Marcel Dekker, Inc., New York, 1996.

This routine is particularly suited for systems where the matrix R is ill-conditioned. Internal scaling is used.

To obtain a stabilizing solution of the algebraic Riccati equations set SORT = 'S'. The routine can also compute the anti-stabilizing solutions of the algebraic Riccati equations, by specifying SORT = 'U'.

**Program Text**

* SB02OD EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER NMAX, MMAX, PMAX PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 ) INTEGER NMAX2M, NMAX2 PARAMETER ( NMAX2M = 2*NMAX+MMAX, NMAX2 = 2*NMAX ) INTEGER LDA, LDB, LDL, LDQ, LDR, LDS, LDT, LDU, LDX PARAMETER ( LDA = NMAX, LDB = NMAX, LDL = NMAX, $ LDQ = MAX(NMAX,PMAX), LDR = MAX(MMAX,PMAX), $ LDS = NMAX2M, LDT = NMAX2M, LDU = NMAX2, $ LDX = NMAX ) INTEGER LIWORK PARAMETER ( LIWORK = MAX(MMAX,NMAX2) ) INTEGER LDWORK PARAMETER ( LDWORK = MAX(14*NMAX+23,16*NMAX) ) INTEGER LBWORK PARAMETER ( LBWORK = NMAX2 ) * .. Local Scalars .. DOUBLE PRECISION RCOND, TOL INTEGER I, INFO, J, M, N, P CHARACTER*1 DICO, FACT, JOBB, JOBL, SORT, UPLO LOGICAL LJOBB * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), ALFAI(NMAX2), ALFAR(NMAX2), $ B(LDB,MMAX), BETA(NMAX2), DWORK(LDWORK), $ L(LDL,MMAX), Q(LDQ,NMAX), R(LDR,MMAX), $ S(LDS,NMAX2M), T(LDT,NMAX2), U(LDU,NMAX2), $ X(LDX,NMAX) INTEGER IWORK(LIWORK) LOGICAL BWORK(LBWORK) C .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL SB02OD * .. 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, P, TOL, DICO, JOBB, FACT, UPLO, JOBL, $ SORT IF ( N.LT.0 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99995 ) N ELSE READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N ) IF ( M.LT.0 .OR. M.GT.MMAX ) THEN WRITE ( NOUT, FMT = 99994 ) M ELSE LJOBB = LSAME( JOBB, 'B' ) IF ( LJOBB ) 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,N ) END IF IF ( P.LT.0 .OR. P.GT.PMAX ) THEN WRITE ( NOUT, FMT = 99993 ) P ELSE IF ( LSAME( FACT, 'N' ) .OR. LSAME( FACT, 'D' ) ) THEN READ ( NIN, FMT = * ) $ ( ( Q(I,J), J = 1,N ), I = 1,N ) ELSE READ ( NIN, FMT = * ) $ ( ( Q(I,J), J = 1,N ), I = 1,P ) END IF IF ( LJOBB ) THEN IF ( LSAME( FACT, 'N' ) .OR. LSAME( FACT, 'C' ) ) THEN READ ( NIN, FMT = * ) $ ( ( R(I,J), J = 1,M ), I = 1,M ) ELSE READ ( NIN, FMT = * ) $ ( ( R(I,J), J = 1,M ), I = 1,P ) END IF IF ( LSAME( JOBL, 'N' ) ) $ READ ( NIN, FMT = * ) $ ( ( L(I,J), J = 1,M ), I = 1,N ) END IF * Find the solution matrix X. CALL SB02OD( DICO, JOBB, FACT, UPLO, JOBL, SORT, N, M, P, $ A, LDA, B, LDB, Q, LDQ, R, LDR, L, LDL, $ RCOND, X, LDX, ALFAR, ALFAI, BETA, S, LDS, $ T, LDT, U, LDU, TOL, IWORK, DWORK, LDWORK, $ BWORK, INFO ) * IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE WRITE ( NOUT, FMT = 99997 ) DO 20 I = 1, N WRITE ( NOUT, FMT = 99996 ) ( X(I,J), J = 1,N ) 20 CONTINUE END IF END IF END IF END IF STOP * 99999 FORMAT (' SB02OD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from SB02OD = ',I2) 99997 FORMAT (' The solution matrix X 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) 99993 FORMAT (/' P is out of range.',/' P = ',I5) END

SB02OD EXAMPLE PROGRAM DATA 2 1 3 0.0 C B B U Z S 0.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0

SB02OD EXAMPLE PROGRAM RESULTS The solution matrix X is 1.7321 1.0000 1.0000 1.7321