**Purpose**

To estimate a forward error bound for the solution X of a real discrete-time Lyapunov matrix equation, op(A)'*X*op(A) - X = C, where op(A) = A or A' (A**T) and C is symmetric (C = C**T). The matrix A, the right hand side C, and the solution X are N-by-N. An absolute residual matrix, which takes into account the rounding errors in forming it, is given in the array R.

SUBROUTINE SB03SX( TRANA, UPLO, LYAPUN, N, XANORM, T, LDT, U, LDU, $ R, LDR, FERR, IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER LYAPUN, TRANA, UPLO INTEGER INFO, LDR, LDT, LDU, LDWORK, N DOUBLE PRECISION FERR, XANORM C .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION DWORK( * ), R( LDR, * ), T( LDT, * ), $ U( LDU, * )

**Mode Parameters**

TRANA CHARACTER*1 Specifies the form of op(A) to be used, as follows: = 'N': op(A) = A (No transpose); = 'T': op(A) = A**T (Transpose); = 'C': op(A) = A**T (Conjugate transpose = Transpose). UPLO CHARACTER*1 Specifies which part of the symmetric matrix R is to be used, as follows: = 'U': Upper triangular part; = 'L': Lower triangular part. LYAPUN CHARACTER*1 Specifies whether or not the original Lyapunov equations should be solved, as follows: = 'O': Solve the original Lyapunov equations, updating the right-hand sides and solutions with the matrix U, e.g., X <-- U'*X*U; = 'R': Solve reduced Lyapunov equations only, without updating the right-hand sides and solutions.

N (input) INTEGER The order of the matrices A and R. N >= 0. XANORM (input) DOUBLE PRECISION The absolute (maximal) norm of the symmetric solution matrix X of the Lyapunov equation. XANORM >= 0. T (input) DOUBLE PRECISION array, dimension (LDT,N) The leading N-by-N upper Hessenberg part of this array must contain the upper quasi-triangular matrix T in Schur canonical form from a Schur factorization of A. LDT INTEGER The leading dimension of array T. LDT >= MAX(1,N). U (input) DOUBLE PRECISION array, dimension (LDU,N) The leading N-by-N part of this array must contain the orthogonal matrix U from a real Schur factorization of A. If LYAPUN = 'R', the array U is not referenced. LDU INTEGER The leading dimension of array U. LDU >= 1, if LYAPUN = 'R'; LDU >= MAX(1,N), if LYAPUN = 'O'. R (input/output) DOUBLE PRECISION array, dimension (LDR,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 absolute residual matrix R, with bounds on rounding errors added. On entry, if UPLO = 'L', the leading N-by-N lower triangular part of this array must contain the lower triangular part of the absolute residual matrix R, with bounds on rounding errors added. On exit, the leading N-by-N part of this array contains the symmetric absolute residual matrix R (with bounds on rounding errors added), fully stored. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,N). FERR (output) DOUBLE PRECISION An estimated forward error bound for the solution X. If XTRUE is the true solution, FERR bounds the magnitude of the largest entry in (X - XTRUE) divided by the magnitude of the largest entry in X. If N = 0 or XANORM = 0, FERR is set to 0, without any calculations.

IWORK INTEGER array, dimension (N*N) DWORK DOUBLE PRECISION array, dimension (LDWORK) LDWORK INTEGER The length of the array DWORK. LDWORK >= 0, if N = 0; LDWORK >= MAX(3,2*N*N), if N > 0.

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = N+1: if T has almost reciprocal eigenvalues; perturbed values were used to solve Lyapunov equations (but the matrix T is unchanged).

The forward error bound is estimated using a practical error bound similar to the one proposed in [1], based on the 1-norm estimator in [2].

[1] Higham, N.J. Perturbation theory and backward error for AX-XB=C. BIT, vol. 33, pp. 124-136, 1993. [2] Higham, N.J. FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation. ACM Trans. Math. Softw., 14, pp. 381-396, 1988.

3 The algorithm requires 0(N ) operations.

The option LYAPUN = 'R' may occasionally produce slightly worse or better estimates, and it is much faster than the option 'O'. The routine can be also used as a final step in estimating a forward error bound for the solution of a discrete-time algebraic matrix Riccati equation.

**Program Text**

None

None

None