001:       SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
002: *     .. Scalar Arguments ..
003:       DOUBLE COMPLEX ALPHA,BETA
004:       INTEGER K,LDA,LDB,LDC,M,N
005:       CHARACTER TRANSA,TRANSB
006: *     ..
007: *     .. Array Arguments ..
008:       DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  ZGEMM  performs one of the matrix-matrix operations
015: *
016: *     C := alpha*op( A )*op( B ) + beta*C,
017: *
018: *  where  op( X ) is one of
019: *
020: *     op( X ) = X   or   op( X ) = X'   or   op( X ) = conjg( X' ),
021: *
022: *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
023: *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
024: *
025: *  Arguments
026: *  ==========
027: *
028: *  TRANSA - CHARACTER*1.
029: *           On entry, TRANSA specifies the form of op( A ) to be used in
030: *           the matrix multiplication as follows:
031: *
032: *              TRANSA = 'N' or 'n',  op( A ) = A.
033: *
034: *              TRANSA = 'T' or 't',  op( A ) = A'.
035: *
036: *              TRANSA = 'C' or 'c',  op( A ) = conjg( A' ).
037: *
038: *           Unchanged on exit.
039: *
040: *  TRANSB - CHARACTER*1.
041: *           On entry, TRANSB specifies the form of op( B ) to be used in
042: *           the matrix multiplication as follows:
043: *
044: *              TRANSB = 'N' or 'n',  op( B ) = B.
045: *
046: *              TRANSB = 'T' or 't',  op( B ) = B'.
047: *
048: *              TRANSB = 'C' or 'c',  op( B ) = conjg( B' ).
049: *
050: *           Unchanged on exit.
051: *
052: *  M      - INTEGER.
053: *           On entry,  M  specifies  the number  of rows  of the  matrix
054: *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
055: *           Unchanged on exit.
056: *
057: *  N      - INTEGER.
058: *           On entry,  N  specifies the number  of columns of the matrix
059: *           op( B ) and the number of columns of the matrix C. N must be
060: *           at least zero.
061: *           Unchanged on exit.
062: *
063: *  K      - INTEGER.
064: *           On entry,  K  specifies  the number of columns of the matrix
065: *           op( A ) and the number of rows of the matrix op( B ). K must
066: *           be at least  zero.
067: *           Unchanged on exit.
068: *
069: *  ALPHA  - COMPLEX*16      .
070: *           On entry, ALPHA specifies the scalar alpha.
071: *           Unchanged on exit.
072: *
073: *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
074: *           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
075: *           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
076: *           part of the array  A  must contain the matrix  A,  otherwise
077: *           the leading  k by m  part of the array  A  must contain  the
078: *           matrix A.
079: *           Unchanged on exit.
080: *
081: *  LDA    - INTEGER.
082: *           On entry, LDA specifies the first dimension of A as declared
083: *           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
084: *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
085: *           least  max( 1, k ).
086: *           Unchanged on exit.
087: *
088: *  B      - COMPLEX*16       array of DIMENSION ( LDB, kb ), where kb is
089: *           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
090: *           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
091: *           part of the array  B  must contain the matrix  B,  otherwise
092: *           the leading  n by k  part of the array  B  must contain  the
093: *           matrix B.
094: *           Unchanged on exit.
095: *
096: *  LDB    - INTEGER.
097: *           On entry, LDB specifies the first dimension of B as declared
098: *           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
099: *           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   - COMPLEX*16      .
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      - COMPLEX*16       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 DCONJG,MAX
143: *     ..
144: *     .. Local Scalars ..
145:       DOUBLE COMPLEX TEMP
146:       INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
147:       LOGICAL CONJA,CONJB,NOTA,NOTB
148: *     ..
149: *     .. Parameters ..
150:       DOUBLE COMPLEX ONE
151:       PARAMETER (ONE= (1.0D+0,0.0D+0))
152:       DOUBLE COMPLEX ZERO
153:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
154: *     ..
155: *
156: *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
157: *     conjugated or transposed, set  CONJA and CONJB  as true if  A  and
158: *     B  respectively are to be  transposed but  not conjugated  and set
159: *     NROWA, NCOLA and  NROWB  as the number of rows and  columns  of  A
160: *     and the number of rows of  B  respectively.
161: *
162:       NOTA = LSAME(TRANSA,'N')
163:       NOTB = LSAME(TRANSB,'N')
164:       CONJA = LSAME(TRANSA,'C')
165:       CONJB = LSAME(TRANSB,'C')
166:       IF (NOTA) THEN
167:           NROWA = M
168:           NCOLA = K
169:       ELSE
170:           NROWA = K
171:           NCOLA = M
172:       END IF
173:       IF (NOTB) THEN
174:           NROWB = K
175:       ELSE
176:           NROWB = N
177:       END IF
178: *
179: *     Test the input parameters.
180: *
181:       INFO = 0
182:       IF ((.NOT.NOTA) .AND. (.NOT.CONJA) .AND.
183:      +    (.NOT.LSAME(TRANSA,'T'))) THEN
184:           INFO = 1
185:       ELSE IF ((.NOT.NOTB) .AND. (.NOT.CONJB) .AND.
186:      +         (.NOT.LSAME(TRANSB,'T'))) THEN
187:           INFO = 2
188:       ELSE IF (M.LT.0) THEN
189:           INFO = 3
190:       ELSE IF (N.LT.0) THEN
191:           INFO = 4
192:       ELSE IF (K.LT.0) THEN
193:           INFO = 5
194:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
195:           INFO = 8
196:       ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
197:           INFO = 10
198:       ELSE IF (LDC.LT.MAX(1,M)) THEN
199:           INFO = 13
200:       END IF
201:       IF (INFO.NE.0) THEN
202:           CALL XERBLA('ZGEMM ',INFO)
203:           RETURN
204:       END IF
205: *
206: *     Quick return if possible.
207: *
208:       IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
209:      +    (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
210: *
211: *     And when  alpha.eq.zero.
212: *
213:       IF (ALPHA.EQ.ZERO) THEN
214:           IF (BETA.EQ.ZERO) THEN
215:               DO 20 J = 1,N
216:                   DO 10 I = 1,M
217:                       C(I,J) = ZERO
218:    10             CONTINUE
219:    20         CONTINUE
220:           ELSE
221:               DO 40 J = 1,N
222:                   DO 30 I = 1,M
223:                       C(I,J) = BETA*C(I,J)
224:    30             CONTINUE
225:    40         CONTINUE
226:           END IF
227:           RETURN
228:       END IF
229: *
230: *     Start the operations.
231: *
232:       IF (NOTB) THEN
233:           IF (NOTA) THEN
234: *
235: *           Form  C := alpha*A*B + beta*C.
236: *
237:               DO 90 J = 1,N
238:                   IF (BETA.EQ.ZERO) THEN
239:                       DO 50 I = 1,M
240:                           C(I,J) = ZERO
241:    50                 CONTINUE
242:                   ELSE IF (BETA.NE.ONE) THEN
243:                       DO 60 I = 1,M
244:                           C(I,J) = BETA*C(I,J)
245:    60                 CONTINUE
246:                   END IF
247:                   DO 80 L = 1,K
248:                       IF (B(L,J).NE.ZERO) THEN
249:                           TEMP = ALPHA*B(L,J)
250:                           DO 70 I = 1,M
251:                               C(I,J) = C(I,J) + TEMP*A(I,L)
252:    70                     CONTINUE
253:                       END IF
254:    80             CONTINUE
255:    90         CONTINUE
256:           ELSE IF (CONJA) THEN
257: *
258: *           Form  C := alpha*conjg( A' )*B + beta*C.
259: *
260:               DO 120 J = 1,N
261:                   DO 110 I = 1,M
262:                       TEMP = ZERO
263:                       DO 100 L = 1,K
264:                           TEMP = TEMP + DCONJG(A(L,I))*B(L,J)
265:   100                 CONTINUE
266:                       IF (BETA.EQ.ZERO) THEN
267:                           C(I,J) = ALPHA*TEMP
268:                       ELSE
269:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
270:                       END IF
271:   110             CONTINUE
272:   120         CONTINUE
273:           ELSE
274: *
275: *           Form  C := alpha*A'*B + beta*C
276: *
277:               DO 150 J = 1,N
278:                   DO 140 I = 1,M
279:                       TEMP = ZERO
280:                       DO 130 L = 1,K
281:                           TEMP = TEMP + A(L,I)*B(L,J)
282:   130                 CONTINUE
283:                       IF (BETA.EQ.ZERO) THEN
284:                           C(I,J) = ALPHA*TEMP
285:                       ELSE
286:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
287:                       END IF
288:   140             CONTINUE
289:   150         CONTINUE
290:           END IF
291:       ELSE IF (NOTA) THEN
292:           IF (CONJB) THEN
293: *
294: *           Form  C := alpha*A*conjg( B' ) + beta*C.
295: *
296:               DO 200 J = 1,N
297:                   IF (BETA.EQ.ZERO) THEN
298:                       DO 160 I = 1,M
299:                           C(I,J) = ZERO
300:   160                 CONTINUE
301:                   ELSE IF (BETA.NE.ONE) THEN
302:                       DO 170 I = 1,M
303:                           C(I,J) = BETA*C(I,J)
304:   170                 CONTINUE
305:                   END IF
306:                   DO 190 L = 1,K
307:                       IF (B(J,L).NE.ZERO) THEN
308:                           TEMP = ALPHA*DCONJG(B(J,L))
309:                           DO 180 I = 1,M
310:                               C(I,J) = C(I,J) + TEMP*A(I,L)
311:   180                     CONTINUE
312:                       END IF
313:   190             CONTINUE
314:   200         CONTINUE
315:           ELSE
316: *
317: *           Form  C := alpha*A*B'          + beta*C
318: *
319:               DO 250 J = 1,N
320:                   IF (BETA.EQ.ZERO) THEN
321:                       DO 210 I = 1,M
322:                           C(I,J) = ZERO
323:   210                 CONTINUE
324:                   ELSE IF (BETA.NE.ONE) THEN
325:                       DO 220 I = 1,M
326:                           C(I,J) = BETA*C(I,J)
327:   220                 CONTINUE
328:                   END IF
329:                   DO 240 L = 1,K
330:                       IF (B(J,L).NE.ZERO) THEN
331:                           TEMP = ALPHA*B(J,L)
332:                           DO 230 I = 1,M
333:                               C(I,J) = C(I,J) + TEMP*A(I,L)
334:   230                     CONTINUE
335:                       END IF
336:   240             CONTINUE
337:   250         CONTINUE
338:           END IF
339:       ELSE IF (CONJA) THEN
340:           IF (CONJB) THEN
341: *
342: *           Form  C := alpha*conjg( A' )*conjg( B' ) + beta*C.
343: *
344:               DO 280 J = 1,N
345:                   DO 270 I = 1,M
346:                       TEMP = ZERO
347:                       DO 260 L = 1,K
348:                           TEMP = TEMP + DCONJG(A(L,I))*DCONJG(B(J,L))
349:   260                 CONTINUE
350:                       IF (BETA.EQ.ZERO) THEN
351:                           C(I,J) = ALPHA*TEMP
352:                       ELSE
353:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
354:                       END IF
355:   270             CONTINUE
356:   280         CONTINUE
357:           ELSE
358: *
359: *           Form  C := alpha*conjg( A' )*B' + beta*C
360: *
361:               DO 310 J = 1,N
362:                   DO 300 I = 1,M
363:                       TEMP = ZERO
364:                       DO 290 L = 1,K
365:                           TEMP = TEMP + DCONJG(A(L,I))*B(J,L)
366:   290                 CONTINUE
367:                       IF (BETA.EQ.ZERO) THEN
368:                           C(I,J) = ALPHA*TEMP
369:                       ELSE
370:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
371:                       END IF
372:   300             CONTINUE
373:   310         CONTINUE
374:           END IF
375:       ELSE
376:           IF (CONJB) THEN
377: *
378: *           Form  C := alpha*A'*conjg( B' ) + beta*C
379: *
380:               DO 340 J = 1,N
381:                   DO 330 I = 1,M
382:                       TEMP = ZERO
383:                       DO 320 L = 1,K
384:                           TEMP = TEMP + A(L,I)*DCONJG(B(J,L))
385:   320                 CONTINUE
386:                       IF (BETA.EQ.ZERO) THEN
387:                           C(I,J) = ALPHA*TEMP
388:                       ELSE
389:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
390:                       END IF
391:   330             CONTINUE
392:   340         CONTINUE
393:           ELSE
394: *
395: *           Form  C := alpha*A'*B' + beta*C
396: *
397:               DO 370 J = 1,N
398:                   DO 360 I = 1,M
399:                       TEMP = ZERO
400:                       DO 350 L = 1,K
401:                           TEMP = TEMP + A(L,I)*B(J,L)
402:   350                 CONTINUE
403:                       IF (BETA.EQ.ZERO) THEN
404:                           C(I,J) = ALPHA*TEMP
405:                       ELSE
406:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
407:                       END IF
408:   360             CONTINUE
409:   370         CONTINUE
410:           END IF
411:       END IF
412: *
413:       RETURN
414: *
415: *     End of ZGEMM .
416: *
417:       END
418: