MB01WD

Residuals of Lyapunov or Stein equations for Cholesky factored solutions

Purpose

  To compute the matrix formula
  _
  R = alpha*( op( A )'*op( T )'*op( T ) + op( T )'*op( T )*op( A ) )
      + beta*R,                                                  (1)

  if DICO = 'C', or
  _
  R = alpha*( op( A )'*op( T )'*op( T )*op( A ) -  op( T )'*op( T ))
      + beta*R,                                                  (2)
                                                          _
  if DICO = 'D', where alpha and beta are scalars, R, and R are
  symmetric matrices, T is a triangular matrix, A is a general or
  Hessenberg matrix, and op( M ) is one of

     op( M ) = M   or   op( M ) = M'.

  The result is overwritten on R.

      SUBROUTINE MB01WD( DICO, UPLO, TRANS, HESS, N, ALPHA, BETA, R,
     $                   LDR, A, LDA, T, LDT, INFO )
C     .. Scalar Arguments ..
      CHARACTER         DICO, HESS, TRANS, UPLO
      INTEGER           INFO, LDA, LDR, LDT, N
      DOUBLE PRECISION  ALPHA, BETA
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), R(LDR,*), T(LDT,*)

Mode Parameters

  DICO    CHARACTER*1
          Specifies the formula to be evaluated, as follows:
          = 'C':  formula (1), "continuous-time" case;
          = 'D':  formula (2), "discrete-time" case.

  UPLO    CHARACTER*1
          Specifies which triangles of the symmetric matrix R and
          triangular matrix T are given, as follows:
          = 'U':  the upper triangular parts of R and T are given;
          = 'L':  the lower triangular parts of R and T are given;

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

  HESS    CHARACTER*1
          Specifies the form of the matrix A, as follows:
          = 'F':  matrix A is full;
          = 'H':  matrix A is Hessenberg (or Schur), either upper
                  (if UPLO = 'U'), or lower (if UPLO = 'L').

  N       (input) INTEGER
          The order of the matrices R, A, and T.  N >= 0.

  ALPHA   (input) DOUBLE PRECISION
          The scalar alpha. When alpha is zero then the arrays A
          and T are not referenced.

  BETA    (input) DOUBLE PRECISION
          The scalar beta. When beta is zero then the array R need
          not be set before entry.

  R       (input/output) DOUBLE PRECISION array, dimension (LDR,N)
          On entry with UPLO = 'U', the leading N-by-N upper
          triangular part of this array must contain the upper
          triangular part of the symmetric matrix R.
          On entry with UPLO = 'L', the leading N-by-N lower
          triangular part of this array must contain the lower
          triangular part of the symmetric matrix R.
          On exit, the leading N-by-N upper triangular part (if
          UPLO = 'U'), or lower triangular part (if UPLO = 'L'), of
          this array contains the corresponding triangular part of
                              _
          the computed matrix R.

  LDR     INTEGER
          The leading dimension of array R.  LDR >= MAX(1,N).

  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 HESS = 'H' the elements below the
          first subdiagonal, if UPLO = 'U', or above the first
          superdiagonal, if UPLO = 'L', need not be set to zero,
          and are not referenced if DICO = 'D'.
          On exit, the leading N-by-N part of this array contains
          the following matrix product
             alpha*T'*T*A, if TRANS = 'N', or
             alpha*A*T*T', otherwise,
          if DICO = 'C', or
             T*A, if TRANS = 'N', or
             A*T, otherwise,
          if DICO = 'D' (and in this case, these products have a
          Hessenberg form, if HESS = 'H').

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

  T       (input) DOUBLE PRECISION array, dimension (LDT,N)
          If UPLO = 'U', the leading N-by-N upper triangular part of
          this array must contain the upper triangular matrix T and
          the strictly lower triangular part need not be set to zero
          (and it is not referenced).
          If UPLO = 'L', the leading N-by-N lower triangular part of
          this array must contain the lower triangular matrix T and
          the strictly upper triangular part need not be set to zero
          (and it is not referenced).

  LDT     INTEGER
          The leading dimension of array T.  LDT >= MAX(1,N).

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -k, the k-th argument had an illegal
                value.

  The matrix expression (1) or (2) is efficiently evaluated taking
  the structure into account. BLAS 3 operations (DTRMM, DSYRK and
  their specializations) are used throughout.

  If A is a full matrix, the algorithm requires approximately
   3
  N  operations, if DICO = 'C';
         3
  7/6 x N  operations, if DICO = 'D'.

  None

Program Text

  None

  None

  None

Return to index