**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.

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( * )

**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.

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.

DWORK DOUBLE PRECISION array, dimension (2*M)

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

The solution matrix X is computed column-wise via a back substitution scheme, an extension and refinement of the algorithm in [1], similar to that used in [2] 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.

[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.

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

None

**Program Text**

None

None

None