Purpose
To solve for X the continuous-time Sylvester equation
AX + XB = C
where A, B, C and X are general N-by-N, M-by-M, N-by-M and
N-by-M matrices respectively.
Specification
SUBROUTINE SB04MD( N, M, A, LDA, B, LDB, C, LDC, Z, LDZ, IWORK,
$ DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
INTEGER INFO, LDA, LDB, LDC, LDWORK, LDZ, M, N
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), Z(LDZ,*)
Arguments
Input/Output Parameters
N (input) INTEGER
The order of the matrix A. N >= 0.
M (input) INTEGER
The order of the matrix B. M >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the coefficient matrix A of the equation.
On exit, the leading N-by-N upper Hessenberg part of this
array contains the matrix H, and the remainder of the
leading N-by-N part, together with the elements 2,3,...,N
of array DWORK, contain the orthogonal transformation
matrix U (stored in factored form).
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
On entry, the leading M-by-M part of this array must
contain the coefficient matrix B of the equation.
On exit, the leading M-by-M part of this array contains
the quasi-triangular Schur factor S of the matrix B'.
LDB INTEGER
The leading dimension of array B. LDB >= MAX(1,M).
C (input/output) DOUBLE PRECISION array, dimension (LDC,M)
On entry, the leading N-by-M part of this array must
contain the coefficient matrix C of the equation.
On exit, the leading N-by-M part of this array contains
the solution matrix X of the problem.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,N).
Z (output) DOUBLE PRECISION array, dimension (LDZ,M)
The leading M-by-M part of this array contains the
orthogonal matrix Z used to transform B' to real upper
Schur form.
LDZ INTEGER
The leading dimension of array Z. LDZ >= MAX(1,M).
Workspace
IWORK INTEGER array, dimension (4*N)
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK, and DWORK(2), DWORK(3),..., DWORK(N) contain
the scalar factors of the elementary reflectors used to
reduce A to upper Hessenberg form, as returned by LAPACK
Library routine DGEHRD.
LDWORK INTEGER
The length of the array DWORK.
LDWORK = MAX(1, 2*N*N + 8*N, 5*M, N + M).
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;
> 0: if INFO = i, 1 <= i <= M, the QR algorithm failed to
compute all the eigenvalues (see LAPACK Library
routine DGEES);
> M: if a singular matrix was encountered whilst solving
for the (INFO-M)-th column of matrix X.
Method
The matrix A is transformed to upper Hessenberg form H = U'AU by
the orthogonal transformation matrix U; matrix B' is transformed
to real upper Schur form S = Z'B'Z using the orthogonal
transformation matrix Z. The matrix C is also multiplied by the
transformations, F = U'CZ, and the solution matrix Y of the
transformed system
HY + YS' = F
is computed by back substitution. Finally, the matrix Y is then
multiplied by the orthogonal transformation matrices, X = UYZ', in
order to obtain the solution matrix X to the original problem.
References
[1] 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.
Numerical Aspects
3 3 2 2 The algorithm requires about (5/3) N + 10 M + 5 N M + 2.5 M N operations and is backward stable.Further Comments
NoneExample
Program Text
* SB04MD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX
PARAMETER ( NMAX = 20, MMAX = 20 )
INTEGER LDA, LDB, LDC, LDZ
PARAMETER ( LDA = NMAX, LDB = MMAX, LDC = NMAX,
$ LDZ = MMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = 4*NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( 1, 2*NMAX*NMAX+8*NMAX, 5*MMAX,
$ NMAX+MMAX ) )
* .. Local Scalars ..
INTEGER I, INFO, J, M, N
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,MMAX),
$ DWORK(LDWORK), Z(LDZ,MMAX)
INTEGER IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL SB04MD
* .. 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, M
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,M )
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,M ), I = 1,N )
* Find the solution matrix X.
CALL SB04MD( N, M, A, LDA, B, LDB, C, LDC, Z, LDZ, IWORK,
$ DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,M )
20 CONTINUE
WRITE ( NOUT, FMT = 99995 )
DO 40 I = 1, M
WRITE ( NOUT, FMT = 99996 ) ( Z(I,J), J = 1,M )
40 CONTINUE
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB04MD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04MD = ',I2)
99997 FORMAT (' The solution matrix X is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The orthogonal matrix Z is ')
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' M is out of range.',/' M = ',I5)
END
Program Data
SB04MD EXAMPLE PROGRAM DATA 3 2 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
SB04MD EXAMPLE PROGRAM RESULTS The solution matrix X is -2.7685 0.5498 -1.0531 0.6865 4.5257 -0.4389 The orthogonal matrix Z is -0.9732 -0.2298 0.2298 -0.9732