3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static integer c__3 = 3;
8 static integer c__2 = 2;
9 static integer c__0 = 0;
11 /* Subroutine */ int dstebz_(char *range, char *order, integer *n, doublereal
12 *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol,
13 doublereal *d__, doublereal *e, integer *m, integer *nsplit,
14 doublereal *w, integer *iblock, integer *isplit, doublereal *work,
15 integer *iwork, integer *info)
17 /* System generated locals */
18 integer i__1, i__2, i__3;
19 doublereal d__1, d__2, d__3, d__4, d__5;
21 /* Builtin functions */
22 double sqrt(doublereal), log(doublereal);
25 integer j, ib, jb, ie, je, nb;
32 doublereal ulp, wlu, wul;
34 doublereal tmp1, tmp2;
35 integer iend, ioff, iout, itmp1, jdisc;
36 extern logical lsame_(char *, char *);
42 doublereal wkill, rtoli, tnorm;
43 extern doublereal dlamch_(char *);
45 extern /* Subroutine */ int dlaebz_(integer *, integer *, integer *,
46 integer *, integer *, integer *, doublereal *, doublereal *,
47 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
48 doublereal *, doublereal *, integer *, integer *, doublereal *,
49 integer *, integer *);
50 integer irange, idiscl;
53 extern /* Subroutine */ int xerbla_(char *, integer *);
54 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
55 integer *, integer *);
56 integer idiscu, iorder;
62 /* -- LAPACK routine (version 3.1) -- */
63 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
65 /* 8-18-00: Increase FUDGE factor for T3E (eca) */
67 /* .. Scalar Arguments .. */
69 /* .. Array Arguments .. */
75 /* DSTEBZ computes the eigenvalues of a symmetric tridiagonal */
76 /* matrix T. The user may ask for all eigenvalues, all eigenvalues */
77 /* in the half-open interval (VL, VU], or the IL-th through IU-th */
80 /* To avoid overflow, the matrix must be scaled so that its */
81 /* largest element is no greater than overflow**(1/2) * */
82 /* underflow**(1/4) in absolute value, and for greatest */
83 /* accuracy, it should not be much smaller than that. */
85 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
86 /* Matrix", Report CS41, Computer Science Dept., Stanford */
87 /* University, July 21, 1966. */
92 /* RANGE (input) CHARACTER*1 */
93 /* = 'A': ("All") all eigenvalues will be found. */
94 /* = 'V': ("Value") all eigenvalues in the half-open interval */
95 /* (VL, VU] will be found. */
96 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
97 /* entire matrix) will be found. */
99 /* ORDER (input) CHARACTER*1 */
100 /* = 'B': ("By Block") the eigenvalues will be grouped by */
101 /* split-off block (see IBLOCK, ISPLIT) and */
102 /* ordered from smallest to largest within */
104 /* = 'E': ("Entire matrix") */
105 /* the eigenvalues for the entire matrix */
106 /* will be ordered from smallest to */
109 /* N (input) INTEGER */
110 /* The order of the tridiagonal matrix T. N >= 0. */
112 /* VL (input) DOUBLE PRECISION */
113 /* VU (input) DOUBLE PRECISION */
114 /* If RANGE='V', the lower and upper bounds of the interval to */
115 /* be searched for eigenvalues. Eigenvalues less than or equal */
116 /* to VL, or greater than VU, will not be returned. VL < VU. */
117 /* Not referenced if RANGE = 'A' or 'I'. */
119 /* IL (input) INTEGER */
120 /* IU (input) INTEGER */
121 /* If RANGE='I', the indices (in ascending order) of the */
122 /* smallest and largest eigenvalues to be returned. */
123 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
124 /* Not referenced if RANGE = 'A' or 'V'. */
126 /* ABSTOL (input) DOUBLE PRECISION */
127 /* The absolute tolerance for the eigenvalues. An eigenvalue */
128 /* (or cluster) is considered to be located if it has been */
129 /* determined to lie in an interval whose width is ABSTOL or */
130 /* less. If ABSTOL is less than or equal to zero, then ULP*|T| */
131 /* will be used, where |T| means the 1-norm of T. */
133 /* Eigenvalues will be computed most accurately when ABSTOL is */
134 /* set to twice the underflow threshold 2*DLAMCH('S'), not zero. */
136 /* D (input) DOUBLE PRECISION array, dimension (N) */
137 /* The n diagonal elements of the tridiagonal matrix T. */
139 /* E (input) DOUBLE PRECISION array, dimension (N-1) */
140 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */
142 /* M (output) INTEGER */
143 /* The actual number of eigenvalues found. 0 <= M <= N. */
144 /* (See also the description of INFO=2,3.) */
146 /* NSPLIT (output) INTEGER */
147 /* The number of diagonal blocks in the matrix T. */
148 /* 1 <= NSPLIT <= N. */
150 /* W (output) DOUBLE PRECISION array, dimension (N) */
151 /* On exit, the first M elements of W will contain the */
152 /* eigenvalues. (DSTEBZ may use the remaining N-M elements as */
155 /* IBLOCK (output) INTEGER array, dimension (N) */
156 /* At each row/column j where E(j) is zero or small, the */
157 /* matrix T is considered to split into a block diagonal */
158 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */
159 /* block (from 1 to the number of blocks) the eigenvalue W(i) */
160 /* belongs. (DSTEBZ may use the remaining N-M elements as */
163 /* ISPLIT (output) INTEGER array, dimension (N) */
164 /* The splitting points, at which T breaks up into submatrices. */
165 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */
166 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
167 /* etc., and the NSPLIT-th consists of rows/columns */
168 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
169 /* (Only the first NSPLIT elements will actually be used, but */
170 /* since the user cannot know a priori what value NSPLIT will */
171 /* have, N words must be reserved for ISPLIT.) */
173 /* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
175 /* IWORK (workspace) INTEGER array, dimension (3*N) */
177 /* INFO (output) INTEGER */
178 /* = 0: successful exit */
179 /* < 0: if INFO = -i, the i-th argument had an illegal value */
180 /* > 0: some or all of the eigenvalues failed to converge or */
181 /* were not computed: */
182 /* =1 or 3: Bisection failed to converge for some */
183 /* eigenvalues; these eigenvalues are flagged by a */
184 /* negative block number. The effect is that the */
185 /* eigenvalues may not be as accurate as the */
186 /* absolute and relative tolerances. This is */
187 /* generally caused by unexpectedly inaccurate */
189 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */
190 /* IL:IU were found. */
191 /* Effect: M < IU+1-IL */
192 /* Cause: non-monotonic arithmetic, causing the */
193 /* Sturm sequence to be non-monotonic. */
194 /* Cure: recalculate, using RANGE='A', and pick */
195 /* out eigenvalues IL:IU. In some cases, */
196 /* increasing the PARAMETER "FUDGE" may */
197 /* make things work. */
198 /* = 4: RANGE='I', and the Gershgorin interval */
199 /* initially used was too small. No eigenvalues */
201 /* Probable cause: your machine has sloppy */
202 /* floating-point arithmetic. */
203 /* Cure: Increase the PARAMETER "FUDGE", */
204 /* recompile, and try again. */
206 /* Internal Parameters */
207 /* =================== */
209 /* RELFAC DOUBLE PRECISION, default = 2.0e0 */
210 /* The relative tolerance. An interval (a,b] lies within */
211 /* "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), */
212 /* where "ulp" is the machine precision (distance from 1 to */
213 /* the next larger floating point number.) */
215 /* FUDGE DOUBLE PRECISION, default = 2 */
216 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */
217 /* a value of 1 should work, but on machines with sloppy */
218 /* arithmetic, this needs to be larger. The default for */
219 /* publicly released versions should be large enough to handle */
220 /* the worst machine around. Note that this has no effect */
221 /* on accuracy of the solution. */
223 /* ===================================================================== */
225 /* .. Parameters .. */
227 /* .. Local Scalars .. */
229 /* .. Local Arrays .. */
231 /* .. External Functions .. */
233 /* .. External Subroutines .. */
235 /* .. Intrinsic Functions .. */
237 /* .. Executable Statements .. */
239 /* Parameter adjustments */
253 if (lsame_(range, "A")) {
255 } else if (lsame_(range, "V")) {
257 } else if (lsame_(range, "I")) {
265 if (lsame_(order, "B")) {
267 } else if (lsame_(order, "E")) {
273 /* Check for Errors */
277 } else if (iorder <= 0) {
281 } else if (irange == 2) {
285 } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) {
287 } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) {
293 xerbla_("DSTEBZ", &i__1);
297 /* Initialize error flags */
303 /* Quick return if possible */
310 /* Simplifications: */
312 if (irange == 3 && *il == 1 && *iu == *n) {
316 /* Get machine constants */
317 /* NB is the minimum vector length for vector bisection, or 0 */
318 /* if only scalar is to be done. */
320 safemn = dlamch_("S");
323 nb = ilaenv_(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
328 /* Special Case when N=1 */
333 if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
343 /* Compute Splitting Points */
351 for (j = 2; j <= i__1; ++j) {
352 /* Computing 2nd power */
355 /* Computing 2nd power */
357 if ((d__1 = d__[j] * d__[j - 1], abs(d__1)) * (d__2 * d__2) + safemn
359 isplit[*nsplit] = j - 1;
364 pivmin = max(pivmin,tmp1);
368 isplit[*nsplit] = *n;
371 /* Compute Interval and ATOLI */
375 /* RANGE='I': Compute the interval containing eigenvalues */
378 /* Compute Gershgorin interval for entire (split) matrix */
379 /* and use it as the initial interval */
386 for (j = 1; j <= i__1; ++j) {
387 tmp2 = sqrt(work[j]);
389 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
392 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
399 d__1 = gu, d__2 = d__[*n] + tmp1;
402 d__1 = gl, d__2 = d__[*n] - tmp1;
405 d__1 = abs(gl), d__2 = abs(gu);
406 tnorm = max(d__1,d__2);
407 gl = gl - tnorm * 2.1 * ulp * *n - pivmin * 4.2000000000000002;
408 gu = gu + tnorm * 2.1 * ulp * *n + pivmin * 2.1;
410 /* Compute Iteration parameters */
412 itmax = (integer) ((log(tnorm + pivmin) - log(pivmin)) / log(2.)) + 2;
432 dlaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
433 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
434 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
436 if (iwork[6] == *iu) {
452 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
458 /* RANGE='A' or 'V' -- Set ATOLI */
461 d__3 = abs(d__[1]) + abs(e[1]), d__4 = (d__1 = d__[*n], abs(d__1)) + (
462 d__2 = e[*n - 1], abs(d__2));
463 tnorm = max(d__3,d__4);
466 for (j = 2; j <= i__1; ++j) {
468 d__4 = tnorm, d__5 = (d__1 = d__[j], abs(d__1)) + (d__2 = e[j - 1]
469 , abs(d__2)) + (d__3 = e[j], abs(d__3));
470 tnorm = max(d__4,d__5);
489 /* Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. */
490 /* NWL accumulates the number of eigenvalues .le. WL, */
491 /* NWU accumulates the number of eigenvalues .le. WU */
500 for (jb = 1; jb <= i__1; ++jb) {
508 /* Special Case -- IN=1 */
510 if (irange == 1 || wl >= d__[ibegin] - pivmin) {
513 if (irange == 1 || wu >= d__[ibegin] - pivmin) {
516 if (irange == 1 || wl < d__[ibegin] - pivmin && wu >= d__[ibegin]
524 /* General Case -- IN > 1 */
526 /* Compute Gershgorin Interval */
527 /* and use it as the initial interval */
534 for (j = ibegin; j <= i__2; ++j) {
535 tmp2 = (d__1 = e[j], abs(d__1));
537 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
540 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
547 d__1 = gu, d__2 = d__[iend] + tmp1;
550 d__1 = gl, d__2 = d__[iend] - tmp1;
553 d__1 = abs(gl), d__2 = abs(gu);
554 bnorm = max(d__1,d__2);
555 gl = gl - bnorm * 2.1 * ulp * in - pivmin * 2.1;
556 gu = gu + bnorm * 2.1 * ulp * in + pivmin * 2.1;
558 /* Compute ATOLI for the current submatrix */
562 d__1 = abs(gl), d__2 = abs(gu);
563 atoli = ulp * max(d__1,d__2);
581 /* Set Up Initial Interval */
584 work[*n + in + 1] = gu;
585 dlaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
586 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
587 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
588 w[*m + 1], &iblock[*m + 1], &iinfo);
591 nwu += iwork[in + 1];
592 iwoff = *m - iwork[1];
594 /* Compute Eigenvalues */
596 itmax = (integer) ((log(gu - gl + pivmin) - log(pivmin)) / log(2.)
598 dlaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
599 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
600 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
601 &w[*m + 1], &iblock[*m + 1], &iinfo);
603 /* Copy Eigenvalues Into W and IBLOCK */
604 /* Use -JB for block number for unconverged eigenvalues. */
607 for (j = 1; j <= i__2; ++j) {
608 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
610 /* Flag non-convergence. */
612 if (j > iout - iinfo) {
618 i__3 = iwork[j + in] + iwoff;
619 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
633 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
634 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
638 idiscl = *il - 1 - nwl;
641 if (idiscl > 0 || idiscu > 0) {
643 for (je = 1; je <= i__1; ++je) {
644 if (w[je] <= wlu && idiscl > 0) {
646 } else if (w[je] >= wul && idiscu > 0) {
651 iblock[im] = iblock[je];
657 if (idiscl > 0 || idiscu > 0) {
659 /* Code to deal with effects of bad arithmetic: */
660 /* Some low eigenvalues to be discarded are not in (WL,WLU], */
661 /* or high eigenvalues to be discarded are not in (WUL,WU] */
662 /* so just kill off the smallest IDISCL/largest IDISCU */
663 /* eigenvalues, by simply finding the smallest/largest */
666 /* (If N(w) is monotone non-decreasing, this should never */
672 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
675 for (je = 1; je <= i__2; ++je) {
676 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
690 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
693 for (je = 1; je <= i__2; ++je) {
694 if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
706 for (je = 1; je <= i__1; ++je) {
707 if (iblock[je] != 0) {
710 iblock[im] = iblock[je];
716 if (idiscl < 0 || idiscu < 0) {
721 /* If ORDER='B', do nothing -- the eigenvalues are already sorted */
723 /* If ORDER='E', sort the eigenvalues from smallest to largest */
725 if (iorder == 1 && *nsplit > 1) {
727 for (je = 1; je <= i__1; ++je) {
731 for (j = je + 1; j <= i__2; ++j) {
742 iblock[ie] = iblock[je];