SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) **************************************************************************** * * * DATA PARALLEL BLAS based on MPL * * * * Version 1.0 16/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 DOUBLE PRECISION ALPHA * .. Array Arguments .. double precision, array(:,:) :: a,b intent(in) :: a intent(inout) :: b * .. * * 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'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - 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. * * 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 right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * 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. * * * .. Local Arrays .. double precision, array(m,n) :: aloc, bloc logical, array(m,n) :: mask1, mask2 * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX INTRINSIC spread INTRINSIC eoshift * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA INTEGER m1, p1, ip DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * 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( N.EQ.0 ) $ RETURN m1 = -1 p1 = 1 * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN B = ZERO RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * if (upper) then do i = m, 1, -1 ip = i - 1 if (nounit) b(i,:) = b(i,:)/a(i,i) b(1:ip,1:n) = b(1:ip,1:n) & - spread(a(1:ip,i),2,n)*spread(b(i,1:n),1,i) enddo else do i = 1 , m ip = i + 1 if (nounit) b(i,:) = b(i,:)/a(i,i) b(ip:m,1:n) = b(ip:m,1:n) & - spread(a(ip:m,i),2,n)*spread(b(i,1:n),1,m-i) enddo endif else * * Form B := alpha*inv( A' )*B. * if (upper) then do i = 1 , m ip = i + 1 if (nounit) b(i,:) = b(i,:)/a(i,i) b(ip:m,1:n) = b(ip:m,1:n) & - spread(a(i,ip:m),2,n)*spread(b(i,1:n),1,m-i) enddo else do i = m, 1, -1 ip = i - 1 if (nounit) b(i,:) = b(i,:)/a(i,i) b(1:ip,1:n) = b(1:ip,1:n) & - spread(a(i,1:ip),2,n)*spread(b(i,1:n),1,i) enddo endif endif ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * if (upper) then do i = 1 , n ip = i + 1 if (nounit) b(:,i) = b(:,i)/a(i,i) b(1:m,ip:n) = b(1:m,ip:n) & - spread(a(i,ip:n),1,m)*spread(b(1:m,i),2,n-i) enddo else do i = n, 1, -1 ip = i - 1 if (nounit) b(:,i) = b(:,i)/a(i,i) b(1:m,1:ip) = b(1:m,1:ip) & - spread(a(i,1:ip),1,m)*spread(b(1:m,i),2,i) enddo endif ELSE * * Form B := alpha*B*inv( A' ). * if (upper) then do i = n, 1, -1 ip = i - 1 if (nounit) b(:,i) = b(:,i)/a(i,i) b(1:m,1:ip) = b(1:m,1:ip) & - spread(a(1:ip,i),1,m)*spread(b(1:m,i),2,i) enddo else do i = 1 , n ip = i + 1 if (nounit) b(:,i) = b(:,i)/a(i,i) b(1:m,ip:n) = b(1:m,ip:n) & - spread(a(ip:n,i),1,m)*spread(b(1:m,i),2,n-i) enddo endif ENDIF ENDIF c IF( alpha .ne. ONE ) b = alpha * b * RETURN * * End of DTRSM . * END