SG03BX

Solving (for Cholesky factor) generalized stable 2-by-2 continuous- or discrete-time Lyapunov equations, with pencil A - lambda E having complex conjugate eigenvalues (E upper triangular)

Purpose

  To solve for X = op(U)**T * op(U) either the generalized c-stable
  continuous-time Lyapunov equation

          T                    T
     op(A)  * X * op(E) + op(E)  * X * op(A)

              2        T
     = - SCALE  * op(B)  * op(B),                                (1)

  or the generalized d-stable discrete-time Lyapunov equation

          T                    T
     op(A)  * X * op(A) - op(E)  * X * op(E)

              2        T
     = - SCALE  * op(B)  * op(B),                                (2)

  where op(K) is either K or K**T for K = A, B, E, U. The Cholesky
  factor U of the solution is computed without first finding X.

  Furthermore, the auxiliary matrices

                                -1        -1
     M1 := op(U) * op(A) * op(E)   * op(U)

                        -1        -1
     M2 := op(B) * op(E)   * op(U)

  are computed in a numerically reliable way.

  The matrices A, B, E, M1, M2, and U are real 2-by-2 matrices. The
  pencil A - lambda * E must have a pair of complex conjugate
  eigenvalues. The eigenvalues must be in the open right half plane
  (in the continuous-time case) or inside the unit circle (in the
  discrete-time case).

  The resulting matrix U is upper triangular. The entries on its
  main diagonal are non-negative. SCALE is an output scale factor
  set to avoid overflow in U.

      SUBROUTINE SG03BX( DICO, TRANS, A, LDA, E, LDE, B, LDB, U, LDU,
     $                   SCALE, M1, LDM1, M2, LDM2, INFO )
C     .. Scalar Arguments ..
      CHARACTER         DICO, TRANS
      DOUBLE PRECISION  SCALE
      INTEGER           INFO, LDA, LDB, LDE, LDM1, LDM2, LDU
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), E(LDE,*), M1(LDM1,*),
     $                  M2(LDM2,*), U(LDU,*)

Mode Parameters

  DICO    CHARACTER*1
          Specifies whether the continuous-time or the discrete-time
          equation is to be solved:
          = 'C':  Solve continuous-time equation (1);
          = 'D':  Solve discrete-time equation (2).

  TRANS   CHARACTER*1
          Specifies whether the transposed equation is to be solved
          or not:
          = 'N':  op(K) = K,     K = A, B, E, U;
          = 'T':  op(K) = K**T,  K = A, B, E, U.

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

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= 2.

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

  LDE     INTEGER
          The leading dimension of the array E.  LDE >= 2.

  B       (input) DOUBLE PRECISION array, dimension (LDB,2)
          The leading 2-by-2 upper triangular part of this array
          must contain the matrix B.

  LDB     INTEGER
          The leading dimension of the array B.  LDB >= 2.

  U       (output) DOUBLE PRECISION array, dimension (LDU,2)
          The leading 2-by-2 part of this array contains the upper
          triangular matrix U.

  LDU     INTEGER
          The leading dimension of the array U.  LDU >= 2.

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

  M1      (output) DOUBLE PRECISION array, dimension (LDM1,2)
          The leading 2-by-2 part of this array contains the
          matrix M1.

  LDM1    INTEGER
          The leading dimension of the array M1.  LDM1 >= 2.

  M2      (output) DOUBLE PRECISION array, dimension (LDM2,2)
          The leading 2-by-2 part of this array contains the
          matrix M2.

  LDM2    INTEGER
          The leading dimension of the array M2.  LDM2 >= 2.

  INFO    INTEGER
          = 0:  successful exit;
          = 2:  the eigenvalues of the pencil A - lambda * E are not
                a pair of complex conjugate numbers;
          = 3:  the eigenvalues of the pencil A - lambda * E are
                not in the open right half plane (in the continuous-
                time case) or inside the unit circle (in the
                discrete-time case).

  The method used by the routine is based on a generalization of the
  method due to Hammarling ([1], section 6) for Lyapunov equations
  of order 2. A more detailed description is given in [2].

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

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

  If the solution matrix U is singular, the matrices M1 and M2 are
  properly set (see [1], equation (6.21)).

Program Text

  None

  None

  None