SG03BZ

Solving (for Cholesky factor) stable generalized continuous- or discrete-time complex Lyapunov equations

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

Purpose

  To compute the Cholesky factor U of the matrix X,

              H
     X = op(U)  * op(U),

  which is the solution of either the generalized c-stable
  continuous-time Lyapunov equation

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

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

  or the generalized d-stable discrete-time Lyapunov equation

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

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

  without first finding X and without the need to form the matrix
  op(B)**H * op(B).

  op(K) is either K or K**H for K = A, B, E, U. A and E are N-by-N
  matrices, op(B) is an M-by-N matrix. The resulting matrix U is an
  N-by-N upper triangular matrix with non-negative entries on its
  main diagonal. SCALE is an output scale factor set to avoid
  overflow in U.

  In the continuous-time case (1) the pencil A - lambda * E must be
  c-stable (that is, all eigenvalues must have negative real parts).
  In the discrete-time case (2) the pencil A - lambda * E must be
  d-stable (that is, the moduli of all eigenvalues must be smaller
  than one).

Specification
      SUBROUTINE SG03BZ( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q,
     $                   LDQ, Z, LDZ, B, LDB, SCALE, ALPHA, BETA, DWORK,
     $                   ZWORK, LZWORK, INFO )
C     .. Scalar Arguments ..
      DOUBLE PRECISION  SCALE
      INTEGER           INFO, LDA, LDB, LDE, LDQ, LDZ, LZWORK, M, N
      CHARACTER         DICO, FACT, TRANS
C     .. Array Arguments ..
      COMPLEX*16        A(LDA,*), ALPHA(*), B(LDB,*), BETA(*), E(LDE,*),
     $                  Q(LDQ,*), Z(LDZ,*), ZWORK(*)
      DOUBLE PRECISION  DWORK(*)

Arguments

Mode Parameters

  DICO    CHARACTER*1
          Specifies which type of the equation is considered:
          = 'C':  Continuous-time equation (1);
          = 'D':  Discrete-time equation (2).

  FACT    CHARACTER*1
          Specifies whether the generalized (complex) Schur
          factorization of the pencil A - lambda * E is supplied on
          entry or not:
          = 'N':  Factorization is not supplied;
          = 'F':  Factorization is supplied.

  TRANS   CHARACTER*1
          Specifies whether the conjugate transposed equation is to
          be solved or not:
          = 'N':  op(A) = A,    op(E) = E;
          = 'C':  op(A) = A**H, op(E) = E**H.

Input/Output Parameters
  N       (input) INTEGER
          The order of the matrix A.  N >= 0.

  M       (input) INTEGER
          The number of rows in the matrix op(B).  M >= 0.
          If M = 0, A and E are unchanged on exit, and Q, Z, ALPHA
          and BETA are not set.

  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
          On entry, if FACT = 'F', then the leading N-by-N upper
          triangular part of this array must contain the generalized
          Schur factor A_s of the matrix A (see definition (3) in
          section METHOD). A_s must be an upper triangular matrix.
          The elements below the upper triangular part of the array
          A are used as workspace.
          If FACT = 'N', then the leading N-by-N part of this array
          must contain the matrix A.
          On exit, if FACT = 'N', the leading N-by-N upper
          triangular part of this array contains the generalized
          Schur factor A_s of the matrix A. (A_s is an upper
          triangular matrix.) If FACT = 'F', the leading N-by-N
          upper triangular part of this array is unchanged.

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

  E       (input/output) COMPLEX*16 array, dimension (LDE,N)
          On entry, if FACT = 'F', then the leading N-by-N upper
          triangular part of this array must contain the generalized
          Schur factor E_s of the matrix E (see definition (4) in
          section METHOD). E_s must be an upper triangular matrix.
          The elements below the upper triangular part of the array
          E are used as workspace.
          If FACT = 'N', then the leading N-by-N part of this array
          must contain the coefficient matrix E of the equation.
          On exit, if FACT = 'N', the leading N-by-N upper
          triangular part of this array contains the generalized
          Schur factor E_s of the matrix E. (E_s is an upper
          triangular matrix.) If FACT = 'F', the leading N-by-N
          upper triangular part of this array is unchanged.

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

  Q       (input/output) COMPLEX*16 array, dimension (LDQ,N)
          On entry, if FACT = 'F', then the leading N-by-N part of
          this array must contain the unitary matrix Q from the
          generalized Schur factorization (see definitions (3) and
          (4) in section METHOD), or an identity matrix (if the
          original equation has upper triangular matrices A and E).
          If FACT = 'N', Q need not be set on entry.
          On exit, if FACT = 'N', the leading N-by-N part of this
          array contains the unitary matrix Q from the generalized
          Schur factorization. If FACT = 'F', this array is
          unchanged.

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

  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
          On entry, if FACT = 'F', then the leading N-by-N part of
          this array must contain the unitary matrix Z from the
          generalized Schur factorization (see definitions (3) and
          (4) in section METHOD), or an identity matrix (if the
          original equation has upper triangular matrices A and E).
          If FACT = 'N', Z need not be set on entry.
          On exit, if FACT = 'N', the leading N-by-N part of this
          array contains the unitary matrix Z from the generalized
          Schur factorization. If FACT = 'F', this array is
          unchanged.

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

  B       (input/output) COMPLEX*16 array, dimension (LDB,N1)
          On entry, if TRANS = 'C', the leading N-by-M part of this
          array must contain the matrix B and N1 >= MAX(M,N).
          If TRANS = 'N', the leading M-by-N part of this array
          must contain the matrix B and N1 >= N.
          On exit, if INFO = 0, the leading N-by-N part of this
          array contains the Cholesky factor U of the solution
          matrix X of the problem, X = op(U)**H * op(U).
          If M = 0 and N > 0, then U is set to zero.

  LDB     INTEGER
          The leading dimension of the array B.
          If TRANS = 'C',  LDB >= MAX(1,N).
          If TRANS = 'N',  LDB >= MAX(1,M,N).

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

  ALPHA   (output) COMPLEX*16 arrays, dimension (N)
  BETA    If INFO = 0, 5, 6, or 7, then ALPHA(j)/BETA(j),
          j = 1, ... , N, are the eigenvalues of the matrix pencil
          A - lambda * E (the diagonals of the complex Schur form).
          All BETA(j) are non-negative real numbers.
          ALPHA will be always less than and usually comparable with
          norm(A) in magnitude, and BETA always less than and
          usually comparable with norm(B).

Workspace
  DWORK   DOUBLE PRECISION array, dimension (LDWORK), where
          LDWORK = 0,           if MIN(M,N) = 0 or
                                   FACT = 'F' and N <= 1; else,
          LDWORK = N-1,         if FACT = 'F' and DICO = 'C';
          LDWORK = MAX(N-1,10), if FACT = 'F' and DICO = 'D';
          LDWORK = 8*N,         if FACT = 'N'.

  ZWORK   COMPLEX*16 array, dimension (LZWORK)
          On exit, if INFO = 0, ZWORK(1) returns the optimal value
          of LZWORK.
          On exit, if INFO = -21, ZWORK(1) returns the minimum value
          of LZWORK.

  LZWORK  INTEGER
          The dimension of the array ZWORK.
          LZWORK >= MAX(1,3*N-3,2*N).
          For good performance, LZWORK should be larger.

          If LZWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal size of the ZWORK
          array, returns this value as the first entry of the ZWORK
          array, and no error message related to LZWORK is issued by
          XERBLA.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 4:  FACT = 'N' and the pencil A - lambda * E cannot be
                reduced to generalized Schur form: LAPACK routine
                ZGGES has failed to converge;
          = 5:  DICO = 'C' and the pencil A - lambda * E is not
                c-stable;
          = 6:  DICO = 'D' and the pencil A - lambda * E is not
                d-stable;
          = 7:  the LAPACK routine ZSTEIN utilized to factorize M3
                failed to converge in the discrete-time case (see
                section METHOD for SLICOT Library routine SG03BS).
                This error is unlikely to occur.

