Annotation of rpl/lapack/blas/dgemm.f, revision 1.17

1.8       bertrand    1: *> \brief \b DGEMM
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.14      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *  Definition:
                      9: *  ===========
                     10: *
                     11: *       SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
1.14      bertrand   12: *
1.8       bertrand   13: *       .. Scalar Arguments ..
                     14: *       DOUBLE PRECISION ALPHA,BETA
                     15: *       INTEGER K,LDA,LDB,LDC,M,N
                     16: *       CHARACTER TRANSA,TRANSB
                     17: *       ..
                     18: *       .. Array Arguments ..
                     19: *       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
                     20: *       ..
1.14      bertrand   21: *
1.8       bertrand   22: *
                     23: *> \par Purpose:
                     24: *  =============
                     25: *>
                     26: *> \verbatim
                     27: *>
                     28: *> DGEMM  performs one of the matrix-matrix operations
                     29: *>
                     30: *>    C := alpha*op( A )*op( B ) + beta*C,
                     31: *>
                     32: *> where  op( X ) is one of
                     33: *>
                     34: *>    op( X ) = X   or   op( X ) = X**T,
                     35: *>
                     36: *> alpha and beta are scalars, and A, B and C are matrices, with op( A )
                     37: *> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
                     38: *> \endverbatim
                     39: *
                     40: *  Arguments:
                     41: *  ==========
                     42: *
                     43: *> \param[in] TRANSA
                     44: *> \verbatim
                     45: *>          TRANSA is CHARACTER*1
                     46: *>           On entry, TRANSA specifies the form of op( A ) to be used in
                     47: *>           the matrix multiplication as follows:
                     48: *>
                     49: *>              TRANSA = 'N' or 'n',  op( A ) = A.
                     50: *>
                     51: *>              TRANSA = 'T' or 't',  op( A ) = A**T.
                     52: *>
                     53: *>              TRANSA = 'C' or 'c',  op( A ) = A**T.
                     54: *> \endverbatim
                     55: *>
                     56: *> \param[in] TRANSB
                     57: *> \verbatim
                     58: *>          TRANSB is CHARACTER*1
                     59: *>           On entry, TRANSB specifies the form of op( B ) to be used in
                     60: *>           the matrix multiplication as follows:
                     61: *>
                     62: *>              TRANSB = 'N' or 'n',  op( B ) = B.
                     63: *>
                     64: *>              TRANSB = 'T' or 't',  op( B ) = B**T.
                     65: *>
                     66: *>              TRANSB = 'C' or 'c',  op( B ) = B**T.
                     67: *> \endverbatim
                     68: *>
                     69: *> \param[in] M
                     70: *> \verbatim
                     71: *>          M is INTEGER
                     72: *>           On entry,  M  specifies  the number  of rows  of the  matrix
                     73: *>           op( A )  and of the  matrix  C.  M  must  be at least  zero.
                     74: *> \endverbatim
                     75: *>
                     76: *> \param[in] N
                     77: *> \verbatim
                     78: *>          N is INTEGER
                     79: *>           On entry,  N  specifies the number  of columns of the matrix
                     80: *>           op( B ) and the number of columns of the matrix C. N must be
                     81: *>           at least zero.
                     82: *> \endverbatim
                     83: *>
                     84: *> \param[in] K
                     85: *> \verbatim
                     86: *>          K is INTEGER
                     87: *>           On entry,  K  specifies  the number of columns of the matrix
                     88: *>           op( A ) and the number of rows of the matrix op( B ). K must
                     89: *>           be at least  zero.
                     90: *> \endverbatim
                     91: *>
                     92: *> \param[in] ALPHA
                     93: *> \verbatim
                     94: *>          ALPHA is DOUBLE PRECISION.
                     95: *>           On entry, ALPHA specifies the scalar alpha.
                     96: *> \endverbatim
                     97: *>
                     98: *> \param[in] A
                     99: *> \verbatim
1.15      bertrand  100: *>          A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
1.8       bertrand  101: *>           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
                    102: *>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
                    103: *>           part of the array  A  must contain the matrix  A,  otherwise
                    104: *>           the leading  k by m  part of the array  A  must contain  the
                    105: *>           matrix A.
                    106: *> \endverbatim
                    107: *>
                    108: *> \param[in] LDA
                    109: *> \verbatim
                    110: *>          LDA is INTEGER
                    111: *>           On entry, LDA specifies the first dimension of A as declared
                    112: *>           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
                    113: *>           LDA must be at least  max( 1, m ), otherwise  LDA must be at
                    114: *>           least  max( 1, k ).
                    115: *> \endverbatim
                    116: *>
                    117: *> \param[in] B
                    118: *> \verbatim
