SB03RD

Solution of continuous-time Lyapunov equations and separation estimation

Purpose

  To solve the real Lyapunov matrix equation

         op(A)'*X + X*op(A) = scale*C

  and/or estimate the separation between the matrices op(A) and
  -op(A)', where op(A) = A or A' (A**T) and C is symmetric (C = C').
  (A' denotes the transpose of the matrix A.) A is N-by-N, the right
  hand side C and the solution X are N-by-N, and scale is an output
  scale factor, set less than or equal to 1 to avoid overflow in X.

      SUBROUTINE SB03RD( JOB, FACT, TRANA, N, A, LDA, U, LDU, C, LDC,
     $                   SCALE, SEP, FERR, WR, WI, IWORK, DWORK,
     $                   LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER          FACT, JOB, TRANA
      INTEGER            INFO, LDA, LDC, LDU, LDWORK, N
      DOUBLE PRECISION   FERR, SCALE, SEP
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), DWORK( * ),
     $                   U( LDU, * ), WI( * ), WR( * )

Mode Parameters

  JOB     CHARACTER*1
          Specifies the computation to be performed, as follows:
          = 'X':  Compute the solution only;
          = 'S':  Compute the separation only;
          = 'B':  Compute both the solution and the separation.

  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 U 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 U.

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

  N       (input) INTEGER
          The order of the matrices A, X, and C.  N >= 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 in Schur canonical form.
          On exit, if INFO = 0 or INFO = N+1, the leading N-by-N
          part of this array contains the upper quasi-triangular
          matrix 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).

  U       (input or output) DOUBLE PRECISION array, dimension
          (LDU,N)
          If FACT = 'F', then U is an input argument and on entry
          it must contain the orthogonal matrix U from the real
          Schur factorization of A.
          If FACT = 'N', then U is an output argument and on exit,
          if INFO = 0 or INFO = N+1, it contains the orthogonal
          N-by-N matrix from the real Schur factorization of A.

  LDU     INTEGER
          The leading dimension of array U.  LDU >= MAX(1,N).

  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
          On entry with JOB = 'X' or 'B', the leading N-by-N part of
          this array must contain the symmetric matrix C.
          On exit with JOB = 'X' or 'B', if INFO = 0 or INFO = N+1,
          the leading N-by-N part of C has been overwritten by the
          symmetric solution matrix X.
          If JOB = 'S', C is not referenced.

  LDC     INTEGER
          The leading dimension of array C.
          LDC >= 1,        if JOB = 'S';
          LDC >= MAX(1,N), otherwise.

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

  SEP     (output) DOUBLE PRECISION
          If JOB = 'S' or JOB = 'B', and INFO = 0 or INFO = N+1, SEP
          contains the estimated separation of the matrices op(A)
          and -op(A)'.
          If JOB = 'X' or N = 0, SEP is not referenced.

  FERR    (output) DOUBLE PRECISION
          If JOB = 'B', and INFO = 0 or INFO = N+1, FERR contains
          an estimated forward error bound for the solution X.
          If XTRUE is the true solution, FERR bounds the relative
          error in the computed solution, measured in the Frobenius
          norm:  norm(X - XTRUE)/norm(XTRUE).
          If JOB = 'X' or JOB = 'S', FERR is not referenced.

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

  IWORK   INTEGER array, dimension (N*N)
          This array is not referenced if JOB = 'X'.

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

  LDWORK  INTEGER
          The length of the array DWORK.  LDWORK >= 1 and
          If JOB = 'X' then
             If FACT = 'F', LDWORK >= N*N;
             If FACT = 'N', LDWORK >= MAX(N*N,3*N).
          If JOB = 'S' or JOB = 'B' then
             If FACT = 'F', LDWORK >= 2*N*N;
             If FACT = 'N', LDWORK >= MAX(2*N*N,3*N).
          For optimum performance LDWORK should be larger.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          > 0:  if INFO = i, the QR algorithm failed to compute all
                the eigenvalues (see LAPACK Library routine DGEES);
                elements i+1:n of WR and WI contain eigenvalues
                which have converged, and A contains the partially
                converged Schur form;
          = N+1:  if the matrices A and -A' have common or very
                close eigenvalues; perturbed values were used to
                solve the equation (but the matrix A is unchanged).

  After reducing matrix A to real Schur canonical form (if needed),
  the Bartels-Stewart algorithm is used. A set of equivalent linear
  algebraic systems of equations of order at most four are formed
  and solved using Gaussian elimination with complete pivoting.

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

                            3
  The algorithm requires 0(N ) operations.

  SEP is defined as the separation of op(A) and -op(A)':

         sep( op(A), -op(A)' ) = sigma_min( T )

  where sigma_min(T) is the smallest singular value of the
  N*N-by-N*N matrix

     T = kprod( I(N), op(A)' ) + kprod( op(A), I(N) ).

  I(N) is an N-by-N identity matrix, and kprod denotes the Kronecker
  product. The program estimates sigma_min(T) by the reciprocal of
  an estimate of the 1-norm of inverse(T). The true reciprocal
  1-norm of inverse(T) cannot differ from sigma_min(T) by more
  than a factor of N.

  When SEP is small, small changes in A, C can cause large changes
  in the solution of the equation. An approximate bound on the
  maximum relative error in the computed solution is

                         EPS * norm(A) / SEP

  where EPS is the machine precision.

Program Text

  None

  None

  None

Return to index