Method
  An extension [2] of Hammarling's method [1] to generalized
  Lyapunov equations is utilized to solve (1) or (2).

  First the pencil A - lambda * E is reduced to complex generalized
  Schur form A_s - lambda * E_s by means of unitary transformations
  (QZ-algorithm):

     A_s = Q**H * A * Z   (upper triangular),                    (3)

     E_s = Q**H * E * Z   (upper triangular).                    (4)

  If the pencil A - lambda * E has already been factorized prior to
  calling the routine, however, then the factors A_s, E_s, Q and Z
  may be supplied and the initial factorization omitted.

  Depending on the parameters TRANS and M, the N-by-N upper
  triangular matrix B_s is defined as follows. In any case Q_B is
  an M-by-M unitary matrix, which need not be accumulated.

  1. If TRANS = 'N' and M < N, B_s is the upper triangular matrix
     from the QR-factorization

        ( Q_B  O )           ( B * Z )
        (        ) * B_s  =  (       ),
        (  O   I )           (   O   )

     where the O's are zero matrices of proper size and I is the
     identity matrix of order N-M.

  2. If TRANS = 'N' and M >= N, B_s is the upper triangular matrix
     from the (rectangular) QR-factorization

              ( B_s )
        Q_B * (     )  =  B * Z,
              (  O  )

     where O is the (M-N)-by-N zero matrix.

  3. If TRANS = 'C' and M < N, B_s is the upper triangular matrix
     from the RQ-factorization

                    ( Q_B  O )
        (B_s  O ) * (        )  =  ( Q**H * B   O ).
                    (  O   I )

  4. If TRANS = 'C' and M >= N, B_s is the upper triangular matrix
     from the (rectangular) RQ-factorization

        ( B_s   O ) * Q_B  =  Q**H * B,

     where O is the N-by-(M-N) zero matrix.

  Assuming SCALE = 1, the transformation of A, E and B described
  above leads to the reduced continuous-time equation

              H        H
       op(A_s)  op(U_s)  op(U_s) op(E_s)

              H        H
     + op(E_s)  op(U_s)  op(U_s) op(A_s)

                 H
     =  - op(B_s)  op(B_s)                                       (5)

  or to the reduced discrete-time equation

              H        H
       op(A_s)  op(U_s)  op(U_s) op(A_s)

              H        H
     - op(E_s)  op(U_s)  op(U_s) op(E_s)

                 H
     =  - op(B_s)  op(B_s).                                      (6)

  For brevity we restrict ourself to equation (5) and the case
  TRANS = 'N'. The other three cases can be treated in a similar
  fashion.

  We use the following partitioning for the matrices A_s, E_s, B_s,
  and U_s

              ( A11   A12 )          ( E11   E12 )
        A_s = (           ),   E_s = (           ),
              (   0   A22 )          (   0   E22 )

              ( B11   B12 )          ( U11   U12 )
        B_s = (           ),   U_s = (           ).              (7)
              (   0   B22 )          (   0   U22 )

  The size of the (1,1)-blocks is 1-by-1.

  We compute U11, U12**H, and U22 in three steps.

  Step I:

     From (5) and (7) we get the 1-by-1 equation

             H      H                   H      H
          A11  * U11  * U11 * E11  + E11  * U11  * U11 * A11

                 H
          = - B11  * B11.

     For brevity, details are omitted here. See [2]. The technique
     for computing U11 is similar to those applied to standard
     Lyapunov equations in Hammarling's algorithm ([1], section 5).

     Furthermore, the auxiliary scalars M1 and M2 defined as follows

        M1 = A11 / E11 ,

        M2 = B11 / E11 / U11 ,

     are computed in a numerically reliable way.

  Step II:

     The generalized Sylvester equation

           H      H      H      H
        A22  * U12  + E22  * U12  * M1  =

             H           H      H      H      H
        - B12  * M2 - A12  * U11  - E12  * U11  * M1

     is solved for U12**H, as a linear system of order N-1.

  Step III:

     It can be shown that

           H      H                  H      H
        A22  * U22  * U22 * E22 + E22  * U22  * U22 * A22  =

             H              H
        - B22  * B22 - y * y                                     (8)

     holds, where y is defined as

               H        H      H      H      H
        y = B12  - ( E12  * U11  + E22  * U12  ) * M2 .

     If B22_tilde is the square triangular matrix arising from the
     (rectangular) QR-factorization

                    ( B22_tilde )     ( B22  )
        Q_B_tilde * (           )  =  (      ),
                    (     O     )     ( y**H )

     where Q_B_tilde is a unitary matrix of order N, then

             H              H                H
        - B22  * B22 - y * y   =  - B22_tilde  * B22_tilde.

     Replacing the right hand side in (8) by the term
     - B22_tilde**H * B22_tilde leads to a reduced generalized
     Lyapunov equation like (5), but of dimension N-1.

  The recursive application of the steps I to III yields the
  solution U_s of the equation (5).

  It remains to compute the solution matrix U of the original
  problem (1) or (2) from the matrix U_s. To this end we transform
  the solution back (with respect to the transformation that led
  from (1) to (5) (from (2) to (6)) and apply the QR-factorization
  (RQ-factorization). The upper triangular solution matrix U is
  obtained by

     Q_U * U  =  U_s * Q**H     (if TRANS = 'N'),

  or

     U * Q_U  =  Z * U_s        (if TRANS = 'C'),

  where Q_U is an N-by-N unitary matrix. Again, the unitary matrix
  Q_U need not be accumulated.

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.

Numerical Aspects
  The number of flops required by the routine is given by the
  following table. Note that we count a single floating point
  arithmetic operation as one flop.

              |           FACT = 'F'                  FACT = 'N'
     ---------+--------------------------------------------------
      M <= N  |     (13*N**3+6*M*N**2         (211*N**3+6*M*N**2
              |   +6*M**2*N-2*M**3)/3        +6*M**2*N-2*M**3)/3
              |
       M > N  | (11*N**3+12*M*N**2)/3     (209*N**3+12*M*N**2)/3

Further Comments
  The Lyapunov equation may be very ill-conditioned. In particular,
  if DICO = 'D' and the pencil A - lambda * E has a pair of almost
  reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost
  degenerate pair of eigenvalues, then the Lyapunov equation will be
  ill-conditioned. Perturbed values were used to solve the equation.
  A condition estimate can be obtained from the routine SG03AD.

Example

Program Text

  None
Program Data
  None
Program Results
  None

Return to index