### Solution of continuous- or discrete-time generalized Lyapunov equations and separation estimation

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

Purpose

```  To solve for X either the generalized continuous-time Lyapunov
equation

T                T
op(A)  X op(E) + op(E)  X op(A) = SCALE * Y,                (1)

or the generalized discrete-time Lyapunov equation

T                T
op(A)  X op(A) - op(E)  X op(E) = SCALE * Y,                (2)

where op(M) is either M or M**T for M = A, E and the right hand
side Y is symmetric. A, E, Y, and the solution X are N-by-N
matrices. SCALE is an output scale factor, set to avoid overflow
in X.

Estimates of the separation and the relative forward error norm
are provided.

```
Specification
```      SUBROUTINE SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E,
\$                   LDE, Q, LDQ, Z, LDZ, X, LDX, SCALE, SEP, FERR,
\$                   ALPHAR, ALPHAI, BETA, IWORK, DWORK, LDWORK,
\$                   INFO )
C     .. Scalar Arguments ..
CHARACTER         DICO, FACT, JOB, TRANS, UPLO
DOUBLE PRECISION  FERR, SCALE, SEP
INTEGER           INFO, LDA, LDE, LDQ, LDWORK, LDX, LDZ, N
C     .. Array Arguments ..
DOUBLE PRECISION  A(LDA,*), ALPHAI(*), ALPHAR(*), BETA(*),
\$                  DWORK(*), E(LDE,*), Q(LDQ,*), X(LDX,*),
\$                  Z(LDZ,*)
INTEGER           IWORK(*)

```
Arguments

Mode Parameters

