[37] | 1 | #include "mex.h" |
---|
| 2 | |
---|
| 3 | #if 1 |
---|
| 4 | |
---|
| 5 | /* Minka version, from Golub and van Loan. |
---|
| 6 | * Does not use blocking, hence not as efficient as the BLAS routine. |
---|
| 7 | */ |
---|
| 8 | |
---|
| 9 | int dtrsm(char *side, char *uplo, char *transa, char *diag, |
---|
| 10 | int *m, int *n, double *alpha, double *a, int *lda, |
---|
| 11 | double *b, int *ldb) |
---|
| 12 | { |
---|
| 13 | int i,j,k; |
---|
| 14 | #define A(I,J) a[(I) + (J)*(*lda)] |
---|
| 15 | #define B(I,J) b[(I) + (J)*(*ldb)] |
---|
| 16 | |
---|
| 17 | if(*uplo == 'U') { |
---|
| 18 | /* Alg 3.1.4 on p90 */ |
---|
| 19 | for(j=0;j<*n;j++) { |
---|
| 20 | for(k=*m-1;k>=0;k--) { |
---|
| 21 | if(B(k,j) != 0.) { |
---|
| 22 | B(k,j) /= A(k,k); |
---|
| 23 | for(i=0;i<k;i++) { |
---|
| 24 | B(i,j) -= B(k,j) * A(i,k); |
---|
| 25 | } |
---|
| 26 | } |
---|
| 27 | } |
---|
| 28 | } |
---|
| 29 | } else { |
---|
| 30 | for(j=0;j<*n;j++) { |
---|
| 31 | for(k=0;k<*m;k++) { |
---|
| 32 | if(B(k,j) != 0.) { |
---|
| 33 | B(k,j) /= A(k,k); |
---|
| 34 | for(i=k+1;i<*m;i++) { |
---|
| 35 | B(i,j) -= B(k,j) * A(i,k); |
---|
| 36 | } |
---|
| 37 | } |
---|
| 38 | } |
---|
| 39 | } |
---|
| 40 | } |
---|
| 41 | } |
---|
| 42 | |
---|
| 43 | #else |
---|
| 44 | |
---|
| 45 | /* BLAS code from netlib.org */ |
---|
| 46 | |
---|
| 47 | typedef int logical; |
---|
| 48 | |
---|
| 49 | logical lsame_(char *ca, char *cb) |
---|
| 50 | { |
---|
| 51 | return(*ca == *cb); |
---|
| 52 | } |
---|
| 53 | |
---|
| 54 | int xerbla_(char *srname, int *info) |
---|
| 55 | { |
---|
| 56 | mexErrMsgTxt(srname); |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | /* Subroutine */ |
---|
| 60 | int dtrsm(char *side, char *uplo, char *transa, char *diag, |
---|
| 61 | int *m, int *n, double *alpha, double *a, int *lda, |
---|
| 62 | double *b, int *ldb) |
---|
| 63 | { |
---|
| 64 | /* System generated locals */ |
---|
| 65 | int a_dim1, a_offset, b_dim1, b_offset; |
---|
| 66 | |
---|
| 67 | /* Local variables */ |
---|
| 68 | static int info; |
---|
| 69 | static double temp; |
---|
| 70 | static int i, j, k; |
---|
| 71 | static logical lside; |
---|
| 72 | static int nrowa; |
---|
| 73 | static logical upper; |
---|
| 74 | static logical nounit; |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | /* Purpose |
---|
| 78 | ======= |
---|
| 79 | |
---|
| 80 | DTRSM solves one of the matrix equations |
---|
| 81 | |
---|
| 82 | op( A )*X = alpha*B, or X*op( A ) = alpha*B, |
---|
| 83 | |
---|
| 84 | where alpha is a scalar, X and B are m by n matrices, A is a unit, or |
---|
| 85 | |
---|
| 86 | non-unit, upper or lower triangular matrix and op( A ) is one of |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | op( A ) = A or op( A ) = A'. |
---|
| 90 | |
---|
| 91 | The matrix X is overwritten on B. |
---|
| 92 | |
---|
| 93 | Parameters |
---|
| 94 | ========== |
---|
| 95 | |
---|
| 96 | SIDE - CHARACTER*1. |
---|
| 97 | On entry, SIDE specifies whether op( A ) appears on the left |
---|
| 98 | |
---|
| 99 | or right of X as follows: |
---|
| 100 | |
---|
| 101 | SIDE = 'L' or 'l' op( A )*X = alpha*B. |
---|
| 102 | |
---|
| 103 | SIDE = 'R' or 'r' X*op( A ) = alpha*B. |
---|
| 104 | |
---|
| 105 | Unchanged on exit. |
---|
| 106 | |
---|
| 107 | UPLO - CHARACTER*1. |
---|
| 108 | On entry, UPLO specifies whether the matrix A is an upper or |
---|
| 109 | |
---|
| 110 | lower triangular matrix as follows: |
---|
| 111 | |
---|
| 112 | UPLO = 'U' or 'u' A is an upper triangular matrix. |
---|
| 113 | |
---|
| 114 | UPLO = 'L' or 'l' A is a lower triangular matrix. |
---|
| 115 | |
---|
| 116 | Unchanged on exit. |
---|
| 117 | |
---|
| 118 | TRANSA - CHARACTER*1. |
---|
| 119 | On entry, TRANSA specifies the form of op( A ) to be used in |
---|
| 120 | |
---|
| 121 | the matrix multiplication as follows: |
---|
| 122 | |
---|
| 123 | TRANSA = 'N' or 'n' op( A ) = A. |
---|
| 124 | |
---|
| 125 | TRANSA = 'T' or 't' op( A ) = A'. |
---|
| 126 | |
---|
| 127 | TRANSA = 'C' or 'c' op( A ) = A'. |
---|
| 128 | |
---|
| 129 | Unchanged on exit. |
---|
| 130 | |
---|
| 131 | DIAG - CHARACTER*1. |
---|
| 132 | On entry, DIAG specifies whether or not A is unit triangular |
---|
| 133 | |
---|
| 134 | as follows: |
---|
| 135 | |
---|
| 136 | DIAG = 'U' or 'u' A is assumed to be unit triangular. |
---|
| 137 | |
---|
| 138 | DIAG = 'N' or 'n' A is not assumed to be unit |
---|
| 139 | triangular. |
---|
| 140 | |
---|
| 141 | Unchanged on exit. |
---|
| 142 | |
---|
| 143 | M - INTEGER. |
---|
| 144 | On entry, M specifies the number of rows of B. M must be at |
---|
| 145 | |
---|
| 146 | least zero. |
---|
| 147 | Unchanged on exit. |
---|
| 148 | |
---|
| 149 | N - INTEGER. |
---|
| 150 | On entry, N specifies the number of columns of B. N must be |
---|
| 151 | |
---|
| 152 | at least zero. |
---|
| 153 | Unchanged on exit. |
---|
| 154 | |
---|
| 155 | ALPHA - DOUBLE PRECISION. |
---|
| 156 | On entry, ALPHA specifies the scalar alpha. When alpha is |
---|
| 157 | |
---|
| 158 | zero then A is not referenced and B need not be set before |
---|
| 159 | |
---|
| 160 | entry. |
---|
| 161 | Unchanged on exit. |
---|
| 162 | |
---|
| 163 | A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m |
---|
| 164 | |
---|
| 165 | when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. |
---|
| 166 | |
---|
| 167 | Before entry with UPLO = 'U' or 'u', the leading k by k |
---|
| 168 | |
---|
| 169 | upper triangular part of the array A must contain the upper |
---|
| 170 | |
---|
| 171 | triangular matrix and the strictly lower triangular part of |
---|
| 172 | |
---|
| 173 | A is not referenced. |
---|
| 174 | Before entry with UPLO = 'L' or 'l', the leading k by k |
---|
| 175 | |
---|
| 176 | lower triangular part of the array A must contain the lower |
---|
| 177 | |
---|
| 178 | triangular matrix and the strictly upper triangular part of |
---|
| 179 | |
---|
| 180 | A is not referenced. |
---|
| 181 | Note that when DIAG = 'U' or 'u', the diagonal elements of |
---|
| 182 | |
---|
| 183 | A are not referenced either, but are assumed to be unity. |
---|
| 184 | |
---|
| 185 | Unchanged on exit. |
---|
| 186 | |
---|
| 187 | LDA - INTEGER. |
---|
| 188 | On entry, LDA specifies the first dimension of A as declared |
---|
| 189 | |
---|
| 190 | in the calling (sub) program. When SIDE = 'L' or 'l' then |
---|
| 191 | |
---|
| 192 | LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' |
---|
| 193 | |
---|
| 194 | then LDA must be at least max( 1, n ). |
---|
| 195 | Unchanged on exit. |
---|
| 196 | |
---|
| 197 | B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). |
---|
| 198 | Before entry, the leading m by n part of the array B must |
---|
| 199 | |
---|
| 200 | contain the right-hand side matrix B, and on exit is |
---|
| 201 | |
---|
| 202 | overwritten by the solution matrix X. |
---|
| 203 | |
---|
| 204 | LDB - INTEGER. |
---|
| 205 | On entry, LDB specifies the first dimension of B as declared |
---|
| 206 | |
---|
| 207 | in the calling (sub) program. LDB must be at least |
---|
| 208 | |
---|
| 209 | max( 1, m ). |
---|
| 210 | Unchanged on exit. |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | Level 3 Blas routine. |
---|
| 214 | |
---|
| 215 | |
---|
| 216 | -- Written on 8-February-1989. |
---|
| 217 | Jack Dongarra, Argonne National Laboratory. |
---|
| 218 | Iain Duff, AERE Harwell. |
---|
| 219 | Jeremy Du Croz, Numerical Algorithms Group Ltd. |
---|
| 220 | Sven Hammarling, Numerical Algorithms Group Ltd. |
---|
| 221 | |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | Test the input parameters. |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | Parameter adjustments |
---|
| 228 | Function Body */ |
---|
| 229 | |
---|
| 230 | #define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)] |
---|
| 231 | #define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)] |
---|
| 232 | |
---|
| 233 | lside = lsame_(side, "L"); |
---|
| 234 | if (lside) { |
---|
| 235 | nrowa = *m; |
---|
| 236 | } else { |
---|
| 237 | nrowa = *n; |
---|
| 238 | } |
---|
| 239 | nounit = lsame_(diag, "N"); |
---|
| 240 | upper = lsame_(uplo, "U"); |
---|
| 241 | |
---|
| 242 | info = 0; |
---|
| 243 | if (! lside && ! lsame_(side, "R")) { |
---|
| 244 | info = 1; |
---|
| 245 | } else if (! upper && ! lsame_(uplo, "L")) { |
---|
| 246 | info = 2; |
---|
| 247 | } else if (! lsame_(transa, "N") && ! lsame_(transa, "T") |
---|
| 248 | && ! lsame_(transa, "C")) { |
---|
| 249 | info = 3; |
---|
| 250 | } else if (! lsame_(diag, "U") && ! lsame_(diag, "N")) { |
---|
| 251 | info = 4; |
---|
| 252 | } else if (*m < 0) { |
---|
| 253 | info = 5; |
---|
| 254 | } else if (*n < 0) { |
---|
| 255 | info = 6; |
---|
| 256 | } else if (*lda < nrowa) { |
---|
| 257 | info = 9; |
---|
| 258 | } else if (*ldb < *m) { |
---|
| 259 | info = 11; |
---|
| 260 | } |
---|
| 261 | if (info != 0) { |
---|
| 262 | xerbla_("DTRSM ", &info); |
---|
| 263 | return 0; |
---|
| 264 | } |
---|
| 265 | |
---|
| 266 | /* Quick return if possible. */ |
---|
| 267 | |
---|
| 268 | if (*n == 0) { |
---|
| 269 | return 0; |
---|
| 270 | } |
---|
| 271 | |
---|
| 272 | /* And when alpha.eq.zero. */ |
---|
| 273 | |
---|
| 274 | if (*alpha == 0.) { |
---|
| 275 | for (j = 1; j <= *n; ++j) { |
---|
| 276 | for (i = 1; i <= *m; ++i) { |
---|
| 277 | B(i,j) = 0.; |
---|
| 278 | /* L10: */ |
---|
| 279 | } |
---|
| 280 | /* L20: */ |
---|
| 281 | } |
---|
| 282 | return 0; |
---|
| 283 | } |
---|
| 284 | |
---|
| 285 | /* Start the operations. */ |
---|
| 286 | |
---|
| 287 | if (lside) { |
---|
| 288 | if (lsame_(transa, "N")) { |
---|
| 289 | |
---|
| 290 | /* Form B := alpha*inv( A )*B. */ |
---|
| 291 | |
---|
| 292 | if (upper) { |
---|
| 293 | for (j = 1; j <= *n; ++j) { |
---|
| 294 | if (*alpha != 1.) { |
---|
| 295 | for (i = 1; i <= *m; ++i) { |
---|
| 296 | B(i,j) = *alpha * B(i,j); |
---|
| 297 | /* L30: */ |
---|
| 298 | } |
---|
| 299 | } |
---|
| 300 | /* Alg 3.1.4 on p90 */ |
---|
| 301 | for (k = *m; k >= 1; --k) { |
---|
| 302 | if (B(k,j) != 0.) { |
---|
| 303 | printf("%d %d %g %g\n",k,j,B(k,j),A(k,k)); |
---|
| 304 | B(k,j) /= A(k,k); |
---|
| 305 | for (i = 1; i <= k-1; ++i) { |
---|
| 306 | B(i,j) -= B(k,j) * A(i,k); |
---|
| 307 | /* L40: */ |
---|
| 308 | } |
---|
| 309 | } |
---|
| 310 | /* L50: */ |
---|
| 311 | } |
---|
| 312 | /* L60: */ |
---|
| 313 | } |
---|
| 314 | } else { |
---|
| 315 | for (j = 1; j <= *n; ++j) { |
---|
| 316 | if (*alpha != 1.) { |
---|
| 317 | for (i = 1; i <= *m; ++i) { |
---|
| 318 | B(i,j) = *alpha * B(i,j); |
---|
| 319 | /* L70: */ |
---|
| 320 | } |
---|
| 321 | } |
---|
| 322 | for (k = 1; k <= *m; ++k) { |
---|
| 323 | if (B(k,j) != 0.) { |
---|
| 324 | if (nounit) { |
---|
| 325 | B(k,j) /= A(k,k); |
---|
| 326 | } |
---|
| 327 | for (i = k + 1; i <= *m; ++i) { |
---|
| 328 | B(i,j) -= B(k,j) * A(i,k); |
---|
| 329 | /* L80: */ |
---|
| 330 | } |
---|
| 331 | } |
---|
| 332 | /* L90: */ |
---|
| 333 | } |
---|
| 334 | /* L100: */ |
---|
| 335 | } |
---|
| 336 | } |
---|
| 337 | } else { |
---|
| 338 | |
---|
| 339 | /* Form B := alpha*inv( A' )*B. */ |
---|
| 340 | |
---|
| 341 | if (upper) { |
---|
| 342 | for (j = 1; j <= *n; ++j) { |
---|
| 343 | for (i = 1; i <= *m; ++i) { |
---|
| 344 | temp = *alpha * B(i,j); |
---|
| 345 | for (k = 1; k <= i-1; ++k) { |
---|
| 346 | temp -= A(k,i) * B(k,j); |
---|
| 347 | /* L110: */ |
---|
| 348 | } |
---|
| 349 | if (nounit) { |
---|
| 350 | temp /= A(i,i); |
---|
| 351 | } |
---|
| 352 | B(i,j) = temp; |
---|
| 353 | /* L120: */ |
---|
| 354 | } |
---|
| 355 | /* L130: */ |
---|
| 356 | } |
---|
| 357 | } else { |
---|
| 358 | for (j = 1; j <= *n; ++j) { |
---|
| 359 | for (i = *m; i >= 1; --i) { |
---|
| 360 | temp = *alpha * B(i,j); |
---|
| 361 | for (k = i + 1; k <= *m; ++k) { |
---|
| 362 | temp -= A(k,i) * B(k,j); |
---|
| 363 | /* L140: */ |
---|
| 364 | } |
---|
| 365 | if (nounit) { |
---|
| 366 | temp /= A(i,i); |
---|
| 367 | } |
---|
| 368 | B(i,j) = temp; |
---|
| 369 | /* L150: */ |
---|
| 370 | } |
---|
| 371 | /* L160: */ |
---|
| 372 | } |
---|
| 373 | } |
---|
| 374 | } |
---|
| 375 | } else { |
---|
| 376 | if (lsame_(transa, "N")) { |
---|
| 377 | |
---|
| 378 | /* Form B := alpha*B*inv( A ). */ |
---|
| 379 | |
---|
| 380 | if (upper) { |
---|
| 381 | for (j = 1; j <= *n; ++j) { |
---|
| 382 | if (*alpha != 1.) { |
---|
| 383 | for (i = 1; i <= *m; ++i) { |
---|
| 384 | B(i,j) = *alpha * B(i,j); |
---|
| 385 | /* L170: */ |
---|
| 386 | } |
---|
| 387 | } |
---|
| 388 | for (k = 1; k <= j-1; ++k) { |
---|
| 389 | if (A(k,j) != 0.) { |
---|
| 390 | for (i = 1; i <= *m; ++i) { |
---|
| 391 | B(i,j) -= A(k,j) * B(i,k); |
---|
| 392 | /* L180: */ |
---|
| 393 | } |
---|
| 394 | } |
---|
| 395 | /* L190: */ |
---|
| 396 | } |
---|
| 397 | if (nounit) { |
---|
| 398 | temp = 1. / A(j,j); |
---|
| 399 | for (i = 1; i <= *m; ++i) { |
---|
| 400 | B(i,j) = temp * B(i,j); |
---|
| 401 | /* L200: */ |
---|
| 402 | } |
---|
| 403 | } |
---|
| 404 | /* L210: */ |
---|
| 405 | } |
---|
| 406 | } else { |
---|
| 407 | for (j = *n; j >= 1; --j) { |
---|
| 408 | if (*alpha != 1.) { |
---|
| 409 | for (i = 1; i <= *m; ++i) { |
---|
| 410 | B(i,j) = *alpha * B(i,j); |
---|
| 411 | /* L220: */ |
---|
| 412 | } |
---|
| 413 | } |
---|
| 414 | for (k = j + 1; k <= *n; ++k) { |
---|
| 415 | if (A(k,j) != 0.) { |
---|
| 416 | for (i = 1; i <= *m; ++i) { |
---|
| 417 | B(i,j) -= A(k,j) * B(i,k); |
---|
| 418 | /* L230: */ |
---|
| 419 | } |
---|
| 420 | } |
---|
| 421 | /* L240: */ |
---|
| 422 | } |
---|
| 423 | if (nounit) { |
---|
| 424 | temp = 1. / A(j,j); |
---|
| 425 | for (i = 1; i <= *m; ++i) { |
---|
| 426 | B(i,j) = temp * B(i,j); |
---|
| 427 | /* L250: */ |
---|
| 428 | } |
---|
| 429 | } |
---|
| 430 | /* L260: */ |
---|
| 431 | } |
---|
| 432 | } |
---|
| 433 | } else { |
---|
| 434 | |
---|
| 435 | /* Form B := alpha*B*inv( A' ). */ |
---|
| 436 | |
---|
| 437 | if (upper) { |
---|
| 438 | for (k = *n; k >= 1; --k) { |
---|
| 439 | if (nounit) { |
---|
| 440 | temp = 1. / A(k,k); |
---|
| 441 | for (i = 1; i <= *m; ++i) { |
---|
| 442 | B(i,k) = temp * B(i,k); |
---|
| 443 | /* L270: */ |
---|
| 444 | } |
---|
| 445 | } |
---|
| 446 | for (j = 1; j <= k-1; ++j) { |
---|
| 447 | if (A(j,k) != 0.) { |
---|
| 448 | temp = A(j,k); |
---|
| 449 | for (i = 1; i <= *m; ++i) { |
---|
| 450 | B(i,j) -= temp * B(i,k); |
---|
| 451 | /* L280: */ |
---|
| 452 | } |
---|
| 453 | } |
---|
| 454 | /* L290: */ |
---|
| 455 | } |
---|
| 456 | if (*alpha != 1.) { |
---|
| 457 | for (i = 1; i <= *m; ++i) { |
---|
| 458 | B(i,k) = *alpha * B(i,k); |
---|
| 459 | /* L300: */ |
---|
| 460 | } |
---|
| 461 | } |
---|
| 462 | /* L310: */ |
---|
| 463 | } |
---|
| 464 | } else { |
---|
| 465 | for (k = 1; k <= *n; ++k) { |
---|
| 466 | if (nounit) { |
---|
| 467 | temp = 1. / A(k,k); |
---|
| 468 | for (i = 1; i <= *m; ++i) { |
---|
| 469 | B(i,k) = temp * B(i,k); |
---|
| 470 | /* L320: */ |
---|
| 471 | } |
---|
| 472 | } |
---|
| 473 | for (j = k + 1; j <= *n; ++j) { |
---|
| 474 | if (A(j,k) != 0.) { |
---|
| 475 | temp = A(j,k); |
---|
| 476 | for (i = 1; i <= *m; ++i) { |
---|
| 477 | B(i,j) -= temp * B(i,k); |
---|
| 478 | /* L330: */ |
---|
| 479 | } |
---|
| 480 | } |
---|
| 481 | /* L340: */ |
---|
| 482 | } |
---|
| 483 | if (*alpha != 1.) { |
---|
| 484 | for (i = 1; i <= *m; ++i) { |
---|
| 485 | B(i,k) = *alpha * B(i,k); |
---|
| 486 | /* L350: */ |
---|
| 487 | } |
---|
| 488 | } |
---|
| 489 | /* L360: */ |
---|
| 490 | } |
---|
| 491 | } |
---|
| 492 | } |
---|
| 493 | } |
---|
| 494 | |
---|
| 495 | return 0; |
---|
| 496 | |
---|
| 497 | /* End of DTRSM . */ |
---|
| 498 | |
---|
| 499 | } /* dtrsm_ */ |
---|
| 500 | #endif |
---|