![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgemm.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 dgemm.f CVS log | Up to [local] / rpl / lapack / blas |
1.1 bertrand 1: SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION ALPHA,BETA
4: INTEGER K,LDA,LDB,LDC,M,N
5: CHARACTER TRANSA,TRANSB
6: * ..
7: * .. Array Arguments ..
8: DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DGEMM performs one of the matrix-matrix operations
15: *
16: * C := alpha*op( A )*op( B ) + beta*C,
17: *
18: * where op( X ) is one of
19: *
20: * op( X ) = X or op( X ) = X',
21: *
22: * alpha and beta are scalars, and A, B and C are matrices, with op( A )
23: * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
24: *
25: * Arguments
26: * ==========
27: *
28: * TRANSA - CHARACTER*1.
29: * On entry, TRANSA specifies the form of op( A ) to be used in
30: * the matrix multiplication as follows:
31: *
32: * TRANSA = 'N' or 'n', op( A ) = A.
33: *
34: * TRANSA = 'T' or 't', op( A ) = A'.
35: *
36: * TRANSA = 'C' or 'c', op( A ) = A'.
37: *
38: * Unchanged on exit.
39: *
40: * TRANSB - CHARACTER*1.
41: * On entry, TRANSB specifies the form of op( B ) to be used in
42: * the matrix multiplication as follows:
43: *
44: * TRANSB = 'N' or 'n', op( B ) = B.
45: *
46: * TRANSB = 'T' or 't', op( B ) = B'.
47: *
48: * TRANSB = 'C' or 'c', op( B ) = B'.
49: *
50: * Unchanged on exit.
51: *
52: * M - INTEGER.
53: * On entry, M specifies the number of rows of the matrix
54: * op( A ) and of the matrix C. M must be at least zero.
55: * Unchanged on exit.
56: *
57: * N - INTEGER.
58: * On entry, N specifies the number of columns of the matrix
59: * op( B ) and the number of columns of the matrix C. N must be
60: * at least zero.
61: * Unchanged on exit.
62: *
63: * K - INTEGER.
64: * On entry, K specifies the number of columns of the matrix
65: * op( A ) and the number of rows of the matrix op( B ). K must
66: * be at least zero.
67: * Unchanged on exit.
68: *
69: * ALPHA - DOUBLE PRECISION.
70: * On entry, ALPHA specifies the scalar alpha.
71: * Unchanged on exit.
72: *
73: * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
74: * k when TRANSA = 'N' or 'n', and is m otherwise.
75: * Before entry with TRANSA = 'N' or 'n', the leading m by k
76: * part of the array A must contain the matrix A, otherwise
77: * the leading k by m part of the array A must contain the
78: * matrix A.
79: * Unchanged on exit.
80: *
81: * LDA - INTEGER.
82: * On entry, LDA specifies the first dimension of A as declared
83: * in the calling (sub) program. When TRANSA = 'N' or 'n' then
84: * LDA must be at least max( 1, m ), otherwise LDA must be at
85: * least max( 1, k ).
86: * Unchanged on exit.
87: *
88: * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
89: * n when TRANSB = 'N' or 'n', and is k otherwise.
90: * Before entry with TRANSB = 'N' or 'n', the leading k by n
91: * part of the array B must contain the matrix B, otherwise
92: * the leading n by k part of the array B must contain the
93: * matrix B.
94: * Unchanged on exit.
95: *
96: * LDB - INTEGER.
97: * On entry, LDB specifies the first dimension of B as declared
98: * in the calling (sub) program. When TRANSB = 'N' or 'n' then
99: * LDB must be at least max( 1, k ), otherwise LDB must be at
100: * least max( 1, n ).
101: * Unchanged on exit.
102: *
103: * BETA - DOUBLE PRECISION.
104: * On entry, BETA specifies the scalar beta. When BETA is
105: * supplied as zero then C need not be set on input.
106: * Unchanged on exit.
107: *
108: * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
109: * Before entry, the leading m by n part of the array C must
110: * contain the matrix C, except when beta is zero, in which
111: * case C need not be set on entry.
112: * On exit, the array C is overwritten by the m by n matrix
113: * ( alpha*op( A )*op( B ) + beta*C ).
114: *
115: * LDC - INTEGER.
116: * On entry, LDC specifies the first dimension of C as declared
117: * in the calling (sub) program. LDC must be at least
118: * max( 1, m ).
119: * Unchanged on exit.
120: *
121: * Further Details
122: * ===============
123: *
124: * Level 3 Blas routine.
125: *
126: * -- Written on 8-February-1989.
127: * Jack Dongarra, Argonne National Laboratory.
128: * Iain Duff, AERE Harwell.
129: * Jeremy Du Croz, Numerical Algorithms Group Ltd.
130: * Sven Hammarling, Numerical Algorithms Group Ltd.
131: *
132: * =====================================================================
133: *
134: * .. External Functions ..
135: LOGICAL LSAME
136: EXTERNAL LSAME
137: * ..
138: * .. External Subroutines ..
139: EXTERNAL XERBLA
140: * ..
141: * .. Intrinsic Functions ..
142: INTRINSIC MAX
143: * ..
144: * .. Local Scalars ..
145: DOUBLE PRECISION TEMP
146: INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
147: LOGICAL NOTA,NOTB
148: * ..
149: * .. Parameters ..
150: DOUBLE PRECISION ONE,ZERO
151: PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
152: * ..
153: *
154: * Set NOTA and NOTB as true if A and B respectively are not
155: * transposed and set NROWA, NCOLA and NROWB as the number of rows
156: * and columns of A and the number of rows of B respectively.
157: *
158: NOTA = LSAME(TRANSA,'N')
159: NOTB = LSAME(TRANSB,'N')
160: IF (NOTA) THEN
161: NROWA = M
162: NCOLA = K
163: ELSE
164: NROWA = K
165: NCOLA = M
166: END IF
167: IF (NOTB) THEN
168: NROWB = K
169: ELSE
170: NROWB = N
171: END IF
172: *
173: * Test the input parameters.
174: *
175: INFO = 0
176: IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
177: + (.NOT.LSAME(TRANSA,'T'))) THEN
178: INFO = 1
179: ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
180: + (.NOT.LSAME(TRANSB,'T'))) THEN
181: INFO = 2
182: ELSE IF (M.LT.0) THEN
183: INFO = 3
184: ELSE IF (N.LT.0) THEN
185: INFO = 4
186: ELSE IF (K.LT.0) THEN
187: INFO = 5
188: ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
189: INFO = 8
190: ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
191: INFO = 10
192: ELSE IF (LDC.LT.MAX(1,M)) THEN
193: INFO = 13
194: END IF
195: IF (INFO.NE.0) THEN
196: CALL XERBLA('DGEMM ',INFO)
197: RETURN
198: END IF
199: *
200: * Quick return if possible.
201: *
202: IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
203: + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
204: *
205: * And if alpha.eq.zero.
206: *
207: IF (ALPHA.EQ.ZERO) THEN
208: IF (BETA.EQ.ZERO) THEN
209: DO 20 J = 1,N
210: DO 10 I = 1,M
211: C(I,J) = ZERO
212: 10 CONTINUE
213: 20 CONTINUE
214: ELSE
215: DO 40 J = 1,N
216: DO 30 I = 1,M
217: C(I,J) = BETA*C(I,J)
218: 30 CONTINUE
219: 40 CONTINUE
220: END IF
221: RETURN
222: END IF
223: *
224: * Start the operations.
225: *
226: IF (NOTB) THEN
227: IF (NOTA) THEN
228: *
229: * Form C := alpha*A*B + beta*C.
230: *
231: DO 90 J = 1,N
232: IF (BETA.EQ.ZERO) THEN
233: DO 50 I = 1,M
234: C(I,J) = ZERO
235: 50 CONTINUE
236: ELSE IF (BETA.NE.ONE) THEN
237: DO 60 I = 1,M
238: C(I,J) = BETA*C(I,J)
239: 60 CONTINUE
240: END IF
241: DO 80 L = 1,K
242: IF (B(L,J).NE.ZERO) THEN
243: TEMP = ALPHA*B(L,J)
244: DO 70 I = 1,M
245: C(I,J) = C(I,J) + TEMP*A(I,L)
246: 70 CONTINUE
247: END IF
248: 80 CONTINUE
249: 90 CONTINUE
250: ELSE
251: *
252: * Form C := alpha*A'*B + beta*C
253: *
254: DO 120 J = 1,N
255: DO 110 I = 1,M
256: TEMP = ZERO
257: DO 100 L = 1,K
258: TEMP = TEMP + A(L,I)*B(L,J)
259: 100 CONTINUE
260: IF (BETA.EQ.ZERO) THEN
261: C(I,J) = ALPHA*TEMP
262: ELSE
263: C(I,J) = ALPHA*TEMP + BETA*C(I,J)
264: END IF
265: 110 CONTINUE
266: 120 CONTINUE
267: END IF
268: ELSE
269: IF (NOTA) THEN
270: *
271: * Form C := alpha*A*B' + beta*C
272: *
273: DO 170 J = 1,N
274: IF (BETA.EQ.ZERO) THEN
275: DO 130 I = 1,M
276: C(I,J) = ZERO
277: 130 CONTINUE
278: ELSE IF (BETA.NE.ONE) THEN
279: DO 140 I = 1,M
280: C(I,J) = BETA*C(I,J)
281: 140 CONTINUE
282: END IF
283: DO 160 L = 1,K
284: IF (B(J,L).NE.ZERO) THEN
285: TEMP = ALPHA*B(J,L)
286: DO 150 I = 1,M
287: C(I,J) = C(I,J) + TEMP*A(I,L)
288: 150 CONTINUE
289: END IF
290: 160 CONTINUE
291: 170 CONTINUE
292: ELSE
293: *
294: * Form C := alpha*A'*B' + beta*C
295: *
296: DO 200 J = 1,N
297: DO 190 I = 1,M
298: TEMP = ZERO
299: DO 180 L = 1,K
300: TEMP = TEMP + A(L,I)*B(J,L)
301: 180 CONTINUE
302: IF (BETA.EQ.ZERO) THEN
303: C(I,J) = ALPHA*TEMP
304: ELSE
305: C(I,J) = ALPHA*TEMP + BETA*C(I,J)
306: END IF
307: 190 CONTINUE
308: 200 CONTINUE
309: END IF
310: END IF
311: *
312: RETURN
313: *
314: * End of DGEMM .
315: *
316: END