1.15      bertrand  119: *>          B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
1.8       bertrand  120: *>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
                    121: *>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
                    122: *>           part of the array  B  must contain the matrix  B,  otherwise
                    123: *>           the leading  n by k  part of the array  B  must contain  the
                    124: *>           matrix B.
                    125: *> \endverbatim
                    126: *>
                    127: *> \param[in] LDB
                    128: *> \verbatim
                    129: *>          LDB is INTEGER
                    130: *>           On entry, LDB specifies the first dimension of B as declared
                    131: *>           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
                    132: *>           LDB must be at least  max( 1, k ), otherwise  LDB must be at
                    133: *>           least  max( 1, n ).
                    134: *> \endverbatim
                    135: *>
                    136: *> \param[in] BETA
                    137: *> \verbatim
                    138: *>          BETA is DOUBLE PRECISION.
                    139: *>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
                    140: *>           supplied as zero then C need not be set on input.
                    141: *> \endverbatim
                    142: *>
                    143: *> \param[in,out] C
                    144: *> \verbatim
1.15      bertrand  145: *>          C is DOUBLE PRECISION array, dimension ( LDC, N )
1.8       bertrand  146: *>           Before entry, the leading  m by n  part of the array  C must
                    147: *>           contain the matrix  C,  except when  beta  is zero, in which
                    148: *>           case C need not be set on entry.
                    149: *>           On exit, the array  C  is overwritten by the  m by n  matrix
                    150: *>           ( alpha*op( A )*op( B ) + beta*C ).
                    151: *> \endverbatim
                    152: *>
                    153: *> \param[in] LDC
                    154: *> \verbatim
                    155: *>          LDC is INTEGER
                    156: *>           On entry, LDC specifies the first dimension of C as declared
                    157: *>           in  the  calling  (sub)  program.   LDC  must  be  at  least
                    158: *>           max( 1, m ).
                    159: *> \endverbatim
                    160: *
                    161: *  Authors:
                    162: *  ========
                    163: *
1.14      bertrand  164: *> \author Univ. of Tennessee
                    165: *> \author Univ. of California Berkeley
                    166: *> \author Univ. of Colorado Denver
                    167: *> \author NAG Ltd.
1.8       bertrand  168: *
                    169: *> \ingroup double_blas_level3
                    170: *
                    171: *> \par Further Details:
                    172: *  =====================
                    173: *>
                    174: *> \verbatim
                    175: *>
                    176: *>  Level 3 Blas routine.
                    177: *>
                    178: *>  -- Written on 8-February-1989.
                    179: *>     Jack Dongarra, Argonne National Laboratory.
                    180: *>     Iain Duff, AERE Harwell.
                    181: *>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
                    182: *>     Sven Hammarling, Numerical Algorithms Group Ltd.
                    183: *> \endverbatim
                    184: *>
                    185: *  =====================================================================
1.1       bertrand  186:       SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
1.8       bertrand  187: *
1.17    ! bertrand  188: *  -- Reference BLAS level3 routine --
1.8       bertrand  189: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
                    190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    191: *
1.1       bertrand  192: *     .. Scalar Arguments ..
                    193:       DOUBLE PRECISION ALPHA,BETA
                    194:       INTEGER K,LDA,LDB,LDC,M,N
                    195:       CHARACTER TRANSA,TRANSB
                    196: *     ..
                    197: *     .. Array Arguments ..
                    198:       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
                    199: *     ..
                    200: *
                    201: *  =====================================================================
                    202: *
                    203: *     .. External Functions ..
                    204:       LOGICAL LSAME
                    205:       EXTERNAL LSAME
                    206: *     ..
                    207: *     .. External Subroutines ..
                    208:       EXTERNAL XERBLA
                    209: *     ..
                    210: *     .. Intrinsic Functions ..
                    211:       INTRINSIC MAX
                    212: *     ..
                    213: *     .. Local Scalars ..
                    214:       DOUBLE PRECISION TEMP
1.17    ! bertrand  215:       INTEGER I,INFO,J,L,NROWA,NROWB
1.1       bertrand  216:       LOGICAL NOTA,NOTB
                    217: *     ..
                    218: *     .. Parameters ..
                    219:       DOUBLE PRECISION ONE,ZERO
                    220:       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
                    221: *     ..
                    222: *
                    223: *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
1.17    ! bertrand  224: *     transposed and set  NROWA and NROWB  as the number of rows of  A
        !           225: *     and  B  respectively.
1.1       bertrand  226: *
                    227:       NOTA = LSAME(TRANSA,'N')
                    228:       NOTB = LSAME(TRANSB,'N')
                    229:       IF (NOTA) THEN
                    230:           NROWA = M
                    231:       ELSE
                    232:           NROWA = K
                    233:       END IF
                    234:       IF (NOTB) THEN
                    235:           NROWB = K
                    236:       ELSE
                    237:           NROWB = N
                    238:       END IF
                    239: *
                    240: *     Test the input parameters.
                    241: *
                    242:       INFO = 0
                    243:       IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
                    244:      +    (.NOT.LSAME(TRANSA,'T'))) THEN
                    245:           INFO = 1
                    246:       ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
                    247:      +         (.NOT.LSAME(TRANSB,'T'))) THEN
                    248:           INFO = 2
                    249:       ELSE IF (M.LT.0) THEN
                    250:           INFO = 3
                    251:       ELSE IF (N.LT.0) THEN
                    252:           INFO = 4
                    253:       ELSE IF (K.LT.0) THEN
                    254:           INFO = 5
                    255:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
                    256:           INFO = 8
                    257:       ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
                    258:           INFO = 10
                    259:       ELSE IF (LDC.LT.MAX(1,M)) THEN
                    260:           INFO = 13
                    261:       END IF
                    262:       IF (INFO.NE.0) THEN
                    263:           CALL XERBLA('DGEMM ',INFO)
                    264:           RETURN
                    265:       END IF
                    266: *
                    267: *     Quick return if possible.
                    268: *
                    269:       IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
                    270:      +    (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
                    271: *
                    272: *     And if  alpha.eq.zero.
                    273: *
                    274:       IF (ALPHA.EQ.ZERO) THEN
                    275:           IF (BETA.EQ.ZERO) THEN
                    276:               DO 20 J = 1,N
                    277:                   DO 10 I = 1,M
                    278:                       C(I,J) = ZERO
                    279:    10             CONTINUE
                    280:    20         CONTINUE
                    281:           ELSE
                    282:               DO 40 J = 1,N
                    283:                   DO 30 I = 1,M
                    284:                       C(I,J) = BETA*C(I,J)
                    285:    30             CONTINUE
                    286:    40         CONTINUE
                    287:           END IF
                    288:           RETURN
                    289:       END IF
                    290: *
                    291: *     Start the operations.
                    292: *
                    293:       IF (NOTB) THEN
                    294:           IF (NOTA) THEN
                    295: *
                    296: *           Form  C := alpha*A*B + beta*C.
                    297: *
                    298:               DO 90 J = 1,N
                    299:                   IF (BETA.EQ.ZERO) THEN
                    300:                       DO 50 I = 1,M
                    301:                           C(I,J) = ZERO
                    302:    50                 CONTINUE
                    303:                   ELSE IF (BETA.NE.ONE) THEN
                    304:                       DO 60 I = 1,M
                    305:                           C(I,J) = BETA*C(I,J)
                    306:    60                 CONTINUE
                    307:                   END IF
                    308:                   DO 80 L = 1,K
1.12      bertrand  309:                       TEMP = ALPHA*B(L,J)
                    310:                       DO 70 I = 1,M
                    311:                           C(I,J) = C(I,J) + TEMP*A(I,L)
                    312:    70                 CONTINUE
1.1       bertrand  313:    80             CONTINUE
                    314:    90         CONTINUE
                    315:           ELSE
                    316: *
1.7       bertrand  317: *           Form  C := alpha*A**T*B + beta*C
1.1       bertrand  318: *
                    319:               DO 120 J = 1,N
                    320:                   DO 110 I = 1,M
                    321:                       TEMP = ZERO
                    322:                       DO 100 L = 1,K
                    323:                           TEMP = TEMP + A(L,I)*B(L,J)
                    324:   100                 CONTINUE
                    325:                       IF (BETA.EQ.ZERO) THEN
                    326:                           C(I,J) = ALPHA*TEMP
                    327:                       ELSE
                    328:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
                    329:                       END IF
                    330:   110             CONTINUE
                    331:   120         CONTINUE
                    332:           END IF
                    333:       ELSE
                    334:           IF (NOTA) THEN
                    335: *
1.7       bertrand  336: *           Form  C := alpha*A*B**T + beta*C
1.1       bertrand  337: *
                    338:               DO 170 J = 1,N
                    339:                   IF (BETA.EQ.ZERO) THEN
                    340:                       DO 130 I = 1,M
                    341:                           C(I,J) = ZERO
                    342:   130                 CONTINUE
                    343:                   ELSE IF (BETA.NE.ONE) THEN
                    344:                       DO 140 I = 1,M
                    345:                           C(I,J) = BETA*C(I,J)
                    346:   140                 CONTINUE
                    347:                   END IF
                    348:                   DO 160 L = 1,K
1.12      bertrand  349:                       TEMP = ALPHA*B(J,L)
                    350:                       DO 150 I = 1,M
                    351:                           C(I,J) = C(I,J) + TEMP*A(I,L)
                    352:   150                 CONTINUE
1.1       bertrand  353:   160             CONTINUE
                    354:   170         CONTINUE
                    355:           ELSE
                    356: *
1.7       bertrand  357: *           Form  C := alpha*A**T*B**T + beta*C
1.1       bertrand  358: *
                    359:               DO 200 J = 1,N
                    360:                   DO 190 I = 1,M
                    361:                       TEMP = ZERO
                    362:                       DO 180 L = 1,K
                    363:                           TEMP = TEMP + A(L,I)*B(J,L)
                    364:   180                 CONTINUE
                    365:                       IF (BETA.EQ.ZERO) THEN
                    366:                           C(I,J) = ALPHA*TEMP
                    367:                       ELSE
                    368:                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
                    369:                       END IF
                    370:   190             CONTINUE
                    371:   200         CONTINUE
                    372:           END IF
                    373:       END IF
                    374: *
                    375:       RETURN
                    376: *
1.17    ! bertrand  377: *     End of DGEMM
1.1       bertrand  378: *
                    379:       END

CVSweb interface <joel.bertrand@systella.fr>