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)

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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 matrices E and B are upper triangular.

  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.

Specification
      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,*)

Arguments

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.

Input/Output Parameters
  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.

Error Indicator
  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);
          = 4:  the LAPACK routine ZSTEIN utilized to factorize M3
                (see SLICOT Library routine SG03BS) failed to
                converge. This error is unlikely to occur.

Method
  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].

References
  [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.

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

Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to Supporting Routines index