SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) **************************************************************************** * * * DATA PARALLEL BLAS based on MPL * * * * Version 1.0 1/9-92 , * * For MasPar MP-1 computers * * * * para//ab, University of Bergen, NORWAY * * * * These programs must be called using F90 style array syntax. * * Note that the F77 style calling sequence has been retained * * in this version for compatibility reasons, be aware that * * parameters related to the array dimensions and shape therefore may * * be redundant and without any influence. * * The calling sequence may be changed in a future version. * * Please report any BUGs, ideas for improvement or other * * comments to * * adm@parallab.uib.no * * * * Future versions may then reflect your suggestions. * * The most current version of this software is available * * from netlib@nac.no , send the message `send index from maspar' * * * * REVISIONS: * * * **************************************************************************** implicit none * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB REAL ALPHA * .. Array Arguments .. double precision, array(:,:) :: a,b intent(in) :: a intent(inout) :: b * .. * * 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'. * * Parameters * ========== * * SIDE - 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 ). * * Unchanged on exit. * * UPLO - 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. * * Unchanged on exit. * * TRANSA - 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'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - 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. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - 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. * Unchanged on exit. * * A - DOUBLE PRECISION array of 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. * Unchanged on exit. * * LDA - 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 ). * Unchanged on exit. * * B - DOUBLE PRECISION array of 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. * * LDB - 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 ). * Unchanged on exit. * * * 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. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL DGEMM * .. Intrinsic Functions INTRINSIC MAX INTRINSIC matmul INTRINSIC transpose * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Local Arrays .. integer, array(max(m,n),max(m,n)) :: maska double precision, array(max(m,n),max(m,n)) :: aloc double precision, array(m,n) :: btmp * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * cts cmpf mpl mpl_itrmsk cmpf ondpu a,b * 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( N.EQ.0 ) RETURN IF( m.EQ.0 ) RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN B (1:m,1:n) = ZERO RETURN END IF * * Start the operations. * IF( LSIDE )THEN * * Form B := alpha*op(A)*B. * forall (i = 1:m, j = 1:m) maska(i,j) = j-i cts CALL MPL_ITRMSK(maska(1:m,1:m),m) c aloc(1:m,1:m) = zero c if (upper) then where ( maska(1:m,1:m) .ge. 0 ) aloc(1:m,1:m) = a(1:m,1:m) else where ( maska(1:m,1:m) .le. 0 ) aloc(1:m,1:m) = a(1:m,1:m) endif c if (.not. nounit) then where ( maska(1:m,1:m) .eq. 0 ) aloc(1:m,1:m) = one endif c IF( LSAME( TRANSA, 'N' ) )THEN * b(1:m,1:n)=alpha*matmul(aloc(1:m,1:m),b(1:m,1:n)) cts CALL DGEMM('N','N',m,n,m,alpha,aloc(1:m,1:m),m cts $ ,b(1:m,1:n),m,ZERO,btmp(1:m,1:n),m) CALL DGEMM('N','N',m,n,m,alpha,aloc,m $ ,b,m,ZERO,btmp,m) b = btmp ELSE * b(1:m,1:n)=alpha*matmul(transpose(aloc(1:m,1:m)),b(1:m,1:n)) cts CALL DGEMM('T','N',m,n,m,alpha,aloc(1:m,1:m),m cts $ ,b(1:m,1:n),m,ZERO,btmp(1:m,1:n),m) CALL DGEMM('T','N',m,n,m,alpha,aloc,m $ ,b,m,ZERO,btmp,m) b = btmp ENDIF ELSE * * Form B := alpha*B*op(A). * forall (i = 1:n, j = 1:n) maska(i,j) = j-i cts CALL MPL_ITRMSK(maska(1:n,1:n),n) c aloc(1:n,1:n) = zero c if (upper) then where ( maska(1:n,1:n) .ge. 0 ) aloc(1:n,1:n) = a(1:n,1:n) else where ( maska(1:n,1:n) .le. 0 ) aloc(1:n,1:n) = a(1:n,1:n) endif c if (.not. nounit) then where ( maska(1:n,1:n) .eq. 0 ) aloc(1:n,1:n) = one endif c IF( LSAME( TRANSA, 'N' ) )THEN * b(1:m,1:n)=alpha*matmul(b(1:m,1:n),aloc(1:n,1:n)) cts CALL DGEMM('N','N',m,n,n,alpha,b(1:m,1:n),m cts $ ,aloc(1:n,1:n),n,ZERO,btmp(1:m,1:n),m) CALL DGEMM('N','N',m,n,n,alpha,b,m $ ,aloc,n,ZERO,btmp,m) b = btmp ELSE * b(1:m,1:n)=alpha*matmul(b(1:m,1:n),transpose(aloc(1:n,1:n))) cts CALL DGEMM('N','T',m,n,n,alpha,b(1:m,1:n),m cts $ ,aloc(1:n,1:n),n,ZERO,btmp(1:m,1:n),m) CALL DGEMM('N','T',m,n,n,alpha,b,m $ ,aloc,n,ZERO,btmp,m) b = btmp ENDIF ENDIF * RETURN * * End of DTRMM . * END