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_STEBZ_HEADER
00036 #define TEMPLATE_LAPACK_STEBZ_HEADER
00037
00038
00039 template<class Treal>
00040 int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal
00041 *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol,
00042 const Treal *d__, const Treal *e, integer *m, integer *nsplit,
00043 Treal *w, integer *iblock, integer *isplit, Treal *work,
00044 integer *iwork, 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 integer c__1 = 1;
00209 integer c_n1 = -1;
00210 integer c__3 = 3;
00211 integer c__2 = 2;
00212 integer c__0 = 0;
00213
00214
00215 integer i__1, i__2, i__3;
00216 Treal d__1, d__2, d__3, d__4, d__5;
00217
00218 integer iend, ioff, iout, itmp1, j, jdisc;
00219 integer iinfo;
00220 Treal atoli;
00221 integer iwoff;
00222 Treal bnorm;
00223 integer itmax;
00224 Treal wkill, rtoli, tnorm;
00225 integer ib, jb, ie, je, nb;
00226 Treal gl;
00227 integer im, in;
00228 integer ibegin;
00229 Treal gu;
00230 integer iw;
00231 Treal wl;
00232 integer irange, idiscl;
00233 Treal safemn, wu;
00234 integer idumma[1];
00235 integer idiscu, iorder;
00236 logical ncnvrg;
00237 Treal pivmin;
00238 logical toofew;
00239 integer nwl;
00240 Treal ulp, wlu, wul;
00241 integer nwu;
00242 Treal tmp1, tmp2;
00243
00244
00245 --iwork;
00246 --work;
00247 --isplit;
00248 --iblock;
00249 --w;
00250 --e;
00251 --d__;
00252
00253
00254 wlu = wul = 0;
00255
00256 *info = 0;
00257
00258
00259
00260 if (template_blas_lsame(range, "A")) {
00261 irange = 1;
00262 } else if (template_blas_lsame(range, "V")) {
00263 irange = 2;
00264 } else if (template_blas_lsame(range, "I")) {
00265 irange = 3;
00266 } else {
00267 irange = 0;
00268 }
00269
00270
00271
00272 if (template_blas_lsame(order, "B")) {
00273 iorder = 2;
00274 } else if (template_blas_lsame(order, "E")) {
00275 iorder = 1;
00276 } else {
00277 iorder = 0;
00278 }
00279
00280
00281
00282 if (irange <= 0) {
00283 *info = -1;
00284 } else if (iorder <= 0) {
00285 *info = -2;
00286 } else if (*n < 0) {
00287 *info = -3;
00288 } else if (irange == 2) {
00289 if (*vl >= *vu) {
00290 *info = -5;
00291 }
00292 } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) {
00293 *info = -6;
00294 } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) {
00295 *info = -7;
00296 }
00297
00298 if (*info != 0) {
00299 i__1 = -(*info);
00300 template_blas_erbla("STEBZ ", &i__1);
00301 return 0;
00302 }
00303
00304
00305
00306 *info = 0;
00307 ncnvrg = FALSE_;
00308 toofew = FALSE_;
00309
00310
00311
00312 *m = 0;
00313 if (*n == 0) {
00314 return 0;
00315 }
00316
00317
00318
00319 if (irange == 3 && *il == 1 && *iu == *n) {
00320 irange = 1;
00321 }
00322
00323
00324
00325
00326
00327 safemn = template_lapack_lamch("S", (Treal)0);
00328 ulp = template_lapack_lamch("P", (Treal)0);
00329 rtoli = ulp * 2.;
00330 nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
00331 ftnlen)1);
00332 if (nb <= 1) {
00333 nb = 0;
00334 }
00335
00336
00337
00338 if (*n == 1) {
00339 *nsplit = 1;
00340 isplit[1] = 1;
00341 if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
00342 *m = 0;
00343 } else {
00344 w[1] = d__[1];
00345 iblock[1] = 1;
00346 *m = 1;
00347 }
00348 return 0;
00349 }
00350
00351
00352
00353 *nsplit = 1;
00354 work[*n] = 0.;
00355 pivmin = 1.;
00356
00357
00358 i__1 = *n;
00359 for (j = 2; j <= i__1; ++j) {
00360
00361 d__1 = e[j - 1];
00362 tmp1 = d__1 * d__1;
00363
00364 d__2 = ulp;
00365 if ((d__1 = d__[j] * d__[j - 1], absMACRO(d__1)) * (d__2 * d__2) + safemn
00366 > tmp1) {
00367 isplit[*nsplit] = j - 1;
00368 ++(*nsplit);
00369 work[j - 1] = 0.;
00370 } else {
00371 work[j - 1] = tmp1;
00372 pivmin = maxMACRO(pivmin,tmp1);
00373 }
00374
00375 }
00376 isplit[*nsplit] = *n;
00377 pivmin *= safemn;
00378
00379
00380
00381 if (irange == 3) {
00382
00383
00384
00385
00386
00387
00388
00389 gu = d__[1];
00390 gl = d__[1];
00391 tmp1 = 0.;
00392
00393 i__1 = *n - 1;
00394 for (j = 1; j <= i__1; ++j) {
00395 tmp2 = template_blas_sqrt(work[j]);
00396
00397 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
00398 gu = maxMACRO(d__1,d__2);
00399
00400 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
00401 gl = minMACRO(d__1,d__2);
00402 tmp1 = tmp2;
00403
00404 }
00405
00406
00407 d__1 = gu, d__2 = d__[*n] + tmp1;
00408 gu = maxMACRO(d__1,d__2);
00409
00410 d__1 = gl, d__2 = d__[*n] - tmp1;
00411 gl = minMACRO(d__1,d__2);
00412
00413 d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00414 tnorm = maxMACRO(d__1,d__2);
00415 gl = gl - tnorm * 2. * ulp * *n - pivmin * 4.;
00416 gu = gu + tnorm * 2. * ulp * *n + pivmin * 2.;
00417
00418
00419
00420 itmax = (integer) ((template_blas_log(tnorm + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)) + 2;
00421 if (*abstol <= 0.) {
00422 atoli = ulp * tnorm;
00423 } else {
00424 atoli = *abstol;
00425 }
00426
00427 work[*n + 1] = gl;
00428 work[*n + 2] = gl;
00429 work[*n + 3] = gu;
00430 work[*n + 4] = gu;
00431 work[*n + 5] = gl;
00432 work[*n + 6] = gu;
00433 iwork[1] = -1;
00434 iwork[2] = -1;
00435 iwork[3] = *n + 1;
00436 iwork[4] = *n + 1;
00437 iwork[5] = *il - 1;
00438 iwork[6] = *iu;
00439
00440 template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
00441 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
00442 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
00443
00444 if (iwork[6] == *iu) {
00445 wl = work[*n + 1];
00446 wlu = work[*n + 3];
00447 nwl = iwork[1];
00448 wu = work[*n + 4];
00449 wul = work[*n + 2];
00450 nwu = iwork[4];
00451 } else {
00452 wl = work[*n + 2];
00453 wlu = work[*n + 4];
00454 nwl = iwork[2];
00455 wu = work[*n + 3];
00456 wul = work[*n + 1];
00457 nwu = iwork[3];
00458 }
00459
00460 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
00461 *info = 4;
00462 return 0;
00463 }
00464 } else {
00465
00466
00467
00468
00469 d__3 = absMACRO(d__[1]) + absMACRO(e[1]), d__4 = (d__1 = d__[*n], absMACRO(d__1)) + (
00470 d__2 = e[*n - 1], absMACRO(d__2));
00471 tnorm = maxMACRO(d__3,d__4);
00472
00473 i__1 = *n - 1;
00474 for (j = 2; j <= i__1; ++j) {
00475
00476 d__4 = tnorm, d__5 = (d__1 = d__[j], absMACRO(d__1)) + (d__2 = e[j - 1]
00477 , absMACRO(d__2)) + (d__3 = e[j], absMACRO(d__3));
00478 tnorm = maxMACRO(d__4,d__5);
00479
00480 }
00481
00482 if (*abstol <= 0.) {
00483 atoli = ulp * tnorm;
00484 } else {
00485 atoli = *abstol;
00486 }
00487
00488 if (irange == 2) {
00489 wl = *vl;
00490 wu = *vu;
00491 } else {
00492 wl = 0.;
00493 wu = 0.;
00494 }
00495 }
00496
00497
00498
00499
00500
00501 *m = 0;
00502 iend = 0;
00503 *info = 0;
00504 nwl = 0;
00505 nwu = 0;
00506
00507 i__1 = *nsplit;
00508 for (jb = 1; jb <= i__1; ++jb) {
00509 ioff = iend;
00510 ibegin = ioff + 1;
00511 iend = isplit[jb];
00512 in = iend - ioff;
00513
00514 if (in == 1) {
00515
00516
00517
00518 if (irange == 1 || wl >= d__[ibegin] - pivmin) {
00519 ++nwl;
00520 }
00521 if (irange == 1 || wu >= d__[ibegin] - pivmin) {
00522 ++nwu;
00523 }
00524 if (irange == 1 || ( wl < d__[ibegin] - pivmin && wu >= d__[ibegin]
00525 - pivmin ) ) {
00526 ++(*m);
00527 w[*m] = d__[ibegin];
00528 iblock[*m] = jb;
00529 }
00530 } else {
00531
00532
00533
00534
00535
00536
00537 gu = d__[ibegin];
00538 gl = d__[ibegin];
00539 tmp1 = 0.;
00540
00541 i__2 = iend - 1;
00542 for (j = ibegin; j <= i__2; ++j) {
00543 tmp2 = (d__1 = e[j], absMACRO(d__1));
00544
00545 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
00546 gu = maxMACRO(d__1,d__2);
00547
00548 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
00549 gl = minMACRO(d__1,d__2);
00550 tmp1 = tmp2;
00551
00552 }
00553
00554
00555 d__1 = gu, d__2 = d__[iend] + tmp1;
00556 gu = maxMACRO(d__1,d__2);
00557
00558 d__1 = gl, d__2 = d__[iend] - tmp1;
00559 gl = minMACRO(d__1,d__2);
00560
00561 d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00562 bnorm = maxMACRO(d__1,d__2);
00563 gl = gl - bnorm * 2. * ulp * in - pivmin * 2.;
00564 gu = gu + bnorm * 2. * ulp * in + pivmin * 2.;
00565
00566
00567
00568 if (*abstol <= 0.) {
00569
00570 d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00571 atoli = ulp * maxMACRO(d__1,d__2);
00572 } else {
00573 atoli = *abstol;
00574 }
00575
00576 if (irange > 1) {
00577 if (gu < wl) {
00578 nwl += in;
00579 nwu += in;
00580 goto L70;
00581 }
00582 gl = maxMACRO(gl,wl);
00583 gu = minMACRO(gu,wu);
00584 if (gl >= gu) {
00585 goto L70;
00586 }
00587 }
00588
00589
00590
00591 work[*n + 1] = gl;
00592 work[*n + in + 1] = gu;
00593 template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
00594 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
00595 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
00596 w[*m + 1], &iblock[*m + 1], &iinfo);
00597
00598 nwl += iwork[1];
00599 nwu += iwork[in + 1];
00600 iwoff = *m - iwork[1];
00601
00602
00603
00604 itmax = (integer) ((template_blas_log(gu - gl + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)
00605 ) + 2;
00606 template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
00607 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
00608 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
00609 &w[*m + 1], &iblock[*m + 1], &iinfo);
00610
00611
00612
00613
00614 i__2 = iout;
00615 for (j = 1; j <= i__2; ++j) {
00616 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
00617
00618
00619
00620 if (j > iout - iinfo) {
00621 ncnvrg = TRUE_;
00622 ib = -jb;
00623 } else {
00624 ib = jb;
00625 }
00626 i__3 = iwork[j + in] + iwoff;
00627 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
00628 w[je] = tmp1;
00629 iblock[je] = ib;
00630
00631 }
00632
00633 }
00634
00635 *m += im;
00636 }
00637 L70:
00638 ;
00639 }
00640
00641
00642
00643
00644 if (irange == 3) {
00645 im = 0;
00646 idiscl = *il - 1 - nwl;
00647 idiscu = nwu - *iu;
00648
00649 if (idiscl > 0 || idiscu > 0) {
00650 i__1 = *m;
00651 for (je = 1; je <= i__1; ++je) {
00652 if (w[je] <= wlu && idiscl > 0) {
00653 --idiscl;
00654 } else if (w[je] >= wul && idiscu > 0) {
00655 --idiscu;
00656 } else {
00657 ++im;
00658 w[im] = w[je];
00659 iblock[im] = iblock[je];
00660 }
00661
00662 }
00663 *m = im;
00664 }
00665 if (idiscl > 0 || idiscu > 0) {
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677 if (idiscl > 0) {
00678 wkill = wu;
00679 i__1 = idiscl;
00680 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
00681 iw = 0;
00682 i__2 = *m;
00683 for (je = 1; je <= i__2; ++je) {
00684 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
00685 iw = je;
00686 wkill = w[je];
00687 }
00688
00689 }
00690 iblock[iw] = 0;
00691
00692 }
00693 }
00694 if (idiscu > 0) {
00695
00696 wkill = wl;
00697 i__1 = idiscu;
00698 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
00699 iw = 0;
00700 i__2 = *m;
00701 for (je = 1; je <= i__2; ++je) {
00702 if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
00703 iw = je;
00704 wkill = w[je];
00705 }
00706
00707 }
00708 iblock[iw] = 0;
00709
00710 }
00711 }
00712 im = 0;
00713 i__1 = *m;
00714 for (je = 1; je <= i__1; ++je) {
00715 if (iblock[je] != 0) {
00716 ++im;
00717 w[im] = w[je];
00718 iblock[im] = iblock[je];
00719 }
00720
00721 }
00722 *m = im;
00723 }
00724 if (idiscl < 0 || idiscu < 0) {
00725 toofew = TRUE_;
00726 }
00727 }
00728
00729
00730
00731
00732
00733 if (iorder == 1 && *nsplit > 1) {
00734 i__1 = *m - 1;
00735 for (je = 1; je <= i__1; ++je) {
00736 ie = 0;
00737 tmp1 = w[je];
00738 i__2 = *m;
00739 for (j = je + 1; j <= i__2; ++j) {
00740 if (w[j] < tmp1) {
00741 ie = j;
00742 tmp1 = w[j];
00743 }
00744
00745 }
00746
00747 if (ie != 0) {
00748 itmp1 = iblock[ie];
00749 w[ie] = w[je];
00750 iblock[ie] = iblock[je];
00751 w[je] = tmp1;
00752 iblock[je] = itmp1;
00753 }
00754
00755 }
00756 }
00757
00758 *info = 0;
00759 if (ncnvrg) {
00760 ++(*info);
00761 }
00762 if (toofew) {
00763 *info += 2;
00764 }
00765 return 0;
00766
00767
00768
00769 }
00770
00771 #endif