**Purpose**

To estimate the separation between the matrices op(A) and -op(A)', sep(op(A),-op(A)') = min norm(op(A)'*X + X*op(A))/norm(X) = 1 / norm(inv(Omega)) and/or the 1-norm of Theta, where op(A) = A or A' (A**T), and Omega and Theta are linear operators associated to the real continuous-time Lyapunov matrix equation op(A)'*X + X*op(A) = C, defined by Omega(W) = op(A)'*W + W*op(A), Theta(W) = inv(Omega(op(W)'*X + X*op(W))). The 1-norm condition estimators are used.

SUBROUTINE SB03QY( JOB, TRANA, LYAPUN, N, T, LDT, U, LDU, X, LDX, $ SEP, THNORM, IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER JOB, LYAPUN, TRANA INTEGER INFO, LDT, LDU, LDWORK, LDX, N DOUBLE PRECISION SEP, THNORM C .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION DWORK( * ), T( LDT, * ), U( LDU, * ), $ X( LDX, * )

**Mode Parameters**

JOB CHARACTER*1 Specifies the computation to be performed, as follows: = 'S': Compute the separation only; = 'T': Compute the norm of Theta only; = 'B': Compute both the separation and the norm of Theta. 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). 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 X. N >= 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'. X (input) DOUBLE PRECISION array, dimension (LDX,N) The leading N-by-N part of this array must contain the solution matrix X of the Lyapunov equation (reduced Lyapunov equation if LYAPUN = 'R'). If JOB = 'S', the array X is not referenced. LDX INTEGER The leading dimension of array X. LDX >= 1, if JOB = 'S'; LDX >= MAX(1,N), if JOB = 'T' or 'B'. SEP (output) DOUBLE PRECISION If JOB = 'S' or JOB = 'B', and INFO >= 0, SEP contains the estimated separation of the matrices op(A) and -op(A)'. If JOB = 'T' or N = 0, SEP is not referenced. THNORM (output) DOUBLE PRECISION If JOB = 'T' or JOB = 'B', and INFO >= 0, THNORM contains the estimated 1-norm of operator Theta. If JOB = 'S' or N = 0, THNORM is not referenced.

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

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = N+1: if the matrices T and -T' have common or very close eigenvalues; perturbed values were used to solve Lyapunov equations (but the matrix T is unchanged).

SEP is defined as the separation of op(A) and -op(A)': sep( op(A), -op(A)' ) = sigma_min( K ) where sigma_min(K) is the smallest singular value of the N*N-by-N*N matrix K = kprod( I(N), op(A)' ) + kprod( op(A)', I(N) ). I(N) is an N-by-N identity matrix, and kprod denotes the Kronecker product. The routine estimates sigma_min(K) by the reciprocal of an estimate of the 1-norm of inverse(K), computed as suggested in [1]. This involves the solution of several continuous-time Lyapunov equations, either direct or transposed. The true reciprocal 1-norm of inverse(K) cannot differ from sigma_min(K) by more than a factor of N. The 1-norm of Theta is estimated similarly.

[1] 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.

When SEP is zero, the routine returns immediately, with THNORM (if requested) not set. In this case, the equation is singular. The option LYAPUN = 'R' may occasionally produce slightly worse or better estimates, and it is much faster than the option 'O'.

**Program Text**

None

None

None