3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
13 * .. Scalar Arguments ..
14 * DOUBLE PRECISION ALPHA
16 * CHARACTER DIAG,SIDE,TRANSA,UPLO
18 * .. Array Arguments ..
19 * DOUBLE PRECISION A(LDA,*),B(LDB,*)
28 *> DTRMM performs one of the matrix-matrix operations
30 *> B := alpha*op( A )*B, or B := alpha*B*op( A ),
32 *> where alpha is a scalar, B is an m by n matrix, A is a unit, or
33 *> non-unit, upper or lower triangular matrix and op( A ) is one of
35 *> op( A ) = A or op( A ) = A**T.
43 *> SIDE is CHARACTER*1
44 *> On entry, SIDE specifies whether op( A ) multiplies B from
45 *> the left or right as follows:
47 *> SIDE = 'L' or 'l' B := alpha*op( A )*B.
49 *> SIDE = 'R' or 'r' B := alpha*B*op( A ).
54 *> UPLO is CHARACTER*1
55 *> On entry, UPLO specifies whether the matrix A is an upper or
56 *> lower triangular matrix as follows:
58 *> UPLO = 'U' or 'u' A is an upper triangular matrix.
60 *> UPLO = 'L' or 'l' A is a lower triangular matrix.
65 *> TRANSA is CHARACTER*1
66 *> On entry, TRANSA specifies the form of op( A ) to be used in
67 *> the matrix multiplication as follows:
69 *> TRANSA = 'N' or 'n' op( A ) = A.
71 *> TRANSA = 'T' or 't' op( A ) = A**T.
73 *> TRANSA = 'C' or 'c' op( A ) = A**T.
78 *> DIAG is CHARACTER*1
79 *> On entry, DIAG specifies whether or not A is unit triangular
82 *> DIAG = 'U' or 'u' A is assumed to be unit triangular.
84 *> DIAG = 'N' or 'n' A is not assumed to be unit
91 *> On entry, M specifies the number of rows of B. M must be at
98 *> On entry, N specifies the number of columns of B. N must be
104 *> ALPHA is DOUBLE PRECISION.
105 *> On entry, ALPHA specifies the scalar alpha. When alpha is
106 *> zero then A is not referenced and B need not be set before
112 *> A is DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
113 *> when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
114 *> Before entry with UPLO = 'U' or 'u', the leading k by k
115 *> upper triangular part of the array A must contain the upper
116 *> triangular matrix and the strictly lower triangular part of
117 *> A is not referenced.
118 *> Before entry with UPLO = 'L' or 'l', the leading k by k
119 *> lower triangular part of the array A must contain the lower
120 *> triangular matrix and the strictly upper triangular part of
121 *> A is not referenced.
122 *> Note that when DIAG = 'U' or 'u', the diagonal elements of
123 *> A are not referenced either, but are assumed to be unity.
129 *> On entry, LDA specifies the first dimension of A as declared
130 *> in the calling (sub) program. When SIDE = 'L' or 'l' then
131 *> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
132 *> then LDA must be at least max( 1, n ).
137 *> B is DOUBLE PRECISION array of DIMENSION ( LDB, n ).
138 *> Before entry, the leading m by n part of the array B must
139 *> contain the matrix B, and on exit is overwritten by the
140 *> transformed matrix.
146 *> On entry, LDB specifies the first dimension of B as declared
147 *> in the calling (sub) program. LDB must be at least
154 *> \author Univ. of Tennessee
155 *> \author Univ. of California Berkeley
156 *> \author Univ. of Colorado Denver
159 *> \date November 2011
161 *> \ingroup double_blas_level3
163 *> \par Further Details:
164 * =====================
168 *> Level 3 Blas routine.
170 *> -- Written on 8-February-1989.
171 *> Jack Dongarra, Argonne National Laboratory.
172 *> Iain Duff, AERE Harwell.
173 *> Jeremy Du Croz, Numerical Algorithms Group Ltd.
174 *> Sven Hammarling, Numerical Algorithms Group Ltd.
177 * =====================================================================
178 SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
180 * -- Reference BLAS level3 routine (version 3.4.0) --
181 * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185 * .. Scalar Arguments ..
186 DOUBLE PRECISION ALPHA
188 CHARACTER DIAG,SIDE,TRANSA,UPLO
190 * .. Array Arguments ..
191 DOUBLE PRECISION A(LDA,*),B(LDB,*)
194 * =====================================================================
196 * .. External Functions ..
200 * .. External Subroutines ..
203 * .. Intrinsic Functions ..
206 * .. Local Scalars ..
207 DOUBLE PRECISION TEMP
208 INTEGER I,INFO,J,K,NROWA
209 LOGICAL LSIDE,NOUNIT,UPPER
212 DOUBLE PRECISION ONE,ZERO
213 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
216 * Test the input parameters.
218 LSIDE = LSAME(SIDE,'L')
224 NOUNIT = LSAME(DIAG,'N')
225 UPPER = LSAME(UPLO,'U')
228 IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
230 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
232 ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
233 + (.NOT.LSAME(TRANSA,'T')) .AND.
234 + (.NOT.LSAME(TRANSA,'C'))) THEN
236 ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
238 ELSE IF (M.LT.0) THEN
240 ELSE IF (N.LT.0) THEN
242 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
244 ELSE IF (LDB.LT.MAX(1,M)) THEN
248 CALL XERBLA('DTRMM ',INFO)
252 * Quick return if possible.
254 IF (M.EQ.0 .OR. N.EQ.0) RETURN
256 * And when alpha.eq.zero.
258 IF (ALPHA.EQ.ZERO) THEN
267 * Start the operations.
270 IF (LSAME(TRANSA,'N')) THEN
272 * Form B := alpha*A*B.
277 IF (B(K,J).NE.ZERO) THEN
280 B(I,J) = B(I,J) + TEMP*A(I,K)
282 IF (NOUNIT) TEMP = TEMP*A(K,K)
290 IF (B(K,J).NE.ZERO) THEN
293 IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
295 B(I,J) = B(I,J) + TEMP*A(I,K)
303 * Form B := alpha*A**T*B.
309 IF (NOUNIT) TEMP = TEMP*A(I,I)
311 TEMP = TEMP + A(K,I)*B(K,J)
320 IF (NOUNIT) TEMP = TEMP*A(I,I)
322 TEMP = TEMP + A(K,I)*B(K,J)
330 IF (LSAME(TRANSA,'N')) THEN
332 * Form B := alpha*B*A.
337 IF (NOUNIT) TEMP = TEMP*A(J,J)
342 IF (A(K,J).NE.ZERO) THEN
345 B(I,J) = B(I,J) + TEMP*B(I,K)
353 IF (NOUNIT) TEMP = TEMP*A(J,J)
358 IF (A(K,J).NE.ZERO) THEN
361 B(I,J) = B(I,J) + TEMP*B(I,K)
369 * Form B := alpha*B*A**T.
374 IF (A(J,K).NE.ZERO) THEN
377 B(I,J) = B(I,J) + TEMP*B(I,K)
382 IF (NOUNIT) TEMP = TEMP*A(K,K)
383 IF (TEMP.NE.ONE) THEN
392 IF (A(J,K).NE.ZERO) THEN
395 B(I,J) = B(I,J) + TEMP*B(I,K)
400 IF (NOUNIT) TEMP = TEMP*A(K,K)
401 IF (TEMP.NE.ONE) THEN