SB03OD

Solution of stable continuous- or discrete-time Lyapunov equations (Cholesky factor)

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

Purpose

  To solve for X = op(U)'*op(U) either the stable non-negative
  definite continuous-time Lyapunov equation
                                2
     op(A)'*X + X*op(A) = -scale *op(B)'*op(B)                   (1)

  or the convergent non-negative definite discrete-time Lyapunov
  equation
                                2
     op(A)'*X*op(A) - X = -scale *op(B)'*op(B)                   (2)

  where op(K) = K or K' (i.e., the transpose of the matrix K), A is
  an N-by-N matrix, op(B) is an M-by-N matrix, U is an upper
  triangular matrix containing the Cholesky factor of the solution
  matrix X, X = op(U)'*op(U), and scale is an output scale factor,
  set less than or equal to 1 to avoid overflow in X. If matrix B
  has full rank then the solution matrix X will be positive-definite
  and hence the Cholesky factor U will be nonsingular, but if B is
  rank deficient then X may be only positive semi-definite and U
  will be singular.

  In the case of equation (1) the matrix A must be stable (that
  is, all the eigenvalues of A must have negative real parts),
  and for equation (2) the matrix A must be convergent (that is,
  all the eigenvalues of A must lie inside the unit circle).

Specification
      SUBROUTINE SB03OD( DICO, FACT, TRANS, N, M, A, LDA, Q, LDQ, B,
     $                   LDB, SCALE, WR, WI, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER         DICO, FACT, TRANS
      INTEGER           INFO, LDA, LDB, LDQ, LDWORK, M, N
      DOUBLE PRECISION  SCALE
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), Q(LDQ,*), WI(*),
     $                  WR(*)

Arguments

Mode Parameters

  DICO    CHARACTER*1
          Specifies the type of Lyapunov equation to be solved as
          follows:
          = 'C':  Equation (1), continuous-time case;
          = 'D':  Equation (2), discrete-time case.

  FACT    CHARACTER*1
          Specifies whether or not the real Schur factorization
          of the matrix A is supplied on entry, as follows:
          = 'F':  On entry, A and Q contain the factors from the
                  real Schur factorization of the matrix A;
          = 'N':  The Schur factorization of A will be computed
                  and the factors will be stored in A and Q.

  TRANS   CHARACTER*1
          Specifies the form of op(K) to be used, as follows:
          = 'N':  op(K) = K    (No transpose);
          = 'T':  op(K) = K**T (Transpose).

Input/Output Parameters
  N       (input) INTEGER
          The order of the matrix A and the number of columns in
          matrix op(B).  N >= 0.

  M       (input) INTEGER
          The number of rows in matrix op(B).  M >= 0.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the leading N-by-N part of this array must
          contain the matrix A. If FACT = 'F', then A contains
          an upper quasi-triangular matrix S in Schur canonical
          form; the elements below the upper Hessenberg part of the
          array A are not referenced.
          On exit, the leading N-by-N upper Hessenberg part of this
          array contains the upper quasi-triangular matrix S in
          Schur canonical form from the Shur factorization of A.
          The contents of array A is not modified if FACT = 'F'.

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

  Q       (input or output) DOUBLE PRECISION array, dimension
          (LDQ,N)
          On entry, if FACT = 'F', then the leading N-by-N part of
          this array must contain the orthogonal matrix Q of the
          Schur factorization of A.
          Otherwise, Q need not be set on entry.
          On exit, the leading N-by-N part of this array contains
          the orthogonal matrix Q of the Schur factorization of A.
          The contents of array Q is not modified if FACT = 'F'.

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

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
          if TRANS = 'N', and dimension (LDB,max(M,N)), if
          TRANS = 'T'.
          On entry, if TRANS = 'N', the leading M-by-N part of this
          array must contain the coefficient matrix B of the
          equation.
          On entry, if TRANS = 'T', the leading N-by-M part of this
          array must contain the coefficient matrix B of the
          equation.
          On exit, the leading N-by-N part of this array contains
          the upper triangular Cholesky factor U of the solution
          matrix X of the problem, X = op(U)'*op(U).
          If M = 0 and N > 0, then U is set to zero.

  LDB     INTEGER
          The leading dimension of array B.
          LDB >= MAX(1,N,M), if TRANS = 'N';
          LDB >= MAX(1,N),   if TRANS = 'T'.

  SCALE   (output) DOUBLE PRECISION
          The scale factor, scale, set less than or equal to 1 to
          prevent the solution overflowing.

  WR      (output) DOUBLE PRECISION array, dimension (N)
  WI      (output) DOUBLE PRECISION array, dimension (N)
          If FACT = 'N', and INFO >= 0 and INFO <= 2, WR and WI
          contain the real and imaginary parts, respectively, of
          the eigenvalues of A.
          If FACT = 'F', WR and WI are not referenced.

