## SB04ND

### Solution of continuous-time Sylvester equations with one matrix in real Schur form and the other matrix in Hessenberg form (Hessenberg-Schur method)

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

Purpose

```  To solve for X the continuous-time Sylvester equation

AX + XB = 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 SB04ND( 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)
This parameter is not referenced if ABSCHU = 'S',
ULA = 'U', and ULB = 'U'.

LDWORK  INTEGER
The length of the array DWORK.
LDWORK = 0, 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
), which involves the solution of triangular systems of
equations that are constructed recursively and which may be nearly
singular if A and -B have close eigenvalues. If near singularity
is detected, then the routine returns with the Error Indicator
(INFO) set to 1.

```
References
```   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
```                                         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
).

```
```  None
```
Example

Program Text

```*     SB04ND 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         SB04ND
*     .. 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 SB04ND( 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 (' SB04ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB04ND = ',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
``` SB04ND EXAMPLE PROGRAM DATA
5     3     0.0     U     U     B
17.0  24.0   1.0   8.0  15.0
23.0   5.0   7.0  14.0  16.0
0.0   6.0  13.0  20.0  22.0
0.0   0.0  19.0  21.0   3.0
0.0   0.0   0.0   2.0   9.0
8.0   1.0   6.0
0.0   5.0   7.0
0.0   9.0   2.0
62.0 -12.0  26.0
59.0 -10.0  31.0
70.0  -6.0   9.0
35.0  31.0  -7.0
36.0 -15.0   7.0
```
Program Results
``` SB04ND EXAMPLE PROGRAM RESULTS

The solution matrix X is
0.0000   0.0000   1.0000
1.0000   0.0000   0.0000
0.0000   1.0000   0.0000
1.0000   1.0000  -1.0000
2.0000  -2.0000   1.0000
```