00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #ifndef TEMPLATE_LAPACK_LATRS_HEADER
00036 #define TEMPLATE_LAPACK_LATRS_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *
00041 normin, const integer *n, const Treal *a, const integer *lda, Treal *x,
00042 Treal *scale, Treal *cnorm, integer *info)
00043 {
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 integer c__1 = 1;
00208 Treal c_b36 = .5;
00209
00210
00211 integer a_dim1, a_offset, i__1, i__2, i__3;
00212 Treal d__1, d__2, d__3;
00213
00214 integer jinc;
00215 Treal xbnd;
00216 integer imax;
00217 Treal tmax, tjjs, xmax, grow, sumj;
00218 integer i__, j;
00219 Treal tscal, uscal;
00220 integer jlast;
00221 logical upper;
00222 Treal xj;
00223 Treal bignum;
00224 logical notran;
00225 integer jfirst;
00226 Treal smlnum;
00227 logical nounit;
00228 Treal rec, tjj;
00229 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00230
00231
00232 a_dim1 = *lda;
00233 a_offset = 1 + a_dim1 * 1;
00234 a -= a_offset;
00235 --x;
00236 --cnorm;
00237
00238
00239 *info = 0;
00240 upper = template_blas_lsame(uplo, "U");
00241 notran = template_blas_lsame(trans, "N");
00242 nounit = template_blas_lsame(diag, "N");
00243
00244
00245
00246 if (! upper && ! template_blas_lsame(uplo, "L")) {
00247 *info = -1;
00248 } else if (! notran && ! template_blas_lsame(trans, "T") && !
00249 template_blas_lsame(trans, "C")) {
00250 *info = -2;
00251 } else if (! nounit && ! template_blas_lsame(diag, "U")) {
00252 *info = -3;
00253 } else if (! template_blas_lsame(normin, "Y") && ! template_blas_lsame(normin,
00254 "N")) {
00255 *info = -4;
00256 } else if (*n < 0) {
00257 *info = -5;
00258 } else if (*lda < maxMACRO(1,*n)) {
00259 *info = -7;
00260 }
00261 if (*info != 0) {
00262 i__1 = -(*info);
00263 template_blas_erbla("LATRS ", &i__1);
00264 return 0;
00265 }
00266
00267
00268
00269 if (*n == 0) {
00270 return 0;
00271 }
00272
00273
00274
00275 smlnum = template_lapack_lamch("Safe minimum", (Treal)0) / template_lapack_lamch("Precision", (Treal)0);
00276 bignum = 1. / smlnum;
00277 *scale = 1.;
00278
00279 if (template_blas_lsame(normin, "N")) {
00280
00281
00282
00283 if (upper) {
00284
00285
00286
00287 i__1 = *n;
00288 for (j = 1; j <= i__1; ++j) {
00289 i__2 = j - 1;
00290 cnorm[j] = template_blas_asum(&i__2, &a_ref(1, j), &c__1);
00291
00292 }
00293 } else {
00294
00295
00296
00297 i__1 = *n - 1;
00298 for (j = 1; j <= i__1; ++j) {
00299 i__2 = *n - j;
00300 cnorm[j] = template_blas_asum(&i__2, &a_ref(j + 1, j), &c__1);
00301
00302 }
00303 cnorm[*n] = 0.;
00304 }
00305 }
00306
00307
00308
00309
00310 imax = template_blas_idamax(n, &cnorm[1], &c__1);
00311 tmax = cnorm[imax];
00312 if (tmax <= bignum) {
00313 tscal = 1.;
00314 } else {
00315 tscal = 1. / (smlnum * tmax);
00316 dscal_(n, &tscal, &cnorm[1], &c__1);
00317 }
00318
00319
00320
00321
00322 j = template_blas_idamax(n, &x[1], &c__1);
00323 xmax = (d__1 = x[j], absMACRO(d__1));
00324 xbnd = xmax;
00325 if (notran) {
00326
00327
00328
00329 if (upper) {
00330 jfirst = *n;
00331 jlast = 1;
00332 jinc = -1;
00333 } else {
00334 jfirst = 1;
00335 jlast = *n;
00336 jinc = 1;
00337 }
00338
00339 if (tscal != 1.) {
00340 grow = 0.;
00341 goto L50;
00342 }
00343
00344 if (nounit) {
00345
00346
00347
00348
00349
00350
00351 grow = 1. / maxMACRO(xbnd,smlnum);
00352 xbnd = grow;
00353 i__1 = jlast;
00354 i__2 = jinc;
00355 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00356
00357
00358
00359 if (grow <= smlnum) {
00360 goto L50;
00361 }
00362
00363
00364
00365 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
00366
00367 d__1 = xbnd, d__2 = minMACRO(1.,tjj) * grow;
00368 xbnd = minMACRO(d__1,d__2);
00369 if (tjj + cnorm[j] >= smlnum) {
00370
00371
00372
00373 grow *= tjj / (tjj + cnorm[j]);
00374 } else {
00375
00376
00377
00378 grow = 0.;
00379 }
00380
00381 }
00382 grow = xbnd;
00383 } else {
00384
00385
00386
00387
00388
00389
00390 d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
00391 grow = minMACRO(d__1,d__2);
00392 i__2 = jlast;
00393 i__1 = jinc;
00394 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00395
00396
00397
00398 if (grow <= smlnum) {
00399 goto L50;
00400 }
00401
00402
00403
00404 grow *= 1. / (cnorm[j] + 1.);
00405
00406 }
00407 }
00408 L50:
00409
00410 ;
00411 } else {
00412
00413
00414
00415 if (upper) {
00416 jfirst = 1;
00417 jlast = *n;
00418 jinc = 1;
00419 } else {
00420 jfirst = *n;
00421 jlast = 1;
00422 jinc = -1;
00423 }
00424
00425 if (tscal != 1.) {
00426 grow = 0.;
00427 goto L80;
00428 }
00429
00430 if (nounit) {
00431
00432
00433
00434
00435
00436
00437 grow = 1. / maxMACRO(xbnd,smlnum);
00438 xbnd = grow;
00439 i__1 = jlast;
00440 i__2 = jinc;
00441 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00442
00443
00444
00445 if (grow <= smlnum) {
00446 goto L80;
00447 }
00448
00449
00450
00451 xj = cnorm[j] + 1.;
00452
00453 d__1 = grow, d__2 = xbnd / xj;
00454 grow = minMACRO(d__1,d__2);
00455
00456
00457
00458 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
00459 if (xj > tjj) {
00460 xbnd *= tjj / xj;
00461 }
00462
00463 }
00464 grow = minMACRO(grow,xbnd);
00465 } else {
00466
00467
00468
00469
00470
00471
00472 d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
00473 grow = minMACRO(d__1,d__2);
00474 i__2 = jlast;
00475 i__1 = jinc;
00476 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00477
00478
00479
00480 if (grow <= smlnum) {
00481 goto L80;
00482 }
00483
00484
00485
00486 xj = cnorm[j] + 1.;
00487 grow /= xj;
00488
00489 }
00490 }
00491 L80:
00492 ;
00493 }
00494
00495 if (grow * tscal > smlnum) {
00496
00497
00498
00499
00500 template_blas_trsv(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1);
00501 } else {
00502
00503
00504
00505 if (xmax > bignum) {
00506
00507
00508
00509
00510 *scale = bignum / xmax;
00511 dscal_(n, scale, &x[1], &c__1);
00512 xmax = bignum;
00513 }
00514
00515 if (notran) {
00516
00517
00518
00519 i__1 = jlast;
00520 i__2 = jinc;
00521 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00522
00523
00524
00525 xj = (d__1 = x[j], absMACRO(d__1));
00526 if (nounit) {
00527 tjjs = a_ref(j, j) * tscal;
00528 } else {
00529 tjjs = tscal;
00530 if (tscal == 1.) {
00531 goto L100;
00532 }
00533 }
00534 tjj = absMACRO(tjjs);
00535 if (tjj > smlnum) {
00536
00537
00538
00539 if (tjj < 1.) {
00540 if (xj > tjj * bignum) {
00541
00542
00543
00544 rec = 1. / xj;
00545 dscal_(n, &rec, &x[1], &c__1);
00546 *scale *= rec;
00547 xmax *= rec;
00548 }
00549 }
00550 x[j] /= tjjs;
00551 xj = (d__1 = x[j], absMACRO(d__1));
00552 } else if (tjj > 0.) {
00553
00554
00555
00556 if (xj > tjj * bignum) {
00557
00558
00559
00560
00561 rec = tjj * bignum / xj;
00562 if (cnorm[j] > 1.) {
00563
00564
00565
00566
00567 rec /= cnorm[j];
00568 }
00569 dscal_(n, &rec, &x[1], &c__1);
00570 *scale *= rec;
00571 xmax *= rec;
00572 }
00573 x[j] /= tjjs;
00574 xj = (d__1 = x[j], absMACRO(d__1));
00575 } else {
00576
00577
00578
00579
00580 i__3 = *n;
00581 for (i__ = 1; i__ <= i__3; ++i__) {
00582 x[i__] = 0.;
00583
00584 }
00585 x[j] = 1.;
00586 xj = 1.;
00587 *scale = 0.;
00588 xmax = 0.;
00589 }
00590 L100:
00591
00592
00593
00594
00595 if (xj > 1.) {
00596 rec = 1. / xj;
00597 if (cnorm[j] > (bignum - xmax) * rec) {
00598
00599
00600
00601 rec *= .5;
00602 dscal_(n, &rec, &x[1], &c__1);
00603 *scale *= rec;
00604 }
00605 } else if (xj * cnorm[j] > bignum - xmax) {
00606
00607
00608
00609 dscal_(n, &c_b36, &x[1], &c__1);
00610 *scale *= .5;
00611 }
00612
00613 if (upper) {
00614 if (j > 1) {
00615
00616
00617
00618
00619 i__3 = j - 1;
00620 d__1 = -x[j] * tscal;
00621 daxpy_(&i__3, &d__1, &a_ref(1, j), &c__1, &x[1], &
00622 c__1);
00623 i__3 = j - 1;
00624 i__ = template_blas_idamax(&i__3, &x[1], &c__1);
00625 xmax = (d__1 = x[i__], absMACRO(d__1));
00626 }
00627 } else {
00628 if (j < *n) {
00629
00630
00631
00632
00633 i__3 = *n - j;
00634 d__1 = -x[j] * tscal;
00635 daxpy_(&i__3, &d__1, &a_ref(j + 1, j), &c__1, &x[j +
00636 1], &c__1);
00637 i__3 = *n - j;
00638 i__ = j + template_blas_idamax(&i__3, &x[j + 1], &c__1);
00639 xmax = (d__1 = x[i__], absMACRO(d__1));
00640 }
00641 }
00642
00643 }
00644
00645 } else {
00646
00647
00648
00649 i__2 = jlast;
00650 i__1 = jinc;
00651 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00652
00653
00654
00655
00656 xj = (d__1 = x[j], absMACRO(d__1));
00657 uscal = tscal;
00658 rec = 1. / maxMACRO(xmax,1.);
00659 if (cnorm[j] > (bignum - xj) * rec) {
00660
00661
00662
00663 rec *= .5;
00664 if (nounit) {
00665 tjjs = a_ref(j, j) * tscal;
00666 } else {
00667 tjjs = tscal;
00668 }
00669 tjj = absMACRO(tjjs);
00670 if (tjj > 1.) {
00671
00672
00673
00674
00675 d__1 = 1., d__2 = rec * tjj;
00676 rec = minMACRO(d__1,d__2);
00677 uscal /= tjjs;
00678 }
00679 if (rec < 1.) {
00680 dscal_(n, &rec, &x[1], &c__1);
00681 *scale *= rec;
00682 xmax *= rec;
00683 }
00684 }
00685
00686 sumj = 0.;
00687 if (uscal == 1.) {
00688
00689
00690
00691
00692 if (upper) {
00693 i__3 = j - 1;
00694 sumj = ddot_(&i__3, &a_ref(1, j), &c__1, &x[1], &c__1)
00695 ;
00696 } else if (j < *n) {
00697 i__3 = *n - j;
00698 sumj = ddot_(&i__3, &a_ref(j + 1, j), &c__1, &x[j + 1]
00699 , &c__1);
00700 }
00701 } else {
00702
00703
00704
00705 if (upper) {
00706 i__3 = j - 1;
00707 for (i__ = 1; i__ <= i__3; ++i__) {
00708 sumj += a_ref(i__, j) * uscal * x[i__];
00709
00710 }
00711 } else if (j < *n) {
00712 i__3 = *n;
00713 for (i__ = j + 1; i__ <= i__3; ++i__) {
00714 sumj += a_ref(i__, j) * uscal * x[i__];
00715
00716 }
00717 }
00718 }
00719
00720 if (uscal == tscal) {
00721
00722
00723
00724
00725 x[j] -= sumj;
00726 xj = (d__1 = x[j], absMACRO(d__1));
00727 if (nounit) {
00728 tjjs = a_ref(j, j) * tscal;
00729 } else {
00730 tjjs = tscal;
00731 if (tscal == 1.) {
00732 goto L150;
00733 }
00734 }
00735
00736
00737
00738 tjj = absMACRO(tjjs);
00739 if (tjj > smlnum) {
00740
00741
00742
00743 if (tjj < 1.) {
00744 if (xj > tjj * bignum) {
00745
00746
00747
00748 rec = 1. / xj;
00749 dscal_(n, &rec, &x[1], &c__1);
00750 *scale *= rec;
00751 xmax *= rec;
00752 }
00753 }
00754 x[j] /= tjjs;
00755 } else if (tjj > 0.) {
00756
00757
00758
00759 if (xj > tjj * bignum) {
00760
00761
00762
00763 rec = tjj * bignum / xj;
00764 dscal_(n, &rec, &x[1], &c__1);
00765 *scale *= rec;
00766 xmax *= rec;
00767 }
00768 x[j] /= tjjs;
00769 } else {
00770
00771
00772
00773
00774 i__3 = *n;
00775 for (i__ = 1; i__ <= i__3; ++i__) {
00776 x[i__] = 0.;
00777
00778 }
00779 x[j] = 1.;
00780 *scale = 0.;
00781 xmax = 0.;
00782 }
00783 L150:
00784 ;
00785 } else {
00786
00787
00788
00789
00790 x[j] = x[j] / tjjs - sumj;
00791 }
00792
00793 d__2 = xmax, d__3 = (d__1 = x[j], absMACRO(d__1));
00794 xmax = maxMACRO(d__2,d__3);
00795
00796 }
00797 }
00798 *scale /= tscal;
00799 }
00800
00801
00802
00803 if (tscal != 1.) {
00804 d__1 = 1. / tscal;
00805 dscal_(n, &d__1, &cnorm[1], &c__1);
00806 }
00807
00808 return 0;
00809
00810
00811
00812 }
00813
00814 #undef a_ref
00815
00816
00817 #endif