**Purpose**

To find the line search parameter alpha minimizing the Frobenius norm (in a MATLAB-style notation) P(alpha) := norm(R(X+alpha*S), 'fro') = norm((1-alpha)*R(X) +/- alpha^2*V, 'fro'), where R(X) is the residual of a (generalized) continuous-time algebraic Riccati equation 0 = op(A)'*X + X*op(A) +/- X*G*X + Q =: R(X), or 0 = op(A)'*X*op(E) + op(E)'*X*op(A) +/- op(E)'*X*G*X*op(E) + Q =: R(X), V = op(E)'*S*G*S*op(E), and op(W) is either W or W'. The matrix S is the Newton step. _-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, R = L'*L, is a weighting matrix of the optimal problem, and L is its (Cholesky) factor. Optionally, V is specified as V = H*K, or V = F*F', but F or H and K must be evaluated in S. See the SLICOT Library routine SG02CW description for more details.

SUBROUTINE SG02CX( JOBE, FLAG, JOBG, UPLO, TRANS, N, M, E, LDE, R, $ LDR, S, LDS, G, LDG, ALPHA, RNORM, DWORK, $ LDWORK, IWARN, INFO ) C .. Scalar Arguments .. CHARACTER FLAG, JOBE, JOBG, TRANS, UPLO INTEGER INFO, IWARN, LDE, LDG, LDR, LDS, LDWORK, M, N DOUBLE PRECISION ALPHA, RNORM C .. Array Arguments .. DOUBLE PRECISION DWORK(*), E(LDE,*), G(LDG,*), R(LDR,*), S(LDS,*)

**Mode Parameters**

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 matrix product V is defined, as follows: = 'G': The matrix G is given: V = op(E)'*S*G*S*op(E); = 'D': The matrix D is given: V = op(E)'*S*D*D'*S*op(E); = 'F': The matrix F is given: V = F*F'; = 'H': The matrices H and K are given: V = H*K. UPLO CHARACTER*1 Specifies which triangles of the symmetric matrices R, G, if JOBG = 'G', and S, if JOBG = 'G' or JOBG = 'D', 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 matrix multiplication, as follows: = 'N': op(W) = W; = 'T': op(W) = W'; = 'C': op(W) = W'.

N (input) INTEGER The order of the matrices E, R, and S. N >= 0. M (input) INTEGER If JOBG <> 'G', the number of columns of the matrices D, F, or K'. M >= 0. If JOBG = 'G', the value of M is meaningless. E (input) DOUBLE PRECISION array, dimension (LDE,*) If JOBE = 'G' and (JOBG = 'G' or JOBG = 'D'), the leading N-by-N part of this array must contain the matrix E. If JOBE = 'I' 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 (JOBG = 'G' or JOBG = 'D'); LDE >= 1, if JOBE = 'I' or JOBG = 'F' or JOBG = 'H'. R (input) DOUBLE PRECISION array, dimension (LDR,N) 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 R(X), the residual of the algebraic Riccati equation. The other strictly triangular part is not referenced. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,N). S (input) DOUBLE PRECISION array, dimension (LDS,*) If JOBG = 'G' or JOBG = 'D', the leading N-by-N part of this array must contain the symmetric Newton step matrix S. If JOBE = 'I', the full matrix must be given. Otherwise, it is sufficient to input only the triangular part specified by UPLO, and the remaining strictly triangular part is not referenced. If JOBG = 'F', this array is not referenced. If JOBG = 'H', the leading M-by-N part of this array must contain the matrix K. LDS INTEGER The leading dimension of array S. LDS >= MAX(1,N), if JOBG = 'G' or JOBG = 'D'; LDS >= 1, if JOBG = 'F'; LDS >= MAX(1,M), if 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. The diagonal elements of this array are modified internally, but are restored on exit. If JOBG = 'D', the leading N-by-M part of this array must contain the matrix D, so that G = D*D'. If JOBG = 'F', the leading N-by-M part of this array must contain the matrix F. If JOBG = 'H', leading N-by-M part of this array must contain the matrix H. LDG INTEGER The leading dimension of array G. LDG >= MAX(1,N). ALPHA (output) DOUBLE PRECISION If INFO = 0, ALPHA contains the real number alpha which minimizes P(alpha) = norm(R(X+alpha*S), 'fro') in the interval [0,2]. If INFO = 1 or IWARN = 2, ALPHA is set equal to 1. RNORM (output) DOUBLE PRECISION On exit, if INFO >= 0, RNORM contains the Frobenius norm of the residual R(X+alpha*S).

DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if LDWORK = -1 on input, then DWORK(1) returns the optimal value of LDWORK. On exit, if LDWORK = -2 on input or INFO = -19, then DWORK(1) returns the minimal value of LDWORK. On exit, if INFO = 0, the leading N-by-N upper or lower triangular part (depending on UPLO) of this array contains the corresponding triangular part of the matrix V. LDWORK The length of the array DWORK. LDWORK >= N*N + MAX( 2*N*N, 51 ), if JOBG = 'G' and JOBE = 'G'; LDWORK >= N*N + MAX( N*N, 51 ), if JOBG = 'G' and JOBE = 'I', or JOBG = 'F', or JOBG = 'H'; LDWORK >= N*N + MAX( MAX( N*N, 51 ), MIN( 2*N*N, N*M ) ), if JOBG = 'D' and JOBE = 'G'; LDWORK >= N*N + MAX( N*N, N*M, 51 ), if JOBG = 'D' and JOBE = 'I'. For M <= N, the last two formulas simplify to LDWORK >= N*N + MAX( N*N, 51). 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.

IWARN INTEGER = 0: no warnings; = 2: no optimal line search parameter t := alpha in [0,2] was found; t = 1 was set.

INFO INTEGER = 0: successful exit; < 0: if INFO = -k, the k-th argument had an illegal value; = 1: an error occurred during the call of the SLICOT Library routine MC01XD: the eigenvalues computation for the 3-by-3 (generalized) eigenproblem failed. ALPHA and RNORM are set as specified above.

The matrix V is computed with the suitable formula, and used to set up a third order polynomial whose roots in [0,2], if any, are candidates for the solution of the minimum residual problem. The roots of the polynomial are computed by solving an equivalent 3-by-3 (generalized) eigenproblem.

[1] Benner, P. Contributions to the Numerical Solution of Algebraic Riccati Equations and Related Eigenvalue Problems. Fakultat fur Mathematik, Technische Universitat Chemnitz- Zwickau, D-09107 Chemnitz, Germany, Feb. 1997.

The calculations are backward stable. The computational effort is of the order of c*N**2*M operations. Here, M = N, if JOBG = 'G', and the coefficient c varies between 0.5 and 2.5, depending on JOBE and JOBG. (An "operation" includes a multiplication, an addition, and some address calculations.) The computed value of norm(R(X+alpha*S),'fro'), returned in RNORM, could be inaccurate if R(X) is small, since then subtraction cancellation could appear in the updating formula which is used. (This can happen, e.g., when solving a Riccati equation by Newton's method with line search, since then the sequence of R(.) tends to zero.) In such a case, it is better to recompute the residual from the data.

The routine does not ckeck if the matrix S is zero, when a quick return with ALPHA = 1 is possible. With suitable input arguments E and G, this routine may also be used for discrete-time algebraic Riccati equations. (Matrix E must be the current closed-loop matrix.)

**Program Text**

None

None

None