4347
|
1 SUBROUTINE SGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, |
|
2 $ BETA, C, LDC ) |
|
3 * .. Scalar Arguments .. |
|
4 CHARACTER*1 TRANSA, TRANSB |
|
5 INTEGER M, N, K, LDA, LDB, LDC |
|
6 REAL ALPHA, BETA |
|
7 * .. Array Arguments .. |
|
8 REAL A( LDA, * ), B( LDB, * ), C( LDC, * ) |
|
9 * .. |
|
10 * |
|
11 * Purpose |
|
12 * ======= |
|
13 * |
|
14 * SGEMM 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 * Parameters |
|
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 - REAL . |
|
70 * On entry, ALPHA specifies the scalar alpha. |
|
71 * Unchanged on exit. |
|
72 * |
|
73 * A - REAL 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 - REAL 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 - REAL . |
|
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 - REAL 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 * |
|
122 * Level 3 Blas routine. |
|
123 * |
|
124 * -- Written on 8-February-1989. |
|
125 * Jack Dongarra, Argonne National Laboratory. |
|
126 * Iain Duff, AERE Harwell. |
|
127 * Jeremy Du Croz, Numerical Algorithms Group Ltd. |
|
128 * Sven Hammarling, Numerical Algorithms Group Ltd. |
|
129 * |
|
130 * |
|
131 * .. External Functions .. |
|
132 LOGICAL LSAME |
|
133 EXTERNAL LSAME |
|
134 * .. External Subroutines .. |
|
135 EXTERNAL XERBLA |
|
136 * .. Intrinsic Functions .. |
|
137 INTRINSIC MAX |
|
138 * .. Local Scalars .. |
|
139 LOGICAL NOTA, NOTB |
|
140 INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB |
|
141 REAL TEMP |
|
142 * .. Parameters .. |
|
143 REAL ONE , ZERO |
|
144 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) |
|
145 * .. |
|
146 * .. Executable Statements .. |
|
147 * |
|
148 * Set NOTA and NOTB as true if A and B respectively are not |
|
149 * transposed and set NROWA, NCOLA and NROWB as the number of rows |
|
150 * and columns of A and the number of rows of B respectively. |
|
151 * |
|
152 NOTA = LSAME( TRANSA, 'N' ) |
|
153 NOTB = LSAME( TRANSB, 'N' ) |
|
154 IF( NOTA )THEN |
|
155 NROWA = M |
|
156 NCOLA = K |
|
157 ELSE |
|
158 NROWA = K |
|
159 NCOLA = M |
|
160 END IF |
|
161 IF( NOTB )THEN |
|
162 NROWB = K |
|
163 ELSE |
|
164 NROWB = N |
|
165 END IF |
|
166 * |
|
167 * Test the input parameters. |
|
168 * |
|
169 INFO = 0 |
|
170 IF( ( .NOT.NOTA ).AND. |
|
171 $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. |
|
172 $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN |
|
173 INFO = 1 |
|
174 ELSE IF( ( .NOT.NOTB ).AND. |
|
175 $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. |
|
176 $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN |
|
177 INFO = 2 |
|
178 ELSE IF( M .LT.0 )THEN |
|
179 INFO = 3 |
|
180 ELSE IF( N .LT.0 )THEN |
|
181 INFO = 4 |
|
182 ELSE IF( K .LT.0 )THEN |
|
183 INFO = 5 |
|
184 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN |
|
185 INFO = 8 |
|
186 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN |
|
187 INFO = 10 |
|
188 ELSE IF( LDC.LT.MAX( 1, M ) )THEN |
|
189 INFO = 13 |
|
190 END IF |
|
191 IF( INFO.NE.0 )THEN |
|
192 CALL XERBLA( 'SGEMM ', INFO ) |
|
193 RETURN |
|
194 END IF |
|
195 * |
|
196 * Quick return if possible. |
|
197 * |
|
198 IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. |
|
199 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) |
|
200 $ RETURN |
|
201 * |
|
202 * And if alpha.eq.zero. |
|
203 * |
|
204 IF( ALPHA.EQ.ZERO )THEN |
|
205 IF( BETA.EQ.ZERO )THEN |
|
206 DO 20, J = 1, N |
|
207 DO 10, I = 1, M |
|
208 C( I, J ) = ZERO |
|
209 10 CONTINUE |
|
210 20 CONTINUE |
|
211 ELSE |
|
212 DO 40, J = 1, N |
|
213 DO 30, I = 1, M |
|
214 C( I, J ) = BETA*C( I, J ) |
|
215 30 CONTINUE |
|
216 40 CONTINUE |
|
217 END IF |
|
218 RETURN |
|
219 END IF |
|
220 * |
|
221 * Start the operations. |
|
222 * |
|
223 IF( NOTB )THEN |
|
224 IF( NOTA )THEN |
|
225 * |
|
226 * Form C := alpha*A*B + beta*C. |
|
227 * |
|
228 DO 90, J = 1, N |
|
229 IF( BETA.EQ.ZERO )THEN |
|
230 DO 50, I = 1, M |
|
231 C( I, J ) = ZERO |
|
232 50 CONTINUE |
|
233 ELSE IF( BETA.NE.ONE )THEN |
|
234 DO 60, I = 1, M |
|
235 C( I, J ) = BETA*C( I, J ) |
|
236 60 CONTINUE |
|
237 END IF |
|
238 DO 80, L = 1, K |
|
239 IF( B( L, J ).NE.ZERO )THEN |
|
240 TEMP = ALPHA*B( L, J ) |
|
241 DO 70, I = 1, M |
|
242 C( I, J ) = C( I, J ) + TEMP*A( I, L ) |
|
243 70 CONTINUE |
|
244 END IF |
|
245 80 CONTINUE |
|
246 90 CONTINUE |
|
247 ELSE |
|
248 * |
|
249 * Form C := alpha*A'*B + beta*C |
|
250 * |
|
251 DO 120, J = 1, N |
|
252 DO 110, I = 1, M |
|
253 TEMP = ZERO |
|
254 DO 100, L = 1, K |
|
255 TEMP = TEMP + A( L, I )*B( L, J ) |
|
256 100 CONTINUE |
|
257 IF( BETA.EQ.ZERO )THEN |
|
258 C( I, J ) = ALPHA*TEMP |
|
259 ELSE |
|
260 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) |
|
261 END IF |
|
262 110 CONTINUE |
|
263 120 CONTINUE |
|
264 END IF |
|
265 ELSE |
|
266 IF( NOTA )THEN |
|
267 * |
|
268 * Form C := alpha*A*B' + beta*C |
|
269 * |
|
270 DO 170, J = 1, N |
|
271 IF( BETA.EQ.ZERO )THEN |
|
272 DO 130, I = 1, M |
|
273 C( I, J ) = ZERO |
|
274 130 CONTINUE |
|
275 ELSE IF( BETA.NE.ONE )THEN |
|
276 DO 140, I = 1, M |
|
277 C( I, J ) = BETA*C( I, J ) |
|
278 140 CONTINUE |
|
279 END IF |
|
280 DO 160, L = 1, K |
|
281 IF( B( J, L ).NE.ZERO )THEN |
|
282 TEMP = ALPHA*B( J, L ) |
|
283 DO 150, I = 1, M |
|
284 C( I, J ) = C( I, J ) + TEMP*A( I, L ) |
|
285 150 CONTINUE |
|
286 END IF |
|
287 160 CONTINUE |
|
288 170 CONTINUE |
|
289 ELSE |
|
290 * |
|
291 * Form C := alpha*A'*B' + beta*C |
|
292 * |
|
293 DO 200, J = 1, N |
|
294 DO 190, I = 1, M |
|
295 TEMP = ZERO |
|
296 DO 180, L = 1, K |
|
297 TEMP = TEMP + A( L, I )*B( J, L ) |
|
298 180 CONTINUE |
|
299 IF( BETA.EQ.ZERO )THEN |
|
300 C( I, J ) = ALPHA*TEMP |
|
301 ELSE |
|
302 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) |
|
303 END IF |
|
304 190 CONTINUE |
|
305 200 CONTINUE |
|
306 END IF |
|
307 END IF |
|
308 * |
|
309 RETURN |
|
310 * |
|
311 * End of SGEMM . |
|
312 * |
|
313 END |