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