LAPACK: BLAS/SRC/zgemm.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
00002 *     .. Scalar Arguments ..
00003       DOUBLE COMPLEX ALPHA,BETA
00004       INTEGER K,LDA,LDB,LDC,M,N
00005       CHARACTER TRANSA,TRANSB
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  ZGEMM  performs one of the matrix-matrix operations
00015 *
00016 *     C := alpha*op( A )*op( B ) + beta*C,
00017 *
00018 *  where  op( X ) is one of
00019 *
00020 *     op( X ) = X   or   op( X ) = X'   or   op( X ) = conjg( X' ),
00021 *
00022 *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
00023 *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
00024 *
00025 *  Arguments
00026 *  ==========
00027 *
00028 *  TRANSA - CHARACTER*1.
00029 *           On entry, TRANSA specifies the form of op( A ) to be used in
00030 *           the matrix multiplication as follows:
00031 *
00032 *              TRANSA = 'N' or 'n',  op( A ) = A.
00033 *
00034 *              TRANSA = 'T' or 't',  op( A ) = A'.
00035 *
00036 *              TRANSA = 'C' or 'c',  op( A ) = conjg( A' ).
00037 *
00038 *           Unchanged on exit.
00039 *
00040 *  TRANSB - CHARACTER*1.
00041 *           On entry, TRANSB specifies the form of op( B ) to be used in
00042 *           the matrix multiplication as follows:
00043 *
00044 *              TRANSB = 'N' or 'n',  op( B ) = B.
00045 *
00046 *              TRANSB = 'T' or 't',  op( B ) = B'.
00047 *
00048 *              TRANSB = 'C' or 'c',  op( B ) = conjg( B' ).
00049 *
00050 *           Unchanged on exit.
00051 *
00052 *  M      - INTEGER.
00053 *           On entry,  M  specifies  the number  of rows  of the  matrix
00054 *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
00055 *           Unchanged on exit.
00056 *
00057 *  N      - INTEGER.
00058 *           On entry,  N  specifies the number  of columns of the matrix
00059 *           op( B ) and the number of columns of the matrix C. N must be
00060 *           at least zero.
00061 *           Unchanged on exit.
00062 *
00063 *  K      - INTEGER.
00064 *           On entry,  K  specifies  the number of columns of the matrix
00065 *           op( A ) and the number of rows of the matrix op( B ). K must
00066 *           be at least  zero.
00067 *           Unchanged on exit.
00068 *
00069 *  ALPHA  - COMPLEX*16      .
00070 *           On entry, ALPHA specifies the scalar alpha.
00071 *           Unchanged on exit.
00072 *
00073 *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
00074 *           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
00075 *           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
00076 *           part of the array  A  must contain the matrix  A,  otherwise
00077 *           the leading  k by m  part of the array  A  must contain  the
00078 *           matrix A.
00079 *           Unchanged on exit.
00080 *
00081 *  LDA    - INTEGER.
00082 *           On entry, LDA specifies the first dimension of A as declared
00083 *           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
00084 *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
00085 *           least  max( 1, k ).
00086 *           Unchanged on exit.
00087 *
00088 *  B      - COMPLEX*16       array of DIMENSION ( LDB, kb ), where kb is
00089 *           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
00090 *           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
00091 *           part of the array  B  must contain the matrix  B,  otherwise
00092 *           the leading  n by k  part of the array  B  must contain  the
00093 *           matrix B.
00094 *           Unchanged on exit.
00095 *
00096 *  LDB    - INTEGER.
00097 *           On entry, LDB specifies the first dimension of B as declared
00098 *           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
00099 *           LDB must be at least  max( 1, k ), otherwise  LDB must be at
00100 *           least  max( 1, n ).
00101 *           Unchanged on exit.
00102 *
00103 *  BETA   - COMPLEX*16      .
00104 *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
00105 *           supplied as zero then C need not be set on input.
00106 *           Unchanged on exit.
00107 *
00108 *  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
00109 *           Before entry, the leading  m by n  part of the array  C must
00110 *           contain the matrix  C,  except when  beta  is zero, in which
00111 *           case C need not be set on entry.
00112 *           On exit, the array  C  is overwritten by the  m by n  matrix
00113 *           ( alpha*op( A )*op( B ) + beta*C ).
00114 *
00115 *  LDC    - INTEGER.
00116 *           On entry, LDC specifies the first dimension of C as declared
00117 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
00118 *           max( 1, m ).
00119 *           Unchanged on exit.
00120 *
00121 *  Further Details
00122 *  ===============
00123 *
00124 *  Level 3 Blas routine.
00125 *
00126 *  -- Written on 8-February-1989.
00127 *     Jack Dongarra, Argonne National Laboratory.
00128 *     Iain Duff, AERE Harwell.
00129 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
00130 *     Sven Hammarling, Numerical Algorithms Group Ltd.
00131 *
00132 *  =====================================================================
00133 *
00134 *     .. External Functions ..
00135       LOGICAL LSAME
00136       EXTERNAL LSAME
00137 *     ..
00138 *     .. External Subroutines ..
00139       EXTERNAL XERBLA
00140 *     ..
00141 *     .. Intrinsic Functions ..
00142       INTRINSIC DCONJG,MAX
00143 *     ..
00144 *     .. Local Scalars ..
00145       DOUBLE COMPLEX TEMP
00146       INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
00147       LOGICAL CONJA,CONJB,NOTA,NOTB
00148 *     ..
00149 *     .. Parameters ..
00150       DOUBLE COMPLEX ONE
00151       PARAMETER (ONE= (1.0D+0,0.0D+0))
00152       DOUBLE COMPLEX ZERO
00153       PARAMETER (ZERO= (0.0D+0,0.0D+0))
00154 *     ..
00155 *
00156 *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
00157 *     conjugated or transposed, set  CONJA and CONJB  as true if  A  and
00158 *     B  respectively are to be  transposed but  not conjugated  and set
00159 *     NROWA, NCOLA and  NROWB  as the number of rows and  columns  of  A
00160 *     and the number of rows of  B  respectively.
00161 *
00162       NOTA = LSAME(TRANSA,'N')
00163       NOTB = LSAME(TRANSB,'N')
00164       CONJA = LSAME(TRANSA,'C')
00165       CONJB = LSAME(TRANSB,'C')
00166       IF (NOTA) THEN
00167           NROWA = M
00168           NCOLA = K
00169       ELSE
00170           NROWA = K
00171           NCOLA = M
00172       END IF
00173       IF (NOTB) THEN
00174           NROWB = K
00175       ELSE
00176           NROWB = N
00177       END IF
00178 *
00179 *     Test the input parameters.
00180 *
00181       INFO = 0
00182       IF ((.NOT.NOTA) .AND. (.NOT.CONJA) .AND.
00183      +    (.NOT.LSAME(TRANSA,'T'))) THEN
00184           INFO = 1
00185       ELSE IF ((.NOT.NOTB) .AND. (.NOT.CONJB) .AND.
00186      +         (.NOT.LSAME(TRANSB,'T'))) THEN
00187           INFO = 2
00188       ELSE IF (M.LT.0) THEN
00189           INFO = 3
00190       ELSE IF (N.LT.0) THEN
00191           INFO = 4
00192       ELSE IF (K.LT.0) THEN
00193           INFO = 5
00194       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00195           INFO = 8
00196       ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
00197           INFO = 10
00198       ELSE IF (LDC.LT.MAX(1,M)) THEN
00199           INFO = 13
00200       END IF
00201       IF (INFO.NE.0) THEN
00202           CALL XERBLA('ZGEMM ',INFO)
00203           RETURN
00204       END IF
00205 *
00206 *     Quick return if possible.
00207 *
00208       IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
00209      +    (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
00210 *
00211 *     And when  alpha.eq.zero.
00212 *
00213       IF (ALPHA.EQ.ZERO) THEN
00214           IF (BETA.EQ.ZERO) THEN
00215               DO 20 J = 1,N
00216                   DO 10 I = 1,M
00217                       C(I,J) = ZERO
00218    10             CONTINUE
00219    20         CONTINUE
00220           ELSE
00221               DO 40 J = 1,N
00222                   DO 30 I = 1,M
00223                       C(I,J) = BETA*C(I,J)
00224    30             CONTINUE
00225    40         CONTINUE
00226           END IF
00227           RETURN
00228       END IF
00229 *
00230 *     Start the operations.
00231 *
00232       IF (NOTB) THEN
00233           IF (NOTA) THEN
00234 *
00235 *           Form  C := alpha*A*B + beta*C.
00236 *
00237               DO 90 J = 1,N
00238                   IF (BETA.EQ.ZERO) THEN
00239                       DO 50 I = 1,M
00240                           C(I,J) = ZERO
00241    50                 CONTINUE
00242                   ELSE IF (BETA.NE.ONE) THEN
00243                       DO 60 I = 1,M
00244                           C(I,J) = BETA*C(I,J)
00245    60                 CONTINUE
00246                   END IF
00247                   DO 80 L = 1,K
00248                       IF (B(L,J).NE.ZERO) THEN
00249                           TEMP = ALPHA*B(L,J)
00250                           DO 70 I = 1,M
00251                               C(I,J) = C(I,J) + TEMP*A(I,L)
00252    70                     CONTINUE
00253                       END IF
00254    80             CONTINUE
00255    90         CONTINUE
00256           ELSE IF (CONJA) THEN
00257 *
00258 *           Form  C := alpha*conjg( A' )*B + beta*C.
00259 *
00260               DO 120 J = 1,N
00261                   DO 110 I = 1,M
00262                       TEMP = ZERO
00263                       DO 100 L = 1,K
00264                           TEMP = TEMP + DCONJG(A(L,I))*B(L,J)
00265   100                 CONTINUE
00266                       IF (BETA.EQ.ZERO) THEN
00267                           C(I,J) = ALPHA*TEMP
00268                       ELSE
00269                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00270                       END IF
00271   110             CONTINUE
00272   120         CONTINUE
00273           ELSE
00274 *
00275 *           Form  C := alpha*A'*B + beta*C
00276 *
00277               DO 150 J = 1,N
00278                   DO 140 I = 1,M
00279                       TEMP = ZERO
00280                       DO 130 L = 1,K
00281                           TEMP = TEMP + A(L,I)*B(L,J)
00282   130                 CONTINUE
00283                       IF (BETA.EQ.ZERO) THEN
00284                           C(I,J) = ALPHA*TEMP
00285                       ELSE
00286                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00287                       END IF
00288   140             CONTINUE
00289   150         CONTINUE
00290           END IF
00291       ELSE IF (NOTA) THEN
00292           IF (CONJB) THEN
00293 *
00294 *           Form  C := alpha*A*conjg( B' ) + beta*C.
00295 *
00296               DO 200 J = 1,N
00297                   IF (BETA.EQ.ZERO) THEN
00298                       DO 160 I = 1,M
00299                           C(I,J) = ZERO
00300   160                 CONTINUE
00301                   ELSE IF (BETA.NE.ONE) THEN
00302                       DO 170 I = 1,M
00303                           C(I,J) = BETA*C(I,J)
00304   170                 CONTINUE
00305                   END IF
00306                   DO 190 L = 1,K
00307                       IF (B(J,L).NE.ZERO) THEN
00308                           TEMP = ALPHA*DCONJG(B(J,L))
00309                           DO 180 I = 1,M
00310                               C(I,J) = C(I,J) + TEMP*A(I,L)
00311   180                     CONTINUE
00312                       END IF
00313   190             CONTINUE
00314   200         CONTINUE
00315           ELSE
00316 *
00317 *           Form  C := alpha*A*B'          + beta*C
00318 *
00319               DO 250 J = 1,N
00320                   IF (BETA.EQ.ZERO) THEN
00321                       DO 210 I = 1,M
00322                           C(I,J) = ZERO
00323   210                 CONTINUE
00324                   ELSE IF (BETA.NE.ONE) THEN
00325                       DO 220 I = 1,M
00326                           C(I,J) = BETA*C(I,J)
00327   220                 CONTINUE
00328                   END IF
00329                   DO 240 L = 1,K
00330                       IF (B(J,L).NE.ZERO) THEN
00331                           TEMP = ALPHA*B(J,L)
00332                           DO 230 I = 1,M
00333                               C(I,J) = C(I,J) + TEMP*A(I,L)
00334   230                     CONTINUE
00335                       END IF
00336   240             CONTINUE
00337   250         CONTINUE
00338           END IF
00339       ELSE IF (CONJA) THEN
00340           IF (CONJB) THEN
00341 *
00342 *           Form  C := alpha*conjg( A' )*conjg( B' ) + beta*C.
00343 *
00344               DO 280 J = 1,N
00345                   DO 270 I = 1,M
00346                       TEMP = ZERO
00347                       DO 260 L = 1,K
00348                           TEMP = TEMP + DCONJG(A(L,I))*DCONJG(B(J,L))
00349   260                 CONTINUE
00350                       IF (BETA.EQ.ZERO) THEN
00351                           C(I,J) = ALPHA*TEMP
00352                       ELSE
00353                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00354                       END IF
00355   270             CONTINUE
00356   280         CONTINUE
00357           ELSE
00358 *
00359 *           Form  C := alpha*conjg( A' )*B' + beta*C
00360 *
00361               DO 310 J = 1,N
00362                   DO 300 I = 1,M
00363                       TEMP = ZERO
00364                       DO 290 L = 1,K
00365                           TEMP = TEMP + DCONJG(A(L,I))*B(J,L)
00366   290                 CONTINUE
00367                       IF (BETA.EQ.ZERO) THEN
00368                           C(I,J) = ALPHA*TEMP
00369                       ELSE
00370                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00371                       END IF
00372   300             CONTINUE
00373   310         CONTINUE
00374           END IF
00375       ELSE
00376           IF (CONJB) THEN
00377 *
00378 *           Form  C := alpha*A'*conjg( B' ) + beta*C
00379 *
00380               DO 340 J = 1,N
00381                   DO 330 I = 1,M
00382                       TEMP = ZERO
00383                       DO 320 L = 1,K
00384                           TEMP = TEMP + A(L,I)*DCONJG(B(J,L))
00385   320                 CONTINUE
00386                       IF (BETA.EQ.ZERO) THEN
00387                           C(I,J) = ALPHA*TEMP
00388                       ELSE
00389                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00390                       END IF
00391   330             CONTINUE
00392   340         CONTINUE
00393           ELSE
00394 *
00395 *           Form  C := alpha*A'*B' + beta*C
00396 *
00397               DO 370 J = 1,N
00398                   DO 360 I = 1,M
00399                       TEMP = ZERO
00400                       DO 350 L = 1,K
00401                           TEMP = TEMP + A(L,I)*B(J,L)
00402   350                 CONTINUE
00403                       IF (BETA.EQ.ZERO) THEN
00404                           C(I,J) = ALPHA*TEMP
00405                       ELSE
00406                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00407                       END IF
00408   360             CONTINUE
00409   370         CONTINUE
00410           END IF
00411       END IF
00412 *
00413       RETURN
00414 *
00415 *     End of ZGEMM .
00416 *
00417       END