double

Functions
subroutine	dgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
	DGEMM
subroutine	dsymm (side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
	DSYMM
subroutine	dsyr2k (uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
	DSYR2K
subroutine	dsyrk (uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
	DSYRK
subroutine	dtrmm (side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
	DTRMM
subroutine	dtrsm (side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
	DTRSM

Detailed Description

This is the group of double LEVEL 3 BLAS routines.

Function Documentation

dgemm()

subroutine dgemm	(	character	transa,
		character	transb,
		integer	m,
		integer	n,
		integer	k,
		double precision	alpha,
		double precision, dimension(lda,*)	a,
		integer	lda,
		double precision, dimension(ldb,*)	b,
		integer	ldb,
		double precision	beta,
		double precision, dimension(ldc,*)	c,
		integer	ldc )

DGEMM

Purpose:

!>
!> DGEMM  performs one of the matrix-matrix operations
!>
!>    C := alpha*op( A )*op( B ) + beta*C,
!>
!> where  op( X ) is one of
!>
!>    op( X ) = X   or   op( X ) = X**T,
!>
!> alpha and beta are scalars, and A, B and C are matrices, with op( A )
!> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
!> 

Parameters

[in]	TRANSA	!> TRANSA is CHARACTER*1 !> On entry, TRANSA specifies the form of op( A ) to be used in !> the matrix multiplication as follows: !> !> TRANSA = 'N' or 'n', op( A ) = A. !> !> TRANSA = 'T' or 't', op( A ) = A**T. !> !> TRANSA = 'C' or 'c', op( A ) = A**T. !>
[in]	TRANSB	!> TRANSB is CHARACTER*1 !> On entry, TRANSB specifies the form of op( B ) to be used in !> the matrix multiplication as follows: !> !> TRANSB = 'N' or 'n', op( B ) = B. !> !> TRANSB = 'T' or 't', op( B ) = B**T. !> !> TRANSB = 'C' or 'c', op( B ) = B**T. !>
[in]	M	!> M is INTEGER !> On entry, M specifies the number of rows of the matrix !> op( A ) and of the matrix C. M must be at least zero. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the number of columns of the matrix !> op( B ) and the number of columns of the matrix C. N must be !> at least zero. !>
[in]	K	!> K is INTEGER !> On entry, K specifies the number of columns of the matrix !> op( A ) and the number of rows of the matrix op( B ). K must !> be at least zero. !>
[in]	ALPHA	!> ALPHA is DOUBLE PRECISION. !> On entry, ALPHA specifies the scalar alpha. !>
[in]	A	!> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is !> k when TRANSA = 'N' or 'n', and is m otherwise. !> Before entry with TRANSA = 'N' or 'n', the leading m by k !> part of the array A must contain the matrix A, otherwise !> the leading k by m part of the array A must contain the !> matrix A. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. When TRANSA = 'N' or 'n' then !> LDA must be at least max( 1, m ), otherwise LDA must be at !> least max( 1, k ). !>
[in]	B	!> B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is !> n when TRANSB = 'N' or 'n', and is k otherwise. !> Before entry with TRANSB = 'N' or 'n', the leading k by n !> part of the array B must contain the matrix B, otherwise !> the leading n by k part of the array B must contain the !> matrix B. !>
[in]	LDB	!> LDB is INTEGER !> On entry, LDB specifies the first dimension of B as declared !> in the calling (sub) program. When TRANSB = 'N' or 'n' then !> LDB must be at least max( 1, k ), otherwise LDB must be at !> least max( 1, n ). !>
[in]	BETA	!> BETA is DOUBLE PRECISION. !> On entry, BETA specifies the scalar beta. When BETA is !> supplied as zero then C need not be set on input. !>
[in,out]	C	!> C is DOUBLE PRECISION array, dimension ( LDC, N ) !> Before entry, the leading m by n part of the array C must !> contain the matrix C, except when beta is zero, in which !> case C need not be set on entry. !> On exit, the array C is overwritten by the m by n matrix !> ( alpha*op( A )*op( B ) + beta*C ). !>
[in]	LDC	!> LDC is INTEGER !> On entry, LDC specifies the first dimension of C as declared !> in the calling (sub) program. LDC must be at least !> max( 1, m ). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  Level 3 Blas routine.
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> 

Definition at line 186 of file dgemm.f.

*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA,BETA
      INTEGER K,LDA,LDB,LDC,M,N
      CHARACTER TRANSA,TRANSB
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER I,INFO,J,L,NROWA,NROWB
      LOGICAL NOTA,NOTB
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      parameter(one=1.0d+0,zero=0.0d+0)
*     ..
*
*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
*     transposed and set  NROWA and NROWB  as the number of rows of  A
*     and  B  respectively.
*
      nota = lsame(transa,'N')
      notb = lsame(transb,'N')
      IF (nota) THEN
          nrowa = m
      ELSE
          nrowa = k
      END IF
      IF (notb) THEN
          nrowb = k
      ELSE
          nrowb = n
      END IF
*
*     Test the input parameters.
*
      info = 0
      IF ((.NOT.nota) .AND. (.NOT.lsame(transa,'C')) .AND.
     +    (.NOT.lsame(transa,'T'))) THEN
          info = 1
      ELSE IF ((.NOT.notb) .AND. (.NOT.lsame(transb,'C')) .AND.
     +         (.NOT.lsame(transb,'T'))) THEN
          info = 2
      ELSE IF (m.LT.0) THEN
          info = 3
      ELSE IF (n.LT.0) THEN
          info = 4
      ELSE IF (k.LT.0) THEN
          info = 5
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 8
      ELSE IF (ldb.LT.max(1,nrowb)) THEN
          info = 10
      ELSE IF (ldc.LT.max(1,m)) THEN
          info = 13
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DGEMM ',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF ((m.EQ.0) .OR. (n.EQ.0) .OR.
     +    (((alpha.EQ.zero).OR. (k.EQ.0)).AND. (beta.EQ.one))) RETURN
*
*     And if  alpha.eq.zero.
*
      IF (alpha.EQ.zero) THEN
          IF (beta.EQ.zero) THEN
              DO 20 j = 1,n
                  DO 10 i = 1,m
                      c(i,j) = zero
   10             CONTINUE
   20         CONTINUE
          ELSE
              DO 40 j = 1,n
                  DO 30 i = 1,m
                      c(i,j) = beta*c(i,j)
   30             CONTINUE
   40         CONTINUE
          END IF
          RETURN
      END IF
*
*     Start the operations.
*
      IF (notb) THEN
          IF (nota) THEN
*
*           Form  C := alpha*A*B + beta*C.
*
              DO 90 j = 1,n
                  IF (beta.EQ.zero) THEN
                      DO 50 i = 1,m
                          c(i,j) = zero
   50                 CONTINUE
                  ELSE IF (beta.NE.one) THEN
                      DO 60 i = 1,m
                          c(i,j) = beta*c(i,j)
   60                 CONTINUE
                  END IF
                  DO 80 l = 1,k
                      temp = alpha*b(l,j)
                      DO 70 i = 1,m
                          c(i,j) = c(i,j) + temp*a(i,l)
   70                 CONTINUE
   80             CONTINUE
   90         CONTINUE
          ELSE
*
*           Form  C := alpha*A**T*B + beta*C
*
              DO 120 j = 1,n
                  DO 110 i = 1,m
                      temp = zero
                      DO 100 l = 1,k
                          temp = temp + a(l,i)*b(l,j)
  100                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = alpha*temp
                      ELSE
                          c(i,j) = alpha*temp + beta*c(i,j)
                      END IF
  110             CONTINUE
  120         CONTINUE
          END IF
      ELSE
          IF (nota) THEN
*
*           Form  C := alpha*A*B**T + beta*C
*
              DO 170 j = 1,n
                  IF (beta.EQ.zero) THEN
                      DO 130 i = 1,m
                          c(i,j) = zero
  130                 CONTINUE
                  ELSE IF (beta.NE.one) THEN
                      DO 140 i = 1,m
                          c(i,j) = beta*c(i,j)
  140                 CONTINUE
                  END IF
                  DO 160 l = 1,k
                      temp = alpha*b(j,l)
                      DO 150 i = 1,m
                          c(i,j) = c(i,j) + temp*a(i,l)
  150                 CONTINUE
  160             CONTINUE
  170         CONTINUE
          ELSE
*
*           Form  C := alpha*A**T*B**T + beta*C
*
              DO 200 j = 1,n
                  DO 190 i = 1,m
                      temp = zero
                      DO 180 l = 1,k
                          temp = temp + a(l,i)*b(j,l)
  180                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = alpha*temp
                      ELSE
                          c(i,j) = alpha*temp + beta*c(i,j)
                      END IF
  190             CONTINUE
  200         CONTINUE
          END IF
      END IF
*
      RETURN
*
*     End of DGEMM
*

dsymm()

subroutine dsymm	(	character	side,
		character	uplo,
		integer	m,
		integer	n,
		double precision	alpha,
		double precision, dimension(lda,*)	a,
		integer	lda,
		double precision, dimension(ldb,*)	b,
		integer	ldb,
		double precision	beta,
		double precision, dimension(ldc,*)	c,
		integer	ldc )

DSYMM

Purpose:

!>
!> DSYMM  performs one of the matrix-matrix operations
!>
!>    C := alpha*A*B + beta*C,
!>
!> or
!>
!>    C := alpha*B*A + beta*C,
!>
!> where alpha and beta are scalars,  A is a symmetric matrix and  B and
!> C are  m by n matrices.
!> 

Parameters

[in]	SIDE	!> SIDE is CHARACTER*1 !> On entry, SIDE specifies whether the symmetric matrix A !> appears on the left or right in the operation as follows: !> !> SIDE = 'L' or 'l' C := alpha*A*B + beta*C, !> !> SIDE = 'R' or 'r' C := alpha*B*A + beta*C, !>
[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the upper or lower !> triangular part of the symmetric matrix A is to be !> referenced as follows: !> !> UPLO = 'U' or 'u' Only the upper triangular part of the !> symmetric matrix is to be referenced. !> !> UPLO = 'L' or 'l' Only the lower triangular part of the !> symmetric matrix is to be referenced. !>
[in]	M	!> M is INTEGER !> On entry, M specifies the number of rows of the matrix C. !> M must be at least zero. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the number of columns of the matrix C. !> N must be at least zero. !>
[in]	ALPHA	!> ALPHA is DOUBLE PRECISION. !> On entry, ALPHA specifies the scalar alpha. !>
[in]	A	!> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is !> m when SIDE = 'L' or 'l' and is n otherwise. !> Before entry with SIDE = 'L' or 'l', the m by m part of !> the array A must contain the symmetric matrix, such that !> when UPLO = 'U' or 'u', the leading m by m upper triangular !> part of the array A must contain the upper triangular part !> of the symmetric matrix and the strictly lower triangular !> part of A is not referenced, and when UPLO = 'L' or 'l', !> the leading m by m lower triangular part of the array A !> must contain the lower triangular part of the symmetric !> matrix and the strictly upper triangular part of A is not !> referenced. !> Before entry with SIDE = 'R' or 'r', the n by n part of !> the array A must contain the symmetric matrix, such that !> when UPLO = 'U' or 'u', the leading n by n upper triangular !> part of the array A must contain the upper triangular part !> of the symmetric matrix and the strictly lower triangular !> part of A is not referenced, and when UPLO = 'L' or 'l', !> the leading n by n lower triangular part of the array A !> must contain the lower triangular part of the symmetric !> matrix and the strictly upper triangular part of A is not !> referenced. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. When SIDE = 'L' or 'l' then !> LDA must be at least max( 1, m ), otherwise LDA must be at !> least max( 1, n ). !>
[in]	B	!> B is DOUBLE PRECISION array, dimension ( LDB, N ) !> Before entry, the leading m by n part of the array B must !> contain the matrix B. !>
[in]	LDB	!> LDB is INTEGER !> On entry, LDB specifies the first dimension of B as declared !> in the calling (sub) program. LDB must be at least !> max( 1, m ). !>
[in]	BETA	!> BETA is DOUBLE PRECISION. !> On entry, BETA specifies the scalar beta. When BETA is !> supplied as zero then C need not be set on input. !>
[in,out]	C	!> C is DOUBLE PRECISION array, dimension ( LDC, N ) !> Before entry, the leading m by n part of the array C must !> contain the matrix C, except when beta is zero, in which !> case C need not be set on entry. !> On exit, the array C is overwritten by the m by n updated !> matrix. !>
[in]	LDC	!> LDC is INTEGER !> On entry, LDC specifies the first dimension of C as declared !> in the calling (sub) program. LDC must be at least !> max( 1, m ). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  Level 3 Blas routine.
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> 

Definition at line 188 of file dsymm.f.

*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA,BETA
      INTEGER LDA,LDB,LDC,M,N
      CHARACTER SIDE,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP1,TEMP2
      INTEGER I,INFO,J,K,NROWA
      LOGICAL UPPER
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      parameter(one=1.0d+0,zero=0.0d+0)
*     ..
*
*     Set NROWA as the number of rows of A.
*
      IF (lsame(side,'L')) THEN
          nrowa = m
      ELSE
          nrowa = n
      END IF
      upper = lsame(uplo,'U')
*
*     Test the input parameters.
*
      info = 0
      IF ((.NOT.lsame(side,'L')) .AND. (.NOT.lsame(side,'R'))) THEN
          info = 1
      ELSE IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
          info = 2
      ELSE IF (m.LT.0) THEN
          info = 3
      ELSE IF (n.LT.0) THEN
          info = 4
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 7
      ELSE IF (ldb.LT.max(1,m)) THEN
          info = 9
      ELSE IF (ldc.LT.max(1,m)) THEN
          info = 12
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DSYMM ',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF ((m.EQ.0) .OR. (n.EQ.0) .OR.
     +    ((alpha.EQ.zero).AND. (beta.EQ.one))) RETURN
*
*     And when  alpha.eq.zero.
*
      IF (alpha.EQ.zero) THEN
          IF (beta.EQ.zero) THEN
              DO 20 j = 1,n
                  DO 10 i = 1,m
                      c(i,j) = zero
   10             CONTINUE
   20         CONTINUE
          ELSE
              DO 40 j = 1,n
                  DO 30 i = 1,m
                      c(i,j) = beta*c(i,j)
   30             CONTINUE
   40         CONTINUE
          END IF
          RETURN
      END IF
*
*     Start the operations.
*
      IF (lsame(side,'L')) THEN
*
*        Form  C := alpha*A*B + beta*C.
*
          IF (upper) THEN
              DO 70 j = 1,n
                  DO 60 i = 1,m
                      temp1 = alpha*b(i,j)
                      temp2 = zero
                      DO 50 k = 1,i - 1
                          c(k,j) = c(k,j) + temp1*a(k,i)
                          temp2 = temp2 + b(k,j)*a(k,i)
   50                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = temp1*a(i,i) + alpha*temp2
                      ELSE
                          c(i,j) = beta*c(i,j) + temp1*a(i,i) +
     +                             alpha*temp2
                      END IF
   60             CONTINUE
   70         CONTINUE
          ELSE
              DO 100 j = 1,n
                  DO 90 i = m,1,-1
                      temp1 = alpha*b(i,j)
                      temp2 = zero
                      DO 80 k = i + 1,m
                          c(k,j) = c(k,j) + temp1*a(k,i)
                          temp2 = temp2 + b(k,j)*a(k,i)
   80                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = temp1*a(i,i) + alpha*temp2
                      ELSE
                          c(i,j) = beta*c(i,j) + temp1*a(i,i) +
     +                             alpha*temp2
                      END IF
   90             CONTINUE
  100         CONTINUE
          END IF
      ELSE
*
*        Form  C := alpha*B*A + beta*C.
*
          DO 170 j = 1,n
              temp1 = alpha*a(j,j)
              IF (beta.EQ.zero) THEN
                  DO 110 i = 1,m
                      c(i,j) = temp1*b(i,j)
  110             CONTINUE
              ELSE
                  DO 120 i = 1,m
                      c(i,j) = beta*c(i,j) + temp1*b(i,j)
  120             CONTINUE
              END IF
              DO 140 k = 1,j - 1
                  IF (upper) THEN
                      temp1 = alpha*a(k,j)
                  ELSE
                      temp1 = alpha*a(j,k)
                  END IF
                  DO 130 i = 1,m
                      c(i,j) = c(i,j) + temp1*b(i,k)
  130             CONTINUE
  140         CONTINUE
              DO 160 k = j + 1,n
                  IF (upper) THEN
                      temp1 = alpha*a(j,k)
                  ELSE
                      temp1 = alpha*a(k,j)
                  END IF
                  DO 150 i = 1,m
                      c(i,j) = c(i,j) + temp1*b(i,k)
  150             CONTINUE
  160         CONTINUE
  170     CONTINUE
      END IF
*
      RETURN
*
*     End of DSYMM
*

dsyr2k()

subroutine dsyr2k	(	character	uplo,
		character	trans,
		integer	n,
		integer	k,
		double precision	alpha,
		double precision, dimension(lda,*)	a,
		integer	lda,
		double precision, dimension(ldb,*)	b,
		integer	ldb,
		double precision	beta,
		double precision, dimension(ldc,*)	c,
		integer	ldc )

DSYR2K

Purpose:

!>
!> DSYR2K  performs one of the symmetric rank 2k operations
!>
!>    C := alpha*A*B**T + alpha*B*A**T + beta*C,
!>
!> or
!>
!>    C := alpha*A**T*B + alpha*B**T*A + beta*C,
!>
!> where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
!> and  A and B  are  n by k  matrices  in the  first  case  and  k by n
!> matrices in the second case.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the upper or lower !> triangular part of the array C is to be referenced as !> follows: !> !> UPLO = 'U' or 'u' Only the upper triangular part of C !> is to be referenced. !> !> UPLO = 'L' or 'l' Only the lower triangular part of C !> is to be referenced. !>
[in]	TRANS	!> TRANS is CHARACTER*1 !> On entry, TRANS specifies the operation to be performed as !> follows: !> !> TRANS = 'N' or 'n' C := alpha*A*B**T + alpha*B*A**T + !> beta*C. !> !> TRANS = 'T' or 't' C := alpha*A**T*B + alpha*B**T*A + !> beta*C. !> !> TRANS = 'C' or 'c' C := alpha*A**T*B + alpha*B**T*A + !> beta*C. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the order of the matrix C. N must be !> at least zero. !>
[in]	K	!> K is INTEGER !> On entry with TRANS = 'N' or 'n', K specifies the number !> of columns of the matrices A and B, and on entry with !> TRANS = 'T' or 't' or 'C' or 'c', K specifies the number !> of rows of the matrices A and B. K must be at least zero. !>
[in]	ALPHA	!> ALPHA is DOUBLE PRECISION. !> On entry, ALPHA specifies the scalar alpha. !>
[in]	A	!> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is !> k when TRANS = 'N' or 'n', and is n otherwise. !> Before entry with TRANS = 'N' or 'n', the leading n by k !> part of the array A must contain the matrix A, otherwise !> the leading k by n part of the array A must contain the !> matrix A. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. When TRANS = 'N' or 'n' !> then LDA must be at least max( 1, n ), otherwise LDA must !> be at least max( 1, k ). !>
[in]	B	!> B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is !> k when TRANS = 'N' or 'n', and is n otherwise. !> Before entry with TRANS = 'N' or 'n', the leading n by k !> part of the array B must contain the matrix B, otherwise !> the leading k by n part of the array B must contain the !> matrix B. !>
[in]	LDB	!> LDB is INTEGER !> On entry, LDB specifies the first dimension of B as declared !> in the calling (sub) program. When TRANS = 'N' or 'n' !> then LDB must be at least max( 1, n ), otherwise LDB must !> be at least max( 1, k ). !>
[in]	BETA	!> BETA is DOUBLE PRECISION. !> On entry, BETA specifies the scalar beta. !>
[in,out]	C	!> C is DOUBLE PRECISION array, dimension ( LDC, N ) !> Before entry with UPLO = 'U' or 'u', the leading n by n !> upper triangular part of the array C must contain the upper !> triangular part of the symmetric matrix and the strictly !> lower triangular part of C is not referenced. On exit, the !> upper triangular part of the array C is overwritten by the !> upper triangular part of the updated matrix. !> Before entry with UPLO = 'L' or 'l', the leading n by n !> lower triangular part of the array C must contain the lower !> triangular part of the symmetric matrix and the strictly !> upper triangular part of C is not referenced. On exit, the !> lower triangular part of the array C is overwritten by the !> lower triangular part of the updated matrix. !>
[in]	LDC	!> LDC is INTEGER !> On entry, LDC specifies the first dimension of C as declared !> in the calling (sub) program. LDC must be at least !> max( 1, n ). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  Level 3 Blas routine.
!>
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> 

Definition at line 191 of file dsyr2k.f.

*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA,BETA
      INTEGER K,LDA,LDB,LDC,N
      CHARACTER TRANS,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP1,TEMP2
      INTEGER I,INFO,J,L,NROWA
      LOGICAL UPPER
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      parameter(one=1.0d+0,zero=0.0d+0)
*     ..
*
*     Test the input parameters.
*
      IF (lsame(trans,'N')) THEN
          nrowa = n
      ELSE
          nrowa = k
      END IF
      upper = lsame(uplo,'U')
*
      info = 0
      IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
          info = 1
      ELSE IF ((.NOT.lsame(trans,'N')) .AND.
     +         (.NOT.lsame(trans,'T')) .AND.
     +         (.NOT.lsame(trans,'C'))) THEN
          info = 2
      ELSE IF (n.LT.0) THEN
          info = 3
      ELSE IF (k.LT.0) THEN
          info = 4
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 7
      ELSE IF (ldb.LT.max(1,nrowa)) THEN
          info = 9
      ELSE IF (ldc.LT.max(1,n)) THEN
          info = 12
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DSYR2K',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF ((n.EQ.0) .OR. (((alpha.EQ.zero).OR.
     +    (k.EQ.0)).AND. (beta.EQ.one))) RETURN
*
*     And when  alpha.eq.zero.
*
      IF (alpha.EQ.zero) THEN
          IF (upper) THEN
              IF (beta.EQ.zero) THEN
                  DO 20 j = 1,n
                      DO 10 i = 1,j
                          c(i,j) = zero
   10                 CONTINUE
   20             CONTINUE
              ELSE
                  DO 40 j = 1,n
                      DO 30 i = 1,j
                          c(i,j) = beta*c(i,j)
   30                 CONTINUE
   40             CONTINUE
              END IF
          ELSE
              IF (beta.EQ.zero) THEN
                  DO 60 j = 1,n
                      DO 50 i = j,n
                          c(i,j) = zero
   50                 CONTINUE
   60             CONTINUE
              ELSE
                  DO 80 j = 1,n
                      DO 70 i = j,n
                          c(i,j) = beta*c(i,j)
   70                 CONTINUE
   80             CONTINUE
              END IF
          END IF
          RETURN
      END IF
*
*     Start the operations.
*
      IF (lsame(trans,'N')) THEN
*
*        Form  C := alpha*A*B**T + alpha*B*A**T + C.
*
          IF (upper) THEN
              DO 130 j = 1,n
                  IF (beta.EQ.zero) THEN
                      DO 90 i = 1,j
                          c(i,j) = zero
   90                 CONTINUE
                  ELSE IF (beta.NE.one) THEN
                      DO 100 i = 1,j
                          c(i,j) = beta*c(i,j)
  100                 CONTINUE
                  END IF
                  DO 120 l = 1,k
                      IF ((a(j,l).NE.zero) .OR. (b(j,l).NE.zero)) THEN
                          temp1 = alpha*b(j,l)
                          temp2 = alpha*a(j,l)
                          DO 110 i = 1,j
                              c(i,j) = c(i,j) + a(i,l)*temp1 +
     +                                 b(i,l)*temp2
  110                     CONTINUE
                      END IF
  120             CONTINUE
  130         CONTINUE
          ELSE
              DO 180 j = 1,n
                  IF (beta.EQ.zero) THEN
                      DO 140 i = j,n
                          c(i,j) = zero
  140                 CONTINUE
                  ELSE IF (beta.NE.one) THEN
                      DO 150 i = j,n
                          c(i,j) = beta*c(i,j)
  150                 CONTINUE
                  END IF
                  DO 170 l = 1,k
                      IF ((a(j,l).NE.zero) .OR. (b(j,l).NE.zero)) THEN
                          temp1 = alpha*b(j,l)
                          temp2 = alpha*a(j,l)
                          DO 160 i = j,n
                              c(i,j) = c(i,j) + a(i,l)*temp1 +
     +                                 b(i,l)*temp2
  160                     CONTINUE
                      END IF
  170             CONTINUE
  180         CONTINUE
          END IF
      ELSE
*
*        Form  C := alpha*A**T*B + alpha*B**T*A + C.
*
          IF (upper) THEN
              DO 210 j = 1,n
                  DO 200 i = 1,j
                      temp1 = zero
                      temp2 = zero
                      DO 190 l = 1,k
                          temp1 = temp1 + a(l,i)*b(l,j)
                          temp2 = temp2 + b(l,i)*a(l,j)
  190                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = alpha*temp1 + alpha*temp2
                      ELSE
                          c(i,j) = beta*c(i,j) + alpha*temp1 +
     +                             alpha*temp2
                      END IF
  200             CONTINUE
  210         CONTINUE
          ELSE
              DO 240 j = 1,n
                  DO 230 i = j,n
                      temp1 = zero
                      temp2 = zero
                      DO 220 l = 1,k
                          temp1 = temp1 + a(l,i)*b(l,j)
                          temp2 = temp2 + b(l,i)*a(l,j)
  220                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = alpha*temp1 + alpha*temp2
                      ELSE
                          c(i,j) = beta*c(i,j) + alpha*temp1 +
     +                             alpha*temp2
                      END IF
  230             CONTINUE
  240         CONTINUE
          END IF
      END IF
*
      RETURN
*
*     End of DSYR2K
*

dsyrk()

subroutine dsyrk	(	character	uplo,
		character	trans,
		integer	n,
		integer	k,
		double precision	alpha,
		double precision, dimension(lda,*)	a,
		integer	lda,
		double precision	beta,
		double precision, dimension(ldc,*)	c,
		integer	ldc )

DSYRK

Purpose:

!>
!> DSYRK  performs one of the symmetric rank k operations
!>
!>    C := alpha*A*A**T + beta*C,
!>
!> or
!>
!>    C := alpha*A**T*A + beta*C,
!>
!> where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
!> and  A  is an  n by k  matrix in the first case and a  k by n  matrix
!> in the second case.
!> 

Parameters

[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the upper or lower !> triangular part of the array C is to be referenced as !> follows: !> !> UPLO = 'U' or 'u' Only the upper triangular part of C !> is to be referenced. !> !> UPLO = 'L' or 'l' Only the lower triangular part of C !> is to be referenced. !>
[in]	TRANS	!> TRANS is CHARACTER*1 !> On entry, TRANS specifies the operation to be performed as !> follows: !> !> TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C. !> !> TRANS = 'T' or 't' C := alpha*A**T*A + beta*C. !> !> TRANS = 'C' or 'c' C := alpha*A**T*A + beta*C. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the order of the matrix C. N must be !> at least zero. !>
[in]	K	!> K is INTEGER !> On entry with TRANS = 'N' or 'n', K specifies the number !> of columns of the matrix A, and on entry with !> TRANS = 'T' or 't' or 'C' or 'c', K specifies the number !> of rows of the matrix A. K must be at least zero. !>
[in]	ALPHA	!> ALPHA is DOUBLE PRECISION. !> On entry, ALPHA specifies the scalar alpha. !>
[in]	A	!> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is !> k when TRANS = 'N' or 'n', and is n otherwise. !> Before entry with TRANS = 'N' or 'n', the leading n by k !> part of the array A must contain the matrix A, otherwise !> the leading k by n part of the array A must contain the !> matrix A. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. When TRANS = 'N' or 'n' !> then LDA must be at least max( 1, n ), otherwise LDA must !> be at least max( 1, k ). !>
[in]	BETA	!> BETA is DOUBLE PRECISION. !> On entry, BETA specifies the scalar beta. !>
[in,out]	C	!> C is DOUBLE PRECISION array, dimension ( LDC, N ) !> Before entry with UPLO = 'U' or 'u', the leading n by n !> upper triangular part of the array C must contain the upper !> triangular part of the symmetric matrix and the strictly !> lower triangular part of C is not referenced. On exit, the !> upper triangular part of the array C is overwritten by the !> upper triangular part of the updated matrix. !> Before entry with UPLO = 'L' or 'l', the leading n by n !> lower triangular part of the array C must contain the lower !> triangular part of the symmetric matrix and the strictly !> upper triangular part of C is not referenced. On exit, the !> lower triangular part of the array C is overwritten by the !> lower triangular part of the updated matrix. !>
[in]	LDC	!> LDC is INTEGER !> On entry, LDC specifies the first dimension of C as declared !> in the calling (sub) program. LDC must be at least !> max( 1, n ). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  Level 3 Blas routine.
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> 

Definition at line 168 of file dsyrk.f.

*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA,BETA
      INTEGER K,LDA,LDC,N
      CHARACTER TRANS,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*),C(LDC,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER I,INFO,J,L,NROWA
      LOGICAL UPPER
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      parameter(one=1.0d+0,zero=0.0d+0)
*     ..
*
*     Test the input parameters.
*
      IF (lsame(trans,'N')) THEN
          nrowa = n
      ELSE
          nrowa = k
      END IF
      upper = lsame(uplo,'U')
*
      info = 0
      IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
          info = 1
      ELSE IF ((.NOT.lsame(trans,'N')) .AND.
     +         (.NOT.lsame(trans,'T')) .AND.
     +         (.NOT.lsame(trans,'C'))) THEN
          info = 2
      ELSE IF (n.LT.0) THEN
          info = 3
      ELSE IF (k.LT.0) THEN
          info = 4
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 7
      ELSE IF (ldc.LT.max(1,n)) THEN
          info = 10
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DSYRK ',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF ((n.EQ.0) .OR. (((alpha.EQ.zero).OR.
     +    (k.EQ.0)).AND. (beta.EQ.one))) RETURN
*
*     And when  alpha.eq.zero.
*
      IF (alpha.EQ.zero) THEN
          IF (upper) THEN
              IF (beta.EQ.zero) THEN
                  DO 20 j = 1,n
                      DO 10 i = 1,j
                          c(i,j) = zero
   10                 CONTINUE
   20             CONTINUE
              ELSE
                  DO 40 j = 1,n
                      DO 30 i = 1,j
                          c(i,j) = beta*c(i,j)
   30                 CONTINUE
   40             CONTINUE
              END IF
          ELSE
              IF (beta.EQ.zero) THEN
                  DO 60 j = 1,n
                      DO 50 i = j,n
                          c(i,j) = zero
   50                 CONTINUE
   60             CONTINUE
              ELSE
                  DO 80 j = 1,n
                      DO 70 i = j,n
                          c(i,j) = beta*c(i,j)
   70                 CONTINUE
   80             CONTINUE
              END IF
          END IF
          RETURN
      END IF
*
*     Start the operations.
*
      IF (lsame(trans,'N')) THEN
*
*        Form  C := alpha*A*A**T + beta*C.
*
          IF (upper) THEN
              DO 130 j = 1,n
                  IF (beta.EQ.zero) THEN
                      DO 90 i = 1,j
                          c(i,j) = zero
   90                 CONTINUE
                  ELSE IF (beta.NE.one) THEN
                      DO 100 i = 1,j
                          c(i,j) = beta*c(i,j)
  100                 CONTINUE
                  END IF
                  DO 120 l = 1,k
                      IF (a(j,l).NE.zero) THEN
                          temp = alpha*a(j,l)
                          DO 110 i = 1,j
                              c(i,j) = c(i,j) + temp*a(i,l)
  110                     CONTINUE
                      END IF
  120             CONTINUE
  130         CONTINUE
          ELSE
              DO 180 j = 1,n
                  IF (beta.EQ.zero) THEN
                      DO 140 i = j,n
                          c(i,j) = zero
  140                 CONTINUE
                  ELSE IF (beta.NE.one) THEN
                      DO 150 i = j,n
                          c(i,j) = beta*c(i,j)
  150                 CONTINUE
                  END IF
                  DO 170 l = 1,k
                      IF (a(j,l).NE.zero) THEN
                          temp = alpha*a(j,l)
                          DO 160 i = j,n
                              c(i,j) = c(i,j) + temp*a(i,l)
  160                     CONTINUE
                      END IF
  170             CONTINUE
  180         CONTINUE
          END IF
      ELSE
