![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrmm.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas|
Return to dtrmm.f CVS log | Up to [local] / rpl / lapack / blas |
1.1 bertrand 1: SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION ALPHA
4: INTEGER LDA,LDB,M,N
5: CHARACTER DIAG,SIDE,TRANSA,UPLO
6: * ..
7: * .. Array Arguments ..
8: DOUBLE PRECISION A(LDA,*),B(LDB,*)
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DTRMM performs one of the matrix-matrix operations
15: *
16: * B := alpha*op( A )*B, or B := alpha*B*op( A ),
17: *
18: * where alpha is a scalar, B is an m by n matrix, A is a unit, or
19: * non-unit, upper or lower triangular matrix and op( A ) is one of
20: *
1.7 ! bertrand 21: * op( A ) = A or op( A ) = A**T.
1.1 bertrand 22: *
23: * Arguments
24: * ==========
25: *
26: * SIDE - CHARACTER*1.
27: * On entry, SIDE specifies whether op( A ) multiplies B from
28: * the left or right as follows:
29: *
30: * SIDE = 'L' or 'l' B := alpha*op( A )*B.
31: *
32: * SIDE = 'R' or 'r' B := alpha*B*op( A ).
33: *
34: * Unchanged on exit.
35: *
36: * UPLO - CHARACTER*1.
37: * On entry, UPLO specifies whether the matrix A is an upper or
38: * lower triangular matrix as follows:
39: *
40: * UPLO = 'U' or 'u' A is an upper triangular matrix.
41: *
42: * UPLO = 'L' or 'l' A is a lower triangular matrix.
43: *
44: * Unchanged on exit.
45: *
46: * TRANSA - CHARACTER*1.
47: * On entry, TRANSA specifies the form of op( A ) to be used in
48: * the matrix multiplication as follows:
49: *
50: * TRANSA = 'N' or 'n' op( A ) = A.
51: *
1.7 ! bertrand 52: * TRANSA = 'T' or 't' op( A ) = A**T.
1.1 bertrand 53: *
1.7 ! bertrand 54: * TRANSA = 'C' or 'c' op( A ) = A**T.
1.1 bertrand 55: *
56: * Unchanged on exit.
57: *
58: * DIAG - CHARACTER*1.
59: * On entry, DIAG specifies whether or not A is unit triangular
60: * as follows:
61: *
62: * DIAG = 'U' or 'u' A is assumed to be unit triangular.
63: *
64: * DIAG = 'N' or 'n' A is not assumed to be unit
65: * triangular.
66: *
67: * Unchanged on exit.
68: *
69: * M - INTEGER.
70: * On entry, M specifies the number of rows of B. M must be at
71: * least zero.
72: * Unchanged on exit.
73: *
74: * N - INTEGER.
75: * On entry, N specifies the number of columns of B. N must be
76: * at least zero.
77: * Unchanged on exit.
78: *
79: * ALPHA - DOUBLE PRECISION.
80: * On entry, ALPHA specifies the scalar alpha. When alpha is
81: * zero then A is not referenced and B need not be set before
82: * entry.
83: * Unchanged on exit.
84: *
85: * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
86: * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
87: * Before entry with UPLO = 'U' or 'u', the leading k by k
88: * upper triangular part of the array A must contain the upper
89: * triangular matrix and the strictly lower triangular part of
90: * A is not referenced.
91: * Before entry with UPLO = 'L' or 'l', the leading k by k
92: * lower triangular part of the array A must contain the lower
93: * triangular matrix and the strictly upper triangular part of
94: * A is not referenced.
95: * Note that when DIAG = 'U' or 'u', the diagonal elements of
96: * A are not referenced either, but are assumed to be unity.
97: * Unchanged on exit.
98: *
99: * 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: *
1.7 ! bertrand 237: * Form B := alpha*A**T*B.
1.1 bertrand 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: *
1.7 ! bertrand 303: * Form B := alpha*B*A**T.
1.1 bertrand 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