SB03MX

Solving a discrete-time Lyapunov equation with matrix A quasi-triangular

Purpose

  To solve the real discrete Lyapunov matrix equation

         op(A)'*X*op(A) - X = scale*C

  where op(A) = A or A' (A**T), A is upper quasi-triangular and C is
  symmetric (C = C'). (A' denotes the transpose of the matrix A.)
  A is N-by-N, the right hand side C and the solution X are N-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 must be in Schur canonical form (as returned by LAPACK routines
  DGEES or 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 SB03MX( TRANA, N, A, LDA, C, LDC, SCALE, DWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER          TRANA
      INTEGER            INFO, LDA, LDC, N
      DOUBLE PRECISION   SCALE
C     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), 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).

  N       (input) INTEGER
          The order of the matrices A, X, and C.  N >= 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N 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,N).

  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
          On entry, the leading N-by-N part of this array must
          contain the symmetric matrix C.
          On exit, if INFO >= 0, the leading N-by-N part of this
          array contains the symmetric solution matrix X.

  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.

  DWORK   DOUBLE PRECISION array, dimension (2*N)

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if A has almost reciprocal eigenvalues; perturbed
                values were used to solve the equation (but the
                matrix A is unchanged).

  A discrete-time version of the Bartels-Stewart algorithm is used.
  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] Barraud, A.Y.                   T
      A numerical algorithm to solve A XA - X = Q.
      IEEE Trans. Auto. Contr., AC-22, pp. 883-885, 1977.

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

                            3
  The algorithm requires 0(N ) operations.

  None

Program Text

  None

  None

  None