Purpose
To solve for X either the real continuous-time Sylvester equation
op(A)*X + ISGN*X*op(B) = scale*C, (1)
or the real discrete-time Sylvester equation
op(A)*X*op(B) + ISGN*X = scale*C, (2)
where op(M) = M or M**T, and ISGN = 1 or -1. A is M-by-M and
B is N-by-N; the right hand side C and the solution X are M-by-N;
and scale is an output scale factor, set less than or equal to 1
to avoid overflow in X. The solution matrix X is overwritten
onto C.
If A and/or B are not (upper) quasi-triangular, that is, block
upper triangular with 1-by-1 and 2-by-2 diagonal blocks, they are
reduced to Schur canonical form, that is, quasi-triangular with
each 2-by-2 diagonal block having its diagonal elements equal and
its off-diagonal elements of opposite sign.
Specification
SUBROUTINE SB04PD( DICO, FACTA, FACTB, TRANA, TRANB, ISGN, M, N,
$ A, LDA, U, LDU, B, LDB, V, LDV, C, LDC, SCALE,
$ DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, FACTA, FACTB, TRANA, TRANB
INTEGER INFO, ISGN, LDA, LDB, LDC, LDU, LDV, LDWORK, M,
$ N
DOUBLE PRECISION SCALE
C .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
$ DWORK( * ), U( LDU, * ), V( LDV, * )
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the equation from which X is to be determined
as follows:
= 'C': Equation (1), continuous-time case;
= 'D': Equation (2), discrete-time case.
FACTA CHARACTER*1
Specifies whether or not the real Schur factorization
of the matrix A is supplied on entry, as follows:
= 'F': On entry, A and U contain the factors from the
real Schur factorization of the matrix A;
= 'N': The Schur factorization of A will be computed
and the factors will be stored in A and U;
= 'S': The matrix A is quasi-triangular (or Schur).
FACTB CHARACTER*1
Specifies whether or not the real Schur factorization
of the matrix B is supplied on entry, as follows:
= 'F': On entry, B and V contain the factors from the
real Schur factorization of the matrix B;
= 'N': The Schur factorization of B will be computed
and the factors will be stored in B and V;
= 'S': The matrix B is quasi-triangular (or Schur).
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).
TRANB CHARACTER*1
Specifies the form of op(B) to be used, as follows:
= 'N': op(B) = B (No transpose);
= 'T': op(B) = B**T (Transpose);
= 'C': op(B) = B**T (Conjugate transpose = Transpose).
ISGN INTEGER
Specifies the sign of the equation as described before.
ISGN may only be 1 or -1.
Input/Output Parameters
M (input) INTEGER
The order of the matrix A, and the number of rows in the
matrices X and C. M >= 0.
N (input) INTEGER
The order of the matrix B, and the number of columns in
the matrices X and C. N >= 0.
A (input or input/output) DOUBLE PRECISION array,
dimension (LDA,M)
On entry, the leading M-by-M part of this array must
contain the matrix A. If FACTA = 'S', then A contains
a quasi-triangular matrix, and if FACTA = 'F', then A
is in Schur canonical form; the elements below the upper
Hessenberg part of the array A are not referenced.
On exit, if FACTA = 'N', and INFO = 0 or INFO >= M+1, the
leading M-by-M upper Hessenberg part of this array
contains the upper quasi-triangular matrix in Schur
canonical form from the Schur factorization of A. The
contents of array A is not modified if FACTA = 'F' or 'S'.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,M).
U (input or output) DOUBLE PRECISION array, dimension
(LDU,M)
If FACTA = 'F', then U is an input argument and on entry
the leading M-by-M part of this array must contain the
orthogonal matrix U of the real Schur factorization of A.
If FACTA = 'N', then U is an output argument and on exit,
if INFO = 0 or INFO >= M+1, it contains the orthogonal
M-by-M matrix from the real Schur factorization of A.
If FACTA = 'S', the array U is not referenced.
LDU INTEGER
The leading dimension of array U.
LDU >= MAX(1,M), if FACTA = 'F' or 'N';
LDU >= 1, if FACTA = 'S'.
B (input or input/output) DOUBLE PRECISION array,
dimension (LDB,N)
On entry, the leading N-by-N part of this array must
contain the matrix B. If FACTB = 'S', then B contains
a quasi-triangular matrix, and if FACTB = 'F', then B
is in Schur canonical form; the elements below the upper
Hessenberg part of the array B are not referenced.
On exit, if FACTB = 'N', and INFO = 0 or INFO = M+N+1,
the leading N-by-N upper Hessenberg part of this array
contains the upper quasi-triangular matrix in Schur
canonical form from the Schur factorization of B. The
contents of array B is not modified if FACTB = 'F' or 'S'.
LDB (input) INTEGER
The leading dimension of the array B. LDB >= max(1,N).
V (input or output) DOUBLE PRECISION array, dimension
(LDV,N)
If FACTB = 'F', then V is an input argument and on entry
the leading N-by-N part of this array must contain the
orthogonal matrix V of the real Schur factorization of B.
If FACTB = 'N', then V is an output argument and on exit,
if INFO = 0 or INFO = M+N+1, it contains the orthogonal
N-by-N matrix from the real Schur factorization of B.
If FACTB = 'S', the array V is not referenced.
LDV INTEGER
The leading dimension of array V.
LDV >= MAX(1,N), if FACTB = 'F' or 'N';
LDV >= 1, if FACTB = 'S'.
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.
On exit, if INFO = 0 or INFO = M+N+1, the leading M-by-N
part of this array contains the solution matrix X.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,M).
SCALE (output) DOUBLE PRECISION
The scale factor, scale, set less than or equal to 1 to
prevent the solution overflowing.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0 or M+N+1, then: DWORK(1) returns the
optimal value of LDWORK; if FACTA = 'N', DWORK(1+i) and
DWORK(1+M+i), i = 1,...,M, contain the real and imaginary
parts, respectively, of the eigenvalues of A; and, if
FACTB = 'N', DWORK(1+f+j) and DWORK(1+f+N+j), j = 1,...,N,
with f = 2*M if FACTA = 'N', and f = 0, otherwise, contain
the real and imaginary parts, respectively, of the
eigenvalues of B.
LDWORK INTEGER
The length of the array DWORK.
LDWORK >= MAX( 1, a+MAX( c, b+d, b+e ) ),
where a = 1+2*M, if FACTA = 'N',
a = 0, if FACTA <> 'N',
b = 2*N, if FACTB = 'N', FACTA = 'N',
b = 1+2*N, if FACTB = 'N', FACTA <> 'N',
b = 0, if FACTB <> 'N',
c = 3*M, if FACTA = 'N',
c = M, if FACTA = 'F',
c = 0, if FACTA = 'S',
d = 3*N, if FACTB = 'N',
d = N, if FACTB = 'F',
d = 0, if FACTB = 'S',
e = M, if DICO = 'C', FACTA <> 'S',
e = 0, if DICO = 'C', FACTA = 'S',
e = 2*M, if DICO = 'D'.
An upper bound is
LDWORK = 1+2*M+MAX( 3*M, 5*N, 2*N+2*M ).
For good performance, LDWORK should be larger, e.g.,
LDWORK = 1+2*M+MAX( 3*M, 5*N, 2*N+2*M, 2*N+M*N ).
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= i: if INFO = i, i = 1,...,M, the QR algorithm failed
to compute all the eigenvalues of the matrix A
(see LAPACK Library routine DGEES); the elements
2+i:1+M and 2+i+M:1+2*M of DWORK contain the real
and imaginary parts, respectively, of the
eigenvalues of A which have converged, and the
array A contains the partially converged Schur form;
= M+j: if INFO = M+j, j = 1,...,N, the QR algorithm
failed to compute all the eigenvalues of the matrix
B (see LAPACK Library routine DGEES); the elements
2+f+j:1+f+N and 2+f+j+N:1+f+2*N of DWORK contain the
real and imaginary parts, respectively, of the
eigenvalues of B which have converged, and the
array B contains the partially converged Schur form;
as defined for the parameter DWORK,
f = 2*M, if FACTA = 'N',
f = 0, if FACTA <> 'N';
= M+N+1: if DICO = 'C', and the matrices A and -ISGN*B
have common or very close eigenvalues, or
if DICO = 'D', and the matrices A and -ISGN*B have
almost reciprocal eigenvalues (that is, if lambda(i)
and mu(j) are eigenvalues of A and -ISGN*B, then
lambda(i) = 1/mu(j) for some i and j);
perturbed values were used to solve the equation
(but the matrices A and B are unchanged).
Method
An extension and refinement of the algorithms in [1,2] is used.
If the matrices A and/or B are not quasi-triangular (see PURPOSE),
they are reduced to Schur canonical form
A = U*S*U', B = V*T*V',
where U, V are orthogonal, and S, T are block upper triangular
with 1-by-1 and 2-by-2 blocks on their diagonal. The right hand
side matrix C is updated accordingly,
C = U'*C*V;
then, the solution matrix X of the "reduced" Sylvester equation
(with A and B in (1) or (2) replaced by S and T, respectively),
is computed column-wise via a back substitution scheme. A set of
equivalent linear algebraic systems of equations of order at most
four are formed and solved using Gaussian elimination with
complete pivoting. Finally, the solution X of the original
equation is obtained from the updating formula
X = U*X*V'.
If A and/or B are already quasi-triangular (or in Schur form), the
initial factorizations and the corresponding updating steps are
omitted.
References
[1] Bartels, R.H. and Stewart, G.W. T
Solution of the matrix equation A X + XB = C.
Comm. A.C.M., 15, pp. 820-826, 1972.
[2] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J.,
Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A.,
Ostrouchov, S., and Sorensen, D.
LAPACK Users' Guide: Second Edition.
SIAM, Philadelphia, 1995.
Numerical Aspects
The algorithm is stable and reliable, since orthogonal transformations and Gaussian elimination with complete pivoting are used. If INFO = M+N+1, the Sylvester equation is numerically singular.Further Comments
NoneExample
Program Text
* SB04PD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER MMAX, NMAX
PARAMETER ( MMAX = 20, NMAX = 20 )
INTEGER LDA, LDB, LDC, LDU, LDV
PARAMETER ( LDA = MMAX, LDB = NMAX, LDC = MMAX,
$ LDU = MMAX, LDV = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = 1 + 2*MMAX + MAX( 3*MMAX, 5*NMAX,
$ 2*( NMAX + MMAX ) ) )
* .. Local Scalars ..
CHARACTER DICO, FACTA, FACTB, TRANA, TRANB
INTEGER I, INFO, ISGN, J, M, N
DOUBLE PRECISION SCALE
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,MMAX), B(LDB,NMAX), C(LDC,NMAX),
$ DWORK(LDWORK), U(LDU,MMAX), V(LDV,NMAX)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL SB04PD
* .. 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, ISGN, DICO, FACTA, FACTB, TRANA, TRANB
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) M
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,M ), I = 1,M )
IF ( LSAME( FACTA, 'F' ) )
$ READ ( NIN, FMT = * ) ( ( U(I,J), J = 1,M ), I = 1,M )
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) N
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
IF ( LSAME( FACTB, 'F' ) )
$ READ ( NIN, FMT = * ) ( ( V(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,M )
* Find the solution matrix X.
CALL SB04PD( DICO, FACTA, FACTB, TRANA, TRANB, ISGN, M, N,
$ A, LDA, U, LDU, B, LDB, V, LDV, C, LDC, SCALE,
$ DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 )
$ WRITE ( NOUT, FMT = 99998 ) INFO
IF ( INFO.EQ.0 .OR. INFO.EQ.M+N+1 ) THEN
WRITE ( NOUT, FMT = 99997 )
DO 20 I = 1, M
WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99995 ) SCALE
IF ( LSAME( FACTA, 'N' ) ) THEN
WRITE ( NOUT, FMT = 99994 )
DO 40 I = 1, M
WRITE ( NOUT, FMT = 99996 ) ( U(I,J), J = 1,M )
40 CONTINUE
END IF
IF ( LSAME( FACTB, 'N' ) ) THEN
WRITE ( NOUT, FMT = 99993 )
DO 60 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( V(I,J), J = 1,N )
60 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB04PD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04PD = ',I2)
99997 FORMAT (' The solution matrix X is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' Scaling factor = ',F8.4)
99994 FORMAT (/' The orthogonal matrix U is ')
99993 FORMAT (/' The orthogonal matrix V is ')
99992 FORMAT (/' M is out of range.',/' M = ',I5)
99991 FORMAT (/' N is out of range.',/' N = ',I5)
END
Program Data
SB04PD EXAMPLE PROGRAM DATA 3 2 1 D N N N N 2.0 1.0 3.0 0.0 2.0 1.0 6.0 1.0 2.0 2.0 1.0 1.0 6.0 2.0 1.0 1.0 4.0 0.0 5.0Program Results
SB04PD EXAMPLE PROGRAM RESULTS The solution matrix X is -0.3430 0.1995 -0.1856 0.4192 0.6922 -0.2952 Scaling factor = 1.0000 The orthogonal matrix U is 0.5396 -0.7797 0.3178 0.1954 -0.2512 -0.9480 -0.8190 -0.5736 -0.0168 The orthogonal matrix V is -0.9732 -0.2298 0.2298 -0.9732