## SB04OD

### Solution of generalized Sylvester equations with separation estimation

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

Purpose

```  To solve for R and L one of the generalized Sylvester equations

A * R - L * B = scale * C )
)                                 (1)
D * R - L * E = scale * F )

or

A' * R + D' * L = scale * C    )
)                            (2)
R * B' + L * E' = scale * (-F) )

where A and D are M-by-M matrices, B and E are N-by-N matrices and
C, F, R and L are M-by-N matrices.

The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an
output scaling factor chosen to avoid overflow.

The routine also optionally computes a Dif estimate, which
measures the separation of the spectrum of the matrix pair (A,D)
from the spectrum of the matrix pair (B,E), Dif[(A,D),(B,E)].

```
Specification
```      SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C,
\$                   LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, P,
\$                   LDP, Q, LDQ, U, LDU, V, LDV, IWORK, DWORK,
\$                   LDWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER         JOBD, REDUCE, TRANS
INTEGER           INFO, LDA, LDB, LDC, LDD, LDE, LDF, LDP, LDQ,
\$                  LDU, LDV, LDWORK, M, N
DOUBLE PRECISION  DIF, SCALE
C     .. Array Arguments ..
INTEGER           IWORK(*)
DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
\$                  DWORK(*), E(LDE,*), F(LDF,*), P(LDP,*),
\$                  Q(LDQ,*), U(LDU,*), V(LDV,*)

```
Arguments

Mode Parameters

