Purpose
To solve for X the discrete-time Sylvester equation
X + AXB = C,
with at least one of the matrices A or B in Schur form and the
other in Hessenberg or Schur form (both either upper or lower);
A, B, C and X are N-by-N, M-by-M, N-by-M, and N-by-M matrices,
respectively.
Specification
SUBROUTINE SB04RD( ABSCHU, ULA, ULB, N, M, A, LDA, B, LDB, C,
$ LDC, TOL, IWORK, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER ABSCHU, ULA, ULB
INTEGER INFO, LDA, LDB, LDC, LDWORK, M, N
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*)
Arguments
Mode Parameters
ABSCHU CHARACTER*1
Indicates whether A and/or B is/are in Schur or
Hessenberg form as follows:
= 'A': A is in Schur form, B is in Hessenberg form;
= 'B': B is in Schur form, A is in Hessenberg form;
= 'S': Both A and B are in Schur form.
ULA CHARACTER*1
Indicates whether A is in upper or lower Schur form or
upper or lower Hessenberg form as follows:
= 'U': A is in upper Hessenberg form if ABSCHU = 'B' and
upper Schur form otherwise;
= 'L': A is in lower Hessenberg form if ABSCHU = 'B' and
lower Schur form otherwise.
ULB CHARACTER*1
Indicates whether B is in upper or lower Schur form or
upper or lower Hessenberg form as follows:
= 'U': B is in upper Hessenberg form if ABSCHU = 'A' and
upper Schur form otherwise;
= 'L': B is in lower Hessenberg form if ABSCHU = 'A' and
lower Schur form otherwise.
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) DOUBLE PRECISION array, dimension (LDA,N)
The leading N-by-N part of this array must contain the
coefficient matrix A of the equation.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,N).
B (input) DOUBLE PRECISION array, dimension (LDB,M)
The leading M-by-M part of this array must contain the
coefficient matrix B of the equation.
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, if INFO = 0, 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).
Tolerances
TOL DOUBLE PRECISION
The tolerance to be used to test for near singularity in
the Sylvester equation. If the user sets TOL > 0, then the
given value of TOL is used as a lower bound for the
reciprocal condition number; a matrix whose estimated
condition number is less than 1/TOL is considered to be
nonsingular. If the user sets TOL <= 0, then a default
tolerance, defined by TOLDEF = EPS, is used instead, where
EPS is the machine precision (see LAPACK Library routine
DLAMCH).
This parameter is not referenced if ABSCHU = 'S',
ULA = 'U', and ULB = 'U'.
Workspace
IWORK INTEGER array, dimension (2*MAX(M,N))
This parameter is not referenced if ABSCHU = 'S',
ULA = 'U', and ULB = 'U'.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
LDWORK INTEGER
The length of the array DWORK.
LDWORK = 2*N, if ABSCHU = 'S', ULA = 'U', and ULB = 'U';
LDWORK = 2*MAX(M,N)*(4 + 2*MAX(M,N)), otherwise.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if a (numerically) singular matrix T was encountered
during the computation of the solution matrix X.
That is, the estimated reciprocal condition number
of T is less than or equal to TOL.
Method
Matrices A and B are assumed to be in (upper or lower) Hessenberg or Schur form (with at least one of them in Schur form). The solution matrix X is then computed by rows or columns via the back substitution scheme proposed by Golub, Nash and Van Loan (see [1]), which involves the solution of triangular systems of equations that are constructed recursively and which may be nearly singular if A and -B have almost reciprocal eigenvalues. If near singularity is detected, then the routine returns with the Error Indicator (INFO) set to 1.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.
[2] Sima, V.
Algorithms for Linear-quadratic Optimization.
Marcel Dekker, Inc., New York, 1996.
Numerical Aspects
2 2
The algorithm requires approximately 5M N + 0.5MN operations in
2 2
the worst case and 2.5M N + 0.5MN operations in the best case
(where M is the order of the matrix in Hessenberg form and N is
the order of the matrix in Schur form) and is mixed stable (see
[1]).
Further Comments
NoneExample
Program Text
* SB04RD 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
PARAMETER ( LDA = NMAX, LDB = MMAX, LDC = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = 2*( MAX( NMAX,MMAX ) )*
$ ( 4+2*( MAX( NMAX,MMAX ) ) ) )
INTEGER LIWORK
PARAMETER ( LIWORK = 2*MAX( NMAX,MMAX ) )
* .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER I, INFO, J, M, N
CHARACTER*1 ABSCHU, ULA, ULB
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,MMAX),
$ DWORK(LDWORK)
INTEGER IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL SB04RD
* .. 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, TOL, ULA, ULB, ABSCHU
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99995 ) 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 = 99994 ) 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 SB04RD( ABSCHU, ULA, ULB, N, M, A, LDA, B, LDB, C,
$ LDC, TOL, 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
END IF
END IF
END IF
STOP
*
99999 FORMAT (' SB04RD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04RD = ',I2)
99997 FORMAT (' The solution matrix X is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' N is out of range.',/' N = ',I5)
99994 FORMAT (/' M is out of range.',/' M = ',I5)
END
Program Data
SB04RD EXAMPLE PROGRAM DATA 5 5 0.0 U U B 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 1.0 0.0 2.0 3.0 4.0 5.0 0.0 0.0 6.0 7.0 8.0 0.0 0.0 0.0 9.0 1.0 1.0 2.0 3.0 4.0 5.0 0.0 1.0 2.0 3.0 4.0 0.0 0.0 1.0 2.0 3.0 0.0 0.0 0.0 1.0 -5.0 0.0 0.0 0.0 4.0 1.0 2.0 4.0 10.0 40.0 7.0 6.0 20.0 40.0 74.0 38.0 0.0 2.0 8.0 36.0 2.0 0.0 0.0 6.0 52.0 -9.0 0.0 0.0 0.0 13.0 -43.0Program Results
SB04RD EXAMPLE PROGRAM RESULTS The solution matrix X is 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000