```  DICO    CHARACTER*1
Specifies which type of the equation is considered:
= 'C':  Continuous-time equation (1);
= 'D':  Discrete-time equation (2).

JOB     CHARACTER*1
Specifies if the solution is to be computed and if the
separation is to be estimated:
= 'X':  Compute the solution only;
= 'S':  Estimate the separation only;
= 'B':  Compute the solution and estimate the separation.

FACT    CHARACTER*1
Specifies whether the generalized real Schur
factorization of the pencil A - lambda * E is supplied
on entry or not:
= 'N':  Factorization is not supplied;
= 'F':  Factorization is supplied.

TRANS   CHARACTER*1
Specifies whether the transposed equation is to be solved
or not:
= 'N':  op(A) = A,    op(E) = E;
= 'T':  op(A) = A**T, op(E) = E**T.

UPLO    CHARACTER*1
Specifies whether the lower or the upper triangle of the
array X is needed on input:
= 'L':  Only the lower triangle is needed on input;
= 'U':  Only the upper triangle is needed on input.

```
Input/Output Parameters
```  N       (input) INTEGER
The order of the matrix A.  N >= 0.

A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, if FACT = 'F', then the leading N-by-N upper
Hessenberg part of this array must contain the
generalized Schur factor A_s of the matrix A (see
definition (3) in section METHOD). A_s must be an upper
quasitriangular matrix. The elements below the upper
Hessenberg part of the array A are not referenced.
If FACT = 'N', then the leading N-by-N part of this
array must contain the matrix A.
On exit, the leading N-by-N part of this array contains
the generalized Schur factor A_s of the matrix A. (A_s is
an upper quasitriangular matrix.)

LDA     INTEGER
The leading dimension of the array A.  LDA >= MAX(1,N).

E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, if FACT = 'F', then the leading N-by-N upper
triangular part of this array must contain the
generalized Schur factor E_s of the matrix E (see
definition (4) in section METHOD). The elements below the
upper triangular part of the array E are not referenced.
If FACT = 'N', then the leading N-by-N part of this
array must contain the coefficient matrix E of the
equation.
On exit, the leading N-by-N part of this array contains
the generalized Schur factor E_s of the matrix E. (E_s is
an upper triangular matrix.)

LDE     INTEGER
The leading dimension of the array E.  LDE >= MAX(1,N).

Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
On entry, if FACT = 'F', then the leading N-by-N part of
this array must contain the orthogonal matrix Q from
the generalized Schur factorization (see definitions (3)
and (4) in section METHOD).
If FACT = 'N', Q need not be set on entry.
On exit, the leading N-by-N part of this array contains
the orthogonal matrix Q from the generalized Schur
factorization.

LDQ     INTEGER
The leading dimension of the array Q.  LDQ >= MAX(1,N).

Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
On entry, if FACT = 'F', then the leading N-by-N part of
this array must contain the orthogonal matrix Z from
the generalized Schur factorization (see definitions (3)
and (4) in section METHOD).
If FACT = 'N', Z need not be set on entry.
On exit, the leading N-by-N part of this array contains
the orthogonal matrix Z from the generalized Schur
factorization.

LDZ     INTEGER
The leading dimension of the array Z.  LDZ >= MAX(1,N).

X       (input/output) DOUBLE PRECISION array, dimension (LDX,N)
On entry, if JOB = 'B' or 'X', then the leading N-by-N
part of this array must contain the right hand side matrix
Y of the equation. Either the lower or the upper
triangular part of this array is needed (see mode
parameter UPLO).
If JOB = 'S', X is not referenced.
On exit, if JOB = 'B' or 'X', and INFO = 0, 3, or 4, then
the leading N-by-N part of this array contains the
solution matrix X of the equation.
If JOB = 'S', X is not referenced.

LDX     INTEGER
The leading dimension of the array X.  LDX >= MAX(1,N).

SCALE   (output) DOUBLE PRECISION
The scale factor set to avoid overflow in X.
(0 < SCALE <= 1)

SEP     (output) DOUBLE PRECISION
If JOB = 'S' or JOB = 'B', and INFO = 0, 3, or 4, then
SEP contains an estimate of the separation of the
Lyapunov operator.

FERR    (output) DOUBLE PRECISION
If JOB = 'B', and INFO = 0, 3, or 4, then FERR contains an
estimated forward error bound for the solution X. If XTRUE
is the true solution, FERR estimates the relative error
in the computed solution, measured in the Frobenius norm:
norm(X - XTRUE) / norm(XTRUE)

ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
BETA    (output) DOUBLE PRECISION array, dimension (N)
If FACT = 'N' and INFO = 0, 3, or 4, then
(ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, are the
eigenvalues of the matrix pencil A - lambda * E.
If FACT = 'F', ALPHAR, ALPHAI, and BETA are not
referenced.

```
Workspace
```  IWORK   INTEGER array, dimension (N**2)
IWORK is not referenced if JOB = 'X'.

DWORK   DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.

LDWORK  INTEGER
The length of the array DWORK. The following table
contains the minimal work space requirements depending
on the choice of JOB and FACT.

JOB        FACT    |  LDWORK
-------------------+-------------------
'X'        'F'     |  MAX(1,N)
'X'        'N'     |  MAX(1,4*N)
'B', 'S'   'F'     |  MAX(1,2*N**2)
'B', 'S'   'N'     |  MAX(1,2*N**2,4*N)

For optimum performance, LDWORK should be larger.

If LDWORK = -1, then a 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 related to LDWORK is issued by
XERBLA.

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  FACT = 'F' and the matrix contained in the upper
Hessenberg part of the array A is not in upper
quasitriangular form;
= 2:  FACT = 'N' and the pencil A - lambda * E cannot be
reduced to generalized Schur form: LAPACK routine
DGEGS (or DGGES) has failed to converge;
= 3:  DICO = 'D' and the pencil A - lambda * E has a
pair of reciprocal eigenvalues. That is, lambda_i =
1/lambda_j for some i and j, where lambda_i and
lambda_j are eigenvalues of A - lambda * E. Hence,
equation (2) is singular;  perturbed values were
used to solve the equation (but the matrices A and
E are unchanged);
= 4:  DICO = 'C' and the pencil A - lambda * E has a
degenerate pair of eigenvalues. That is, lambda_i =
-lambda_j for some i and j, where lambda_i and
lambda_j are eigenvalues of A - lambda * E. Hence,
equation (1) is singular;  perturbed values were
used to solve the equation (but the matrices A and
E are unchanged).

```
Method
```  A straightforward generalization  of the method proposed by
Bartels and Stewart  is utilized to solve (1) or (2).

First the pencil A - lambda * E is reduced to real generalized
Schur form A_s - lambda * E_s by means of orthogonal
transformations (QZ-algorithm):

A_s = Q**T * A * Z   (upper quasitriangular)                (3)

E_s = Q**T * E * Z   (upper triangular).                    (4)

If FACT = 'F', this step is omitted. Assuming SCALE = 1 and
defining

( Z**T * Y * Z   :   TRANS = 'N'
Y_s = <
( Q**T * Y * Q   :   TRANS = 'T'

( Q**T * X * Q    if TRANS = 'N'
X_s = <                                                     (5)
( Z**T * X * Z    if TRANS = 'T'

leads to the reduced Lyapunov equation

T                      T
op(A_s)  X_s op(E_s) + op(E_s)  X_s op(A_s) = Y_s,          (6)

or
T                      T
op(A_s)  X_s op(A_s) - op(E_s)  X_s op(E_s) = Y_s,          (7)

which are equivalent to (1) or (2), respectively. The solution X_s
of (6) or (7) is computed via block back substitution (if TRANS =
'N') or block forward substitution (if TRANS = 'T'), where the
block order is at most 2. (See  and  for details.)
Equation (5) yields the solution matrix X.

For fast computation the estimates of the separation and the
forward error are gained from (6) or (7) rather than (1) or
(2), respectively. We consider (6) and (7) as special cases of the
generalized Sylvester equation

R * X * S + U * X * V = Y,                                  (8)

whose separation is defined as follows

sep = sep(R,S,U,V) =   min   || R * X * S + U * X * V || .
||X|| = 1                           F
F

Equation (8) is equivalent to the system of linear equations

K * vec(X) = (kron(S**T,R) + kron(V**T,U)) * vec(X) = vec(Y),

where kron is the Kronecker product of two matrices and vec
is the mapping that stacks the columns of a matrix. If K is
nonsingular then

sep = 1 / ||K**(-1)|| .
2

We estimate ||K**(-1)|| by a method devised by Higham . Note
that this method yields an estimation for the 1-norm but we use it
as an approximation for the 2-norm. Estimates for the forward
error norm are provided by

FERR = 2 * EPS * ||A_s||  * ||E_s||  / sep
F          F

in the continuous-time case (1) and

FERR = EPS * ( ||A_s|| **2 + ||E_s|| **2 ) / sep
F             F

in the discrete-time case (2).
The reciprocal condition number, RCOND, of the Lyapunov equation
can be estimated by FERR/EPS.

```
References
```   Bartels, R.H., Stewart, G.W.
Solution of the equation A X + X B = C.
Comm. A.C.M., 15, pp. 820-826, 1972.

 Higham, N.J.
FORTRAN codes for estimating the one-norm of a real or complex
matrix, with applications to condition estimation.
A.C.M. Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, 1988.

 Penzl, T.
Numerical solution of generalized Lyapunov equations.
Advances in Comp. Math., vol. 8, pp. 33-48, 1998.

```
Numerical Aspects
```  The number of flops required by the routine is given by the
following table. Note that we count a single floating point
arithmetic operation as one flop. c is an integer number of modest
size (say 4 or 5).

|  FACT = 'F'            FACT = 'N'
-----------+------------------------------------------
JOB = 'B'  |  (26+8*c)/3 * N**3     (224+8*c)/3 * N**3
JOB = 'S'  |  8*c/3 * N**3          (198+8*c)/3 * N**3
JOB = 'X'  |  26/3 * N**3           224/3 * N**3

The algorithm is backward stable if the eigenvalues of the pencil
A - lambda * E are real. Otherwise, linear systems of order at
most 4 are involved into the computation. These systems are solved
by Gauss elimination with complete pivoting. The loss of stability
of the Gauss elimination with complete pivoting is rarely
encountered in practice.

The Lyapunov equation may be very ill-conditioned. In particular,
if DICO = 'D' and the pencil A - lambda * E has a pair of almost
reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost
degenerate pair of eigenvalues, then the Lyapunov equation will be
ill-conditioned. Perturbed values were used to solve the equation.
Ill-conditioning can be detected by a very small value of the
reciprocal condition number RCOND.

```
```  None
```
Example