```  REDUCE  CHARACTER*1
Indicates whether the matrix pairs (A,D) and/or (B,E) are
to be reduced to generalized Schur form as follows:
= 'R':  The matrix pairs (A,D) and (B,E) are to be reduced
to generalized (real) Schur canonical form;
= 'A':  The matrix pair (A,D) only is to be reduced
to generalized (real) Schur canonical form,
and the matrix pair (B,E) already is in this form;
= 'B':  The matrix pair (B,E) only is to be reduced
to generalized (real) Schur canonical form,
and the matrix pair (A,D) already is in this form;
= 'N':  The matrix pairs (A,D) and (B,E) are already in
generalized (real) Schur canonical form, as
produced by LAPACK routine DGGES.

TRANS   CHARACTER*1
Indicates which of the equations, (1) or (2), is to be
solved as follows:
= 'N':  The generalized Sylvester equation (1) is to be
solved;
= 'T':  The "transposed" generalized Sylvester equation
(2) is to be solved.

JOBD    CHARACTER*1
Indicates whether the Dif estimator is to be computed as
follows:
= '1':  Only the one-norm-based Dif estimate is computed
and stored in DIF;
= '2':  Only the Frobenius norm-based Dif estimate is
computed and stored in DIF;
= 'D':  The equation (1) is solved and the one-norm-based
Dif estimate is computed and stored in DIF;
= 'F':  The equation (1) is solved and the Frobenius norm-
based Dif estimate is computed and stored in DIF;
= 'N':  The Dif estimator is not required and hence DIF is
not referenced. (Solve either (1) or (2) only.)
JOBD is not referenced if TRANS = 'T'.

```
Input/Output Parameters
```  M       (input) INTEGER
The order of the matrices A and D and the number of rows
of the matrices C, F, R and L.  M >= 0.

N       (input) INTEGER
The order of the matrices B and E and the number of
columns of the matrices C, F, R and L.  N >= 0.
No computations are performed if N = 0 or M = 0, but SCALE
and DIF (if JOB <> 'N') are set to 1.

A       (input/output) DOUBLE PRECISION array, dimension (LDA,M)
On entry, the leading M-by-M part of this array must
contain the coefficient matrix A of the equation; A must
be in upper quasi-triangular form if REDUCE = 'B' or 'N'.
On exit, the leading M-by-M part of this array contains
the upper quasi-triangular form of A.

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

B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
On entry, the leading N-by-N part of this array must
contain the coefficient matrix B of the equation; B must
be in upper quasi-triangular form if REDUCE = 'A' or 'N'.
On exit, the leading N-by-N part of this array contains
the upper quasi-triangular form of B.

LDB     INTEGER
The leading dimension of array B.  LDB >= MAX(1,N).

C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading M-by-N part of this array must
contain the right-hand side matrix C of the first equation
in (1) or (2).
On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N
part of this array contains the solution matrix R of the
problem; if JOBD = '1' or '2' and TRANS = 'N', the leading
M-by-N part of this array contains the solution matrix R
achieved during the computation of the Dif estimate.

LDC     INTEGER
The leading dimension of array C.  LDC >= MAX(1,M).

D       (input/output) DOUBLE PRECISION array, dimension (LDD,M)
On entry, the leading M-by-M part of this array must
contain the coefficient matrix D of the equation; D must
be in upper triangular form if REDUCE = 'B' or 'N'.
On exit, the leading M-by-M part of this array contains
the upper triangular form of D.

LDD     INTEGER
The leading dimension of array D.  LDD >= MAX(1,M).

E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading N-by-N part of this array must
contain the coefficient matrix E of the equation; E must
be in upper triangular form if REDUCE = 'A' or 'N'.
On exit, the leading N-by-N part of this array contains
the upper triangular form of E.

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

F       (input/output) DOUBLE PRECISION array, dimension (LDF,N)
On entry, the leading M-by-N part of this array must
contain the right-hand side matrix F of the second
equation in (1) or (2).
On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N
part of this array contains the solution matrix L of the
problem; if JOBD = '1' or '2' and TRANS = 'N', the leading
M-by-N part of this array contains the solution matrix L
achieved during the computation of the Dif estimate.

LDF     INTEGER
The leading dimension of array F.  LDF >= MAX(1,M).

SCALE   (output) DOUBLE PRECISION
The scaling factor in (1) or (2). If 0 < SCALE < 1, C and
F hold the solutions R and L, respectively, to a slightly
perturbed system, but the computed generalized (real)
Schur canonical form matrices P'*A*Q, U'*B*V, P'*D*Q and
U'*E*V (see METHOD), or input matrices A, B, D, and E (if
already reduced to these forms), have not been changed.
If SCALE = 0, C and F hold the solutions R and L,
respectively, to the homogeneous system with C = F = 0.
Normally, SCALE = 1.

DIF     (output) DOUBLE PRECISION
If TRANS = 'N' and JOBD <> 'N', then DIF contains the
value of the Dif estimator, which is an upper bound of
-1
Dif[(A,D),(B,E)] = sigma_min(Z) = 1/||Z  ||, in either the
one-norm, or Frobenius norm, respectively (see METHOD).
Otherwise, DIF is not referenced.

P       (output) DOUBLE PRECISION array, dimension (LDP,*)
If REDUCE = 'R' or 'A', then the leading M-by-M part of
this array contains the (left) transformation matrix used
to reduce (A,D) to generalized Schur form.
Otherwise, P is not referenced and can be supplied as a
dummy array (i.e. set parameter LDP = 1 and declare this
array to be P(1,1) in the calling program).

LDP     INTEGER
The leading dimension of array P.
LDP >= MAX(1,M) if REDUCE = 'R' or 'A',
LDP >= 1        if REDUCE = 'B' or 'N'.

Q       (output) DOUBLE PRECISION array, dimension (LDQ,*)
If REDUCE = 'R' or 'A', then the leading M-by-M part of
this array contains the (right) transformation matrix used
to reduce (A,D) to generalized Schur form.
Otherwise, Q is not referenced and can be supplied as a
dummy array (i.e. set parameter LDQ = 1 and declare this
array to be Q(1,1) in the calling program).

LDQ     INTEGER
The leading dimension of array Q.
LDQ >= MAX(1,M) if REDUCE = 'R' or 'A',
LDQ >= 1        if REDUCE = 'B' or 'N'.

U       (output) DOUBLE PRECISION array, dimension (LDU,*)
If REDUCE = 'R' or 'B', then the leading N-by-N part of
this array contains the (left) transformation matrix used
to reduce (B,E) to generalized Schur form.
Otherwise, U is not referenced and can be supplied as a
dummy array (i.e. set parameter LDU = 1 and declare this
array to be U(1,1) in the calling program).

LDU     INTEGER
The leading dimension of array U.
LDU >= MAX(1,N) if REDUCE = 'R' or 'B',
LDU >= 1        if REDUCE = 'A' or 'N'.

V       (output) DOUBLE PRECISION array, dimension (LDV,*)
If REDUCE = 'R' or 'B', then the leading N-by-N part of
this array contains the (right) transformation matrix used
to reduce (B,E) to generalized Schur form.
Otherwise, V is not referenced and can be supplied as a
dummy array (i.e. set parameter LDV = 1 and declare this
array to be V(1,1) in the calling program).

LDV     INTEGER
The leading dimension of array V.
LDV >= MAX(1,N) if REDUCE = 'R' or 'B',
LDV >= 1        if REDUCE = 'A' or 'N'.

```
Workspace
```  IWORK   INTEGER array, dimension (M+N+6)

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.
If TRANS = 'N' and JOBD = 'D' or 'F', then
LDWORK = MAX(1,11*MN,10*MN+23,2*M*N) if REDUCE = 'R';
LDWORK = MAX(1,11*M, 10*M+23, 2*M*N) if REDUCE = 'A';
LDWORK = MAX(1,11*N, 10*N+23, 2*M*N) if REDUCE = 'B';
LDWORK = MAX(1,2*M*N)                if REDUCE = 'N',
where MN = max(M,N).
Otherwise, the term 2*M*N above should be omitted.
For optimum performance LDWORK should be larger.

If LDWORK = -1  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 is issued by XERBLA.

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  if REDUCE <> 'N' and either (A,D) and/or (B,E)
cannot be reduced to generalized Schur form;
= 2:  if REDUCE = 'N' and either A or B is not in
upper quasi-triangular form;
= 3:  if a singular matrix was encountered during the
computation of the solution matrices R and L, that
is (A,D) and (B,E) have common or close eigenvalues.

```
Method
```  For the case TRANS = 'N', and REDUCE = 'R' or 'N', the algorithm
used by the routine consists of four steps (see [1] and [2]) as
follows:

(a) if REDUCE = 'R', then the matrix pairs (A,D) and (B,E) are
transformed to generalized Schur form, i.e. orthogonal
matrices P, Q, U and V are computed such that P' * A * Q
and U' * B * V are in upper quasi-triangular form and
P' * D * Q and U' * E * V are in upper triangular form;
(b) if REDUCE = 'R', then the matrices C and F are transformed
to give P' * C * V and P' * F * V respectively;
(c) if REDUCE = 'R', then the transformed system

P' * A * Q * R1 - L1 * U' * B * V = scale * P' * C * V
P' * D * Q * R1 - L1 * U' * E * V = scale * P' * F * V

is solved to give R1 and L1; otherwise, equation (1) is
solved to give R and L directly. The Dif estimator
is also computed if JOBD <> 'N'.
(d) if REDUCE = 'R', then the solution is transformed back
to give R = Q * R1 * V' and L = P * L1 * U'.

By using Kronecker products, equation (1) can also be written as
the system of linear equations Z * x = scale*y (see [1]), where

| I*A    I*D  |
Z = |             |.
|-B'*I  -E'*I |

-1
If JOBD <> 'N', then a lower bound on ||Z  ||, in either the one-
norm or Frobenius norm, is computed, which in most cases is
a reliable estimate of the true value. Notice that since Z is a
matrix of order 2 * M * N, the exact value of Dif (i.e., in the
Frobenius norm case, the smallest singular value of Z) may be very
expensive to compute.

The case TRANS = 'N', and REDUCE = 'A' or 'B', is similar, but
only one of the matrix pairs should be reduced and the
calculations simplify.

For the case TRANS = 'T', and REDUCE = 'R' or 'N', the algorithm
is similar, but the steps (b), (c), and (d) are as follows:

(b) if REDUCE = 'R', then the matrices C and F are transformed
to give Q' * C * V and P' * F * U respectively;
(c) if REDUCE = 'R', then the transformed system

Q' * A' * P * R1 + Q' * D' * P * L1 =  scale * Q' * C * V
R1 * V' * B' * U + L1 * V' * E' * U = -scale * P' * F * U

is solved to give R1 and L1; otherwise, equation (2) is
solved to give R and L directly.
(d) if REDUCE = 'R', then the solution is transformed back
to give R = P * R1 * V' and L = P * L1 * V'.

```
References
```  [1] Kagstrom, B. and Westin, L.
Generalized Schur Methods with Condition Estimators for
Solving the Generalized Sylvester Equation.
IEEE Trans. Auto. Contr., 34, pp. 745-751, 1989.
[2] Kagstrom, B. and Westin, L.
GSYLV - Fortran Routines for the Generalized Schur Method with
Dif Estimators for Solving the Generalized Sylvester
Equation.
Report UMINF-132.86, Institute of Information Processing,
Univ. of Umea, Sweden, July 1987.
[3] Golub, G.H., Nash, S. and Van Loan, C.F.
A Hessenberg-Schur Method for the Problem AX + XB = C.
IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979.
[4] Kagstrom, B. and Van Dooren, P.
Additive Decomposition of a Transfer Function with respect to
a Specified Region.
In: "Signal Processing, Scattering and Operator Theory, and
Numerical Methods" (Eds. M.A. Kaashoek et al.).
Proceedings of MTNS-89, Vol. 3, pp. 469-477, Birkhauser Boston
Inc., 1990.
[5] Kagstrom, B. and Van Dooren, P.
A Generalized State-space Approach for the Additive
Decomposition of a Transfer Matrix.
Report UMINF-91.12, Institute of Information Processing, Univ.
of Umea, Sweden, April 1991.

```
Numerical Aspects
```  The algorithm is backward stable. A reliable estimate for the
condition number of Z in the Frobenius norm, is (see [1])

K(Z) = SQRT(  ||A||**2 + ||B||**2 + ||C||**2 + ||D||**2 )/DIF.

If mu is an upper bound on the relative error of the elements of
the matrices A, B, C, D, E and F, then the relative error in the
actual solution is approximately mu * K(Z).

The relative error in the computed solution (due to rounding
errors) is approximately EPS * K(Z), where EPS is the machine
precision (see LAPACK Library routine DLAMCH).

```
```  For applications of the generalized Sylvester equation in control
theory, see [4] and [5].

```
Example

