SB03OR

Solving continuous- or discrete-time Sylvester equations, with matrix S quasi-triangular, for an n-by-m matrix X, 1 <= m <= 2

Purpose

  To compute the solution of the Sylvester equations

     op(S)'*X + X*op(A) = scale*C, if DISCR = .FALSE.  or

     op(S)'*X*op(A) - X = scale*C, if DISCR = .TRUE.

  where op(K) = K or K' (i.e., the transpose of the matrix K), S is
  an N-by-N block upper triangular matrix with one-by-one and
  two-by-two blocks on the diagonal, A is an M-by-M matrix (M = 1 or
  M = 2), X and C are each N-by-M matrices, and scale is an output
  scale factor, set less than or equal to 1 to avoid overflow in X.
  The solution X is overwritten on C.

  SB03OR  is a service routine for the Lyapunov solver  SB03OT.

      SUBROUTINE SB03OR( DISCR, LTRANS, N, M, S, LDS, A, LDA, C, LDC,
     $                   SCALE, INFO )
C     .. Scalar Arguments ..
      LOGICAL            DISCR, LTRANS
      INTEGER            INFO, LDA, LDS, LDC, M, N
      DOUBLE PRECISION   SCALE
C     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), S( LDS, * )

Mode Parameters

  DISCR   LOGICAL
          Specifies the equation to be solved:
          = .FALSE.:  op(S)'*X + X*op(A) = scale*C;
          = .TRUE. :  op(S)'*X*op(A) - X = scale*C.

  LTRANS  LOGICAL
          Specifies the form of op(K) to be used, as follows:
          = .FALSE.:  op(K) = K    (No transpose);
          = .TRUE. :  op(K) = K**T (Transpose).

  N       (input) INTEGER
          The order of the matrix  S  and also the number of rows of
          matrices  X and C.  N >= 0.

  M       (input) INTEGER
          The order of the matrix  A  and also the number of columns
          of matrices  X and C.  M = 1 or M = 2.

  S       (input) DOUBLE PRECISION array, dimension (LDS,N)
          The leading  N-by-N  upper Hessenberg part of the array  S
          must contain the block upper triangular matrix. The
          elements below the upper Hessenberg part of the array  S
          are not referenced.  The array  S  must not contain
          diagonal blocks larger than two-by-two and the two-by-two
          blocks must only correspond to complex conjugate pairs of
          eigenvalues, not to real eigenvalues.

  LDS     INTEGER
          The leading dimension of array S.  LDS >= MAX(1,N).

  A       (input) DOUBLE PRECISION array, dimension (LDS,M)
          The leading  M-by-M  part of this array must contain a
          given matrix, where M = 1 or M = 2.

  LDA     INTEGER
          The leading dimension of array A.  LDA >= M.

  C       (input/output) DOUBLE PRECISION array, dimension (LDC,M)
          On entry, C must contain an N-by-M matrix, where M = 1 or
          M = 2.
          On exit, C contains the N-by-M matrix X, the solution of
          the Sylvester equation.

  LDC     INTEGER
          The leading dimension of array C.  LDC >= MAX(1,N).

  SCALE   (output) DOUBLE PRECISION
          The scale factor, scale, set less than or equal to 1 to
          prevent the solution overflowing.

  INFO    INTEGER
          = 0:  successful exit;
          = 1:  if DISCR = .FALSE., and S and -A have common
                eigenvalues, or if DISCR = .TRUE., and S and A have
                eigenvalues whose product is equal to unity;
                a solution has been computed using slightly
                perturbed values.

  The LAPACK scheme for solving Sylvester equations is adapted.

  [1] Hammarling, S.J.
      Numerical solution of the stable, non-negative definite
      Lyapunov equation.
      IMA J. Num. Anal., 2, pp. 303-325, 1982.

                            2
  The algorithm requires 0(N M) operations and is backward stable.

  None

Program Text

  None

  None

  None