*
*        Form  C := alpha*A**T*A + beta*C.
*
          IF (upper) THEN
              DO 210 j = 1,n
                  DO 200 i = 1,j
                      temp = zero
                      DO 190 l = 1,k
                          temp = temp + a(l,i)*a(l,j)
  190                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = alpha*temp
                      ELSE
                          c(i,j) = alpha*temp + beta*c(i,j)
                      END IF
  200             CONTINUE
  210         CONTINUE
          ELSE
              DO 240 j = 1,n
                  DO 230 i = j,n
                      temp = zero
                      DO 220 l = 1,k
                          temp = temp + a(l,i)*a(l,j)
  220                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = alpha*temp
                      ELSE
                          c(i,j) = alpha*temp + beta*c(i,j)
                      END IF
  230             CONTINUE
  240         CONTINUE
          END IF
      END IF
*
      RETURN
*
*     End of DSYRK
*

dtrmm()

subroutine dtrmm	(	character	side,
		character	uplo,
		character	transa,
		character	diag,
		integer	m,
		integer	n,
		double precision	alpha,
		double precision, dimension(lda,*)	a,
		integer	lda,
		double precision, dimension(ldb,*)	b,
		integer	ldb )

DTRMM

Purpose:

!>
!> DTRMM  performs one of the matrix-matrix operations
!>
!>    B := alpha*op( A )*B,   or   B := alpha*B*op( A ),
!>
!> where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
!> non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
!>
!>    op( A ) = A   or   op( A ) = A**T.
!> 

