## SB03SX

### Estimating a forward error bound for the solution of a discrete-time Lyapunov equation

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

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.

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

```
Arguments

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.

```
Input/Output Parameters
```  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
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
On exit, the leading N-by-N part of this array contains
the symmetric absolute residual matrix R (with bounds on

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.

```
Workspace
```  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.

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

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

```
References
```   Higham, N.J.
Perturbation theory and backward error for AX-XB=C.
BIT, vol. 33, pp. 124-136, 1993.

 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.

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

```
Example

Program Text

```  None
```
Program Data
```  None
```
Program Results
```  None
```