1       SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  2 *     .. Scalar Arguments ..
  3       DOUBLE COMPLEX ALPHA,BETA
  4       INTEGER K,LDA,LDB,LDC,M,N
  5       CHARACTER TRANSA,TRANSB
  6 *     ..
  7 *     .. Array Arguments ..
  8       DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
  9 *     ..
 10 *
 11 *  Purpose
 12 *  =======
 13 *
 14 *  ZGEMM  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**T   or   op( X ) = X**H,
 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**T.
 35 *
 36 *              TRANSA = 'C' or 'c',  op( A ) = A**H.
 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**T.
 47 *
 48 *              TRANSB = 'C' or 'c',  op( B ) = B**H.
 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  - COMPLEX*16      .
 70 *           On entry, ALPHA specifies the scalar alpha.
 71 *           Unchanged on exit.
 72 *
 73 *  A      - COMPLEX*16       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      - COMPLEX*16       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   - 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.0THEN
189           INFO = 3
190       ELSE IF (N.LT.0THEN
191           INFO = 4
192       ELSE IF (K.LT.0THEN
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.0THEN
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*A**H*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**T*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*B**H + 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**T          + 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*A**H*B**H + 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*A**H*B**T + 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**T*B**H + 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**T*B**T + 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