MB01LD

Computation of matrix expression alpha*R + beta*A*X*trans(A) with skew-symmetric matrices R and X

Purpose

  To compute the matrix formula
     _
     R = alpha*R + beta*op( A )*X*op( A )',
                                              _
  where alpha and beta are scalars, R, X, and R are skew-symmetric
  matrices, A is a general matrix, and op( A ) is one of

     op( A ) = A   or   op( A ) = A'.

  The result is overwritten on R.

      SUBROUTINE MB01LD( UPLO, TRANS, M, N, ALPHA, BETA, R, LDR, A, LDA,
     $                   X, LDX, DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER         TRANS, UPLO
      INTEGER           INFO, LDA, LDR, LDWORK, LDX, M, N
      DOUBLE PRECISION  ALPHA, BETA
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), DWORK(*), R(LDR,*), X(LDX,*)

Mode Parameters

  UPLO    CHARACTER*1
          Specifies which triangles of the skew-symmetric matrices R
          and X are given, as follows:
          = 'U':  the strictly upper triangular part is given;
          = 'L':  the strictly lower triangular part is given.

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

  M       (input) INTEGER           _
          The order of the matrices R and R and the number of rows
          of the matrix op( A ).  M >= 0.

  N       (input) INTEGER
          The order of the matrix X and the number of columns of the
          matrix op( A ).  N >= 0.

  ALPHA   (input) DOUBLE PRECISION
          The scalar alpha. When alpha is zero then R need not be
          set before entry, except when R is identified with X in
          the call.

  BETA    (input) DOUBLE PRECISION
          The scalar beta. When beta is zero or N <= 1, or M <= 1,
          then A and X are not referenced.

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

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

  A       (input) DOUBLE PRECISION array, dimension (LDA,k)
          where k is N when TRANS = 'N' and is M when TRANS = 'T' or
          TRANS = 'C'.
          On entry with TRANS = 'N', the leading M-by-N part of this
          array must contain the matrix A.
          On entry with TRANS = 'T' or TRANS = 'C', the leading
          N-by-M part of this array must contain the matrix A.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= MAX(1,k),
          where k is M when TRANS = 'N' and is N when TRANS = 'T' or
          TRANS = 'C'.

  X       (input or input/output) DOUBLE PRECISION array, dimension
          (LDX,K), where K = N, if UPLO = 'U' or  LDWORK >= M*(N-1),
               or K = MAX(N,M), if UPLO = 'L' and LDWORK <  M*(N-1).
          On entry, if UPLO = 'U', the leading N-by-N strictly upper
          triangular part of this array must contain the strictly
          upper triangular part of the skew-symmetric matrix X and
          the lower triangular part of the array is not referenced.
          On entry, if UPLO = 'L', the leading N-by-N strictly lower
          triangular part of this array must contain the strictly
          lower triangular part of the skew-symmetric matrix X and
          the upper triangular part of the array is not referenced.
          If LDWORK < M*(N-1), this array is overwritten with the
          matrix op(A)*X, if UPLO = 'U', or X*op(A)', if UPLO = 'L'.

  LDX     INTEGER
          The leading dimension of the array X.
          LDX >= MAX(1,N),   if UPLO = 'L' or  LDWORK >= M*(N-1);
          LDX >= MAX(1,N,M), if UPLO = 'U' and LDWORK <  M*(N-1).

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          This array is not referenced when beta = 0, or M <= 1, or
          N <= 1.

  LDWORK  The length of the array DWORK.
          LDWORK >= N, if  beta <> 0, and M > 0, and N >  1;
          LDWORK >= 0, if  beta =  0, or  M = 0, or  N <= 1.
          For optimum performance, LDWORK >= M*(N-1), if  beta <> 0,
          M > 1, and N > 1.

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

  The matrix expression is efficiently evaluated taking the skew-
  symmetry into account. If LDWORK >= M*(N-1), a BLAS 3 like
  implementation is used. Specifically, let X = T - T', with T a
  strictly upper or strictly lower triangular matrix, defined by

     T = striu( X ),  if UPLO = 'U',
     T = stril( X ),  if UPLO = 'L',

  where striu and stril denote the strictly upper triangular part
  and strictly lower triangular part of X, respectively. Then,

     A*X*A' = ( A*T )*A' - A*( A*T )',  for TRANS = 'N',
     A'*X*A = A'*( T*A ) - ( T*A )'*A,  for TRANS = 'T', or 'C',

  which involve BLAS 3 operations DTRMM and the skew-symmetric
  correspondent of DSYR2K (with a Fortran implementation available
  in the SLICOT Library routine MB01KD).
  If LDWORK < M*(N-1), a BLAS 2 implementation is used.

  The algorithm requires approximately

                2         2
     3/2 x M x N + 1/2 x M

  operations.

  None

Program Text

  None

  None

  None