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).
Further Comments
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