Program Text

```*     SG03AD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
INTEGER           NIN, NOUT
PARAMETER         ( NIN = 5, NOUT = 6 )
INTEGER           NMAX
PARAMETER         ( NMAX = 20 )
INTEGER           LDA, LDE, LDQ, LDX, LDZ
PARAMETER         ( LDA = NMAX, LDE = NMAX, LDQ = NMAX,
\$                    LDX = NMAX, LDZ = NMAX )
INTEGER           LIWORK, LDWORK
PARAMETER         ( LIWORK = NMAX**2,
\$                    LDWORK = MAX( 2*NMAX**2, 4*NMAX ) )
*     .. Local Scalars ..
CHARACTER*1       DICO, FACT, JOB, TRANS, UPLO
DOUBLE PRECISION  FERR, SCALE, SEP
INTEGER           I, INFO, J, N
*     .. Local Arrays ..
INTEGER           IWORK(LIWORK)
DOUBLE PRECISION  A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX),
\$                  BETA(NMAX), DWORK(LDWORK), E(LDE,NMAX),
\$                  Q(LDQ,NMAX), X(LDX,NMAX), Z(LDZ,NMAX)
*     .. External Functions ..
LOGICAL           LSAME
EXTERNAL          LSAME
*     .. External Subroutines ..
*     .. Intrinsic Functions ..
INTRINSIC         MAX
*     .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, JOB, DICO, FACT, TRANS, UPLO
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
IF ( LSAME ( FACT, 'F' ) ) THEN
READ ( NIN, FMT = * ) ( ( Q(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( Z(I,J), J = 1,N ), I = 1,N )
END IF
IF ( .NOT.LSAME ( JOB, 'S' ) )
\$      READ ( NIN, FMT = * ) ( ( X(I,J), J = 1,N ), I = 1,N )
*        Find the solution matrix X and the scalar SEP.
CALL SG03AD( DICO, JOB, FACT, TRANS, UPLO, N, A, LDA, E, LDE,
\$                Q, LDQ, Z, LDZ, X, LDX, SCALE, SEP, FERR, ALPHAR,
\$                ALPHAI, BETA, IWORK, DWORK, LDWORK, INFO )
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF ( LSAME ( JOB, 'B' ) .OR. LSAME ( JOB, 'S' ) ) THEN
WRITE ( NOUT, FMT = 99997 ) SEP
WRITE ( NOUT, FMT = 99996 ) FERR
END IF
IF ( LSAME ( JOB, 'B' ) .OR. LSAME ( JOB, 'X' ) ) THEN
WRITE ( NOUT, FMT = 99995 ) SCALE
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99994 ) ( X(I,J), J = 1,N )
20          CONTINUE
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SG03AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SG03AD = ',I2)
99997 FORMAT (' SEP =   ',D8.2)
99996 FORMAT (' FERR =  ',D8.2)
99995 FORMAT (' SCALE = ',D8.2,//' The solution matrix X is ')
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' N is out of range.',/' N = ',I5)
END
```
Program Data
``` SG03AD EXAMPLE PROGRAM DATA
3       B       C       N       N       U
3.0     1.0     1.0
1.0     3.0     0.0
1.0     0.0     2.0
1.0     3.0     0.0
3.0     2.0     1.0
1.0     0.0     1.0
-64.0   -73.0   -28.0
0.0   -70.0   -25.0
0.0     0.0   -18.0
```
Program Results
``` SG03AD EXAMPLE PROGRAM RESULTS

SEP =   0.29D+00
FERR =  0.40D-13
SCALE = 0.10D+01

The solution matrix X is
-2.0000  -1.0000   0.0000
-1.0000  -3.0000  -1.0000
0.0000  -1.0000  -3.0000
```