SG03AX

Solving a discrete-time generalized Lyapunov equation with matrix A quasi-triangular and matrix E upper triangular

Purpose

  To solve for X either the reduced generalized discrete-time
  Lyapunov equation

      T            T
     A  * X * A - E  * X * E  =  SCALE * Y                       (1)

  or

              T            T
     A * X * A  - E * X * E   =  SCALE * Y                       (2)

  where the right hand side Y is symmetric. A, E, Y, and the
  solution X are N-by-N matrices. The pencil A - lambda * E must be
  in generalized Schur form (A upper quasitriangular, E upper
  triangular). SCALE is an output scale factor, set to avoid
  overflow in X.

      SUBROUTINE SG03AX( TRANS, N, A, LDA, E, LDE, X, LDX, SCALE, INFO )
C     .. Scalar Arguments ..
      CHARACTER         TRANS
      DOUBLE PRECISION  SCALE
      INTEGER           INFO, LDA, LDE, LDX, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), E(LDE,*), X(LDX,*)

Mode Parameters

  TRANS   CHARACTER*1
          Specifies whether the transposed equation is to be solved
          or not:
          = 'N':  Solve equation (1);
          = 'T':  Solve equation (2).

  N       (input) INTEGER
          The order of the matrix A.  N >= 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N upper Hessenberg part of this array
          must contain the quasitriangular matrix A.

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

  E       (input) DOUBLE PRECISION array, dimension (LDE,N)
          The leading N-by-N upper triangular part of this array
          must contain the matrix E.

  LDE     INTEGER
          The leading dimension of the array E.  LDE >= MAX(1,N).

  X       (input/output) DOUBLE PRECISION array, dimension (LDX,N)
          On entry, the leading N-by-N part of this array must
          contain the right hand side matrix Y of the equation. Only
          the upper triangular part of this matrix need be given.
          On exit, the leading N-by-N part of this array contains
          the solution matrix X of the equation.

  LDX     INTEGER
          The leading dimension of the array X.  LDX >= MAX(1,N).

  SCALE   (output) DOUBLE PRECISION
          The scale factor set to avoid overflow in X.
          (0 < SCALE <= 1)

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  equation is (almost) singular to working precision;
                perturbed values were used to solve the equation
                (but the matrices A and E are unchanged).

  The solution X of (1) or (2) is computed via block back
  substitution or block forward substitution, respectively. (See
  [1] and [2] for details.)

  [1] Bartels, R.H., Stewart, G.W.
      Solution of the equation A X + X B = C.
      Comm. A.C.M., 15, pp. 820-826, 1972.

  [2] Penzl, T.
      Numerical solution of generalized Lyapunov equations.
      Advances in Comp. Math., vol. 8, pp. 33-48, 1998.

  8/3 * N**3 flops are required by the routine. Note that we count a
  single floating point arithmetic operation as one flop.

  The algorithm is backward stable if the eigenvalues of the pencil
  A - lambda * E are real. Otherwise, linear systems of order at
  most 4 are involved into the computation. These systems are solved
  by Gauss elimination with complete pivoting. The loss of stability
  of the Gauss elimination with complete pivoting is rarely
  encountered in practice.

  None

Program Text

  None

  None

  None