001:       SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
002: *     .. Scalar Arguments ..
003:       DOUBLE PRECISION ALPHA
004:       INTEGER LDA,LDB,M,N
005:       CHARACTER DIAG,SIDE,TRANSA,UPLO
006: *     ..
007: *     .. Array Arguments ..
008:       DOUBLE PRECISION A(LDA,*),B(LDB,*)
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  DTRMM  performs one of the matrix-matrix operations
015: *
016: *     B := alpha*op( A )*B,   or   B := alpha*B*op( A ),
017: *
018: *  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
019: *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
020: *
021: *     op( A ) = A   or   op( A ) = A'.
022: *
023: *  Arguments
024: *  ==========
025: *
026: *  SIDE   - CHARACTER*1.
027: *           On entry,  SIDE specifies whether  op( A ) multiplies B from
028: *           the left or right as follows:
029: *
030: *              SIDE = 'L' or 'l'   B := alpha*op( A )*B.
031: *
032: *              SIDE = 'R' or 'r'   B := alpha*B*op( A ).
033: *
034: *           Unchanged on exit.
035: *
036: *  UPLO   - CHARACTER*1.
037: *           On entry, UPLO specifies whether the matrix A is an upper or
038: *           lower triangular matrix as follows:
039: *
040: *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
041: *
042: *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
043: *
044: *           Unchanged on exit.
045: *
046: *  TRANSA - CHARACTER*1.
047: *           On entry, TRANSA specifies the form of op( A ) to be used in
048: *           the matrix multiplication as follows:
049: *
050: *              TRANSA = 'N' or 'n'   op( A ) = A.
051: *
052: *              TRANSA = 'T' or 't'   op( A ) = A'.
053: *
054: *              TRANSA = 'C' or 'c'   op( A ) = A'.
055: *
056: *           Unchanged on exit.
057: *
058: *  DIAG   - CHARACTER*1.
059: *           On entry, DIAG specifies whether or not A is unit triangular
060: *           as follows:
061: *
062: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
063: *
064: *              DIAG = 'N' or 'n'   A is not assumed to be unit
065: *                                  triangular.
066: *
067: *           Unchanged on exit.
068: *
069: *  M      - INTEGER.
070: *           On entry, M specifies the number of rows of B. M must be at
071: *           least zero.
072: *           Unchanged on exit.
073: *
074: *  N      - INTEGER.
075: *           On entry, N specifies the number of columns of B.  N must be
076: *           at least zero.
077: *           Unchanged on exit.
078: *
079: *  ALPHA  - DOUBLE PRECISION.
080: *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
081: *           zero then  A is not referenced and  B need not be set before
082: *           entry.
083: *           Unchanged on exit.
084: *
085: *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
086: *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
087: *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
088: *           upper triangular part of the array  A must contain the upper
089: *           triangular matrix  and the strictly lower triangular part of
090: *           A is not referenced.
091: *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
092: *           lower triangular part of the array  A must contain the lower
093: *           triangular matrix  and the strictly upper triangular part of
094: *           A is not referenced.
095: *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
096: *           A  are not referenced either,  but are assumed to be  unity.
097: *           Unchanged on exit.
098: *
099: *  LDA    - INTEGER.
100: *           On entry, LDA specifies the first dimension of A as declared
101: *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
102: *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
103: *           then LDA must be at least max( 1, n ).
104: *           Unchanged on exit.
105: *
106: *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
107: *           Before entry,  the leading  m by n part of the array  B must
108: *           contain the matrix  B,  and  on exit  is overwritten  by the
109: *           transformed matrix.
110: *
111: *  LDB    - INTEGER.
112: *           On entry, LDB specifies the first dimension of B as declared
113: *           in  the  calling  (sub)  program.   LDB  must  be  at  least
114: *           max( 1, m ).
115: *           Unchanged on exit.
116: *
117: *  Further Details
118: *  ===============
119: *
120: *  Level 3 Blas routine.
121: *
122: *  -- Written on 8-February-1989.
123: *     Jack Dongarra, Argonne National Laboratory.
124: *     Iain Duff, AERE Harwell.
125: *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
126: *     Sven Hammarling, Numerical Algorithms Group Ltd.
127: *
128: *  =====================================================================
129: *
130: *     .. External Functions ..
131:       LOGICAL LSAME
132:       EXTERNAL LSAME
133: *     ..
134: *     .. External Subroutines ..
135:       EXTERNAL XERBLA
136: *     ..
137: *     .. Intrinsic Functions ..
138:       INTRINSIC MAX
139: *     ..
140: *     .. Local Scalars ..
141:       DOUBLE PRECISION TEMP
142:       INTEGER I,INFO,J,K,NROWA
143:       LOGICAL LSIDE,NOUNIT,UPPER
144: *     ..
145: *     .. Parameters ..
146:       DOUBLE PRECISION ONE,ZERO
147:       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
148: *     ..
149: *
150: *     Test the input parameters.
151: *
152:       LSIDE = LSAME(SIDE,'L')
153:       IF (LSIDE) THEN
154:           NROWA = M
155:       ELSE
156:           NROWA = N
157:       END IF
158:       NOUNIT = LSAME(DIAG,'N')
159:       UPPER = LSAME(UPLO,'U')
160: *
161:       INFO = 0
162:       IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
163:           INFO = 1
164:       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
165:           INFO = 2
166:       ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
167:      +         (.NOT.LSAME(TRANSA,'T')) .AND.
168:      +         (.NOT.LSAME(TRANSA,'C'))) THEN
169:           INFO = 3
170:       ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
171:           INFO = 4
172:       ELSE IF (M.LT.0) THEN
173:           INFO = 5
174:       ELSE IF (N.LT.0) THEN
175:           INFO = 6
176:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
177:           INFO = 9
178:       ELSE IF (LDB.LT.MAX(1,M)) THEN
179:           INFO = 11
180:       END IF
181:       IF (INFO.NE.0) THEN
182:           CALL XERBLA('DTRMM ',INFO)
183:           RETURN
184:       END IF
185: *
186: *     Quick return if possible.
187: *
188:       IF (M.EQ.0 .OR. N.EQ.0) RETURN
189: *
190: *     And when  alpha.eq.zero.
191: *
192:       IF (ALPHA.EQ.ZERO) THEN
193:           DO 20 J = 1,N
194:               DO 10 I = 1,M
195:                   B(I,J) = ZERO
196:    10         CONTINUE
197:    20     CONTINUE
198:           RETURN
199:       END IF
200: *
201: *     Start the operations.
202: *
203:       IF (LSIDE) THEN
204:           IF (LSAME(TRANSA,'N')) THEN
205: *
206: *           Form  B := alpha*A*B.
207: *
208:               IF (UPPER) THEN
209:                   DO 50 J = 1,N
210:                       DO 40 K = 1,M
211:                           IF (B(K,J).NE.ZERO) THEN
212:                               TEMP = ALPHA*B(K,J)
213:                               DO 30 I = 1,K - 1
214:                                   B(I,J) = B(I,J) + TEMP*A(I,K)
215:    30                         CONTINUE
216:                               IF (NOUNIT) TEMP = TEMP*A(K,K)
217:                               B(K,J) = TEMP
218:                           END IF
219:    40                 CONTINUE
220:    50             CONTINUE
221:               ELSE
222:                   DO 80 J = 1,N
223:                       DO 70 K = M,1,-1
224:                           IF (B(K,J).NE.ZERO) THEN
225:                               TEMP = ALPHA*B(K,J)
226:                               B(K,J) = TEMP
227:                               IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
228:                               DO 60 I = K + 1,M
229:                                   B(I,J) = B(I,J) + TEMP*A(I,K)
230:    60                         CONTINUE
231:                           END IF
232:    70                 CONTINUE
233:    80             CONTINUE
234:               END IF
235:           ELSE
236: *
237: *           Form  B := alpha*A'*B.
238: *
239:               IF (UPPER) THEN
240:                   DO 110 J = 1,N
241:                       DO 100 I = M,1,-1
242:                           TEMP = B(I,J)
243:                           IF (NOUNIT) TEMP = TEMP*A(I,I)
244:                           DO 90 K = 1,I - 1
245:                               TEMP = TEMP + A(K,I)*B(K,J)
246:    90                     CONTINUE
247:                           B(I,J) = ALPHA*TEMP
248:   100                 CONTINUE
249:   110             CONTINUE
250:               ELSE
251:                   DO 140 J = 1,N
252:                       DO 130 I = 1,M
253:                           TEMP = B(I,J)
254:                           IF (NOUNIT) TEMP = TEMP*A(I,I)
255:                           DO 120 K = I + 1,M
256:                               TEMP = TEMP + A(K,I)*B(K,J)
257:   120                     CONTINUE
258:                           B(I,J) = ALPHA*TEMP
259:   130                 CONTINUE
260:   140             CONTINUE
261:               END IF
262:           END IF
263:       ELSE
264:           IF (LSAME(TRANSA,'N')) THEN
265: *
266: *           Form  B := alpha*B*A.
267: *
268:               IF (UPPER) THEN
269:                   DO 180 J = N,1,-1
270:                       TEMP = ALPHA
271:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
272:                       DO 150 I = 1,M
273:                           B(I,J) = TEMP*B(I,J)
274:   150                 CONTINUE
275:                       DO 170 K = 1,J - 1
276:                           IF (A(K,J).NE.ZERO) THEN
277:                               TEMP = ALPHA*A(K,J)
278:                               DO 160 I = 1,M
279:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
280:   160                         CONTINUE
281:                           END IF
282:   170                 CONTINUE
283:   180             CONTINUE
284:               ELSE
285:                   DO 220 J = 1,N
286:                       TEMP = ALPHA
287:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
288:                       DO 190 I = 1,M
289:                           B(I,J) = TEMP*B(I,J)
290:   190                 CONTINUE
291:                       DO 210 K = J + 1,N
292:                           IF (A(K,J).NE.ZERO) THEN
293:                               TEMP = ALPHA*A(K,J)
294:                               DO 200 I = 1,M
295:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
296:   200                         CONTINUE
297:                           END IF
298:   210                 CONTINUE
299:   220             CONTINUE
300:               END IF
301:           ELSE
302: *
303: *           Form  B := alpha*B*A'.
304: *
305:               IF (UPPER) THEN
306:                   DO 260 K = 1,N
307:                       DO 240 J = 1,K - 1
308:                           IF (A(J,K).NE.ZERO) THEN
309:                               TEMP = ALPHA*A(J,K)
310:                               DO 230 I = 1,M
311:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
312:   230                         CONTINUE
313:                           END IF
314:   240                 CONTINUE
315:                       TEMP = ALPHA
316:                       IF (NOUNIT) TEMP = TEMP*A(K,K)
317:                       IF (TEMP.NE.ONE) THEN
318:                           DO 250 I = 1,M
319:                               B(I,K) = TEMP*B(I,K)
320:   250                     CONTINUE
321:                       END IF
322:   260             CONTINUE
323:               ELSE
324:                   DO 300 K = N,1,-1
325:                       DO 280 J = K + 1,N
326:                           IF (A(J,K).NE.ZERO) THEN
327:                               TEMP = ALPHA*A(J,K)
328:                               DO 270 I = 1,M
329:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
330:   270                         CONTINUE
331:                           END IF
332:   280                 CONTINUE
333:                       TEMP = ALPHA
334:                       IF (NOUNIT) TEMP = TEMP*A(K,K)
335:                       IF (TEMP.NE.ONE) THEN
336:                           DO 290 I = 1,M
337:                               B(I,K) = TEMP*B(I,K)
338:   290                     CONTINUE
339:                       END IF
340:   300             CONTINUE
341:               END IF
342:           END IF
343:       END IF
344: *
345:       RETURN
346: *
347: *     End of DTRMM .
348: *
349:       END
350: