1*> \brief \b DGEMM 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 12* 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* .. 21* 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 100*> A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 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 119*> B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 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 145*> C is DOUBLE PRECISION array of DIMENSION ( LDC, n ). 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* 164*> \author Univ. of Tennessee 165*> \author Univ. of California Berkeley 166*> \author Univ. of Colorado Denver 167*> \author NAG Ltd. 168* 169*> \date November 2015 170* 171*> \ingroup double_blas_level3 172* 173*> \par Further Details: 174* ===================== 175*> 176*> \verbatim 177*> 178*> Level 3 Blas routine. 179*> 180*> -- Written on 8-February-1989. 181*> Jack Dongarra, Argonne National Laboratory. 182*> Iain Duff, AERE Harwell. 183*> Jeremy Du Croz, Numerical Algorithms Group Ltd. 184*> Sven Hammarling, Numerical Algorithms Group Ltd. 185*> \endverbatim 186*> 187* ===================================================================== 188 SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 189* 190* -- Reference BLAS level3 routine (version 3.6.0) -- 191* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- 192* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 193* November 2015 194* 195* .. Scalar Arguments .. 196 DOUBLE PRECISION ALPHA,BETA 197 INTEGER K,LDA,LDB,LDC,M,N 198 CHARACTER TRANSA,TRANSB 199* .. 200* .. Array Arguments .. 201 DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) 202* .. 203* 204* ===================================================================== 205* 206* .. External Functions .. 207 LOGICAL LSAME 208 EXTERNAL LSAME 209* .. 210* .. External Subroutines .. 211 EXTERNAL XERBLA 212* .. 213* .. Intrinsic Functions .. 214 INTRINSIC MAX 215* .. 216* .. Local Scalars .. 217 DOUBLE PRECISION TEMP 218 INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB 219 LOGICAL NOTA,NOTB 220* .. 221* .. Parameters .. 222 DOUBLE PRECISION ONE,ZERO 223 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) 224* .. 225* 226* Set NOTA and NOTB as true if A and B respectively are not 227* transposed and set NROWA, NCOLA and NROWB as the number of rows 228* and columns of A and the number of rows of B respectively. 229* 230 NOTA = LSAME(TRANSA,'N') 231 NOTB = LSAME(TRANSB,'N') 232 IF (NOTA) THEN 233 NROWA = M 234 NCOLA = K 235 ELSE 236 NROWA = K 237 NCOLA = M 238 END IF 239 IF (NOTB) THEN 240 NROWB = K 241 ELSE 242 NROWB = N 243 END IF 244* 245* Test the input parameters. 246* 247 INFO = 0 248 IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND. 249 + (.NOT.LSAME(TRANSA,'T'))) THEN 250 INFO = 1 251 ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND. 252 + (.NOT.LSAME(TRANSB,'T'))) THEN 253 INFO = 2 254 ELSE IF (M.LT.0) THEN 255 INFO = 3 256 ELSE IF (N.LT.0) THEN 257 INFO = 4 258 ELSE IF (K.LT.0) THEN 259 INFO = 5 260 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 261 INFO = 8 262 ELSE IF (LDB.LT.MAX(1,NROWB)) THEN 263 INFO = 10 264 ELSE IF (LDC.LT.MAX(1,M)) THEN 265 INFO = 13 266 END IF 267 IF (INFO.NE.0) THEN 268 CALL XERBLA('DGEMM ',INFO) 269 RETURN 270 END IF 271* 272* Quick return if possible. 273* 274 IF ((M.EQ.0) .OR. (N.EQ.0) .OR. 275 + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN 276* 277* And if alpha.eq.zero. 278* 279 IF (ALPHA.EQ.ZERO) THEN 280 IF (BETA.EQ.ZERO) THEN 281 DO 20 J = 1,N 282 DO 10 I = 1,M 283 C(I,J) = ZERO 284 10 CONTINUE 285 20 CONTINUE 286 ELSE 287 DO 40 J = 1,N 288 DO 30 I = 1,M 289 C(I,J) = BETA*C(I,J) 290 30 CONTINUE 291 40 CONTINUE 292 END IF 293 RETURN 294 END IF 295* 296* Start the operations. 297* 298 IF (NOTB) THEN 299 IF (NOTA) THEN 300* 301* Form C := alpha*A*B + beta*C. 302* 303 DO 90 J = 1,N 304 IF (BETA.EQ.ZERO) THEN 305 DO 50 I = 1,M 306 C(I,J) = ZERO 307 50 CONTINUE 308 ELSE IF (BETA.NE.ONE) THEN 309 DO 60 I = 1,M 310 C(I,J) = BETA*C(I,J) 311 60 CONTINUE 312 END IF 313 DO 80 L = 1,K 314 TEMP = ALPHA*B(L,J) 315 DO 70 I = 1,M 316 C(I,J) = C(I,J) + TEMP*A(I,L) 317 70 CONTINUE 318 80 CONTINUE 319 90 CONTINUE 320 ELSE 321* 322* Form C := alpha*A**T*B + beta*C 323* 324 DO 120 J = 1,N 325 DO 110 I = 1,M 326 TEMP = ZERO 327 DO 100 L = 1,K 328 TEMP = TEMP + A(L,I)*B(L,J) 329 100 CONTINUE 330 IF (BETA.EQ.ZERO) THEN 331 C(I,J) = ALPHA*TEMP 332 ELSE 333 C(I,J) = ALPHA*TEMP + BETA*C(I,J) 334 END IF 335 110 CONTINUE 336 120 CONTINUE 337 END IF 338 ELSE 339 IF (NOTA) THEN 340* 341* Form C := alpha*A*B**T + beta*C 342* 343 DO 170 J = 1,N 344 IF (BETA.EQ.ZERO) THEN 345 DO 130 I = 1,M 346 C(I,J) = ZERO 347 130 CONTINUE 348 ELSE IF (BETA.NE.ONE) THEN 349 DO 140 I = 1,M 350 C(I,J) = BETA*C(I,J) 351 140 CONTINUE 352 END IF 353 DO 160 L = 1,K 354 TEMP = ALPHA*B(J,L) 355 DO 150 I = 1,M 356 C(I,J) = C(I,J) + TEMP*A(I,L) 357 150 CONTINUE 358 160 CONTINUE 359 170 CONTINUE 360 ELSE 361* 362* Form C := alpha*A**T*B**T + beta*C 363* 364 DO 200 J = 1,N 365 DO 190 I = 1,M 366 TEMP = ZERO 367 DO 180 L = 1,K 368 TEMP = TEMP + A(L,I)*B(J,L) 369 180 CONTINUE 370 IF (BETA.EQ.ZERO) THEN 371 C(I,J) = ALPHA*TEMP 372 ELSE 373 C(I,J) = ALPHA*TEMP + BETA*C(I,J) 374 END IF 375 190 CONTINUE 376 200 CONTINUE 377 END IF 378 END IF 379* 380 RETURN 381* 382* End of DGEMM . 383* 384 END 385