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
00036
00037 #ifndef TEMPLATE_LAPACK_LATRS_HEADER
00038 #define TEMPLATE_LAPACK_LATRS_HEADER
00039
00040
00041 template<class Treal>
00042 int template_lapack_latrs(const char *uplo, const char *trans, const char *diag, const char *
00043 normin, const integer *n, const Treal *a, const integer *lda, Treal *x,
00044 Treal *scale, Treal *cnorm, integer *info)
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
00208
00209 integer c__1 = 1;
00210 Treal c_b36 = .5;
00211
00212
00213 integer a_dim1, a_offset, i__1, i__2, i__3;
00214 Treal d__1, d__2, d__3;
00215
00216 integer jinc;
00217 Treal xbnd;
00218 integer imax;
00219 Treal tmax, tjjs, xmax, grow, sumj;
00220 integer i__, j;
00221 Treal tscal, uscal;
00222 integer jlast;
00223 logical upper;
00224 Treal xj;
00225 Treal bignum;
00226 logical notran;
00227 integer jfirst;
00228 Treal smlnum;
00229 logical nounit;
00230 Treal rec, tjj;
00231 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00232
00233
00234 a_dim1 = *lda;
00235 a_offset = 1 + a_dim1 * 1;
00236 a -= a_offset;
00237 --x;
00238 --cnorm;
00239
00240
00241 *info = 0;
00242 upper = template_blas_lsame(uplo, "U");
00243 notran = template_blas_lsame(trans, "N");
00244 nounit = template_blas_lsame(diag, "N");
00245
00246
00247
00248 if (! upper && ! template_blas_lsame(uplo, "L")) {
00249 *info = -1;
00250 } else if (! notran && ! template_blas_lsame(trans, "T") && !
00251 template_blas_lsame(trans, "C")) {
00252 *info = -2;
00253 } else if (! nounit && ! template_blas_lsame(diag, "U")) {
00254 *info = -3;
00255 } else if (! template_blas_lsame(normin, "Y") && ! template_blas_lsame(normin,
00256 "N")) {
00257 *info = -4;
00258 } else if (*n < 0) {
00259 *info = -5;
00260 } else if (*lda < maxMACRO(1,*n)) {
00261 *info = -7;
00262 }
00263 if (*info != 0) {
00264 i__1 = -(*info);
00265 template_blas_erbla("LATRS ", &i__1);
00266 return 0;
00267 }
00268
00269
00270
00271 if (*n == 0) {
00272 return 0;
00273 }
00274
00275
00276
00277 smlnum = template_lapack_lamch("Safe minimum", (Treal)0) / template_lapack_lamch("Precision", (Treal)0);
00278 bignum = 1. / smlnum;
00279 *scale = 1.;
00280
00281 if (template_blas_lsame(normin, "N")) {
00282
00283
00284
00285 if (upper) {
00286
00287
00288
00289 i__1 = *n;
00290 for (j = 1; j <= i__1; ++j) {
00291 i__2 = j - 1;
00292 cnorm[j] = template_blas_asum(&i__2, &a_ref(1, j), &c__1);
00293
00294 }
00295 } else {
00296
00297
00298
00299 i__1 = *n - 1;
00300 for (j = 1; j <= i__1; ++j) {
00301 i__2 = *n - j;
00302 cnorm[j] = template_blas_asum(&i__2, &a_ref(j + 1, j), &c__1);
00303
00304 }
00305 cnorm[*n] = 0.;
00306 }
00307 }
00308
00309
00310
00311
00312 imax = template_blas_idamax(n, &cnorm[1], &c__1);
00313 tmax = cnorm[imax];
00314 if (tmax <= bignum) {
00315 tscal = 1.;
00316 } else {
00317 tscal = 1. / (smlnum * tmax);
00318 dscal_(n, &tscal, &cnorm[1], &c__1);
00319 }
00320
00321
00322
00323
00324 j = template_blas_idamax(n, &x[1], &c__1);
00325 xmax = (d__1 = x[j], absMACRO(d__1));
00326 xbnd = xmax;
00327 if (notran) {
00328
00329
00330
00331 if (upper) {
00332 jfirst = *n;
00333 jlast = 1;
00334 jinc = -1;
00335 } else {
00336 jfirst = 1;
00337 jlast = *n;
00338 jinc = 1;
00339 }
00340
00341 if (tscal != 1.) {
00342 grow = 0.;
00343 goto L50;
00344 }
00345
00346 if (nounit) {
00347
00348
00349
00350
00351
00352
00353 grow = 1. / maxMACRO(xbnd,smlnum);
00354 xbnd = grow;
00355 i__1 = jlast;
00356 i__2 = jinc;
00357 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00358
00359
00360
00361 if (grow <= smlnum) {
00362 goto L50;
00363 }
00364
00365
00366
00367 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
00368
00369 d__1 = xbnd, d__2 = minMACRO(1.,tjj) * grow;
00370 xbnd = minMACRO(d__1,d__2);
00371 if (tjj + cnorm[j] >= smlnum) {
00372
00373
00374
00375 grow *= tjj / (tjj + cnorm[j]);
00376 } else {
00377
00378
00379
00380 grow = 0.;
00381 }
00382
00383 }
00384 grow = xbnd;
00385 } else {
00386
00387
00388
00389
00390
00391
00392 d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
00393 grow = minMACRO(d__1,d__2);
00394 i__2 = jlast;
00395 i__1 = jinc;
00396 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00397
00398
00399
00400 if (grow <= smlnum) {
00401 goto L50;
00402 }
00403
00404
00405
00406 grow *= 1. / (cnorm[j] + 1.);
00407
00408 }
00409 }
00410 L50:
00411
00412 ;
00413 } else {
00414
00415
00416
00417 if (upper) {
00418 jfirst = 1;
00419 jlast = *n;
00420 jinc = 1;
00421 } else {
00422 jfirst = *n;
00423 jlast = 1;
00424 jinc = -1;
00425 }
00426
00427 if (tscal != 1.) {
00428 grow = 0.;
00429 goto L80;
00430 }
00431
00432 if (nounit) {
00433
00434
00435
00436
00437
00438
00439 grow = 1. / maxMACRO(xbnd,smlnum);
00440 xbnd = grow;
00441 i__1 = jlast;
00442 i__2 = jinc;
00443 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00444
00445
00446
00447 if (grow <= smlnum) {
00448 goto L80;
00449 }
00450
00451
00452
00453 xj = cnorm[j] + 1.;
00454
00455 d__1 = grow, d__2 = xbnd / xj;
00456 grow = minMACRO(d__1,d__2);
00457
00458
00459
00460 tjj = (d__1 = a_ref(j, j), absMACRO(d__1));
00461 if (xj > tjj) {
00462 xbnd *= tjj / xj;
00463 }
00464
00465 }
00466 grow = minMACRO(grow,xbnd);
00467 } else {
00468
00469
00470
00471
00472
00473
00474 d__1 = 1., d__2 = 1. / maxMACRO(xbnd,smlnum);
00475 grow = minMACRO(d__1,d__2);
00476 i__2 = jlast;
00477 i__1 = jinc;
00478 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00479
00480
00481
00482 if (grow <= smlnum) {
00483 goto L80;
00484 }
00485
00486
00487
00488 xj = cnorm[j] + 1.;
00489 grow /= xj;
00490
00491 }
00492 }
00493 L80:
00494 ;
00495 }
00496
00497 if (grow * tscal > smlnum) {
00498
00499
00500
00501
00502 template_blas_trsv(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1);
00503 } else {
00504
00505
00506
00507 if (xmax > bignum) {
00508
00509
00510
00511
00512 *scale = bignum / xmax;
00513 dscal_(n, scale, &x[1], &c__1);
00514 xmax = bignum;
00515 }
00516
00517 if (notran) {
00518
00519
00520
00521 i__1 = jlast;
00522 i__2 = jinc;
00523 for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00524
00525
00526
00527 xj = (d__1 = x[j], absMACRO(d__1));
00528 if (nounit) {
00529 tjjs = a_ref(j, j) * tscal;
00530 } else {
00531 tjjs = tscal;
00532 if (tscal == 1.) {
00533 goto L100;
00534 }
00535 }
00536 tjj = absMACRO(tjjs);
00537 if (tjj > smlnum) {
00538
00539
00540
00541 if (tjj < 1.) {
00542 if (xj > tjj * bignum) {
00543
00544
00545
00546 rec = 1. / xj;
00547 dscal_(n, &rec, &x[1], &c__1);
00548 *scale *= rec;
00549 xmax *= rec;
00550 }
00551 }
00552 x[j] /= tjjs;
00553 xj = (d__1 = x[j], absMACRO(d__1));
00554 } else if (tjj > 0.) {
00555
00556
00557
00558 if (xj > tjj * bignum) {
00559
00560
00561
00562
00563 rec = tjj * bignum / xj;
00564 if (cnorm[j] > 1.) {
00565
00566
00567
00568
00569 rec /= cnorm[j];
00570 }
00571 dscal_(n, &rec, &x[1], &c__1);
00572 *scale *= rec;
00573 xmax *= rec;
00574 }
00575 x[j] /= tjjs;
00576 xj = (d__1 = x[j], absMACRO(d__1));
00577 } else {
00578
00579
00580
00581
00582 i__3 = *n;
00583 for (i__ = 1; i__ <= i__3; ++i__) {
00584 x[i__] = 0.;
00585
00586 }
00587 x[j] = 1.;
00588 xj = 1.;
00589 *scale = 0.;
00590 xmax = 0.;
00591 }
00592 L100:
00593
00594
00595
00596
00597 if (xj > 1.) {
00598 rec = 1. / xj;
00599 if (cnorm[j] > (bignum - xmax) * rec) {
00600
00601
00602
00603 rec *= .5;
00604 dscal_(n, &rec, &x[1], &c__1);
00605 *scale *= rec;
00606 }
00607 } else if (xj * cnorm[j] > bignum - xmax) {
00608
00609
00610
00611 dscal_(n, &c_b36, &x[1], &c__1);
00612 *scale *= .5;
00613 }
00614
00615 if (upper) {
00616 if (j > 1) {
00617
00618
00619
00620
00621 i__3 = j - 1;
00622 d__1 = -x[j] * tscal;
00623 daxpy_(&i__3, &d__1, &a_ref(1, j), &c__1, &x[1], &
00624 c__1);
00625 i__3 = j - 1;
00626 i__ = template_blas_idamax(&i__3, &x[1], &c__1);
00627 xmax = (d__1 = x[i__], absMACRO(d__1));
00628 }
00629 } else {
00630 if (j < *n) {
00631
00632
00633
00634
00635 i__3 = *n - j;
00636 d__1 = -x[j] * tscal;
00637 daxpy_(&i__3, &d__1, &a_ref(j + 1, j), &c__1, &x[j +
00638 1], &c__1);
00639 i__3 = *n - j;
00640 i__ = j + template_blas_idamax(&i__3, &x[j + 1], &c__1);
00641 xmax = (d__1 = x[i__], absMACRO(d__1));
00642 }
00643 }
00644
00645 }
00646
00647 } else {
00648
00649
00650
00651 i__2 = jlast;
00652 i__1 = jinc;
00653 for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00654
00655
00656
00657
00658 xj = (d__1 = x[j], absMACRO(d__1));
00659 uscal = tscal;
00660 rec = 1. / maxMACRO(xmax,1.);
00661 if (cnorm[j] > (bignum - xj) * rec) {
00662
00663
00664
00665 rec *= .5;
00666 if (nounit) {
00667 tjjs = a_ref(j, j) * tscal;
00668 } else {
00669 tjjs = tscal;
00670 }
00671 tjj = absMACRO(tjjs);
00672 if (tjj > 1.) {
00673
00674
00675
00676
00677 d__1 = 1., d__2 = rec * tjj;
00678 rec = minMACRO(d__1,d__2);
00679 uscal /= tjjs;
00680 }
00681 if (rec < 1.) {
00682 dscal_(n, &rec, &x[1], &c__1);
00683 *scale *= rec;
00684 xmax *= rec;
00685 }
00686 }
00687
00688 sumj = 0.;
00689 if (uscal == 1.) {
00690
00691
00692
00693
00694 if (upper) {
00695 i__3 = j - 1;
00696 sumj = ddot_(&i__3, &a_ref(1, j), &c__1, &x[1], &c__1)
00697 ;
00698 } else if (j < *n) {
00699 i__3 = *n - j;
00700 sumj = ddot_(&i__3, &a_ref(j + 1, j), &c__1, &x[j + 1]
00701 , &c__1);
00702 }
00703 } else {
00704
00705
00706
00707 if (upper) {
00708 i__3 = j - 1;
00709 for (i__ = 1; i__ <= i__3; ++i__) {
00710 sumj += a_ref(i__, j) * uscal * x[i__];
00711
00712 }
00713 } else if (j < *n) {
00714 i__3 = *n;
00715 for (i__ = j + 1; i__ <= i__3; ++i__) {
00716 sumj += a_ref(i__, j) * uscal * x[i__];
00717
00718 }
00719 }
00720 }
00721
00722 if (uscal == tscal) {
00723
00724
00725
00726
00727 x[j] -= sumj;
00728 xj = (d__1 = x[j], absMACRO(d__1));
00729 if (nounit) {
00730 tjjs = a_ref(j, j) * tscal;
00731 } else {
00732 tjjs = tscal;
00733 if (tscal == 1.) {
00734 goto L150;
00735 }
00736 }
00737
00738
00739
00740 tjj = absMACRO(tjjs);
00741 if (tjj > smlnum) {
00742
00743
00744
00745 if (tjj < 1.) {
00746 if (xj > tjj * bignum) {
00747
00748
00749
00750 rec = 1. / xj;
00751 dscal_(n, &rec, &x[1], &c__1);
00752 *scale *= rec;
00753 xmax *= rec;
00754 }
00755 }
00756 x[j] /= tjjs;
00757 } else if (tjj > 0.) {
00758
00759
00760
00761 if (xj > tjj * bignum) {
00762
00763
00764
00765 rec = tjj * bignum / xj;
00766 dscal_(n, &rec, &x[1], &c__1);
00767 *scale *= rec;
00768 xmax *= rec;
00769 }
00770 x[j] /= tjjs;
00771 } else {
00772
00773
00774
00775
00776 i__3 = *n;
00777 for (i__ = 1; i__ <= i__3; ++i__) {
00778 x[i__] = 0.;
00779
00780 }
00781 x[j] = 1.;
00782 *scale = 0.;
00783 xmax = 0.;
00784 }
00785 L150:
00786 ;
00787 } else {
00788
00789
00790
00791
00792 x[j] = x[j] / tjjs - sumj;
00793 }
00794
00795 d__2 = xmax, d__3 = (d__1 = x[j], absMACRO(d__1));
00796 xmax = maxMACRO(d__2,d__3);
00797
00798 }
00799 }
00800 *scale /= tscal;
00801 }
00802
00803
00804
00805 if (tscal != 1.) {
00806 d__1 = 1. / tscal;
00807 dscal_(n, &d__1, &cnorm[1], &c__1);
00808 }
00809
00810 return 0;
00811
00812
00813
00814 }
00815
00816 #undef a_ref
00817
00818
00819 #endif