1
SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
3
* .. Scalar Arguments ..
4
CHARACTER*1 TRANSA, TRANSB
5
INTEGER M, N, K, LDA, LDB, LDC
7
* .. Array Arguments ..
8
COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
14
* ZGEMM performs one of the matrix-matrix operations
16
* C := alpha*op( A )*op( B ) + beta*C,
18
* where op( X ) is one of
20
* op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ),
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.
28
* TRANSA - CHARACTER*1.
29
* On entry, TRANSA specifies the form of op( A ) to be used in
30
* the matrix multiplication as follows:
32
* TRANSA = 'N' or 'n', op( A ) = A.
34
* TRANSA = 'T' or 't', op( A ) = A'.
36
* TRANSA = 'C' or 'c', op( A ) = conjg( A' ).
40
* TRANSB - CHARACTER*1.
41
* On entry, TRANSB specifies the form of op( B ) to be used in
42
* the matrix multiplication as follows:
44
* TRANSB = 'N' or 'n', op( B ) = B.
46
* TRANSB = 'T' or 't', op( B ) = B'.
48
* TRANSB = 'C' or 'c', op( B ) = conjg( B' ).
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.
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
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
69
* ALPHA - COMPLEX*16 .
70
* On entry, ALPHA specifies the scalar alpha.
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
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
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
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
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.
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 ).
116
* On entry, LDC specifies the first dimension of C as declared
117
* in the calling (sub) program. LDC must be at least
122
* Level 3 Blas routine.
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.
131
* .. External Functions ..
134
* .. External Subroutines ..
136
* .. Intrinsic Functions ..
137
INTRINSIC DCONJG, MAX
138
* .. Local Scalars ..
139
LOGICAL CONJA, CONJB, NOTA, NOTB
140
INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB
144
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
146
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
148
* .. Executable Statements ..
150
* Set NOTA and NOTB as true if A and B respectively are not
151
* conjugated or transposed, set CONJA and CONJB as true if A and
152
* B respectively are to be transposed but not conjugated and set
153
* NROWA, NCOLA and NROWB as the number of rows and columns of A
154
* and the number of rows of B respectively.
156
NOTA = LSAME( TRANSA, 'N' )
157
NOTB = LSAME( TRANSB, 'N' )
158
CONJA = LSAME( TRANSA, 'C' )
159
CONJB = LSAME( TRANSB, 'C' )
173
* Test the input parameters.
176
IF( ( .NOT.NOTA ).AND.
177
$ ( .NOT.CONJA ).AND.
178
$ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN
180
ELSE IF( ( .NOT.NOTB ).AND.
181
$ ( .NOT.CONJB ).AND.
182
$ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN
184
ELSE IF( M .LT.0 )THEN
186
ELSE IF( N .LT.0 )THEN
188
ELSE IF( K .LT.0 )THEN
190
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
192
ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
194
ELSE IF( LDC.LT.MAX( 1, M ) )THEN
198
CALL XERBLA( 'ZGEMM ', INFO )
202
* Quick return if possible.
204
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
205
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
208
* And when alpha.eq.zero.
210
IF( ALPHA.EQ.ZERO )THEN
211
IF( BETA.EQ.ZERO )THEN
220
C( I, J ) = BETA*C( I, J )
227
* Start the operations.
232
* Form C := alpha*A*B + beta*C.
235
IF( BETA.EQ.ZERO )THEN
239
ELSE IF( BETA.NE.ONE )THEN
241
C( I, J ) = BETA*C( I, J )
245
IF( B( L, J ).NE.ZERO )THEN
246
TEMP = ALPHA*B( L, J )
248
C( I, J ) = C( I, J ) + TEMP*A( I, L )
255
* Form C := alpha*conjg( A' )*B + beta*C.
261
TEMP = TEMP + DCONJG( A( L, I ) )*B( L, J )
263
IF( BETA.EQ.ZERO )THEN
264
C( I, J ) = ALPHA*TEMP
266
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
272
* Form C := alpha*A'*B + beta*C
278
TEMP = TEMP + A( L, I )*B( L, J )
280
IF( BETA.EQ.ZERO )THEN
281
C( I, J ) = ALPHA*TEMP
283
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
291
* Form C := alpha*A*conjg( B' ) + beta*C.
294
IF( BETA.EQ.ZERO )THEN
298
ELSE IF( BETA.NE.ONE )THEN
300
C( I, J ) = BETA*C( I, J )
304
IF( B( J, L ).NE.ZERO )THEN
305
TEMP = ALPHA*DCONJG( B( J, L ) )
307
C( I, J ) = C( I, J ) + TEMP*A( I, L )
314
* Form C := alpha*A*B' + beta*C
317
IF( BETA.EQ.ZERO )THEN
321
ELSE IF( BETA.NE.ONE )THEN
323
C( I, J ) = BETA*C( I, J )
327
IF( B( J, L ).NE.ZERO )THEN
328
TEMP = ALPHA*B( J, L )
330
C( I, J ) = C( I, J ) + TEMP*A( I, L )
339
* Form C := alpha*conjg( A' )*conjg( B' ) + beta*C.
346
$ DCONJG( A( L, I ) )*DCONJG( B( J, L ) )
348
IF( BETA.EQ.ZERO )THEN
349
C( I, J ) = ALPHA*TEMP
351
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
357
* Form C := alpha*conjg( A' )*B' + beta*C
363
TEMP = TEMP + DCONJG( A( L, I ) )*B( J, L )
365
IF( BETA.EQ.ZERO )THEN
366
C( I, J ) = ALPHA*TEMP
368
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
376
* Form C := alpha*A'*conjg( B' ) + beta*C
382
TEMP = TEMP + A( L, I )*DCONJG( B( J, L ) )
384
IF( BETA.EQ.ZERO )THEN
385
C( I, J ) = ALPHA*TEMP
387
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
393
* Form C := alpha*A'*B' + beta*C
399
TEMP = TEMP + A( L, I )*B( J, L )
401
IF( BETA.EQ.ZERO )THEN
402
C( I, J ) = ALPHA*TEMP
404
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )