github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/netlib/dgemm.f (about) 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