annotate ext/cblas/src/dtrmm.c @ 209:ccd2019190bf msvc

Some MSVC fixes, including (temporarily, probably) renaming the FFT source file to avoid getting it mixed up with the Vamp SDK one in our object dir
author Chris Cannam
date Thu, 01 Feb 2018 16:34:08 +0000
parents 45330e0d2819
children
rev   line source
Chris@202 1 /* dtrmm.f -- translated by f2c (version 20061008).
Chris@202 2 You must link the resulting object file with libf2c:
Chris@202 3 on Microsoft Windows system, link with libf2c.lib;
Chris@202 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
Chris@202 5 or, if you install libf2c.a in a standard place, with -lf2c -lm
Chris@202 6 -- in that order, at the end of the command line, as in
Chris@202 7 cc *.o -lf2c -lm
Chris@202 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
Chris@202 9
Chris@202 10 http://www.netlib.org/f2c/libf2c.zip
Chris@202 11 */
Chris@202 12
Chris@202 13 #include "f2c.h"
Chris@202 14 #include "blaswrap.h"
Chris@202 15
Chris@202 16 /* Subroutine */ int dtrmm_(char *side, char *uplo, char *transa, char *diag,
Chris@202 17 integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
Chris@202 18 lda, doublereal *b, integer *ldb)
Chris@202 19 {
Chris@202 20 /* System generated locals */
Chris@202 21 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
Chris@202 22
Chris@202 23 /* Local variables */
Chris@202 24 integer i__, j, k, info;
Chris@202 25 doublereal temp;
Chris@202 26 logical lside;
Chris@202 27 extern logical lsame_(char *, char *);
Chris@202 28 integer nrowa;
Chris@202 29 logical upper;
Chris@202 30 extern /* Subroutine */ int xerbla_(char *, integer *);
Chris@202 31 logical nounit;
Chris@202 32
Chris@202 33 /* .. Scalar Arguments .. */
Chris@202 34 /* .. */
Chris@202 35 /* .. Array Arguments .. */
Chris@202 36 /* .. */
Chris@202 37
Chris@202 38 /* Purpose */
Chris@202 39 /* ======= */
Chris@202 40
Chris@202 41 /* DTRMM performs one of the matrix-matrix operations */
Chris@202 42
Chris@202 43 /* B := alpha*op( A )*B, or B := alpha*B*op( A ), */
Chris@202 44
Chris@202 45 /* where alpha is a scalar, B is an m by n matrix, A is a unit, or */
Chris@202 46 /* non-unit, upper or lower triangular matrix and op( A ) is one of */
Chris@202 47
Chris@202 48 /* op( A ) = A or op( A ) = A'. */
Chris@202 49
Chris@202 50 /* Arguments */
Chris@202 51 /* ========== */
Chris@202 52
Chris@202 53 /* SIDE - CHARACTER*1. */
Chris@202 54 /* On entry, SIDE specifies whether op( A ) multiplies B from */
Chris@202 55 /* the left or right as follows: */
Chris@202 56
Chris@202 57 /* SIDE = 'L' or 'l' B := alpha*op( A )*B. */
Chris@202 58
Chris@202 59 /* SIDE = 'R' or 'r' B := alpha*B*op( A ). */
Chris@202 60
Chris@202 61 /* Unchanged on exit. */
Chris@202 62
Chris@202 63 /* UPLO - CHARACTER*1. */
Chris@202 64 /* On entry, UPLO specifies whether the matrix A is an upper or */
Chris@202 65 /* lower triangular matrix as follows: */
Chris@202 66
Chris@202 67 /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
Chris@202 68
Chris@202 69 /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
Chris@202 70
Chris@202 71 /* Unchanged on exit. */
Chris@202 72
Chris@202 73 /* TRANSA - CHARACTER*1. */
Chris@202 74 /* On entry, TRANSA specifies the form of op( A ) to be used in */
Chris@202 75 /* the matrix multiplication as follows: */
Chris@202 76
Chris@202 77 /* TRANSA = 'N' or 'n' op( A ) = A. */
Chris@202 78
Chris@202 79 /* TRANSA = 'T' or 't' op( A ) = A'. */
Chris@202 80
Chris@202 81 /* TRANSA = 'C' or 'c' op( A ) = A'. */
Chris@202 82
Chris@202 83 /* Unchanged on exit. */
Chris@202 84
Chris@202 85 /* DIAG - CHARACTER*1. */
Chris@202 86 /* On entry, DIAG specifies whether or not A is unit triangular */
Chris@202 87 /* as follows: */
Chris@202 88
Chris@202 89 /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
Chris@202 90
Chris@202 91 /* DIAG = 'N' or 'n' A is not assumed to be unit */
Chris@202 92 /* triangular. */
Chris@202 93
Chris@202 94 /* Unchanged on exit. */
Chris@202 95
Chris@202 96 /* M - INTEGER. */
Chris@202 97 /* On entry, M specifies the number of rows of B. M must be at */
Chris@202 98 /* least zero. */
Chris@202 99 /* Unchanged on exit. */
Chris@202 100
Chris@202 101 /* N - INTEGER. */
Chris@202 102 /* On entry, N specifies the number of columns of B. N must be */
Chris@202 103 /* at least zero. */
Chris@202 104 /* Unchanged on exit. */
Chris@202 105
Chris@202 106 /* ALPHA - DOUBLE PRECISION. */
Chris@202 107 /* On entry, ALPHA specifies the scalar alpha. When alpha is */
Chris@202 108 /* zero then A is not referenced and B need not be set before */
Chris@202 109 /* entry. */
Chris@202 110 /* Unchanged on exit. */
Chris@202 111
Chris@202 112 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m */
Chris@202 113 /* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. */
Chris@202 114 /* Before entry with UPLO = 'U' or 'u', the leading k by k */
Chris@202 115 /* upper triangular part of the array A must contain the upper */
Chris@202 116 /* triangular matrix and the strictly lower triangular part of */
Chris@202 117 /* A is not referenced. */
Chris@202 118 /* Before entry with UPLO = 'L' or 'l', the leading k by k */
Chris@202 119 /* lower triangular part of the array A must contain the lower */
Chris@202 120 /* triangular matrix and the strictly upper triangular part of */
Chris@202 121 /* A is not referenced. */
Chris@202 122 /* Note that when DIAG = 'U' or 'u', the diagonal elements of */
Chris@202 123 /* A are not referenced either, but are assumed to be unity. */
Chris@202 124 /* Unchanged on exit. */
Chris@202 125
Chris@202 126 /* LDA - INTEGER. */
Chris@202 127 /* On entry, LDA specifies the first dimension of A as declared */
Chris@202 128 /* in the calling (sub) program. When SIDE = 'L' or 'l' then */
Chris@202 129 /* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' */
Chris@202 130 /* then LDA must be at least max( 1, n ). */
Chris@202 131 /* Unchanged on exit. */
Chris@202 132
Chris@202 133 /* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). */
Chris@202 134 /* Before entry, the leading m by n part of the array B must */
Chris@202 135 /* contain the matrix B, and on exit is overwritten by the */
Chris@202 136 /* transformed matrix. */
Chris@202 137
Chris@202 138 /* LDB - INTEGER. */
Chris@202 139 /* On entry, LDB specifies the first dimension of B as declared */
Chris@202 140 /* in the calling (sub) program. LDB must be at least */
Chris@202 141 /* max( 1, m ). */
Chris@202 142 /* Unchanged on exit. */
Chris@202 143
Chris@202 144
Chris@202 145 /* Level 3 Blas routine. */
Chris@202 146
Chris@202 147 /* -- Written on 8-February-1989. */
Chris@202 148 /* Jack Dongarra, Argonne National Laboratory. */
Chris@202 149 /* Iain Duff, AERE Harwell. */
Chris@202 150 /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
Chris@202 151 /* Sven Hammarling, Numerical Algorithms Group Ltd. */
Chris@202 152
Chris@202 153
Chris@202 154 /* .. External Functions .. */
Chris@202 155 /* .. */
Chris@202 156 /* .. External Subroutines .. */
Chris@202 157 /* .. */
Chris@202 158 /* .. Intrinsic Functions .. */
Chris@202 159 /* .. */
Chris@202 160 /* .. Local Scalars .. */
Chris@202 161 /* .. */
Chris@202 162 /* .. Parameters .. */
Chris@202 163 /* .. */
Chris@202 164
Chris@202 165 /* Test the input parameters. */
Chris@202 166
Chris@202 167 /* Parameter adjustments */
Chris@202 168 a_dim1 = *lda;
Chris@202 169 a_offset = 1 + a_dim1;
Chris@202 170 a -= a_offset;
Chris@202 171 b_dim1 = *ldb;
Chris@202 172 b_offset = 1 + b_dim1;
Chris@202 173 b -= b_offset;
Chris@202 174
Chris@202 175 /* Function Body */
Chris@202 176 lside = lsame_(side, "L");
Chris@202 177 if (lside) {
Chris@202 178 nrowa = *m;
Chris@202 179 } else {
Chris@202 180 nrowa = *n;
Chris@202 181 }
Chris@202 182 nounit = lsame_(diag, "N");
Chris@202 183 upper = lsame_(uplo, "U");
Chris@202 184
Chris@202 185 info = 0;
Chris@202 186 if (! lside && ! lsame_(side, "R")) {
Chris@202 187 info = 1;
Chris@202 188 } else if (! upper && ! lsame_(uplo, "L")) {
Chris@202 189 info = 2;
Chris@202 190 } else if (! lsame_(transa, "N") && ! lsame_(transa,
Chris@202 191 "T") && ! lsame_(transa, "C")) {
Chris@202 192 info = 3;
Chris@202 193 } else if (! lsame_(diag, "U") && ! lsame_(diag,
Chris@202 194 "N")) {
Chris@202 195 info = 4;
Chris@202 196 } else if (*m < 0) {
Chris@202 197 info = 5;
Chris@202 198 } else if (*n < 0) {
Chris@202 199 info = 6;
Chris@202 200 } else if (*lda < max(1,nrowa)) {
Chris@202 201 info = 9;
Chris@202 202 } else if (*ldb < max(1,*m)) {
Chris@202 203 info = 11;
Chris@202 204 }
Chris@202 205 if (info != 0) {
Chris@202 206 xerbla_("DTRMM ", &info);
Chris@202 207 return 0;
Chris@202 208 }
Chris@202 209
Chris@202 210 /* Quick return if possible. */
Chris@202 211
Chris@202 212 if (*m == 0 || *n == 0) {
Chris@202 213 return 0;
Chris@202 214 }
Chris@202 215
Chris@202 216 /* And when alpha.eq.zero. */
Chris@202 217
Chris@202 218 if (*alpha == 0.) {
Chris@202 219 i__1 = *n;
Chris@202 220 for (j = 1; j <= i__1; ++j) {
Chris@202 221 i__2 = *m;
Chris@202 222 for (i__ = 1; i__ <= i__2; ++i__) {
Chris@202 223 b[i__ + j * b_dim1] = 0.;
Chris@202 224 /* L10: */
Chris@202 225 }
Chris@202 226 /* L20: */
Chris@202 227 }
Chris@202 228 return 0;
Chris@202 229 }
Chris@202 230
Chris@202 231 /* Start the operations. */
Chris@202 232
Chris@202 233 if (lside) {
Chris@202 234 if (lsame_(transa, "N")) {
Chris@202 235
Chris@202 236 /* Form B := alpha*A*B. */
Chris@202 237
Chris@202 238 if (upper) {
Chris@202 239 i__1 = *n;
Chris@202 240 for (j = 1; j <= i__1; ++j) {
Chris@202 241 i__2 = *m;
Chris@202 242 for (k = 1; k <= i__2; ++k) {
Chris@202 243 if (b[k + j * b_dim1] != 0.) {
Chris@202 244 temp = *alpha * b[k + j * b_dim1];
Chris@202 245 i__3 = k - 1;
Chris@202 246 for (i__ = 1; i__ <= i__3; ++i__) {
Chris@202 247 b[i__ + j * b_dim1] += temp * a[i__ + k *
Chris@202 248 a_dim1];
Chris@202 249 /* L30: */
Chris@202 250 }
Chris@202 251 if (nounit) {
Chris@202 252 temp *= a[k + k * a_dim1];
Chris@202 253 }
Chris@202 254 b[k + j * b_dim1] = temp;
Chris@202 255 }
Chris@202 256 /* L40: */
Chris@202 257 }
Chris@202 258 /* L50: */
Chris@202 259 }
Chris@202 260 } else {
Chris@202 261 i__1 = *n;
Chris@202 262 for (j = 1; j <= i__1; ++j) {
Chris@202 263 for (k = *m; k >= 1; --k) {
Chris@202 264 if (b[k + j * b_dim1] != 0.) {
Chris@202 265 temp = *alpha * b[k + j * b_dim1];
Chris@202 266 b[k + j * b_dim1] = temp;
Chris@202 267 if (nounit) {
Chris@202 268 b[k + j * b_dim1] *= a[k + k * a_dim1];
Chris@202 269 }
Chris@202 270 i__2 = *m;
Chris@202 271 for (i__ = k + 1; i__ <= i__2; ++i__) {
Chris@202 272 b[i__ + j * b_dim1] += temp * a[i__ + k *
Chris@202 273 a_dim1];
Chris@202 274 /* L60: */
Chris@202 275 }
Chris@202 276 }
Chris@202 277 /* L70: */
Chris@202 278 }
Chris@202 279 /* L80: */
Chris@202 280 }
Chris@202 281 }
Chris@202 282 } else {
Chris@202 283
Chris@202 284 /* Form B := alpha*A'*B. */
Chris@202 285
Chris@202 286 if (upper) {
Chris@202 287 i__1 = *n;
Chris@202 288 for (j = 1; j <= i__1; ++j) {
Chris@202 289 for (i__ = *m; i__ >= 1; --i__) {
Chris@202 290 temp = b[i__ + j * b_dim1];
Chris@202 291 if (nounit) {
Chris@202 292 temp *= a[i__ + i__ * a_dim1];
Chris@202 293 }
Chris@202 294 i__2 = i__ - 1;
Chris@202 295 for (k = 1; k <= i__2; ++k) {
Chris@202 296 temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
Chris@202 297 /* L90: */
Chris@202 298 }
Chris@202 299 b[i__ + j * b_dim1] = *alpha * temp;
Chris@202 300 /* L100: */
Chris@202 301 }
Chris@202 302 /* L110: */
Chris@202 303 }
Chris@202 304 } else {
Chris@202 305 i__1 = *n;
Chris@202 306 for (j = 1; j <= i__1; ++j) {
Chris@202 307 i__2 = *m;
Chris@202 308 for (i__ = 1; i__ <= i__2; ++i__) {
Chris@202 309 temp = b[i__ + j * b_dim1];
Chris@202 310 if (nounit) {
Chris@202 311 temp *= a[i__ + i__ * a_dim1];
Chris@202 312 }
Chris@202 313 i__3 = *m;
Chris@202 314 for (k = i__ + 1; k <= i__3; ++k) {
Chris@202 315 temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
Chris@202 316 /* L120: */
Chris@202 317 }
Chris@202 318 b[i__ + j * b_dim1] = *alpha * temp;
Chris@202 319 /* L130: */
Chris@202 320 }
Chris@202 321 /* L140: */
Chris@202 322 }
Chris@202 323 }
Chris@202 324 }
Chris@202 325 } else {
Chris@202 326 if (lsame_(transa, "N")) {
Chris@202 327
Chris@202 328 /* Form B := alpha*B*A. */
Chris@202 329
Chris@202 330 if (upper) {
Chris@202 331 for (j = *n; j >= 1; --j) {
Chris@202 332 temp = *alpha;
Chris@202 333 if (nounit) {
Chris@202 334 temp *= a[j + j * a_dim1];
Chris@202 335 }
Chris@202 336 i__1 = *m;
Chris@202 337 for (i__ = 1; i__ <= i__1; ++i__) {
Chris@202 338 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
Chris@202 339 /* L150: */
Chris@202 340 }
Chris@202 341 i__1 = j - 1;
Chris@202 342 for (k = 1; k <= i__1; ++k) {
Chris@202 343 if (a[k + j * a_dim1] != 0.) {
Chris@202 344 temp = *alpha * a[k + j * a_dim1];
Chris@202 345 i__2 = *m;
Chris@202 346 for (i__ = 1; i__ <= i__2; ++i__) {
Chris@202 347 b[i__ + j * b_dim1] += temp * b[i__ + k *
Chris@202 348 b_dim1];
Chris@202 349 /* L160: */
Chris@202 350 }
Chris@202 351 }
Chris@202 352 /* L170: */
Chris@202 353 }
Chris@202 354 /* L180: */
Chris@202 355 }
Chris@202 356 } else {
Chris@202 357 i__1 = *n;
Chris@202 358 for (j = 1; j <= i__1; ++j) {
Chris@202 359 temp = *alpha;
Chris@202 360 if (nounit) {
Chris@202 361 temp *= a[j + j * a_dim1];
Chris@202 362 }
Chris@202 363 i__2 = *m;
Chris@202 364 for (i__ = 1; i__ <= i__2; ++i__) {
Chris@202 365 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
Chris@202 366 /* L190: */
Chris@202 367 }
Chris@202 368 i__2 = *n;
Chris@202 369 for (k = j + 1; k <= i__2; ++k) {
Chris@202 370 if (a[k + j * a_dim1] != 0.) {
Chris@202 371 temp = *alpha * a[k + j * a_dim1];
Chris@202 372 i__3 = *m;
Chris@202 373 for (i__ = 1; i__ <= i__3; ++i__) {
Chris@202 374 b[i__ + j * b_dim1] += temp * b[i__ + k *
Chris@202 375 b_dim1];
Chris@202 376 /* L200: */
Chris@202 377 }
Chris@202 378 }
Chris@202 379 /* L210: */
Chris@202 380 }
Chris@202 381 /* L220: */
Chris@202 382 }
Chris@202 383 }
Chris@202 384 } else {
Chris@202 385
Chris@202 386 /* Form B := alpha*B*A'. */
Chris@202 387
Chris@202 388 if (upper) {
Chris@202 389 i__1 = *n;
Chris@202 390 for (k = 1; k <= i__1; ++k) {
Chris@202 391 i__2 = k - 1;
Chris@202 392 for (j = 1; j <= i__2; ++j) {
Chris@202 393 if (a[j + k * a_dim1] != 0.) {
Chris@202 394 temp = *alpha * a[j + k * a_dim1];
Chris@202 395 i__3 = *m;
Chris@202 396 for (i__ = 1; i__ <= i__3; ++i__) {
Chris@202 397 b[i__ + j * b_dim1] += temp * b[i__ + k *
Chris@202 398 b_dim1];
Chris@202 399 /* L230: */
Chris@202 400 }
Chris@202 401 }
Chris@202 402 /* L240: */
Chris@202 403 }
Chris@202 404 temp = *alpha;
Chris@202 405 if (nounit) {
Chris@202 406 temp *= a[k + k * a_dim1];
Chris@202 407 }
Chris@202 408 if (temp != 1.) {
Chris@202 409 i__2 = *m;
Chris@202 410 for (i__ = 1; i__ <= i__2; ++i__) {
Chris@202 411 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
Chris@202 412 /* L250: */
Chris@202 413 }
Chris@202 414 }
Chris@202 415 /* L260: */
Chris@202 416 }
Chris@202 417 } else {
Chris@202 418 for (k = *n; k >= 1; --k) {
Chris@202 419 i__1 = *n;
Chris@202 420 for (j = k + 1; j <= i__1; ++j) {
Chris@202 421 if (a[j + k * a_dim1] != 0.) {
Chris@202 422 temp = *alpha * a[j + k * a_dim1];
Chris@202 423 i__2 = *m;
Chris@202 424 for (i__ = 1; i__ <= i__2; ++i__) {
Chris@202 425 b[i__ + j * b_dim1] += temp * b[i__ + k *
Chris@202 426 b_dim1];
Chris@202 427 /* L270: */
Chris@202 428 }
Chris@202 429 }
Chris@202 430 /* L280: */
Chris@202 431 }
Chris@202 432 temp = *alpha;
Chris@202 433 if (nounit) {
Chris@202 434 temp *= a[k + k * a_dim1];
Chris@202 435 }
Chris@202 436 if (temp != 1.) {
Chris@202 437 i__1 = *m;
Chris@202 438 for (i__ = 1; i__ <= i__1; ++i__) {
Chris@202 439 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
Chris@202 440 /* L290: */
Chris@202 441 }
Chris@202 442 }
Chris@202 443 /* L300: */
Chris@202 444 }
Chris@202 445 }
Chris@202 446 }
Chris@202 447 }
Chris@202 448
Chris@202 449 return 0;
Chris@202 450
Chris@202 451 /* End of DTRMM . */
Chris@202 452
Chris@202 453 } /* dtrmm_ */