SB04PY

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

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