Workspace
  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, or INFO = 1, DWORK(1) returns the
          optimal value of LDWORK.

  LDWORK  INTEGER
          The length of the array DWORK.
          If M > 0, LDWORK >= MAX(1,4*N + MIN(M,N));
          If M = 0, LDWORK >= 1.
          For optimum performance LDWORK should sometimes be larger.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if the Lyapunov equation is (nearly) singular
                (warning indicator);
                if DICO = 'C' this means that while the matrix A
                (or the factor S) has computed eigenvalues with
                negative real parts, it is only just stable in the
                sense that small perturbations in A can make one or
                more of the eigenvalues have a non-negative real
                part;
                if DICO = 'D' this means that while the matrix A
                (or the factor S) has computed eigenvalues inside
                the unit circle, it is nevertheless only just
                convergent, in the sense that small perturbations
                in A can make one or more of the eigenvalues lie
                outside the unit circle;
                perturbed values were used to solve the equation;
          = 2:  if FACT = 'N' and DICO = 'C', but the matrix A is
                not stable (that is, one or more of the eigenvalues
                of A has a non-negative real part), or DICO = 'D',
                but the matrix A is not convergent (that is, one or
                more of the eigenvalues of A lies outside the unit
                circle); however, A will still have been factored
                and the eigenvalues of A returned in WR and WI.
          = 3:  if FACT = 'F' and DICO = 'C', but the Schur factor S
                supplied in the array A is not stable (that is, one
                or more of the eigenvalues of S has a non-negative
                real part), or DICO = 'D', but the Schur factor S
                supplied in the array A is not convergent (that is,
                one or more of the eigenvalues of S lies outside the
                unit circle);
          = 4:  if FACT = 'F' and the Schur factor S supplied in
                the array A has two or more consecutive non-zero
                elements on the first sub-diagonal, so that there is
                a block larger than 2-by-2 on the diagonal;
          = 5:  if FACT = 'F' and the Schur factor S supplied in
                the array A has a 2-by-2 diagonal block with real
                eigenvalues instead of a complex conjugate pair;
          = 6:  if FACT = 'N' and the LAPACK Library routine DGEES
                has failed to converge. This failure is not likely
                to occur. The matrix B will be unaltered but A will
                be destroyed.