Parameters

[in]	SIDE	!> SIDE is CHARACTER*1 !> On entry, SIDE specifies whether op( A ) multiplies B from !> the left or right as follows: !> !> SIDE = 'L' or 'l' B := alpha*op( A )*B. !> !> SIDE = 'R' or 'r' B := alpha*B*op( A ). !>
[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the matrix A is an upper or !> lower triangular matrix as follows: !> !> UPLO = 'U' or 'u' A is an upper triangular matrix. !> !> UPLO = 'L' or 'l' A is a lower triangular matrix. !>
[in]	TRANSA	!> TRANSA is CHARACTER*1 !> On entry, TRANSA specifies the form of op( A ) to be used in !> the matrix multiplication as follows: !> !> TRANSA = 'N' or 'n' op( A ) = A. !> !> TRANSA = 'T' or 't' op( A ) = A**T. !> !> TRANSA = 'C' or 'c' op( A ) = A**T. !>
[in]	DIAG	!> DIAG is CHARACTER*1 !> On entry, DIAG specifies whether or not A is unit triangular !> as follows: !> !> DIAG = 'U' or 'u' A is assumed to be unit triangular. !> !> DIAG = 'N' or 'n' A is not assumed to be unit !> triangular. !>
[in]	M	!> M is INTEGER !> On entry, M specifies the number of rows of B. M must be at !> least zero. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the number of columns of B. N must be !> at least zero. !>
[in]	ALPHA	!> ALPHA is DOUBLE PRECISION. !> On entry, ALPHA specifies the scalar alpha. When alpha is !> zero then A is not referenced and B need not be set before !> entry. !>
[in]	A	!> A is DOUBLE PRECISION array, dimension ( LDA, k ), where k is m !> when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. !> Before entry with UPLO = 'U' or 'u', the leading k by k !> upper triangular part of the array A must contain the upper !> triangular matrix and the strictly lower triangular part of !> A is not referenced. !> Before entry with UPLO = 'L' or 'l', the leading k by k !> lower triangular part of the array A must contain the lower !> triangular matrix and the strictly upper triangular part of !> A is not referenced. !> Note that when DIAG = 'U' or 'u', the diagonal elements of !> A are not referenced either, but are assumed to be unity. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. When SIDE = 'L' or 'l' then !> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' !> then LDA must be at least max( 1, n ). !>
[in,out]	B	!> B is DOUBLE PRECISION array, dimension ( LDB, N ) !> Before entry, the leading m by n part of the array B must !> contain the matrix B, and on exit is overwritten by the !> transformed matrix. !>
[in]	LDB	!> LDB is INTEGER !> On entry, LDB specifies the first dimension of B as declared !> in the calling (sub) program. LDB must be at least !> max( 1, m ). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  Level 3 Blas routine.
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> 

Definition at line 176 of file dtrmm.f.

*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*),B(LDB,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER I,INFO,J,K,NROWA
      LOGICAL LSIDE,NOUNIT,UPPER
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      parameter(one=1.0d+0,zero=0.0d+0)
*     ..
*
*     Test the input parameters.
*
      lside = lsame(side,'L')
      IF (lside) THEN
          nrowa = m
      ELSE
          nrowa = n
      END IF
      nounit = lsame(diag,'N')
      upper = lsame(uplo,'U')
*
      info = 0
      IF ((.NOT.lside) .AND. (.NOT.lsame(side,'R'))) THEN
          info = 1
      ELSE IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
          info = 2
      ELSE IF ((.NOT.lsame(transa,'N')) .AND.
     +         (.NOT.lsame(transa,'T')) .AND.
     +         (.NOT.lsame(transa,'C'))) THEN
          info = 3
      ELSE IF ((.NOT.lsame(diag,'U')) .AND. (.NOT.lsame(diag,'N'))) THEN
          info = 4
      ELSE IF (m.LT.0) THEN
          info = 5
      ELSE IF (n.LT.0) THEN
          info = 6
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 9
      ELSE IF (ldb.LT.max(1,m)) THEN
          info = 11
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DTRMM ',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF (m.EQ.0 .OR. n.EQ.0) RETURN
*
*     And when  alpha.eq.zero.
*
      IF (alpha.EQ.zero) THEN
          DO 20 j = 1,n
              DO 10 i = 1,m
                  b(i,j) = zero
   10         CONTINUE
   20     CONTINUE
          RETURN
      END IF
*
*     Start the operations.
*
      IF (lside) THEN
          IF (lsame(transa,'N')) THEN
*
*           Form  B := alpha*A*B.
*
              IF (upper) THEN
                  DO 50 j = 1,n
                      DO 40 k = 1,m
                          IF (b(k,j).NE.zero) THEN
                              temp = alpha*b(k,j)
                              DO 30 i = 1,k - 1
                                  b(i,j) = b(i,j) + temp*a(i,k)
   30                         CONTINUE
                              IF (nounit) temp = temp*a(k,k)
                              b(k,j) = temp
                          END IF
   40                 CONTINUE
   50             CONTINUE
              ELSE
                  DO 80 j = 1,n
                      DO 70 k = m,1,-1
                          IF (b(k,j).NE.zero) THEN
                              temp = alpha*b(k,j)
                              b(k,j) = temp
                              IF (nounit) b(k,j) = b(k,j)*a(k,k)
                              DO 60 i = k + 1,m
                                  b(i,j) = b(i,j) + temp*a(i,k)
   60                         CONTINUE
                          END IF
   70                 CONTINUE
   80             CONTINUE
              END IF
          ELSE
*
*           Form  B := alpha*A**T*B.
*
              IF (upper) THEN
                  DO 110 j = 1,n
                      DO 100 i = m,1,-1
                          temp = b(i,j)
                          IF (nounit) temp = temp*a(i,i)
                          DO 90 k = 1,i - 1
                              temp = temp + a(k,i)*b(k,j)
   90                     CONTINUE
                          b(i,j) = alpha*temp
  100                 CONTINUE
  110             CONTINUE
              ELSE
                  DO 140 j = 1,n
                      DO 130 i = 1,m
                          temp = b(i,j)
                          IF (nounit) temp = temp*a(i,i)
                          DO 120 k = i + 1,m
                              temp = temp + a(k,i)*b(k,j)
  120                     CONTINUE
                          b(i,j) = alpha*temp
  130                 CONTINUE
  140             CONTINUE
              END IF
          END IF
      ELSE
          IF (lsame(transa,'N')) THEN
*
*           Form  B := alpha*B*A.
*
              IF (upper) THEN
                  DO 180 j = n,1,-1
                      temp = alpha
                      IF (nounit) temp = temp*a(j,j)
                      DO 150 i = 1,m
                          b(i,j) = temp*b(i,j)
  150                 CONTINUE
                      DO 170 k = 1,j - 1
                          IF (a(k,j).NE.zero) THEN
                              temp = alpha*a(k,j)
                              DO 160 i = 1,m
                                  b(i,j) = b(i,j) + temp*b(i,k)
  160                         CONTINUE
                          END IF
  170                 CONTINUE
  180             CONTINUE
              ELSE
                  DO 220 j = 1,n
                      temp = alpha
                      IF (nounit) temp = temp*a(j,j)
                      DO 190 i = 1,m
                          b(i,j) = temp*b(i,j)
  190                 CONTINUE
                      DO 210 k = j + 1,n
                          IF (a(k,j).NE.zero) THEN
                              temp = alpha*a(k,j)
                              DO 200 i = 1,m
                                  b(i,j) = b(i,j) + temp*b(i,k)
  200                         CONTINUE
                          END IF
  210                 CONTINUE
  220             CONTINUE
              END IF
          ELSE
*
*           Form  B := alpha*B*A**T.
*
              IF (upper) THEN
                  DO 260 k = 1,n
                      DO 240 j = 1,k - 1
                          IF (a(j,k).NE.zero) THEN
                              temp = alpha*a(j,k)
                              DO 230 i = 1,m
                                  b(i,j) = b(i,j) + temp*b(i,k)
  230                         CONTINUE
                          END IF
  240                 CONTINUE
                      temp = alpha
                      IF (nounit) temp = temp*a(k,k)
                      IF (temp.NE.one) THEN
                          DO 250 i = 1,m
                              b(i,k) = temp*b(i,k)
  250                     CONTINUE
                      END IF
  260             CONTINUE
              ELSE
                  DO 300 k = n,1,-1
                      DO 280 j = k + 1,n
                          IF (a(j,k).NE.zero) THEN
                              temp = alpha*a(j,k)
                              DO 270 i = 1,m
                                  b(i,j) = b(i,j) + temp*b(i,k)
  270                         CONTINUE
                          END IF
  280                 CONTINUE
                      temp = alpha
                      IF (nounit) temp = temp*a(k,k)
                      IF (temp.NE.one) THEN
                          DO 290 i = 1,m
                              b(i,k) = temp*b(i,k)
  290                     CONTINUE
                      END IF
  300             CONTINUE
              END IF
          END IF
      END IF
*
      RETURN
*
*     End of DTRMM
*

dtrsm()

subroutine dtrsm	(	character	side,
		character	uplo,
		character	transa,
		character	diag,
		integer	m,
		integer	n,
		double precision	alpha,
		double precision, dimension(lda,*)	a,
		integer	lda,
		double precision, dimension(ldb,*)	b,
		integer	ldb )

DTRSM

Purpose:

!>
!> DTRSM  solves one of the matrix equations
!>
!>    op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
!>
!> where alpha is a scalar, X and B are m by n matrices, A is a unit, or
!> non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
!>
!>    op( A ) = A   or   op( A ) = A**T.
!>
!> The matrix X is overwritten on B.
!> 

Parameters

[in]	SIDE	!> SIDE is CHARACTER*1 !> On entry, SIDE specifies whether op( A ) appears on the left !> or right of X as follows: !> !> SIDE = 'L' or 'l' op( A )*X = alpha*B. !> !> SIDE = 'R' or 'r' X*op( A ) = alpha*B. !>
[in]	UPLO	!> UPLO is CHARACTER*1 !> On entry, UPLO specifies whether the matrix A is an upper or !> lower triangular matrix as follows: !> !> UPLO = 'U' or 'u' A is an upper triangular matrix. !> !> UPLO = 'L' or 'l' A is a lower triangular matrix. !>
[in]	TRANSA	!> TRANSA is CHARACTER*1 !> On entry, TRANSA specifies the form of op( A ) to be used in !> the matrix multiplication as follows: !> !> TRANSA = 'N' or 'n' op( A ) = A. !> !> TRANSA = 'T' or 't' op( A ) = A**T. !> !> TRANSA = 'C' or 'c' op( A ) = A**T. !>
[in]	DIAG	!> DIAG is CHARACTER*1 !> On entry, DIAG specifies whether or not A is unit triangular !> as follows: !> !> DIAG = 'U' or 'u' A is assumed to be unit triangular. !> !> DIAG = 'N' or 'n' A is not assumed to be unit !> triangular. !>
[in]	M	!> M is INTEGER !> On entry, M specifies the number of rows of B. M must be at !> least zero. !>
[in]	N	!> N is INTEGER !> On entry, N specifies the number of columns of B. N must be !> at least zero. !>
[in]	ALPHA	!> ALPHA is DOUBLE PRECISION. !> On entry, ALPHA specifies the scalar alpha. When alpha is !> zero then A is not referenced and B need not be set before !> entry. !>
[in]	A	!> A is DOUBLE PRECISION array, dimension ( LDA, k ), !> where k is m when SIDE = 'L' or 'l' !> and k is n when SIDE = 'R' or 'r'. !> Before entry with UPLO = 'U' or 'u', the leading k by k !> upper triangular part of the array A must contain the upper !> triangular matrix and the strictly lower triangular part of !> A is not referenced. !> Before entry with UPLO = 'L' or 'l', the leading k by k !> lower triangular part of the array A must contain the lower !> triangular matrix and the strictly upper triangular part of !> A is not referenced. !> Note that when DIAG = 'U' or 'u', the diagonal elements of !> A are not referenced either, but are assumed to be unity. !>
[in]	LDA	!> LDA is INTEGER !> On entry, LDA specifies the first dimension of A as declared !> in the calling (sub) program. When SIDE = 'L' or 'l' then !> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' !> then LDA must be at least max( 1, n ). !>
[in,out]	B	!> B is DOUBLE PRECISION array, dimension ( LDB, N ) !> Before entry, the leading m by n part of the array B must !> contain the right-hand side matrix B, and on exit is !> overwritten by the solution matrix X. !>
[in]	LDB	!> LDB is INTEGER !> On entry, LDB specifies the first dimension of B as declared !> in the calling (sub) program. LDB must be at least !> max( 1, m ). !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  Level 3 Blas routine.
!>
!>
!>  -- Written on 8-February-1989.
!>     Jack Dongarra, Argonne National Laboratory.
!>     Iain Duff, AERE Harwell.
!>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
!>     Sven Hammarling, Numerical Algorithms Group Ltd.
!> 

Definition at line 180 of file dtrsm.f.

*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*),B(LDB,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC max
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER I,INFO,J,K,NROWA
      LOGICAL LSIDE,NOUNIT,UPPER
*     ..
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      parameter(one=1.0d+0,zero=0.0d+0)
*     ..
*
*     Test the input parameters.
*
      lside = lsame(side,'L')
      IF (lside) THEN
          nrowa = m
      ELSE
          nrowa = n
      END IF
      nounit = lsame(diag,'N')
      upper = lsame(uplo,'U')
*
      info = 0
      IF ((.NOT.lside) .AND. (.NOT.lsame(side,'R'))) THEN
          info = 1
      ELSE IF ((.NOT.upper) .AND. (.NOT.lsame(uplo,'L'))) THEN
          info = 2
      ELSE IF ((.NOT.lsame(transa,'N')) .AND.
     +         (.NOT.lsame(transa,'T')) .AND.
     +         (.NOT.lsame(transa,'C'))) THEN
          info = 3
      ELSE IF ((.NOT.lsame(diag,'U')) .AND. (.NOT.lsame(diag,'N'))) THEN
          info = 4
      ELSE IF (m.LT.0) THEN
          info = 5
      ELSE IF (n.LT.0) THEN
          info = 6
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 9
      ELSE IF (ldb.LT.max(1,m)) THEN
          info = 11
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('DTRSM ',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF (m.EQ.0 .OR. n.EQ.0) RETURN
*
*     And when  alpha.eq.zero.
*
      IF (alpha.EQ.zero) THEN
          DO 20 j = 1,n
              DO 10 i = 1,m
                  b(i,j) = zero
   10         CONTINUE
   20     CONTINUE
          RETURN
      END IF
*
*     Start the operations.
*
      IF (lside) THEN
          IF (lsame(transa,'N')) THEN
*
*           Form  B := alpha*inv( A )*B.
*
              IF (upper) THEN
                  DO 60 j = 1,n
                      IF (alpha.NE.one) THEN
                          DO 30 i = 1,m
                              b(i,j) = alpha*b(i,j)
   30                     CONTINUE
                      END IF
                      DO 50 k = m,1,-1
                          IF (b(k,j).NE.zero) THEN
                              IF (nounit) b(k,j) = b(k,j)/a(k,k)
                              DO 40 i = 1,k - 1
                                  b(i,j) = b(i,j) - b(k,j)*a(i,k)
   40                         CONTINUE
                          END IF
   50                 CONTINUE
   60             CONTINUE
              ELSE
                  DO 100 j = 1,n
                      IF (alpha.NE.one) THEN
                          DO 70 i = 1,m
                              b(i,j) = alpha*b(i,j)
   70                     CONTINUE
                      END IF
                      DO 90 k = 1,m
                          IF (b(k,j).NE.zero) THEN
                              IF (nounit) b(k,j) = b(k,j)/a(k,k)
                              DO 80 i = k + 1,m
                                  b(i,j) = b(i,j) - b(k,j)*a(i,k)
   80                         CONTINUE
                          END IF
   90                 CONTINUE
  100             CONTINUE
              END IF
          ELSE
*
*           Form  B := alpha*inv( A**T )*B.
*
              IF (upper) THEN
                  DO 130 j = 1,n
                      DO 120 i = 1,m
                          temp = alpha*b(i,j)
                          DO 110 k = 1,i - 1
                              temp = temp - a(k,i)*b(k,j)
  110                     CONTINUE
                          IF (nounit) temp = temp/a(i,i)
                          b(i,j) = temp
  120                 CONTINUE
  130             CONTINUE
              ELSE
                  DO 160 j = 1,n
                      DO 150 i = m,1,-1
                          temp = alpha*b(i,j)
                          DO 140 k = i + 1,m
                              temp = temp - a(k,i)*b(k,j)
  140                     CONTINUE
                          IF (nounit) temp = temp/a(i,i)
                          b(i,j) = temp
  150                 CONTINUE
  160             CONTINUE
              END IF
          END IF
      ELSE
          IF (lsame(transa,'N')) THEN
*
*           Form  B := alpha*B*inv( A ).
*
              IF (upper) THEN
                  DO 210 j = 1,n
                      IF (alpha.NE.one) THEN
                          DO 170 i = 1,m
                              b(i,j) = alpha*b(i,j)
  170                     CONTINUE
                      END IF
                      DO 190 k = 1,j - 1
                          IF (a(k,j).NE.zero) THEN
                              DO 180 i = 1,m
                                  b(i,j) = b(i,j) - a(k,j)*b(i,k)
  180                         CONTINUE
                          END IF
  190                 CONTINUE
                      IF (nounit) THEN
                          temp = one/a(j,j)
                          DO 200 i = 1,m
                              b(i,j) = temp*b(i,j)
  200                     CONTINUE
                      END IF
  210             CONTINUE
              ELSE
                  DO 260 j = n,1,-1
                      IF (alpha.NE.one) THEN
                          DO 220 i = 1,m
                              b(i,j) = alpha*b(i,j)
  220                     CONTINUE
                      END IF
                      DO 240 k = j + 1,n
                          IF (a(k,j).NE.zero) THEN
                              DO 230 i = 1,m
                                  b(i,j) = b(i,j) - a(k,j)*b(i,k)
  230                         CONTINUE
                          END IF
  240                 CONTINUE
                      IF (nounit) THEN
                          temp = one/a(j,j)
                          DO 250 i = 1,m
                              b(i,j) = temp*b(i,j)
  250                     CONTINUE
                      END IF
  260             CONTINUE
              END IF
          ELSE
*
*           Form  B := alpha*B*inv( A**T ).
*
              IF (upper) THEN
                  DO 310 k = n,1,-1
                      IF (nounit) THEN
                          temp = one/a(k,k)
                          DO 270 i = 1,m
                              b(i,k) = temp*b(i,k)
  270                     CONTINUE
                      END IF
                      DO 290 j = 1,k - 1
                          IF (a(j,k).NE.zero) THEN
                              temp = a(j,k)
                              DO 280 i = 1,m
                                  b(i,j) = b(i,j) - temp*b(i,k)
  280                         CONTINUE
                          END IF
  290                 CONTINUE
                      IF (alpha.NE.one) THEN
                          DO 300 i = 1,m
                              b(i,k) = alpha*b(i,k)
  300                     CONTINUE
                      END IF
  310             CONTINUE
              ELSE
                  DO 360 k = 1,n
                      IF (nounit) THEN
                          temp = one/a(k,k)
                          DO 320 i = 1,m
                              b(i,k) = temp*b(i,k)
  320                     CONTINUE
                      END IF
                      DO 340 j = k + 1,n
                          IF (a(j,k).NE.zero) THEN
                              temp = a(j,k)
                              DO 330 i = 1,m
                                  b(i,j) = b(i,j) - temp*b(i,k)
  330                         CONTINUE
                          END IF
  340                 CONTINUE
                      IF (alpha.NE.one) THEN
                          DO 350 i = 1,m
                              b(i,k) = alpha*b(i,k)
  350                     CONTINUE
                      END IF
  360             CONTINUE
              END IF
          END IF
      END IF
*
      RETURN
*
*     End of DTRSM
*