| @@ -0,0 +1,395 @@ |
| 1 | + |
| 2 | +/* -- translated by f2c (version 19940927). |
| 3 | + You must link the resulting object file with the libraries: |
| 4 | + -lf2c -lm (in that order) |
| 5 | +*/ |
| 6 | + |
| 7 | +#include "f2c.h" |
| 8 | + |
| 9 | +/* Subroutine */ int dgemm_(char *transa, char *transb, integer *m, integer * |
| 10 | + n, integer *k, doublereal *alpha, doublereal *a, integer *lda, |
| 11 | + doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer |
| 12 | + *ldc) |
| 13 | +{ |
| 14 | + |
| 15 | + |
| 16 | + /* System generated locals */ |
| 17 | + integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, |
| 18 | + i__3; |
| 19 | + |
| 20 | + /* Local variables */ |
| 21 | + static integer info; |
| 22 | + static logical nota, notb; |
| 23 | + static doublereal temp; |
| 24 | + static integer i, j, l, ncola; |
| 25 | + extern logical lsame_(char *, char *); |
| 26 | + static integer nrowa, nrowb; |
| 27 | + extern /* Subroutine */ int xerbla_(char *, integer *); |
| 28 | + |
| 29 | + |
| 30 | +/* Purpose |
| 31 | + ======= |
| 32 | + |
| 33 | + DGEMM performs one of the matrix-matrix operations |
| 34 | + |
| 35 | + C := alpha*op( A )*op( B ) + beta*C, |
| 36 | + |
| 37 | + where op( X ) is one of |
| 38 | + |
| 39 | + op( X ) = X or op( X ) = X', |
| 40 | + |
| 41 | + alpha and beta are scalars, and A, B and C are matrices, with op( A ) |
| 42 | + |
| 43 | + an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. |
| 44 | + |
| 45 | + |
| 46 | + Parameters |
| 47 | + ========== |
| 48 | + |
| 49 | + TRANSA - CHARACTER*1. |
| 50 | + On entry, TRANSA specifies the form of op( A ) to be used in |
| 51 | + |
| 52 | + the matrix multiplication as follows: |
| 53 | + |
| 54 | + TRANSA = 'N' or 'n', op( A ) = A. |
| 55 | + |
| 56 | + TRANSA = 'T' or 't', op( A ) = A'. |
| 57 | + |
| 58 | + TRANSA = 'C' or 'c', op( A ) = A'. |
| 59 | + |
| 60 | + Unchanged on exit. |
| 61 | + |
| 62 | + TRANSB - CHARACTER*1. |
| 63 | + On entry, TRANSB specifies the form of op( B ) to be used in |
| 64 | + |
| 65 | + the matrix multiplication as follows: |
| 66 | + |
| 67 | + TRANSB = 'N' or 'n', op( B ) = B. |
| 68 | + |
| 69 | + TRANSB = 'T' or 't', op( B ) = B'. |
| 70 | + |
| 71 | + TRANSB = 'C' or 'c', op( B ) = B'. |
| 72 | + |
| 73 | + Unchanged on exit. |
| 74 | + |
| 75 | + M - INTEGER. |
| 76 | + On entry, M specifies the number of rows of the matrix |
| 77 | + |
| 78 | + op( A ) and of the matrix C. M must be at least zero. |
| 79 | + |
| 80 | + Unchanged on exit. |
| 81 | + |
| 82 | + N - INTEGER. |
| 83 | + On entry, N specifies the number of columns of the matrix |
| 84 | + |
| 85 | + op( B ) and the number of columns of the matrix C. N must be |
| 86 | + |
| 87 | + at least zero. |
| 88 | + Unchanged on exit. |
| 89 | + |
| 90 | + K - INTEGER. |
| 91 | + On entry, K specifies the number of columns of the matrix |
| 92 | + |
| 93 | + op( A ) and the number of rows of the matrix op( B ). K must |
| 94 | + |
| 95 | + be at least zero. |
| 96 | + Unchanged on exit. |
| 97 | + |
| 98 | + ALPHA - DOUBLE PRECISION. |
| 99 | + On entry, ALPHA specifies the scalar alpha. |
| 100 | + Unchanged on exit. |
| 101 | + |
| 102 | + A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is |
| 103 | + |
| 104 | + k when TRANSA = 'N' or 'n', and is m otherwise. |
| 105 | + Before entry with TRANSA = 'N' or 'n', the leading m by k |
| 106 | + |
| 107 | + part of the array A must contain the matrix A, otherwise |
| 108 | + |
| 109 | + the leading k by m part of the array A must contain the |
| 110 | + |
| 111 | + matrix A. |
| 112 | + Unchanged on exit. |
| 113 | + |
| 114 | + LDA - INTEGER. |
| 115 | + On entry, LDA specifies the first dimension of A as declared |
| 116 | + |
| 117 | + in the calling (sub) program. When TRANSA = 'N' or 'n' then |
| 118 | + |
| 119 | + LDA must be at least max( 1, m ), otherwise LDA must be at |
| 120 | + |
| 121 | + least max( 1, k ). |
| 122 | + Unchanged on exit. |
| 123 | + |
| 124 | + B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is |
| 125 | + |
| 126 | + n when TRANSB = 'N' or 'n', and is k otherwise. |
| 127 | + Before entry with TRANSB = 'N' or 'n', the leading k by n |
| 128 | + |
| 129 | + part of the array B must contain the matrix B, otherwise |
| 130 | + |
| 131 | + the leading n by k part of the array B must contain the |
| 132 | + |
| 133 | + matrix B. |
| 134 | + Unchanged on exit. |
| 135 | + |
| 136 | + LDB - INTEGER. |
| 137 | + On entry, LDB specifies the first dimension of B as declared |
| 138 | + |
| 139 | + in the calling (sub) program. When TRANSB = 'N' or 'n' then |
| 140 | + |
| 141 | + LDB must be at least max( 1, k ), otherwise LDB must be at |
| 142 | + |
| 143 | + least max( 1, n ). |
| 144 | + Unchanged on exit. |
| 145 | + |
| 146 | + BETA - DOUBLE PRECISION. |
| 147 | + On entry, BETA specifies the scalar beta. When BETA is |
| 148 | + |
| 149 | + supplied as zero then C need not be set on input. |
| 150 | + Unchanged on exit. |
| 151 | + |
| 152 | + C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). |
| 153 | + Before entry, the leading m by n part of the array C must |
| 154 | + |
| 155 | + contain the matrix C, except when beta is zero, in which |
| 156 | + |
| 157 | + case C need not be set on entry. |
| 158 | + On exit, the array C is overwritten by the m by n matrix |
| 159 | + |
| 160 | + ( alpha*op( A )*op( B ) + beta*C ). |
| 161 | + |
| 162 | + LDC - INTEGER. |
| 163 | + On entry, LDC specifies the first dimension of C as declared |
| 164 | + |
| 165 | + in the calling (sub) program. LDC must be at least |
| 166 | + |
| 167 | + max( 1, m ). |
| 168 | + Unchanged on exit. |
| 169 | + |
| 170 | + |
| 171 | + Level 3 Blas routine. |
| 172 | + |
| 173 | + -- Written on 8-February-1989. |
| 174 | + Jack Dongarra, Argonne National Laboratory. |
| 175 | + Iain Duff, AERE Harwell. |
| 176 | + Jeremy Du Croz, Numerical Algorithms Group Ltd. |
| 177 | + Sven Hammarling, Numerical Algorithms Group Ltd. |
| 178 | + |
| 179 | + |
| 180 | + |
| 181 | + Set NOTA and NOTB as true if A and B respectively are not |
| 182 | + |
| 183 | + transposed and set NROWA, NCOLA and NROWB as the number of rows |
| 184 | + |
| 185 | + and columns of A and the number of rows of B respectively. |
| 186 | + |
| 187 | + |
| 188 | + |
| 189 | + Parameter adjustments |
| 190 | + Function Body */ |
| 191 | + |
| 192 | +#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)] |
| 193 | +#define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)] |
| 194 | +#define C(I,J) c[(I)-1 + ((J)-1)* ( *ldc)] |
| 195 | + |
| 196 | + nota = lsame_(transa, "N"); |
| 197 | + notb = lsame_(transb, "N"); |
| 198 | + if (nota) { |
| 199 | + nrowa = *m; |
| 200 | + ncola = *k; |
| 201 | + } else { |
| 202 | + nrowa = *k; |
| 203 | + ncola = *m; |
| 204 | + } |
| 205 | + if (notb) { |
| 206 | + nrowb = *k; |
| 207 | + } else { |
| 208 | + nrowb = *n; |
| 209 | + } |
| 210 | + |
| 211 | +/* Test the input parameters. */ |
| 212 | + |
| 213 | + info = 0; |
| 214 | + if (! nota && ! lsame_(transa, "C") && ! lsame_(transa, "T")) { |
| 215 | + info = 1; |
| 216 | + } else if (! notb && ! lsame_(transb, "C") && ! lsame_(transb, |
| 217 | + "T")) { |
| 218 | + info = 2; |
| 219 | + } else if (*m < 0) { |
| 220 | + info = 3; |
| 221 | + } else if (*n < 0) { |
| 222 | + info = 4; |
| 223 | + } else if (*k < 0) { |
| 224 | + info = 5; |
| 225 | + } else if (*lda < max(1,nrowa)) { |
| 226 | + info = 8; |
| 227 | + } else if (*ldb < max(1,nrowb)) { |
| 228 | + info = 10; |
| 229 | + } else if (*ldc < max(1,*m)) { |
| 230 | + info = 13; |
| 231 | + } |
| 232 | + if (info != 0) { |
| 233 | + xerbla_("DGEMM ", &info); |
| 234 | + return 0; |
| 235 | + } |
| 236 | + |
| 237 | +/* Quick return if possible. */ |
| 238 | + |
| 239 | + if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) { |
| 240 | + return 0; |
| 241 | + } |
| 242 | + |
| 243 | +/* And if alpha.eq.zero. */ |
| 244 | + |
| 245 | + if (*alpha == 0.) { |
| 246 | + if (*beta == 0.) { |
| 247 | + i__1 = *n; |
| 248 | + for (j = 1; j <= *n; ++j) { |
| 249 | + i__2 = *m; |
| 250 | + for (i = 1; i <= *m; ++i) { |
| 251 | + C(i,j) = 0.; |
| 252 | +/* L10: */ |
| 253 | + } |
| 254 | +/* L20: */ |
| 255 | + } |
| 256 | + } else { |
| 257 | + i__1 = *n; |
| 258 | + for (j = 1; j <= *n; ++j) { |
| 259 | + i__2 = *m; |
| 260 | + for (i = 1; i <= *m; ++i) { |
| 261 | + C(i,j) = *beta * C(i,j); |
| 262 | +/* L30: */ |
| 263 | + } |
| 264 | +/* L40: */ |
| 265 | + } |
| 266 | + } |
| 267 | + return 0; |
| 268 | + } |
| 269 | + |
| 270 | +/* Start the operations. */ |
| 271 | + |
| 272 | + if (notb) { |
| 273 | + if (nota) { |
| 274 | + |
| 275 | +/* Form C := alpha*A*B + beta*C. */ |
| 276 | + |
| 277 | + i__1 = *n; |
| 278 | + for (j = 1; j <= *n; ++j) { |
| 279 | + if (*beta == 0.) { |
| 280 | + i__2 = *m; |
| 281 | + for (i = 1; i <= *m; ++i) { |
| 282 | + C(i,j) = 0.; |
| 283 | +/* L50: */ |
| 284 | + } |
| 285 | + } else if (*beta != 1.) { |
| 286 | + i__2 = *m; |
| 287 | + for (i = 1; i <= *m; ++i) { |
| 288 | + C(i,j) = *beta * C(i,j); |
| 289 | +/* L60: */ |
| 290 | + } |
| 291 | + } |
| 292 | + i__2 = *k; |
| 293 | + for (l = 1; l <= *k; ++l) { |
| 294 | + if (B(l,j) != 0.) { |
| 295 | + temp = *alpha * B(l,j); |
| 296 | + i__3 = *m; |
| 297 | + for (i = 1; i <= *m; ++i) { |
| 298 | + C(i,j) += temp * A(i,l); |
| 299 | +/* L70: */ |
| 300 | + } |
| 301 | + } |
| 302 | +/* L80: */ |
| 303 | + } |
| 304 | +/* L90: */ |
| 305 | + } |
| 306 | + } else { |
| 307 | + |
| 308 | +/* Form C := alpha*A'*B + beta*C */ |
| 309 | + |
| 310 | + i__1 = *n; |
| 311 | + for (j = 1; j <= *n; ++j) { |
| 312 | + i__2 = *m; |
| 313 | + for (i = 1; i <= *m; ++i) { |
| 314 | + temp = 0.; |
| 315 | + i__3 = *k; |
| 316 | + for (l = 1; l <= *k; ++l) { |
| 317 | + temp += A(l,i) * B(l,j); |
| 318 | +/* L100: */ |
| 319 | + } |
| 320 | + if (*beta == 0.) { |
| 321 | + C(i,j) = *alpha * temp; |
| 322 | + } else { |
| 323 | + C(i,j) = *alpha * temp + *beta * C(i,j); |
| 324 | + } |
| 325 | +/* L110: */ |
| 326 | + } |
| 327 | +/* L120: */ |
| 328 | + } |
| 329 | + } |
| 330 | + } else { |
| 331 | + if (nota) { |
| 332 | + |
| 333 | +/* Form C := alpha*A*B' + beta*C */ |
| 334 | + |
| 335 | + i__1 = *n; |
| 336 | + for (j = 1; j <= *n; ++j) { |
| 337 | + if (*beta == 0.) { |
| 338 | + i__2 = *m; |
| 339 | + for (i = 1; i <= *m; ++i) { |
| 340 | + C(i,j) = 0.; |
| 341 | +/* L130: */ |
| 342 | + } |
| 343 | + } else if (*beta != 1.) { |
| 344 | + i__2 = *m; |
| 345 | + for (i = 1; i <= *m; ++i) { |
| 346 | + C(i,j) = *beta * C(i,j); |
| 347 | +/* L140: */ |
| 348 | + } |
| 349 | + } |
| 350 | + i__2 = *k; |
| 351 | + for (l = 1; l <= *k; ++l) { |
| 352 | + if (B(j,l) != 0.) { |
| 353 | + temp = *alpha * B(j,l); |
| 354 | + i__3 = *m; |
| 355 | + for (i = 1; i <= *m; ++i) { |
| 356 | + C(i,j) += temp * A(i,l); |
| 357 | +/* L150: */ |
| 358 | + } |
| 359 | + } |
| 360 | +/* L160: */ |
| 361 | + } |
| 362 | +/* L170: */ |
| 363 | + } |
| 364 | + } else { |
| 365 | + |
| 366 | +/* Form C := alpha*A'*B' + beta*C */ |
| 367 | + |
| 368 | + i__1 = *n; |
| 369 | + for (j = 1; j <= *n; ++j) { |
| 370 | + i__2 = *m; |
| 371 | + for (i = 1; i <= *m; ++i) { |
| 372 | + temp = 0.; |
| 373 | + i__3 = *k; |
| 374 | + for (l = 1; l <= *k; ++l) { |
| 375 | + temp += A(l,i) * B(j,l); |
| 376 | +/* L180: */ |
| 377 | + } |
| 378 | + if (*beta == 0.) { |
| 379 | + C(i,j) = *alpha * temp; |
| 380 | + } else { |
| 381 | + C(i,j) = *alpha * temp + *beta * C(i,j); |
| 382 | + } |
| 383 | +/* L190: */ |
| 384 | + } |
| 385 | +/* L200: */ |
| 386 | + } |
| 387 | + } |
| 388 | + } |
| 389 | + |
| 390 | + return 0; |
| 391 | + |
| 392 | +/* End of DGEMM . */ |
| 393 | + |
| 394 | +} /* dgemm_ */ |
| 395 | + |