LAPACK: BLAS/SRC/dgemm.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
00002 *     .. Scalar Arguments ..
00003       DOUBLE PRECISION ALPHA,BETA
00004       INTEGER K,LDA,LDB,LDC,M,N
00005       CHARACTER TRANSA,TRANSB
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  DGEMM  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',
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 ) = 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 ) = 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  - DOUBLE PRECISION.
00070 *           On entry, ALPHA specifies the scalar alpha.
00071 *           Unchanged on exit.
00072 *
00073 *  A      - DOUBLE PRECISION 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      - DOUBLE PRECISION 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   - DOUBLE PRECISION.
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      - DOUBLE PRECISION 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 MAX
00143 *     ..
00144 *     .. Local Scalars ..
00145       DOUBLE PRECISION TEMP
00146       INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
00147       LOGICAL NOTA,NOTB
00148 *     ..
00149 *     .. Parameters ..
00150       DOUBLE PRECISION ONE,ZERO
00151       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
00152 *     ..
00153 *
00154 *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
00155 *     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
00156 *     and  columns of  A  and the  number of  rows  of  B  respectively.
00157 *
00158       NOTA = LSAME(TRANSA,'N')
00159       NOTB = LSAME(TRANSB,'N')
00160       IF (NOTA) THEN
00161           NROWA = M
00162           NCOLA = K
00163       ELSE
00164           NROWA = K
00165           NCOLA = M
00166       END IF
00167       IF (NOTB) THEN
00168           NROWB = K
00169       ELSE
00170           NROWB = N
00171       END IF
00172 *
00173 *     Test the input parameters.
00174 *
00175       INFO = 0
00176       IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
00177      +    (.NOT.LSAME(TRANSA,'T'))) THEN
00178           INFO = 1
00179       ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
00180      +         (.NOT.LSAME(TRANSB,'T'))) THEN
00181           INFO = 2
00182       ELSE IF (M.LT.0) THEN
00183           INFO = 3
00184       ELSE IF (N.LT.0) THEN
00185           INFO = 4
00186       ELSE IF (K.LT.0) THEN
00187           INFO = 5
00188       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00189           INFO = 8
00190       ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
00191           INFO = 10
00192       ELSE IF (LDC.LT.MAX(1,M)) THEN
00193           INFO = 13
00194       END IF
00195       IF (INFO.NE.0) THEN
00196           CALL XERBLA('DGEMM ',INFO)
00197           RETURN
00198       END IF
00199 *
00200 *     Quick return if possible.
00201 *
00202       IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
00203      +    (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
00204 *
00205 *     And if  alpha.eq.zero.
00206 *
00207       IF (ALPHA.EQ.ZERO) THEN
00208           IF (BETA.EQ.ZERO) THEN
00209               DO 20 J = 1,N
00210                   DO 10 I = 1,M
00211                       C(I,J) = ZERO
00212    10             CONTINUE
00213    20         CONTINUE
00214           ELSE
00215               DO 40 J = 1,N
00216                   DO 30 I = 1,M
00217                       C(I,J) = BETA*C(I,J)
00218    30             CONTINUE
00219    40         CONTINUE
00220           END IF
00221           RETURN
00222       END IF
00223 *
00224 *     Start the operations.
00225 *
00226       IF (NOTB) THEN
00227           IF (NOTA) THEN
00228 *
00229 *           Form  C := alpha*A*B + beta*C.
00230 *
00231               DO 90 J = 1,N
00232                   IF (BETA.EQ.ZERO) THEN
00233                       DO 50 I = 1,M
00234                           C(I,J) = ZERO
00235    50                 CONTINUE
00236                   ELSE IF (BETA.NE.ONE) THEN
00237                       DO 60 I = 1,M
00238                           C(I,J) = BETA*C(I,J)
00239    60                 CONTINUE
00240                   END IF
00241                   DO 80 L = 1,K
00242                       IF (B(L,J).NE.ZERO) THEN
00243                           TEMP = ALPHA*B(L,J)
00244                           DO 70 I = 1,M
00245                               C(I,J) = C(I,J) + TEMP*A(I,L)
00246    70                     CONTINUE
00247                       END IF
00248    80             CONTINUE
00249    90         CONTINUE
00250           ELSE
00251 *
00252 *           Form  C := alpha*A'*B + beta*C
00253 *
00254               DO 120 J = 1,N
00255                   DO 110 I = 1,M
00256                       TEMP = ZERO
00257                       DO 100 L = 1,K
00258                           TEMP = TEMP + A(L,I)*B(L,J)
00259   100                 CONTINUE
00260                       IF (BETA.EQ.ZERO) THEN
00261                           C(I,J) = ALPHA*TEMP
00262                       ELSE
00263                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00264                       END IF
00265   110             CONTINUE
00266   120         CONTINUE
00267           END IF
00268       ELSE
00269           IF (NOTA) THEN
00270 *
00271 *           Form  C := alpha*A*B' + beta*C
00272 *
00273               DO 170 J = 1,N
00274                   IF (BETA.EQ.ZERO) THEN
00275                       DO 130 I = 1,M
00276                           C(I,J) = ZERO
00277   130                 CONTINUE
00278                   ELSE IF (BETA.NE.ONE) THEN
00279                       DO 140 I = 1,M
00280                           C(I,J) = BETA*C(I,J)
00281   140                 CONTINUE
00282                   END IF
00283                   DO 160 L = 1,K
00284                       IF (B(J,L).NE.ZERO) THEN
00285                           TEMP = ALPHA*B(J,L)
00286                           DO 150 I = 1,M
00287                               C(I,J) = C(I,J) + TEMP*A(I,L)
00288   150                     CONTINUE
00289                       END IF
00290   160             CONTINUE
00291   170         CONTINUE
00292           ELSE
00293 *
00294 *           Form  C := alpha*A'*B' + beta*C
00295 *
00296               DO 200 J = 1,N
00297                   DO 190 I = 1,M
00298                       TEMP = ZERO
00299                       DO 180 L = 1,K
00300                           TEMP = TEMP + A(L,I)*B(J,L)
00301   180                 CONTINUE
00302                       IF (BETA.EQ.ZERO) THEN
00303                           C(I,J) = ALPHA*TEMP
00304                       ELSE
00305                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
00306                       END IF
00307   190             CONTINUE
00308   200         CONTINUE
00309           END IF
00310       END IF
00311 *
00312       RETURN
00313 *
00314 *     End of DGEMM .
00315 *
00316       END