Method
  The method used by the routine is based on the Bartels and Stewart
  method [1], except that it finds the upper triangular matrix U
  directly without first finding X and without the need to form the
  normal matrix op(B)'*op(B).

  The Schur factorization of a square matrix A is given by

     A = QSQ',

  where Q is orthogonal and S is an N-by-N block upper triangular
  matrix with 1-by-1 and 2-by-2 blocks on its diagonal (which
  correspond to the eigenvalues of A). If A has already been
  factored prior to calling the routine however, then the factors
  Q and S may be supplied and the initial factorization omitted.

  If TRANS = 'N', the matrix B is factored as (QR factorization)
         _   _                   _   _  _
     B = P ( R ),  M >= N,   B = P ( R  Z ),  M < N,
           ( 0 )
        _                                    _
  where P is an M-by-M orthogonal matrix and R is a square upper
                                      _   _      _     _  _
  triangular matrix. Then, the matrix B = RQ, or B = ( R  Z )Q (if
  M < N) is factored as
     _                       _
     B = P ( R ),  M >= N,   B = P ( R  Z ),  M < N.

  If TRANS = 'T', the matrix B is factored as (RQ factorization)
                                      _
              _   _                 ( Z ) _
     B = ( 0  R ) P,  M >= N,   B = ( _ ) P,  M < N,
                                    ( R )
        _                                    _
  where P is an M-by-M orthogonal matrix and R is a square upper
                                      _     _     _       _   _
  triangular matrix. Then, the matrix B = Q'R, or B = Q'( Z'  R' )'
  (if M < N) is factored as
     _                       _
     B = ( R ) P,  M >= N,   B = ( Z ) P,  M < N.
                                 ( R )

  These factorizations are utilised to either transform the
  continuous-time Lyapunov equation to the canonical form
                                                     2
    op(S)'*op(V)'*op(V) + op(V)'*op(V)*op(S) = -scale *op(F)'*op(F),

  or the discrete-time Lyapunov equation to the canonical form
                                                     2
    op(S)'*op(V)'*op(V)*op(S) - op(V)'*op(V) = -scale *op(F)'*op(F),

  where V and F are upper triangular, and

     F = R,  M >= N,   F = ( R  Z ),  M < N, if TRANS = 'N';
                           ( 0  0 )

     F = R,  M >= N,   F = ( 0  Z ),  M < N, if TRANS = 'T'.
                           ( 0  R )

  The transformed equation is then solved for V, from which U is
  obtained via the QR factorization of V*Q', if TRANS = 'N', or
  via the RQ factorization of Q*V, if TRANS = 'T'.

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

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

Numerical Aspects
                            3
  The algorithm requires 0(N ) operations and is backward stable.

Further Comments
  The Lyapunov equation may be very ill-conditioned. In particular,
  if A is only just stable (or convergent) then the Lyapunov
  equation will be ill-conditioned.  A symptom of ill-conditioning
  is "large" elements in U relative to those of A and B, or a
  "small" value for scale. A condition estimate can be computed
  using SLICOT Library routine SB03MD.

  SB03OD routine can be also used for solving "unstable" Lyapunov
  equations, i.e., when matrix A has all eigenvalues with positive
  real parts, if DICO = 'C', or with moduli greater than one,
  if DICO = 'D'. Specifically, one may solve for X = op(U)'*op(U)
  either the continuous-time Lyapunov equation
                               2
     op(A)'*X + X*op(A) = scale *op(B)'*op(B),                   (3)

  or the discrete-time Lyapunov equation
                               2
     op(A)'*X*op(A) - X = scale *op(B)'*op(B),                   (4)

  provided, for equation (3), the given matrix A is replaced by -A,
  or, for equation (4), the given matrices A and B are replaced by
  inv(A) and B*inv(A), if TRANS = 'N' (or inv(A)*B, if TRANS = 'T'),
  respectively. Although the inversion generally can rise numerical
  problems, in case of equation (4) it is expected that the matrix A
  is enough well-conditioned, having only eigenvalues with moduli
  greater than 1. However, if A is ill-conditioned, it could be
  preferable to use the more general SLICOT Lyapunov solver SB03MD.

Example

Program Text

*     SB03OD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2010 NICONET e.V.
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER        ( ZERO = 0.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX
      PARAMETER        ( NMAX = 20, MMAX = 20 )
      INTEGER          LDA, LDB, LDQ, LDX, LDWORK
      PARAMETER        ( LDA = NMAX, LDB = MAX( MMAX,NMAX ),
     $                   LDQ = NMAX, LDX = NMAX )
      PARAMETER        ( LDWORK = 4*NMAX+MIN(MMAX,NMAX) )
*     .. Local Scalars ..
      DOUBLE PRECISION SCALE, TEMP
      INTEGER          I, INFO, J, K, M, N
      CHARACTER*1      DICO, FACT, TRANS
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,LDB), DWORK(LDWORK),
     $                 Q(LDQ,NMAX), WR(NMAX), WI(NMAX), X(LDX,NMAX)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         SB03OD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, M, DICO, FACT, TRANS
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99994 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         IF ( LSAME( FACT, 'F' ) ) READ ( NIN, FMT = * )
     $                         ( ( Q(I,J), J = 1,N ), I = 1,N )
         IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
            WRITE ( NOUT, FMT = 99993 ) M
         ELSE
            IF ( LSAME( TRANS, 'N' ) ) THEN
               READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,M )
            ELSE
               READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
            END IF
*           Find the Cholesky factor U.
            CALL SB03OD( DICO, FACT, TRANS, N, M, A, LDA, Q, LDQ, B,
     $                   LDB, SCALE, WR, WI, DWORK, LDWORK, INFO )
*
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
            ELSE
               WRITE ( NOUT, FMT = 99997 )
               DO 20 J = 1, N
                  WRITE ( NOUT, FMT = 99996 ) ( B(I,J), I = 1,J )
   20          CONTINUE
*              Form the solution matrix X = op(U)'*op(U).
               IF ( LSAME( TRANS, 'N' ) ) THEN
                  DO 80 I = 1, N
                     DO 60 J = I, N
                        TEMP = ZERO
                        DO 40 K = 1, I
                           TEMP = TEMP + B(K,I)*B(K,J)
   40                   CONTINUE
                        X(I,J) = TEMP
                        X(J,I) = TEMP
   60                CONTINUE
   80             CONTINUE
               ELSE
                  DO 140 I = 1, N
                     DO 120 J = I, N
                        TEMP = ZERO
                        DO 100 K = J, N
                           TEMP = TEMP + B(I,K)*B(J,K)
  100                   CONTINUE
                        X(I,J) = TEMP
                        X(J,I) = TEMP
  120                CONTINUE
  140             CONTINUE
               END IF
               WRITE ( NOUT, FMT = 99995 )
               DO 160 J = 1, N
                  WRITE ( NOUT, FMT = 99996 ) ( X(I,J), I = 1,N )
  160          CONTINUE
               WRITE ( NOUT, FMT = 99992 ) SCALE
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' SB03OD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from SB03OD = ',I2)
99997 FORMAT (' The transpose of the Cholesky factor U is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The solution matrix X = op(U)''*op(U) is ')
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' M is out of range.',/' M = ',I5)
99992 FORMAT (/' Scaling factor = ',F8.4)
      END
Program Data
 SB03OD EXAMPLE PROGRAM DATA
   4     5      C      N      N
  -1.0  37.0 -12.0 -12.0
  -1.0 -10.0   0.0   4.0
   2.0  -4.0   7.0  -6.0
   2.0   2.0   7.0  -9.0
   1.0   2.5   1.0   3.5
   0.0   1.0   0.0   1.0
  -1.0  -2.5  -1.0  -1.5
   1.0   2.5   4.0  -5.5
  -1.0  -2.5  -4.0   3.5
Program Results
 SB03OD EXAMPLE PROGRAM RESULTS

 The transpose of the Cholesky factor U is 
   1.0000
   3.0000   1.0000
   2.0000  -1.0000   1.0000
  -1.0000   1.0000  -2.0000   1.0000

 The solution matrix X = op(U)'*op(U) is 
   1.0000   3.0000   2.0000  -1.0000
   3.0000  10.0000   5.0000  -2.0000
   2.0000   5.0000   6.0000  -5.0000
  -1.0000  -2.0000  -5.0000   7.0000

 Scaling factor =   1.0000

Click here to get a compressed (gzip) tar file containing the source code of the routine, the example program, data, documentation, and related files.

Return to index