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 dlarrd_(char *range, char *order, integer *n, doublereal
12 *vl, doublereal *vu, integer *il, integer *iu, doublereal *gers,
13 doublereal *reltol, doublereal *d__, doublereal *e, doublereal *e2,
14 doublereal *pivmin, integer *nsplit, integer *isplit, integer *m,
15 doublereal *w, doublereal *werr, doublereal *wl, doublereal *wu,
16 integer *iblock, integer *indexw, doublereal *work, integer *iwork,
19 /* System generated locals */
20 integer i__1, i__2, i__3;
21 doublereal d__1, d__2;
23 /* Builtin functions */
24 double log(doublereal);
27 integer i__, j, ib, ie, je, nb;
36 doublereal tmp1, tmp2;
37 integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc;
38 extern logical lsame_(char *, char *);
42 doublereal wkill, rtoli, uflow, 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, idumma[1];
52 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
53 integer *, integer *);
55 logical ncnvrg, toofew;
58 /* -- LAPACK auxiliary routine (version 3.1) -- */
59 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
62 /* .. Scalar Arguments .. */
64 /* .. Array Arguments .. */
70 /* DLARRD computes the eigenvalues of a symmetric tridiagonal */
71 /* matrix T to suitable accuracy. This is an auxiliary code to be */
72 /* called from DSTEMR. */
73 /* The user may ask for all eigenvalues, all eigenvalues */
74 /* in the half-open interval (VL, VU], or the IL-th through IU-th */
77 /* To avoid overflow, the matrix must be scaled so that its */
78 /* largest element is no greater than overflow**(1/2) * */
79 /* underflow**(1/4) in absolute value, and for greatest */
80 /* accuracy, it should not be much smaller than that. */
82 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
83 /* Matrix", Report CS41, Computer Science Dept., Stanford */
84 /* University, July 21, 1966. */
89 /* RANGE (input) CHARACTER */
90 /* = 'A': ("All") all eigenvalues will be found. */
91 /* = 'V': ("Value") all eigenvalues in the half-open interval */
92 /* (VL, VU] will be found. */
93 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
94 /* entire matrix) will be found. */
96 /* ORDER (input) CHARACTER */
97 /* = 'B': ("By Block") the eigenvalues will be grouped by */
98 /* split-off block (see IBLOCK, ISPLIT) and */
99 /* ordered from smallest to largest within */
101 /* = 'E': ("Entire matrix") */
102 /* the eigenvalues for the entire matrix */
103 /* will be ordered from smallest to */
106 /* N (input) INTEGER */
107 /* The order of the tridiagonal matrix T. N >= 0. */
109 /* VL (input) DOUBLE PRECISION */
110 /* VU (input) DOUBLE PRECISION */
111 /* If RANGE='V', the lower and upper bounds of the interval to */
112 /* be searched for eigenvalues. Eigenvalues less than or equal */
113 /* to VL, or greater than VU, will not be returned. VL < VU. */
114 /* Not referenced if RANGE = 'A' or 'I'. */
116 /* IL (input) INTEGER */
117 /* IU (input) INTEGER */
118 /* If RANGE='I', the indices (in ascending order) of the */
119 /* smallest and largest eigenvalues to be returned. */
120 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
121 /* Not referenced if RANGE = 'A' or 'V'. */
123 /* GERS (input) DOUBLE PRECISION array, dimension (2*N) */
124 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */
125 /* is (GERS(2*i-1), GERS(2*i)). */
127 /* RELTOL (input) DOUBLE PRECISION */
128 /* The minimum relative width of an interval. When an interval */
129 /* is narrower than RELTOL times the larger (in */
130 /* magnitude) endpoint, then it is considered to be */
131 /* sufficiently small, i.e., converged. Note: this should */
132 /* always be at least radix*machine epsilon. */
134 /* D (input) DOUBLE PRECISION array, dimension (N) */
135 /* The n diagonal elements of the tridiagonal matrix T. */
137 /* E (input) DOUBLE PRECISION array, dimension (N-1) */
138 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */
140 /* E2 (input) DOUBLE PRECISION array, dimension (N-1) */
141 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
143 /* PIVMIN (input) DOUBLE PRECISION */
144 /* The minimum pivot allowed in the Sturm sequence for T. */
146 /* NSPLIT (input) INTEGER */
147 /* The number of diagonal blocks in the matrix T. */
148 /* 1 <= NSPLIT <= N. */
150 /* ISPLIT (input) INTEGER array, dimension (N) */
151 /* The splitting points, at which T breaks up into submatrices. */
152 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */
153 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
154 /* etc., and the NSPLIT-th consists of rows/columns */
155 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
156 /* (Only the first NSPLIT elements will actually be used, but */
157 /* since the user cannot know a priori what value NSPLIT will */
158 /* have, N words must be reserved for ISPLIT.) */
160 /* M (output) INTEGER */
161 /* The actual number of eigenvalues found. 0 <= M <= N. */
162 /* (See also the description of INFO=2,3.) */
164 /* W (output) DOUBLE PRECISION array, dimension (N) */
165 /* On exit, the first M elements of W will contain the */
166 /* eigenvalue approximations. DLARRD computes an interval */
167 /* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */
168 /* approximation is given as the interval midpoint */
169 /* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */
170 /* WERR(j) = abs( a_j - b_j)/2 */
172 /* WERR (output) DOUBLE PRECISION array, dimension (N) */
173 /* The error bound on the corresponding eigenvalue approximation */
176 /* WL (output) DOUBLE PRECISION */
177 /* WU (output) DOUBLE PRECISION */
178 /* The interval (WL, WU] contains all the wanted eigenvalues. */
179 /* If RANGE='V', then WL=VL and WU=VU. */
180 /* If RANGE='A', then WL and WU are the global Gerschgorin bounds */
181 /* on the spectrum. */
182 /* If RANGE='I', then WL and WU are computed by DLAEBZ from the */
183 /* index range specified. */
185 /* IBLOCK (output) INTEGER array, dimension (N) */
186 /* At each row/column j where E(j) is zero or small, the */
187 /* matrix T is considered to split into a block diagonal */
188 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */
189 /* block (from 1 to the number of blocks) the eigenvalue W(i) */
190 /* belongs. (DLARRD may use the remaining N-M elements as */
193 /* INDEXW (output) INTEGER array, dimension (N) */
194 /* The indices of the eigenvalues within each block (submatrix); */
195 /* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */
196 /* i-th eigenvalue W(i) is the j-th eigenvalue in block k. */
198 /* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
200 /* IWORK (workspace) INTEGER array, dimension (3*N) */
202 /* INFO (output) INTEGER */
203 /* = 0: successful exit */
204 /* < 0: if INFO = -i, the i-th argument had an illegal value */
205 /* > 0: some or all of the eigenvalues failed to converge or */
206 /* were not computed: */
207 /* =1 or 3: Bisection failed to converge for some */
208 /* eigenvalues; these eigenvalues are flagged by a */
209 /* negative block number. The effect is that the */
210 /* eigenvalues may not be as accurate as the */
211 /* absolute and relative tolerances. This is */
212 /* generally caused by unexpectedly inaccurate */
214 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */
215 /* IL:IU were found. */
216 /* Effect: M < IU+1-IL */
217 /* Cause: non-monotonic arithmetic, causing the */
218 /* Sturm sequence to be non-monotonic. */
219 /* Cure: recalculate, using RANGE='A', and pick */
220 /* out eigenvalues IL:IU. In some cases, */
221 /* increasing the PARAMETER "FUDGE" may */
222 /* make things work. */
223 /* = 4: RANGE='I', and the Gershgorin interval */
224 /* initially used was too small. No eigenvalues */
226 /* Probable cause: your machine has sloppy */
227 /* floating-point arithmetic. */
228 /* Cure: Increase the PARAMETER "FUDGE", */
229 /* recompile, and try again. */
231 /* Internal Parameters */
232 /* =================== */
234 /* FUDGE DOUBLE PRECISION, default = 2 */
235 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */
236 /* a value of 1 should work, but on machines with sloppy */
237 /* arithmetic, this needs to be larger. The default for */
238 /* publicly released versions should be large enough to handle */
239 /* the worst machine around. Note that this has no effect */
240 /* on accuracy of the solution. */
242 /* Based on contributions by */
243 /* W. Kahan, University of California, Berkeley, USA */
244 /* Beresford Parlett, University of California, Berkeley, USA */
245 /* Jim Demmel, University of California, Berkeley, USA */
246 /* Inderjit Dhillon, University of Texas, Austin, USA */
247 /* Osni Marques, LBNL/NERSC, USA */
248 /* Christof Voemel, University of California, Berkeley, USA */
250 /* ===================================================================== */
252 /* .. Parameters .. */
254 /* .. Local Scalars .. */
256 /* .. Local Arrays .. */
258 /* .. External Functions .. */
260 /* .. External Subroutines .. */
262 /* .. Intrinsic Functions .. */
264 /* .. Executable Statements .. */
266 /* Parameter adjustments */
284 if (lsame_(range, "A")) {
286 } else if (lsame_(range, "V")) {
288 } else if (lsame_(range, "I")) {
294 /* Check for Errors */
298 } else if (! (lsame_(order, "B") || lsame_(order,
303 } else if (irange == 2) {
307 } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) {
309 } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) {
316 /* Initialize error flags */
320 /* Quick return if possible */
325 /* Simplification: */
326 if (irange == 3 && *il == 1 && *iu == *n) {
329 /* Get machine constants */
331 uflow = dlamch_("U");
332 /* Special Case when N=1 */
333 /* Treat case of 1x1 matrix for quick return */
335 if (irange == 1 || irange == 2 && d__[1] > *vl && d__[1] <= *vu ||
336 irange == 3 && *il == 1 && *iu == 1) {
339 /* The computation error of the eigenvalue is zero */
346 /* NB is the minimum vector length for vector bisection, or 0 */
347 /* if only scalar is to be done. */
348 nb = ilaenv_(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
352 /* Find global spectral radius */
356 for (i__ = 1; i__ <= i__1; ++i__) {
358 d__1 = gl, d__2 = gers[(i__ << 1) - 1];
361 d__1 = gu, d__2 = gers[i__ * 2];
365 /* Compute global Gerschgorin bounds and spectral diameter */
367 d__1 = abs(gl), d__2 = abs(gu);
368 tnorm = max(d__1,d__2);
369 gl = gl - tnorm * 2. * eps * *n - *pivmin * 4.;
370 gu = gu + tnorm * 2. * eps * *n + *pivmin * 4.;
372 /* Input arguments for DLAEBZ: */
373 /* The relative tolerance. An interval (a,b] lies within */
374 /* "relative tolerance" if b-a < RELTOL*max(|a|,|b|), */
376 /* Set the absolute tolerance for interval convergence to zero to force */
377 /* interval convergence based on relative size of the interval. */
378 /* This is dangerous because intervals might not converge when RELTOL is */
379 /* small. But at least a very small number should be selected so that for */
380 /* strongly graded matrices, the code can get relatively accurate */
382 atoli = uflow * 4. + *pivmin * 4.;
384 /* RANGE='I': Compute an interval containing eigenvalues */
385 /* IL through IU. The initial interval [GL,GU] from the global */
386 /* Gerschgorin bounds GL and GU is refined by DLAEBZ. */
387 itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.)) +
402 dlaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, &
403 d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5]
404 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
409 /* On exit, output intervals may not be ordered by ascending negcount */
410 if (iwork[6] == *iu) {
425 /* On exit, the interval [WL, WLU] contains a value with negcount NWL, */
426 /* and [WUL, WU] contains a value with negcount NWU. */
427 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
431 } else if (irange == 2) {
434 } else if (irange == 1) {
438 /* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */
439 /* NWL accumulates the number of eigenvalues .le. WL, */
440 /* NWU accumulates the number of eigenvalues .le. WU */
448 for (jblk = 1; jblk <= i__1; ++jblk) {
456 if (*wl >= d__[ibegin] - *pivmin) {
459 if (*wu >= d__[ibegin] - *pivmin) {
462 if (irange == 1 || *wl < d__[ibegin] - *pivmin && *wu >= d__[
467 /* The gap for a single block doesn't matter for the later */
468 /* algorithm and is assigned an arbitrary large value */
472 /* Disabled 2x2 case because of a failure on the following matrix */
473 /* RANGE = 'I', IL = IU = 4 */
474 /* Original Tridiagonal, d = [ */
475 /* -0.150102010615740E+00 */
476 /* -0.849897989384260E+00 */
477 /* -0.128208148052635E-15 */
478 /* 0.128257718286320E-15 */
481 /* -0.357171383266986E+00 */
482 /* -0.180411241501588E-15 */
483 /* -0.175152352710251E-15 */
486 /* ELSE IF( IN.EQ.2 ) THEN */
488 /* DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */
489 /* TMP1 = HALF*(D(IBEGIN)+D(IEND)) */
490 /* L1 = TMP1 - DISC */
491 /* IF( WL.GE. L1-PIVMIN ) */
492 /* $ NWL = NWL + 1 */
493 /* IF( WU.GE. L1-PIVMIN ) */
494 /* $ NWU = NWU + 1 */
495 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */
496 /* $ L1-PIVMIN ) ) THEN */
499 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */
500 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */
501 /* IBLOCK( M ) = JBLK */
502 /* INDEXW( M ) = 1 */
504 /* L2 = TMP1 + DISC */
505 /* IF( WL.GE. L2-PIVMIN ) */
506 /* $ NWL = NWL + 1 */
507 /* IF( WU.GE. L2-PIVMIN ) */
508 /* $ NWU = NWU + 1 */
509 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */
510 /* $ L2-PIVMIN ) ) THEN */
513 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */
514 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */
515 /* IBLOCK( M ) = JBLK */
516 /* INDEXW( M ) = 2 */
519 /* General Case - block of size IN >= 2 */
520 /* Compute local Gerschgorin interval and use it as the initial */
521 /* interval for DLAEBZ */
526 for (j = ibegin; j <= i__2; ++j) {
528 d__1 = gl, d__2 = gers[(j << 1) - 1];
531 d__1 = gu, d__2 = gers[j * 2];
536 gl = gl - spdiam * 2. * eps * in - *pivmin * 2.;
537 gu = gu + spdiam * 2. * eps * in + *pivmin * 2.;
541 /* the local block contains none of the wanted eigenvalues */
546 /* refine search interval if possible, only range (WL,WU] matters */
553 /* Find negcount of initial interval boundaries GL and GU */
555 work[*n + in + 1] = gu;
556 dlaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli,
557 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
558 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
559 w[*m + 1], &iblock[*m + 1], &iinfo);
566 nwu += iwork[in + 1];
567 iwoff = *m - iwork[1];
568 /* Compute Eigenvalues */
569 itmax = (integer) ((log(gu - gl + *pivmin) - log(*pivmin)) / log(
571 dlaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli,
572 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
573 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
574 &w[*m + 1], &iblock[*m + 1], &iinfo);
580 /* Copy eigenvalues into W and IBLOCK */
581 /* Use -JBLK for block number for unconverged eigenvalues. */
582 /* Loop over the number of output intervals from DLAEBZ */
584 for (j = 1; j <= i__2; ++j) {
585 /* eigenvalue approximation is middle point of interval */
586 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
587 /* semi length of error interval */
588 tmp2 = (d__1 = work[j + *n] - work[j + in + *n], abs(d__1)) *
590 if (j > iout - iinfo) {
591 /* Flag non-convergence. */
597 i__3 = iwork[j + in] + iwoff;
598 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
601 indexw[je] = je - iwoff;
613 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
614 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
616 idiscl = *il - 1 - nwl;
622 for (je = 1; je <= i__1; ++je) {
623 /* Remove some of the smallest eigenvalues from the left so that */
624 /* at the end IDISCL =0. Move all eigenvalues up to the left. */
625 if (w[je] <= wlu && idiscl > 0) {
631 indexw[im] = indexw[je];
632 iblock[im] = iblock[je];
639 /* Remove some of the largest eigenvalues from the right so that */
640 /* at the end IDISCU =0. Move all eigenvalues up to the left. */
642 for (je = *m; je >= 1; --je) {
643 if (w[je] >= wul && idiscu > 0) {
649 indexw[im] = indexw[je];
650 iblock[im] = iblock[je];
656 for (je = im; je <= i__1; ++je) {
659 werr[jee] = werr[je];
660 indexw[jee] = indexw[je];
661 iblock[jee] = iblock[je];
666 if (idiscl > 0 || idiscu > 0) {
667 /* Code to deal with effects of bad arithmetic. (If N(w) is */
668 /* monotone non-decreasing, this should never happen.) */
669 /* Some low eigenvalues to be discarded are not in (WL,WLU], */
670 /* or high eigenvalues to be discarded are not in (WUL,WU] */
671 /* so just kill off the smallest IDISCL/largest IDISCU */
672 /* eigenvalues, by marking the corresponding IBLOCK = 0 */
676 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
679 for (je = 1; je <= i__2; ++je) {
680 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
693 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
696 for (je = 1; je <= i__2; ++je) {
697 if (iblock[je] != 0 && (w[je] >= wkill || iw == 0)) {
707 /* Now erase all eigenvalues with IBLOCK set to zero */
710 for (je = 1; je <= i__1; ++je) {
711 if (iblock[je] != 0) {
715 indexw[im] = indexw[je];
716 iblock[im] = iblock[je];
722 if (idiscl < 0 || idiscu < 0) {
727 if (irange == 1 && *m != *n || irange == 3 && *m != *iu - *il + 1) {
730 /* If ORDER='B', do nothing the eigenvalues are already sorted by */
732 /* If ORDER='E', sort the eigenvalues from smallest to largest */
733 if (lsame_(order, "E") && *nsplit > 1) {
735 for (je = 1; je <= i__1; ++je) {
739 for (j = je + 1; j <= i__2; ++j) {
752 iblock[ie] = iblock[je];
753 indexw[ie] = indexw[je];