Program Text

```*     SB04OD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER        ( ONE = 1.0D0 )
INTEGER          NIN, NOUT
PARAMETER        ( NIN = 5, NOUT = 6 )
INTEGER          MMAX, NMAX
PARAMETER        ( MMAX = 10, NMAX = 10 )
INTEGER          LDA, LDB, LDC, LDD, LDE, LDF, LDP, LDQ, LDU, LDV
PARAMETER        ( LDA = MMAX, LDB = NMAX, LDC = MMAX, LDD = MMAX,
\$                   LDE = NMAX, LDF = MMAX, LDP = MMAX, LDQ = MMAX,
\$                   LDU = NMAX, LDV = NMAX )
INTEGER          LDWORK, LIWORK
PARAMETER        ( LDWORK = MAX(11*MAX(MMAX,NMAX),
\$                                10*MAX(MMAX,NMAX)+23,2*MMAX*NMAX),
\$                   LIWORK = MMAX+NMAX+6 )
*     .. Local Scalars ..
DOUBLE PRECISION DIF, SCALE
INTEGER          I, INFO, J, M, N
CHARACTER*1      JOBD, REDUCE, TRANS
*     .. Local Arrays ..
DOUBLE PRECISION A(LDA,MMAX), B(LDB,NMAX), C(LDC,NMAX),
\$                 D(LDD,MMAX), DWORK(LDWORK), E(LDE,NMAX),
\$                 F(LDF,NMAX), P(LDP,MMAX), Q(LDQ,MMAX),
\$                 U(LDU,NMAX), V(LDV,NMAX)
INTEGER          IWORK(LIWORK)
*     .. External Functions ..
LOGICAL          LSAME
EXTERNAL         LSAME
*     .. External Subroutines ..
EXTERNAL         SB04OD
*     .. 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 = * ) M, N, REDUCE, TRANS, JOBD
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) M
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,M ), I = 1,M )
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99988 ) N
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,M )
READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,M )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( F(I,J), J = 1,N ), I = 1,M )
*           Find the solution matrices L and R.
CALL SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C,
\$                   LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, P,
\$                   LDP, Q, LDQ, U, LDU, V, LDV, IWORK, DWORK,
\$                   LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 20 I = 1, M
WRITE ( NOUT, FMT = 99991 ) ( F(I,J), J = 1,N )
20          CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 40 I = 1, M
WRITE ( NOUT, FMT = 99991 ) ( C(I,J), J = 1,N )
40          CONTINUE
IF ( LSAME( REDUCE, 'R' ).OR.LSAME( REDUCE, 'A' ) ) THEN
WRITE ( NOUT, FMT = 99995 )
DO 60 I = 1, M
WRITE ( NOUT, FMT = 99991 ) ( P(I,J), J = 1,M )
60             CONTINUE
WRITE ( NOUT, FMT = 99994 )
DO 80 I = 1, M
WRITE ( NOUT, FMT = 99991 ) ( Q(I,J), J = 1,M )
80             CONTINUE
END IF
IF ( LSAME( REDUCE, 'R' ).OR.LSAME( REDUCE, 'B' ) ) THEN
WRITE ( NOUT, FMT = 99993 )
DO 100 I = 1, N
WRITE ( NOUT, FMT = 99991 ) ( U(I,J), J = 1,N )
100             CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 120 I = 1, N
WRITE ( NOUT, FMT = 99991 ) ( V(I,J), J = 1,N )
120             CONTINUE
END IF
IF ( SCALE.NE.ONE ) WRITE ( NOUT, FMT = 99987 ) SCALE
IF ( .NOT.LSAME( JOBD, 'N' ) )
\$            WRITE ( NOUT, FMT = 99990 ) DIF
END IF
END IF
END IF
*
STOP
*
99999 FORMAT (' SB04OD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04OD = ',I2)
99997 FORMAT (' The solution matrix L is ')
99996 FORMAT (/' The solution matrix R is ')
99995 FORMAT (/' The left transformation matrix P is ')
99994 FORMAT (/' The right transformation matrix Q is ')
99993 FORMAT (/' The left transformation matrix U is ')
99992 FORMAT (/' The right transformation matrix V is ')
99991 FORMAT (20(1X,F8.4))
99990 FORMAT (/' DIF = ',F8.4)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
99988 FORMAT (/' N is out of range.',/' N = ',I5)
99987 FORMAT (/' SCALE = ',F8.4)
END
```
Program Data
``` SB04OD EXAMPLE PROGRAM DATA
3     2     R     N     D
1.6   -3.1    1.9
-3.8    4.2    2.4
0.5    2.2   -4.5
1.1    0.1
-1.3   -3.1
-2.0   28.9
-5.7  -11.8
12.9  -31.7
2.5    0.1    1.7
-2.5    0.0    0.9
0.1    5.1   -7.3
6.0    2.4
-3.6    2.5
0.5   23.8
-11.0  -10.4
39.5  -74.8
```
Program Results
``` SB04OD EXAMPLE PROGRAM RESULTS

The solution matrix L is
-0.7538  -1.6210
2.1778   1.7005
-3.5029   2.7961

The solution matrix R is
1.3064   2.7989
0.3698  -5.3376
-0.8767   6.7500

The left transformation matrix P is
-0.3093  -0.9502   0.0383
0.9366  -0.2974   0.1851
-0.1645   0.0932   0.9820

The right transformation matrix Q is
-0.6097  -0.7920  -0.0314
0.6310  -0.5090   0.5854
0.4796  -0.3371  -0.8102

The left transformation matrix U is
-0.8121   0.5835
0.5835   0.8121

The right transformation matrix V is
-0.9861   0.1660
0.1660   0.9861

DIF =   0.1147
```