## SB04PY

### Solving discrete-time Sylvester equations with matrices in real Schur form

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

Purpose

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

op(A)*X*op(B) + ISGN*X = scale*C,

where op(A) = A or A**T, A and B are both upper quasi-triangular,
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.

A and B must be in Schur canonical form (as returned by LAPACK
Library routine DHSEQR), that is, block upper triangular with
1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block has
its diagonal elements equal and its off-diagonal elements of
opposite sign.

```
Specification
```      SUBROUTINE SB04PY( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
\$                   LDC, SCALE, DWORK, INFO )
C     .. Scalar Arguments ..
CHARACTER          TRANA, TRANB
INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
DOUBLE PRECISION   SCALE
C     .. Array Arguments ..
DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
\$                   DWORK( * )

```
Arguments

Mode Parameters

```  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) DOUBLE PRECISION array, dimension (LDA,M)
The leading M-by-M part of this array must contain the
upper quasi-triangular matrix A, in Schur canonical form.
The part of A below the first sub-diagonal is not
referenced.

LDA     INTEGER
The leading dimension of array A.  LDA >= MAX(1,M).

B       (input) DOUBLE PRECISION array, dimension (LDB,N)
The leading N-by-N part of this array must contain the
upper quasi-triangular matrix B, in Schur canonical form.
The part of B below the first sub-diagonal is not
referenced.

LDB     (input) INTEGER
The leading dimension of the 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.
On exit, if INFO >= 0, 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 (2*M)

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
= 1:  A and -ISGN*B have almost reciprocal eigenvalues;
perturbed values were used to solve the equation
(but the matrices A and B are unchanged).

```
Method
```  The solution matrix X is computed column-wise via a back
substitution scheme, an extension and refinement of the algorithm
in , similar to that used in  for continuous-time Sylvester
equations. A set of equivalent linear algebraic systems of
equations of order at most four are formed and solved using
Gaussian elimination with complete pivoting.

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

 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.

```
Numerical Aspects
```  The algorithm is stable and reliable, since Gaussian elimination
with complete pivoting is used.

```
```  None
```
Example

Program Text

```  None
```
Program Data
```  None
```
Program Results
```  None
```