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