001: SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 002: * .. Scalar Arguments .. 003: DOUBLE PRECISION ALPHA,BETA 004: INTEGER K,LDA,LDB,LDC,M,N 005: CHARACTER TRANSA,TRANSB 006: * .. 007: * .. Array Arguments .. 008: DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) 009: * .. 010: * 011: * Purpose 012: * ======= 013: * 014: * DGEMM performs one of the matrix-matrix operations 015: * 016: * C := alpha*op( A )*op( B ) + beta*C, 017: * 018: * where op( X ) is one of 019: * 020: * op( X ) = X or op( X ) = X', 021: * 022: * alpha and beta are scalars, and A, B and C are matrices, with op( A ) 023: * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. 024: * 025: * Arguments 026: * ========== 027: * 028: * TRANSA - CHARACTER*1. 029: * On entry, TRANSA specifies the form of op( A ) to be used in 030: * the matrix multiplication as follows: 031: * 032: * TRANSA = 'N' or 'n', op( A ) = A. 033: * 034: * TRANSA = 'T' or 't', op( A ) = A'. 035: * 036: * TRANSA = 'C' or 'c', op( A ) = A'. 037: * 038: * Unchanged on exit. 039: * 040: * TRANSB - CHARACTER*1. 041: * On entry, TRANSB specifies the form of op( B ) to be used in 042: * the matrix multiplication as follows: 043: * 044: * TRANSB = 'N' or 'n', op( B ) = B. 045: * 046: * TRANSB = 'T' or 't', op( B ) = B'. 047: * 048: * TRANSB = 'C' or 'c', op( B ) = B'. 049: * 050: * Unchanged on exit. 051: * 052: * M - INTEGER. 053: * On entry, M specifies the number of rows of the matrix 054: * op( A ) and of the matrix C. M must be at least zero. 055: * Unchanged on exit. 056: * 057: * N - INTEGER. 058: * On entry, N specifies the number of columns of the matrix 059: * op( B ) and the number of columns of the matrix C. N must be 060: * at least zero. 061: * Unchanged on exit. 062: * 063: * K - INTEGER. 064: * On entry, K specifies the number of columns of the matrix 065: * op( A ) and the number of rows of the matrix op( B ). K must 066: * be at least zero. 067: * Unchanged on exit. 068: * 069: * ALPHA - DOUBLE PRECISION. 070: * On entry, ALPHA specifies the scalar alpha. 071: * Unchanged on exit. 072: * 073: * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 074: * k when TRANSA = 'N' or 'n', and is m otherwise. 075: * Before entry with TRANSA = 'N' or 'n', the leading m by k 076: * part of the array A must contain the matrix A, otherwise 077: * the leading k by m part of the array A must contain the 078: * matrix A. 079: * Unchanged on exit. 080: * 081: * LDA - INTEGER. 082: * On entry, LDA specifies the first dimension of A as declared 083: * in the calling (sub) program. When TRANSA = 'N' or 'n' then 084: * LDA must be at least max( 1, m ), otherwise LDA must be at 085: * least max( 1, k ). 086: * Unchanged on exit. 087: * 088: * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 089: * n when TRANSB = 'N' or 'n', and is k otherwise. 090: * Before entry with TRANSB = 'N' or 'n', the leading k by n 091: * part of the array B must contain the matrix B, otherwise 092: * the leading n by k part of the array B must contain the 093: * matrix B. 094: * Unchanged on exit. 095: * 096: * LDB - INTEGER. 097: * On entry, LDB specifies the first dimension of B as declared 098: * in the calling (sub) program. When TRANSB = 'N' or 'n' then 099: * LDB must be at least max( 1, k ), otherwise LDB must be at 100: * least max( 1, n ). 101: * Unchanged on exit. 102: * 103: * BETA - DOUBLE PRECISION. 104: * On entry, BETA specifies the scalar beta. When BETA is 105: * supplied as zero then C need not be set on input. 106: * Unchanged on exit. 107: * 108: * C - DOUBLE PRECISION 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 ). 114: * 115: * LDC - INTEGER. 116: * On entry, LDC specifies the first dimension of C as declared 117: * in the calling (sub) program. LDC must be at least 118: * max( 1, m ). 119: * Unchanged on exit. 120: * 121: * Further Details 122: * =============== 123: * 124: * Level 3 Blas routine. 125: * 126: * -- Written on 8-February-1989. 127: * Jack Dongarra, Argonne National Laboratory. 128: * Iain Duff, AERE Harwell. 129: * Jeremy Du Croz, Numerical Algorithms Group Ltd. 130: * Sven Hammarling, Numerical Algorithms Group Ltd. 131: * 132: * ===================================================================== 133: * 134: * .. External Functions .. 135: LOGICAL LSAME 136: EXTERNAL LSAME 137: * .. 138: * .. External Subroutines .. 139: EXTERNAL XERBLA 140: * .. 141: * .. Intrinsic Functions .. 142: INTRINSIC MAX 143: * .. 144: * .. Local Scalars .. 145: DOUBLE PRECISION TEMP 146: INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB 147: LOGICAL NOTA,NOTB 148: * .. 149: * .. Parameters .. 150: DOUBLE PRECISION ONE,ZERO 151: PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) 152: * .. 153: * 154: * Set NOTA and NOTB as true if A and B respectively are not 155: * transposed and set NROWA, NCOLA and NROWB as the number of rows 156: * and columns of A and the number of rows of B respectively. 157: * 158: NOTA = LSAME(TRANSA,'N') 159: NOTB = LSAME(TRANSB,'N') 160: IF (NOTA) THEN 161: NROWA = M 162: NCOLA = K 163: ELSE 164: NROWA = K 165: NCOLA = M 166: END IF 167: IF (NOTB) THEN 168: NROWB = K 169: ELSE 170: NROWB = N 171: END IF 172: * 173: * Test the input parameters. 174: * 175: INFO = 0 176: IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND. 177: + (.NOT.LSAME(TRANSA,'T'))) THEN 178: INFO = 1 179: ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND. 180: + (.NOT.LSAME(TRANSB,'T'))) THEN 181: INFO = 2 182: ELSE IF (M.LT.0) THEN 183: INFO = 3 184: ELSE IF (N.LT.0) THEN 185: INFO = 4 186: ELSE IF (K.LT.0) THEN 187: INFO = 5 188: ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 189: INFO = 8 190: ELSE IF (LDB.LT.MAX(1,NROWB)) THEN 191: INFO = 10 192: ELSE IF (LDC.LT.MAX(1,M)) THEN 193: INFO = 13 194: END IF 195: IF (INFO.NE.0) THEN 196: CALL XERBLA('DGEMM ',INFO) 197: RETURN 198: END IF 199: * 200: * Quick return if possible. 201: * 202: IF ((M.EQ.0) .OR. (N.EQ.0) .OR. 203: + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN 204: * 205: * And if alpha.eq.zero. 206: * 207: IF (ALPHA.EQ.ZERO) THEN 208: IF (BETA.EQ.ZERO) THEN 209: DO 20 J = 1,N 210: DO 10 I = 1,M 211: C(I,J) = ZERO 212: 10 CONTINUE 213: 20 CONTINUE 214: ELSE 215: DO 40 J = 1,N 216: DO 30 I = 1,M 217: C(I,J) = BETA*C(I,J) 218: 30 CONTINUE 219: 40 CONTINUE 220: END IF 221: RETURN 222: END IF 223: * 224: * Start the operations. 225: * 226: IF (NOTB) THEN 227: IF (NOTA) THEN 228: * 229: * Form C := alpha*A*B + beta*C. 230: * 231: DO 90 J = 1,N 232: IF (BETA.EQ.ZERO) THEN 233: DO 50 I = 1,M 234: C(I,J) = ZERO 235: 50 CONTINUE 236: ELSE IF (BETA.NE.ONE) THEN 237: DO 60 I = 1,M 238: C(I,J) = BETA*C(I,J) 239: 60 CONTINUE 240: END IF 241: DO 80 L = 1,K 242: IF (B(L,J).NE.ZERO) THEN 243: TEMP = ALPHA*B(L,J) 244: DO 70 I = 1,M 245: C(I,J) = C(I,J) + TEMP*A(I,L) 246: 70 CONTINUE 247: END IF 248: 80 CONTINUE 249: 90 CONTINUE 250: ELSE 251: * 252: * Form C := alpha*A'*B + beta*C 253: * 254: DO 120 J = 1,N 255: DO 110 I = 1,M 256: TEMP = ZERO 257: DO 100 L = 1,K 258: TEMP = TEMP + A(L,I)*B(L,J) 259: 100 CONTINUE 260: IF (BETA.EQ.ZERO) THEN 261: C(I,J) = ALPHA*TEMP 262: ELSE 263: C(I,J) = ALPHA*TEMP + BETA*C(I,J) 264: END IF 265: 110 CONTINUE 266: 120 CONTINUE 267: END IF 268: ELSE 269: IF (NOTA) THEN 270: * 271: * Form C := alpha*A*B' + beta*C 272: * 273: DO 170 J = 1,N 274: IF (BETA.EQ.ZERO) THEN 275: DO 130 I = 1,M 276: C(I,J) = ZERO 277: 130 CONTINUE 278: ELSE IF (BETA.NE.ONE) THEN 279: DO 140 I = 1,M 280: C(I,J) = BETA*C(I,J) 281: 140 CONTINUE 282: END IF 283: DO 160 L = 1,K 284: IF (B(J,L).NE.ZERO) THEN 285: TEMP = ALPHA*B(J,L) 286: DO 150 I = 1,M 287: C(I,J) = C(I,J) + TEMP*A(I,L) 288: 150 CONTINUE 289: END IF 290: 160 CONTINUE 291: 170 CONTINUE 292: ELSE 293: * 294: * Form C := alpha*A'*B' + beta*C 295: * 296: DO 200 J = 1,N 297: DO 190 I = 1,M 298: TEMP = ZERO 299: DO 180 L = 1,K 300: TEMP = TEMP + A(L,I)*B(J,L) 301: 180 CONTINUE 302: IF (BETA.EQ.ZERO) THEN 303: C(I,J) = ALPHA*TEMP 304: ELSE 305: C(I,J) = ALPHA*TEMP + BETA*C(I,J) 306: END IF 307: 190 CONTINUE 308: 200 CONTINUE 309: END IF 310: END IF 311: * 312: RETURN 313: * 314: * End of DGEMM . 315: * 316: END 317: