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 slarrd_(char *range, char *order, integer *n, real *vl,
12 real *vu, integer *il, integer *iu, real *gers, real *reltol, real *
13 d__, real *e, real *e2, real *pivmin, integer *nsplit, integer *
14 isplit, integer *m, real *w, real *werr, real *wl, real *wu, integer *
15 iblock, integer *indexw, real *work, integer *iwork, integer *info)
17 /* System generated locals */
18 integer i__1, i__2, i__3;
21 /* Builtin functions */
22 double log(doublereal);
25 integer i__, j, ib, ie, je, nb;
35 integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc;
36 extern logical lsame_(char *, char *);
40 real wkill, rtoli, uflow, tnorm;
41 integer ibegin, irange, idiscl;
42 extern doublereal slamch_(char *);
45 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
46 integer *, integer *);
48 extern /* Subroutine */ int slaebz_(integer *, integer *, integer *,
49 integer *, integer *, integer *, real *, real *, real *, real *,
50 real *, real *, integer *, real *, real *, integer *, integer *,
51 real *, integer *, integer *);
52 logical ncnvrg, toofew;
55 /* -- LAPACK auxiliary routine (version 3.1) -- */
56 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
59 /* .. Scalar Arguments .. */
61 /* .. Array Arguments .. */
67 /* SLARRD computes the eigenvalues of a symmetric tridiagonal */
68 /* matrix T to suitable accuracy. This is an auxiliary code to be */
69 /* called from SSTEMR. */
70 /* The user may ask for all eigenvalues, all eigenvalues */
71 /* in the half-open interval (VL, VU], or the IL-th through IU-th */
74 /* To avoid overflow, the matrix must be scaled so that its */
75 /* largest element is no greater than overflow**(1/2) * */
76 /* underflow**(1/4) in absolute value, and for greatest */
77 /* accuracy, it should not be much smaller than that. */
79 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
80 /* Matrix", Report CS41, Computer Science Dept., Stanford */
81 /* University, July 21, 1966. */
86 /* RANGE (input) CHARACTER */
87 /* = 'A': ("All") all eigenvalues will be found. */
88 /* = 'V': ("Value") all eigenvalues in the half-open interval */
89 /* (VL, VU] will be found. */
90 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
91 /* entire matrix) will be found. */
93 /* ORDER (input) CHARACTER */
94 /* = 'B': ("By Block") the eigenvalues will be grouped by */
95 /* split-off block (see IBLOCK, ISPLIT) and */
96 /* ordered from smallest to largest within */
98 /* = 'E': ("Entire matrix") */
99 /* the eigenvalues for the entire matrix */
100 /* will be ordered from smallest to */
103 /* N (input) INTEGER */
104 /* The order of the tridiagonal matrix T. N >= 0. */
106 /* VL (input) REAL */
107 /* VU (input) REAL */
108 /* If RANGE='V', the lower and upper bounds of the interval to */
109 /* be searched for eigenvalues. Eigenvalues less than or equal */
110 /* to VL, or greater than VU, will not be returned. VL < VU. */
111 /* Not referenced if RANGE = 'A' or 'I'. */
113 /* IL (input) INTEGER */
114 /* IU (input) INTEGER */
115 /* If RANGE='I', the indices (in ascending order) of the */
116 /* smallest and largest eigenvalues to be returned. */
117 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
118 /* Not referenced if RANGE = 'A' or 'V'. */
120 /* GERS (input) REAL array, dimension (2*N) */
121 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */
122 /* is (GERS(2*i-1), GERS(2*i)). */
124 /* RELTOL (input) REAL */
125 /* The minimum relative width of an interval. When an interval */
126 /* is narrower than RELTOL times the larger (in */
127 /* magnitude) endpoint, then it is considered to be */
128 /* sufficiently small, i.e., converged. Note: this should */
129 /* always be at least radix*machine epsilon. */
131 /* D (input) REAL array, dimension (N) */
132 /* The n diagonal elements of the tridiagonal matrix T. */
134 /* E (input) REAL array, dimension (N-1) */
135 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */
137 /* E2 (input) REAL array, dimension (N-1) */
138 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
140 /* PIVMIN (input) REAL */
141 /* The minimum pivot allowed in the Sturm sequence for T. */
143 /* NSPLIT (input) INTEGER */
144 /* The number of diagonal blocks in the matrix T. */
145 /* 1 <= NSPLIT <= N. */
147 /* ISPLIT (input) INTEGER array, dimension (N) */
148 /* The splitting points, at which T breaks up into submatrices. */
149 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */
150 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
151 /* etc., and the NSPLIT-th consists of rows/columns */
152 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
153 /* (Only the first NSPLIT elements will actually be used, but */
154 /* since the user cannot know a priori what value NSPLIT will */
155 /* have, N words must be reserved for ISPLIT.) */
157 /* M (output) INTEGER */
158 /* The actual number of eigenvalues found. 0 <= M <= N. */
159 /* (See also the description of INFO=2,3.) */
161 /* W (output) REAL array, dimension (N) */
162 /* On exit, the first M elements of W will contain the */
163 /* eigenvalue approximations. SLARRD computes an interval */
164 /* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */
165 /* approximation is given as the interval midpoint */
166 /* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */
167 /* WERR(j) = abs( a_j - b_j)/2 */
169 /* WERR (output) REAL array, dimension (N) */
170 /* The error bound on the corresponding eigenvalue approximation */
173 /* WL (output) REAL */
174 /* WU (output) REAL */
175 /* The interval (WL, WU] contains all the wanted eigenvalues. */
176 /* If RANGE='V', then WL=VL and WU=VU. */
177 /* If RANGE='A', then WL and WU are the global Gerschgorin bounds */
178 /* on the spectrum. */
179 /* If RANGE='I', then WL and WU are computed by SLAEBZ from the */
180 /* index range specified. */
182 /* IBLOCK (output) INTEGER array, dimension (N) */
183 /* At each row/column j where E(j) is zero or small, the */
184 /* matrix T is considered to split into a block diagonal */
185 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */
186 /* block (from 1 to the number of blocks) the eigenvalue W(i) */
187 /* belongs. (SLARRD may use the remaining N-M elements as */
190 /* INDEXW (output) INTEGER array, dimension (N) */
191 /* The indices of the eigenvalues within each block (submatrix); */
192 /* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */
193 /* i-th eigenvalue W(i) is the j-th eigenvalue in block k. */
195 /* WORK (workspace) REAL array, dimension (4*N) */
197 /* IWORK (workspace) INTEGER array, dimension (3*N) */
199 /* INFO (output) INTEGER */
200 /* = 0: successful exit */
201 /* < 0: if INFO = -i, the i-th argument had an illegal value */
202 /* > 0: some or all of the eigenvalues failed to converge or */
203 /* were not computed: */
204 /* =1 or 3: Bisection failed to converge for some */
205 /* eigenvalues; these eigenvalues are flagged by a */
206 /* negative block number. The effect is that the */
207 /* eigenvalues may not be as accurate as the */
208 /* absolute and relative tolerances. This is */
209 /* generally caused by unexpectedly inaccurate */
211 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */
212 /* IL:IU were found. */
213 /* Effect: M < IU+1-IL */
214 /* Cause: non-monotonic arithmetic, causing the */
215 /* Sturm sequence to be non-monotonic. */
216 /* Cure: recalculate, using RANGE='A', and pick */
217 /* out eigenvalues IL:IU. In some cases, */
218 /* increasing the PARAMETER "FUDGE" may */
219 /* make things work. */
220 /* = 4: RANGE='I', and the Gershgorin interval */
221 /* initially used was too small. No eigenvalues */
223 /* Probable cause: your machine has sloppy */
224 /* floating-point arithmetic. */
225 /* Cure: Increase the PARAMETER "FUDGE", */
226 /* recompile, and try again. */
228 /* Internal Parameters */
229 /* =================== */
231 /* FUDGE REAL , default = 2 */
232 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */
233 /* a value of 1 should work, but on machines with sloppy */
234 /* arithmetic, this needs to be larger. The default for */
235 /* publicly released versions should be large enough to handle */
236 /* the worst machine around. Note that this has no effect */
237 /* on accuracy of the solution. */
239 /* Based on contributions by */
240 /* W. Kahan, University of California, Berkeley, USA */
241 /* Beresford Parlett, University of California, Berkeley, USA */
242 /* Jim Demmel, University of California, Berkeley, USA */
243 /* Inderjit Dhillon, University of Texas, Austin, USA */
244 /* Osni Marques, LBNL/NERSC, USA */
245 /* Christof Voemel, University of California, Berkeley, USA */
247 /* ===================================================================== */
249 /* .. Parameters .. */
251 /* .. Local Scalars .. */
253 /* .. Local Arrays .. */
255 /* .. External Functions .. */
257 /* .. External Subroutines .. */
259 /* .. Intrinsic Functions .. */
261 /* .. Executable Statements .. */
263 /* Parameter adjustments */
281 if (lsame_(range, "A")) {
283 } else if (lsame_(range, "V")) {
285 } else if (lsame_(range, "I")) {
291 /* Check for Errors */
295 } else if (! (lsame_(order, "B") || lsame_(order,
300 } else if (irange == 2) {
304 } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) {
306 } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) {
313 /* Initialize error flags */
317 /* Quick return if possible */
322 /* Simplification: */
323 if (irange == 3 && *il == 1 && *iu == *n) {
326 /* Get machine constants */
328 uflow = slamch_("U");
329 /* Special Case when N=1 */
330 /* Treat case of 1x1 matrix for quick return */
332 if (irange == 1 || irange == 2 && d__[1] > *vl && d__[1] <= *vu ||
333 irange == 3 && *il == 1 && *iu == 1) {
336 /* The computation error of the eigenvalue is zero */
343 /* NB is the minimum vector length for vector bisection, or 0 */
344 /* if only scalar is to be done. */
345 nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
349 /* Find global spectral radius */
353 for (i__ = 1; i__ <= i__1; ++i__) {
355 r__1 = gl, r__2 = gers[(i__ << 1) - 1];
356 gl = dmin(r__1,r__2);
358 r__1 = gu, r__2 = gers[i__ * 2];
359 gu = dmax(r__1,r__2);
362 /* Compute global Gerschgorin bounds and spectral diameter */
364 r__1 = dabs(gl), r__2 = dabs(gu);
365 tnorm = dmax(r__1,r__2);
366 gl = gl - tnorm * 2.f * eps * *n - *pivmin * 4.f;
367 gu = gu + tnorm * 2.f * eps * *n + *pivmin * 4.f;
369 /* Input arguments for SLAEBZ: */
370 /* The relative tolerance. An interval (a,b] lies within */
371 /* "relative tolerance" if b-a < RELTOL*max(|a|,|b|), */
373 /* Set the absolute tolerance for interval convergence to zero to force */
374 /* interval convergence based on relative size of the interval. */
375 /* This is dangerous because intervals might not converge when RELTOL is */
376 /* small. But at least a very small number should be selected so that for */
377 /* strongly graded matrices, the code can get relatively accurate */
379 atoli = uflow * 4.f + *pivmin * 4.f;
381 /* RANGE='I': Compute an interval containing eigenvalues */
382 /* IL through IU. The initial interval [GL,GU] from the global */
383 /* Gerschgorin bounds GL and GU is refined by SLAEBZ. */
384 itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.f))
399 slaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, &
400 d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5]
401 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
406 /* On exit, output intervals may not be ordered by ascending negcount */
407 if (iwork[6] == *iu) {
422 /* On exit, the interval [WL, WLU] contains a value with negcount NWL, */
423 /* and [WUL, WU] contains a value with negcount NWU. */
424 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
428 } else if (irange == 2) {
431 } else if (irange == 1) {
435 /* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */
436 /* NWL accumulates the number of eigenvalues .le. WL, */
437 /* NWU accumulates the number of eigenvalues .le. WU */
445 for (jblk = 1; jblk <= i__1; ++jblk) {
453 if (*wl >= d__[ibegin] - *pivmin) {
456 if (*wu >= d__[ibegin] - *pivmin) {
459 if (irange == 1 || *wl < d__[ibegin] - *pivmin && *wu >= d__[
464 /* The gap for a single block doesn't matter for the later */
465 /* algorithm and is assigned an arbitrary large value */
469 /* Disabled 2x2 case because of a failure on the following matrix */
470 /* RANGE = 'I', IL = IU = 4 */
471 /* Original Tridiagonal, d = [ */
472 /* -0.150102010615740E+00 */
473 /* -0.849897989384260E+00 */
474 /* -0.128208148052635E-15 */
475 /* 0.128257718286320E-15 */
478 /* -0.357171383266986E+00 */
479 /* -0.180411241501588E-15 */
480 /* -0.175152352710251E-15 */
483 /* ELSE IF( IN.EQ.2 ) THEN */
485 /* DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */
486 /* TMP1 = HALF*(D(IBEGIN)+D(IEND)) */
487 /* L1 = TMP1 - DISC */
488 /* IF( WL.GE. L1-PIVMIN ) */
489 /* $ NWL = NWL + 1 */
490 /* IF( WU.GE. L1-PIVMIN ) */
491 /* $ NWU = NWU + 1 */
492 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */
493 /* $ L1-PIVMIN ) ) THEN */
496 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */
497 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */
498 /* IBLOCK( M ) = JBLK */
499 /* INDEXW( M ) = 1 */
501 /* L2 = TMP1 + DISC */
502 /* IF( WL.GE. L2-PIVMIN ) */
503 /* $ NWL = NWL + 1 */
504 /* IF( WU.GE. L2-PIVMIN ) */
505 /* $ NWU = NWU + 1 */
506 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */
507 /* $ L2-PIVMIN ) ) THEN */
510 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */
511 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */
512 /* IBLOCK( M ) = JBLK */
513 /* INDEXW( M ) = 2 */
516 /* General Case - block of size IN >= 2 */
517 /* Compute local Gerschgorin interval and use it as the initial */
518 /* interval for SLAEBZ */
523 for (j = ibegin; j <= i__2; ++j) {
525 r__1 = gl, r__2 = gers[(j << 1) - 1];
526 gl = dmin(r__1,r__2);
528 r__1 = gu, r__2 = gers[j * 2];
529 gu = dmax(r__1,r__2);
533 gl = gl - spdiam * 2.f * eps * in - *pivmin * 2.f;
534 gu = gu + spdiam * 2.f * eps * in + *pivmin * 2.f;
538 /* the local block contains none of the wanted eigenvalues */
543 /* refine search interval if possible, only range (WL,WU] matters */
550 /* Find negcount of initial interval boundaries GL and GU */
552 work[*n + in + 1] = gu;
553 slaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli,
554 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
555 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
556 w[*m + 1], &iblock[*m + 1], &iinfo);
563 nwu += iwork[in + 1];
564 iwoff = *m - iwork[1];
565 /* Compute Eigenvalues */
566 itmax = (integer) ((log(gu - gl + *pivmin) - log(*pivmin)) / log(
568 slaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli,
569 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
570 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
571 &w[*m + 1], &iblock[*m + 1], &iinfo);
577 /* Copy eigenvalues into W and IBLOCK */
578 /* Use -JBLK for block number for unconverged eigenvalues. */
579 /* Loop over the number of output intervals from SLAEBZ */
581 for (j = 1; j <= i__2; ++j) {
582 /* eigenvalue approximation is middle point of interval */
583 tmp1 = (work[j + *n] + work[j + in + *n]) * .5f;
584 /* semi length of error interval */
585 tmp2 = (r__1 = work[j + *n] - work[j + in + *n], dabs(r__1)) *
587 if (j > iout - iinfo) {
588 /* Flag non-convergence. */
594 i__3 = iwork[j + in] + iwoff;
595 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
598 indexw[je] = je - iwoff;
610 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
611 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
613 idiscl = *il - 1 - nwl;
619 for (je = 1; je <= i__1; ++je) {
620 /* Remove some of the smallest eigenvalues from the left so that */
621 /* at the end IDISCL =0. Move all eigenvalues up to the left. */
622 if (w[je] <= wlu && idiscl > 0) {
628 indexw[im] = indexw[je];
629 iblock[im] = iblock[je];
636 /* Remove some of the largest eigenvalues from the right so that */
637 /* at the end IDISCU =0. Move all eigenvalues up to the left. */
639 for (je = *m; je >= 1; --je) {
640 if (w[je] >= wul && idiscu > 0) {
646 indexw[im] = indexw[je];
647 iblock[im] = iblock[je];
653 for (je = im; je <= i__1; ++je) {
656 werr[jee] = werr[je];
657 indexw[jee] = indexw[je];
658 iblock[jee] = iblock[je];
663 if (idiscl > 0 || idiscu > 0) {
664 /* Code to deal with effects of bad arithmetic. (If N(w) is */
665 /* monotone non-decreasing, this should never happen.) */
666 /* Some low eigenvalues to be discarded are not in (WL,WLU], */
667 /* or high eigenvalues to be discarded are not in (WUL,WU] */
668 /* so just kill off the smallest IDISCL/largest IDISCU */
669 /* eigenvalues, by marking the corresponding IBLOCK = 0 */
673 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
676 for (je = 1; je <= i__2; ++je) {
677 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)) {
704 /* Now erase all eigenvalues with IBLOCK set to zero */
707 for (je = 1; je <= i__1; ++je) {
708 if (iblock[je] != 0) {
712 indexw[im] = indexw[je];
713 iblock[im] = iblock[je];
719 if (idiscl < 0 || idiscu < 0) {
724 if (irange == 1 && *m != *n || irange == 3 && *m != *iu - *il + 1) {
727 /* If ORDER='B', do nothing the eigenvalues are already sorted by */
729 /* If ORDER='E', sort the eigenvalues from smallest to largest */
730 if (lsame_(order, "E") && *nsplit > 1) {
732 for (je = 1; je <= i__1; ++je) {
736 for (j = je + 1; j <= i__2; ++j) {
749 iblock[ie] = iblock[je];
750 indexw[ie] = indexw[je];