## SB03QY

### Estimating separation between op(A) and -op(A)' and 1-norm of Theta operator for a continuous-time Lyapunov equation

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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.

Specification
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, * )

Arguments

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.

Input/Output Parameters
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.

Workspace
IWORK   INTEGER array, dimension (N*N)

DWORK   DOUBLE PRECISION array, dimension (LDWORK)

LDWORK  INTEGER
The length of the array DWORK.  LDWORK >= 2*N*N.

Error Indicator
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).

Method
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.

References
[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.

Numerical Aspects
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'.

Example

Program Text

None
Program Data
None
Program Results
None