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 sstebz_(char *range, char *order, integer *n, real *vl,
12 real *vu, integer *il, integer *iu, real *abstol, real *d__, real *e,
13 integer *m, integer *nsplit, real *w, integer *iblock, integer *
14 isplit, real *work, integer *iwork, integer *info)
16 /* System generated locals */
17 integer i__1, i__2, i__3;
18 real r__1, r__2, r__3, r__4, r__5;
20 /* Builtin functions */
21 double sqrt(doublereal), log(doublereal);
24 integer j, ib, jb, ie, je, nb;
34 integer iend, ioff, iout, itmp1, jdisc;
35 extern logical lsame_(char *, char *);
41 real wkill, rtoli, tnorm;
42 integer ibegin, irange, idiscl;
43 extern doublereal slamch_(char *);
46 extern /* Subroutine */ int xerbla_(char *, integer *);
47 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
48 integer *, integer *);
50 extern /* Subroutine */ int slaebz_(integer *, integer *, integer *,
51 integer *, integer *, integer *, real *, real *, real *, real *,
52 real *, real *, integer *, real *, real *, integer *, integer *,
53 real *, integer *, integer *);
60 /* -- LAPACK routine (version 3.1) -- */
61 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
63 /* 8-18-00: Increase FUDGE factor for T3E (eca) */
65 /* .. Scalar Arguments .. */
67 /* .. Array Arguments .. */
73 /* SSTEBZ computes the eigenvalues of a symmetric tridiagonal */
74 /* matrix T. The user may ask for all eigenvalues, all eigenvalues */
75 /* in the half-open interval (VL, VU], or the IL-th through IU-th */
78 /* To avoid overflow, the matrix must be scaled so that its */
79 /* largest element is no greater than overflow**(1/2) * */
80 /* underflow**(1/4) in absolute value, and for greatest */
81 /* accuracy, it should not be much smaller than that. */
83 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
84 /* Matrix", Report CS41, Computer Science Dept., Stanford */
85 /* University, July 21, 1966. */
90 /* RANGE (input) CHARACTER*1 */
91 /* = 'A': ("All") all eigenvalues will be found. */
92 /* = 'V': ("Value") all eigenvalues in the half-open interval */
93 /* (VL, VU] will be found. */
94 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
95 /* entire matrix) will be found. */
97 /* ORDER (input) CHARACTER*1 */
98 /* = 'B': ("By Block") the eigenvalues will be grouped by */
99 /* split-off block (see IBLOCK, ISPLIT) and */
100 /* ordered from smallest to largest within */
102 /* = 'E': ("Entire matrix") */
103 /* the eigenvalues for the entire matrix */
104 /* will be ordered from smallest to */
107 /* N (input) INTEGER */
108 /* The order of the tridiagonal matrix T. N >= 0. */
110 /* VL (input) REAL */
111 /* VU (input) REAL */
112 /* If RANGE='V', the lower and upper bounds of the interval to */
113 /* be searched for eigenvalues. Eigenvalues less than or equal */
114 /* to VL, or greater than VU, will not be returned. VL < VU. */
115 /* Not referenced if RANGE = 'A' or 'I'. */
117 /* IL (input) INTEGER */
118 /* IU (input) INTEGER */
119 /* If RANGE='I', the indices (in ascending order) of the */
120 /* smallest and largest eigenvalues to be returned. */
121 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
122 /* Not referenced if RANGE = 'A' or 'V'. */
124 /* ABSTOL (input) REAL */
125 /* The absolute tolerance for the eigenvalues. An eigenvalue */
126 /* (or cluster) is considered to be located if it has been */
127 /* determined to lie in an interval whose width is ABSTOL or */
128 /* less. If ABSTOL is less than or equal to zero, then ULP*|T| */
129 /* will be used, where |T| means the 1-norm of T. */
131 /* Eigenvalues will be computed most accurately when ABSTOL is */
132 /* set to twice the underflow threshold 2*SLAMCH('S'), not zero. */
134 /* D (input) REAL array, dimension (N) */
135 /* The n diagonal elements of the tridiagonal matrix T. */
137 /* E (input) REAL array, dimension (N-1) */
138 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */
140 /* M (output) INTEGER */
141 /* The actual number of eigenvalues found. 0 <= M <= N. */
142 /* (See also the description of INFO=2,3.) */
144 /* NSPLIT (output) INTEGER */
145 /* The number of diagonal blocks in the matrix T. */
146 /* 1 <= NSPLIT <= N. */
148 /* W (output) REAL array, dimension (N) */
149 /* On exit, the first M elements of W will contain the */
150 /* eigenvalues. (SSTEBZ may use the remaining N-M elements as */
153 /* IBLOCK (output) INTEGER array, dimension (N) */
154 /* At each row/column j where E(j) is zero or small, the */
155 /* matrix T is considered to split into a block diagonal */
156 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */
157 /* block (from 1 to the number of blocks) the eigenvalue W(i) */
158 /* belongs. (SSTEBZ may use the remaining N-M elements as */
161 /* ISPLIT (output) INTEGER array, dimension (N) */
162 /* The splitting points, at which T breaks up into submatrices. */
163 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */
164 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
165 /* etc., and the NSPLIT-th consists of rows/columns */
166 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
167 /* (Only the first NSPLIT elements will actually be used, but */
168 /* since the user cannot know a priori what value NSPLIT will */
169 /* have, N words must be reserved for ISPLIT.) */
171 /* WORK (workspace) REAL array, dimension (4*N) */
173 /* IWORK (workspace) INTEGER array, dimension (3*N) */
175 /* INFO (output) INTEGER */
176 /* = 0: successful exit */
177 /* < 0: if INFO = -i, the i-th argument had an illegal value */
178 /* > 0: some or all of the eigenvalues failed to converge or */
179 /* were not computed: */
180 /* =1 or 3: Bisection failed to converge for some */
181 /* eigenvalues; these eigenvalues are flagged by a */
182 /* negative block number. The effect is that the */
183 /* eigenvalues may not be as accurate as the */
184 /* absolute and relative tolerances. This is */
185 /* generally caused by unexpectedly inaccurate */
187 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */
188 /* IL:IU were found. */
189 /* Effect: M < IU+1-IL */
190 /* Cause: non-monotonic arithmetic, causing the */
191 /* Sturm sequence to be non-monotonic. */
192 /* Cure: recalculate, using RANGE='A', and pick */
193 /* out eigenvalues IL:IU. In some cases, */
194 /* increasing the PARAMETER "FUDGE" may */
195 /* make things work. */
196 /* = 4: RANGE='I', and the Gershgorin interval */
197 /* initially used was too small. No eigenvalues */
199 /* Probable cause: your machine has sloppy */
200 /* floating-point arithmetic. */
201 /* Cure: Increase the PARAMETER "FUDGE", */
202 /* recompile, and try again. */
204 /* Internal Parameters */
205 /* =================== */
207 /* RELFAC REAL, default = 2.0e0 */
208 /* The relative tolerance. An interval (a,b] lies within */
209 /* "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), */
210 /* where "ulp" is the machine precision (distance from 1 to */
211 /* the next larger floating point number.) */
213 /* FUDGE REAL, default = 2 */
214 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */
215 /* a value of 1 should work, but on machines with sloppy */
216 /* arithmetic, this needs to be larger. The default for */
217 /* publicly released versions should be large enough to handle */
218 /* the worst machine around. Note that this has no effect */
219 /* on accuracy of the solution. */
221 /* ===================================================================== */
223 /* .. Parameters .. */
225 /* .. Local Scalars .. */
227 /* .. Local Arrays .. */
229 /* .. External Functions .. */
231 /* .. External Subroutines .. */
233 /* .. Intrinsic Functions .. */
235 /* .. Executable Statements .. */
237 /* Parameter adjustments */
251 if (lsame_(range, "A")) {
253 } else if (lsame_(range, "V")) {
255 } else if (lsame_(range, "I")) {
263 if (lsame_(order, "B")) {
265 } else if (lsame_(order, "E")) {
271 /* Check for Errors */
275 } else if (iorder <= 0) {
279 } else if (irange == 2) {
283 } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) {
285 } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) {
291 xerbla_("SSTEBZ", &i__1);
295 /* Initialize error flags */
301 /* Quick return if possible */
308 /* Simplifications: */
310 if (irange == 3 && *il == 1 && *iu == *n) {
314 /* Get machine constants */
315 /* NB is the minimum vector length for vector bisection, or 0 */
316 /* if only scalar is to be done. */
318 safemn = slamch_("S");
321 nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
326 /* Special Case when N=1 */
331 if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
341 /* Compute Splitting Points */
349 for (j = 2; j <= i__1; ++j) {
350 /* Computing 2nd power */
353 /* Computing 2nd power */
355 if ((r__1 = d__[j] * d__[j - 1], dabs(r__1)) * (r__2 * r__2) + safemn
357 isplit[*nsplit] = j - 1;
362 pivmin = dmax(pivmin,tmp1);
366 isplit[*nsplit] = *n;
369 /* Compute Interval and ATOLI */
373 /* RANGE='I': Compute the interval containing eigenvalues */
376 /* Compute Gershgorin interval for entire (split) matrix */
377 /* and use it as the initial interval */
384 for (j = 1; j <= i__1; ++j) {
385 tmp2 = sqrt(work[j]);
387 r__1 = gu, r__2 = d__[j] + tmp1 + tmp2;
388 gu = dmax(r__1,r__2);
390 r__1 = gl, r__2 = d__[j] - tmp1 - tmp2;
391 gl = dmin(r__1,r__2);
397 r__1 = gu, r__2 = d__[*n] + tmp1;
398 gu = dmax(r__1,r__2);
400 r__1 = gl, r__2 = d__[*n] - tmp1;
401 gl = dmin(r__1,r__2);
403 r__1 = dabs(gl), r__2 = dabs(gu);
404 tnorm = dmax(r__1,r__2);
405 gl = gl - tnorm * 2.1f * ulp * *n - pivmin * 4.2000000000000002f;
406 gu = gu + tnorm * 2.1f * ulp * *n + pivmin * 2.1f;
408 /* Compute Iteration parameters */
410 itmax = (integer) ((log(tnorm + pivmin) - log(pivmin)) / log(2.f)) +
412 if (*abstol <= 0.f) {
431 slaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
432 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
433 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
435 if (iwork[6] == *iu) {
451 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
457 /* RANGE='A' or 'V' -- Set ATOLI */
460 r__3 = dabs(d__[1]) + dabs(e[1]), r__4 = (r__1 = d__[*n], dabs(r__1))
461 + (r__2 = e[*n - 1], dabs(r__2));
462 tnorm = dmax(r__3,r__4);
465 for (j = 2; j <= i__1; ++j) {
467 r__4 = tnorm, r__5 = (r__1 = d__[j], dabs(r__1)) + (r__2 = e[j -
468 1], dabs(r__2)) + (r__3 = e[j], dabs(r__3));
469 tnorm = dmax(r__4,r__5);
473 if (*abstol <= 0.f) {
488 /* Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. */
489 /* NWL accumulates the number of eigenvalues .le. WL, */
490 /* NWU accumulates the number of eigenvalues .le. WU */
499 for (jb = 1; jb <= i__1; ++jb) {
507 /* Special Case -- IN=1 */
509 if (irange == 1 || wl >= d__[ibegin] - pivmin) {
512 if (irange == 1 || wu >= d__[ibegin] - pivmin) {
515 if (irange == 1 || wl < d__[ibegin] - pivmin && wu >= d__[ibegin]
523 /* General Case -- IN > 1 */
525 /* Compute Gershgorin Interval */
526 /* and use it as the initial interval */
533 for (j = ibegin; j <= i__2; ++j) {
534 tmp2 = (r__1 = e[j], dabs(r__1));
536 r__1 = gu, r__2 = d__[j] + tmp1 + tmp2;
537 gu = dmax(r__1,r__2);
539 r__1 = gl, r__2 = d__[j] - tmp1 - tmp2;
540 gl = dmin(r__1,r__2);
546 r__1 = gu, r__2 = d__[iend] + tmp1;
547 gu = dmax(r__1,r__2);
549 r__1 = gl, r__2 = d__[iend] - tmp1;
550 gl = dmin(r__1,r__2);
552 r__1 = dabs(gl), r__2 = dabs(gu);
553 bnorm = dmax(r__1,r__2);
554 gl = gl - bnorm * 2.1f * ulp * in - pivmin * 2.1f;
555 gu = gu + bnorm * 2.1f * ulp * in + pivmin * 2.1f;
557 /* Compute ATOLI for the current submatrix */
559 if (*abstol <= 0.f) {
561 r__1 = dabs(gl), r__2 = dabs(gu);
562 atoli = ulp * dmax(r__1,r__2);
580 /* Set Up Initial Interval */
583 work[*n + in + 1] = gu;
584 slaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
585 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
586 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
587 w[*m + 1], &iblock[*m + 1], &iinfo);
590 nwu += iwork[in + 1];
591 iwoff = *m - iwork[1];
593 /* Compute Eigenvalues */
595 itmax = (integer) ((log(gu - gl + pivmin) - log(pivmin)) / log(
597 slaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
598 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
599 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
600 &w[*m + 1], &iblock[*m + 1], &iinfo);
602 /* Copy Eigenvalues Into W and IBLOCK */
603 /* Use -JB for block number for unconverged eigenvalues. */
606 for (j = 1; j <= i__2; ++j) {
607 tmp1 = (work[j + *n] + work[j + in + *n]) * .5f;
609 /* Flag non-convergence. */
611 if (j > iout - iinfo) {
617 i__3 = iwork[j + in] + iwoff;
618 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
632 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
633 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
637 idiscl = *il - 1 - nwl;
640 if (idiscl > 0 || idiscu > 0) {
642 for (je = 1; je <= i__1; ++je) {
643 if (w[je] <= wlu && idiscl > 0) {
645 } else if (w[je] >= wul && idiscu > 0) {
650 iblock[im] = iblock[je];
656 if (idiscl > 0 || idiscu > 0) {
658 /* Code to deal with effects of bad arithmetic: */
659 /* Some low eigenvalues to be discarded are not in (WL,WLU], */
660 /* or high eigenvalues to be discarded are not in (WUL,WU] */
661 /* so just kill off the smallest IDISCL/largest IDISCU */
662 /* eigenvalues, by simply finding the smallest/largest */
665 /* (If N(w) is monotone non-decreasing, this should never */
671 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
674 for (je = 1; je <= i__2; ++je) {
675 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
689 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
692 for (je = 1; je <= i__2; ++je) {
693 if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
705 for (je = 1; je <= i__1; ++je) {
706 if (iblock[je] != 0) {
709 iblock[im] = iblock[je];
715 if (idiscl < 0 || idiscu < 0) {
720 /* If ORDER='B', do nothing -- the eigenvalues are already sorted */
722 /* If ORDER='E', sort the eigenvalues from smallest to largest */
724 if (iorder == 1 && *nsplit > 1) {
726 for (je = 1; je <= i__1; ++je) {
730 for (j = je + 1; j <= i__2; ++j) {
741 iblock[ie] = iblock[je];