3 /* Table of constant values */
5 static doublereal c_b5 = 0.;
6 static integer c__1 = 1;
7 static integer c__2 = 2;
9 /* Subroutine */ int dlarrv_(integer *n, doublereal *vl, doublereal *vu,
10 doublereal *d__, doublereal *l, doublereal *pivmin, integer *isplit,
11 integer *m, integer *dol, integer *dou, doublereal *minrgp,
12 doublereal *rtol1, doublereal *rtol2, doublereal *w, doublereal *werr,
13 doublereal *wgap, integer *iblock, integer *indexw, doublereal *gers,
14 doublereal *z__, integer *ldz, integer *isuppz, doublereal *work,
15 integer *iwork, integer *info)
17 /* System generated locals */
18 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
19 doublereal d__1, d__2;
22 /* Builtin functions */
23 double log(doublereal);
26 integer minwsize, i__, j, k, p, q, miniwsize, ii;
29 doublereal gu, gap, eps, tau, tol, tmp;
35 doublereal rgap, left;
39 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
49 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
50 doublereal *, integer *);
53 integer iindc1, iindc2;
54 extern /* Subroutine */ int dlar1v_(integer *, integer *, integer *,
55 doublereal *, doublereal *, doublereal *, doublereal *,
56 doublereal *, doublereal *, doublereal *, doublereal *, logical *,
57 integer *, doublereal *, doublereal *, integer *, integer *,
58 doublereal *, doublereal *, doublereal *, doublereal *);
61 extern doublereal dlamch_(char *);
62 integer ibegin, indeig;
65 doublereal sgndef, mingma;
66 extern /* Subroutine */ int dlarrb_(integer *, doublereal *, doublereal *,
67 integer *, integer *, doublereal *, doublereal *, integer *,
68 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
69 doublereal *, doublereal *, integer *, integer *);
70 integer oldien, oldncl, wbegin;
73 extern /* Subroutine */ int dlarrf_(integer *, doublereal *, doublereal *,
74 doublereal *, integer *, integer *, doublereal *, doublereal *,
75 doublereal *, doublereal *, doublereal *, doublereal *,
76 doublereal *, doublereal *, doublereal *, doublereal *,
77 doublereal *, integer *);
82 extern /* Subroutine */ int dlaset_(char *, integer *, integer *,
83 doublereal *, doublereal *, doublereal *, integer *);
85 integer iindwk, offset;
87 integer newcls, oldfst, indwrk, windex, oldlst;
89 integer newfst, newftt, parity, windmn, windpl, isupmn, newlst, zusedl;
91 integer newsiz, zusedu, zusedw;
92 doublereal nrminv, rqcorr;
97 /* -- LAPACK auxiliary routine (version 3.1.1) -- */
98 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
101 /* .. Scalar Arguments .. */
103 /* .. Array Arguments .. */
109 /* DLARRV computes the eigenvectors of the tridiagonal matrix */
110 /* T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T. */
111 /* The input eigenvalues should have been computed by DLARRE. */
116 /* N (input) INTEGER */
117 /* The order of the matrix. N >= 0. */
119 /* VL (input) DOUBLE PRECISION */
120 /* VU (input) DOUBLE PRECISION */
121 /* Lower and upper bounds of the interval that contains the desired */
122 /* eigenvalues. VL < VU. Needed to compute gaps on the left or right */
123 /* end of the extremal eigenvalues in the desired RANGE. */
125 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
126 /* On entry, the N diagonal elements of the diagonal matrix D. */
127 /* On exit, D may be overwritten. */
129 /* L (input/output) DOUBLE PRECISION array, dimension (N) */
130 /* On entry, the (N-1) subdiagonal elements of the unit */
131 /* bidiagonal matrix L are in elements 1 to N-1 of L */
132 /* (if the matrix is not splitted.) At the end of each block */
133 /* is stored the corresponding shift as given by DLARRE. */
134 /* On exit, L is overwritten. */
136 /* PIVMIN (in) DOUBLE PRECISION */
137 /* The minimum pivot allowed in the Sturm sequence. */
139 /* ISPLIT (input) INTEGER array, dimension (N) */
140 /* The splitting points, at which T breaks up into blocks. */
141 /* The first block consists of rows/columns 1 to */
142 /* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
143 /* through ISPLIT( 2 ), etc. */
145 /* M (input) INTEGER */
146 /* The total number of input eigenvalues. 0 <= M <= N. */
148 /* DOL (input) INTEGER */
149 /* DOU (input) INTEGER */
150 /* If the user wants to compute only selected eigenvectors from all */
151 /* the eigenvalues supplied, he can specify an index range DOL:DOU. */
152 /* Or else the setting DOL=1, DOU=M should be applied. */
153 /* Note that DOL and DOU refer to the order in which the eigenvalues */
154 /* are stored in W. */
155 /* If the user wants to compute only selected eigenpairs, then */
156 /* the columns DOL-1 to DOU+1 of the eigenvector space Z contain the */
157 /* computed eigenvectors. All other columns of Z are set to zero. */
159 /* MINRGP (input) DOUBLE PRECISION */
161 /* RTOL1 (input) DOUBLE PRECISION */
162 /* RTOL2 (input) DOUBLE PRECISION */
163 /* Parameters for bisection. */
164 /* An interval [LEFT,RIGHT] has converged if */
165 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
167 /* W (input/output) DOUBLE PRECISION array, dimension (N) */
168 /* The first M elements of W contain the APPROXIMATE eigenvalues for */
169 /* which eigenvectors are to be computed. The eigenvalues */
170 /* should be grouped by split-off block and ordered from */
171 /* smallest to largest within the block ( The output array */
172 /* W from DLARRE is expected here ). Furthermore, they are with */
173 /* respect to the shift of the corresponding root representation */
174 /* for their block. On exit, W holds the eigenvalues of the */
175 /* UNshifted matrix. */
177 /* WERR (input/output) DOUBLE PRECISION array, dimension (N) */
178 /* The first M elements contain the semiwidth of the uncertainty */
179 /* interval of the corresponding eigenvalue in W */
181 /* WGAP (input/output) DOUBLE PRECISION array, dimension (N) */
182 /* The separation from the right neighbor eigenvalue in W. */
184 /* IBLOCK (input) INTEGER array, dimension (N) */
185 /* The indices of the blocks (submatrices) associated with the */
186 /* corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue */
187 /* W(i) belongs to the first block from the top, =2 if W(i) */
188 /* belongs to the second block, etc. */
190 /* INDEXW (input) INTEGER array, dimension (N) */
191 /* The indices of the eigenvalues within each block (submatrix); */
192 /* for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the */
193 /* i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. */
195 /* GERS (input) DOUBLE PRECISION array, dimension (2*N) */
196 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */
197 /* is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should */
198 /* be computed from the original UNshifted matrix. */
200 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
201 /* If INFO = 0, the first M columns of Z contain the */
202 /* orthonormal eigenvectors of the matrix T */
203 /* corresponding to the input eigenvalues, with the i-th */
204 /* column of Z holding the eigenvector associated with W(i). */
205 /* Note: the user must ensure that at least max(1,M) columns are */
206 /* supplied in the array Z. */
208 /* LDZ (input) INTEGER */
209 /* The leading dimension of the array Z. LDZ >= 1, and if */
210 /* JOBZ = 'V', LDZ >= max(1,N). */
212 /* ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) */
213 /* The support of the eigenvectors in Z, i.e., the indices */
214 /* indicating the nonzero elements in Z. The I-th eigenvector */
215 /* is nonzero only in elements ISUPPZ( 2*I-1 ) through */
218 /* WORK (workspace) DOUBLE PRECISION array, dimension (12*N) */
220 /* IWORK (workspace) INTEGER array, dimension (7*N) */
222 /* INFO (output) INTEGER */
223 /* = 0: successful exit */
225 /* > 0: A problem occured in DLARRV. */
226 /* < 0: One of the called subroutines signaled an internal problem. */
227 /* Needs inspection of the corresponding parameter IINFO */
228 /* for further information. */
230 /* =-1: Problem in DLARRB when refining a child's eigenvalues. */
231 /* =-2: Problem in DLARRF when computing the RRR of a child. */
232 /* When a child is inside a tight cluster, it can be difficult */
233 /* to find an RRR. A partial remedy from the user's point of */
234 /* view is to make the parameter MINRGP smaller and recompile. */
235 /* However, as the orthogonality of the computed vectors is */
236 /* proportional to 1/MINRGP, the user should be aware that */
237 /* he might be trading in precision when he decreases MINRGP. */
238 /* =-3: Problem in DLARRB when refining a single eigenvalue */
239 /* after the Rayleigh correction was rejected. */
240 /* = 5: The Rayleigh Quotient Iteration failed to converge to */
241 /* full accuracy in MAXITR steps. */
243 /* Further Details */
244 /* =============== */
246 /* Based on contributions by */
247 /* Beresford Parlett, University of California, Berkeley, USA */
248 /* Jim Demmel, University of California, Berkeley, USA */
249 /* Inderjit Dhillon, University of Texas, Austin, USA */
250 /* Osni Marques, LBNL/NERSC, USA */
251 /* Christof Voemel, University of California, Berkeley, USA */
253 /* ===================================================================== */
255 /* .. Parameters .. */
257 /* .. Local Scalars .. */
259 /* .. External Functions .. */
261 /* .. External Subroutines .. */
263 /* .. Intrinsic Functions .. */
265 /* .. Executable Statements .. */
267 /* The first N entries of WORK are reserved for the eigenvalues */
268 /* Parameter adjustments */
279 z_offset = 1 + z_dim1;
287 indlld = (*n << 1) + 1;
291 for (i__ = 1; i__ <= i__1; ++i__) {
295 /* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the */
296 /* factorization used to compute the FP vector */
298 /* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current */
299 /* layer and the one above. */
305 for (i__ = 1; i__ <= i__1; ++i__) {
311 /* Set lower bound for use of Z */
316 /* Set lower bound for use of Z */
319 /* The width of the part of Z that is used */
320 zusedw = zusedu - zusedl + 1;
321 dlaset_("Full", n, &zusedw, &c_b5, &c_b5, &z__[zusedl * z_dim1 + 1], ldz);
322 eps = dlamch_("Precision");
325 /* Set expert flags for standard code. */
327 if (*dol == 1 && *dou == *m) {
329 /* Only selected eigenpairs are computed. Since the other evalues */
330 /* are not refined by RQ iteration, bisection has to compute to full */
335 /* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the */
336 /* desired eigenvalues. The support of the nonzero eigenvector */
337 /* entries is contained in the interval IBEGIN:IEND. */
338 /* Remark that if k eigenpairs are desired, then the eigenvectors */
339 /* are stored in k contiguous columns of Z. */
340 /* DONE is the number of eigenvectors already computed */
345 for (jblk = 1; jblk <= i__1; ++jblk) {
348 /* Find the eigenvectors of the submatrix indexed IBEGIN */
353 if (iblock[wend + 1] == jblk) {
361 } else if (wend < *dol || wbegin > *dou) {
366 /* Find local spectral diameter of the block */
367 gl = gers[(ibegin << 1) - 1];
368 gu = gers[ibegin * 2];
370 for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
372 d__1 = gers[(i__ << 1) - 1];
375 d__1 = gers[i__ * 2];
380 /* OLDIEN is the last index of the previous block */
382 /* Calculate the size of the current block */
383 in = iend - ibegin + 1;
384 /* The number of eigenvalues in the current block */
385 im = wend - wbegin + 1;
386 /* This is for a 1x1 block */
387 if (ibegin == iend) {
389 z__[ibegin + wbegin * z_dim1] = 1.;
390 isuppz[(wbegin << 1) - 1] = ibegin;
391 isuppz[wbegin * 2] = ibegin;
393 work[wbegin] = w[wbegin];
398 /* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) */
399 /* Note that these can be approximations, in this case, the corresp. */
400 /* entries of WERR give the size of the uncertainty interval. */
401 /* The eigenvalue approximations will be refined when necessary as */
402 /* high relative accuracy is required for the computation of the */
403 /* corresponding eigenvectors. */
404 dcopy_(&im, &w[wbegin], &c__1, &work[wbegin], &c__1);
405 /* We store in W the eigenvalue approximations w.r.t. the original */
408 for (i__ = 1; i__ <= i__2; ++i__) {
409 w[wbegin + i__ - 1] += sigma;
412 /* NDEPTH is the current depth of the representation tree */
414 /* PARITY is either 1 or 0 */
416 /* NCLUS is the number of clusters for the next level of the */
417 /* representation tree, we start with NCLUS = 1 for the root */
419 iwork[iindc1 + 1] = 1;
420 iwork[iindc1 + 2] = im;
421 /* IDONE is the number of eigenvectors already computed in the current */
424 /* loop while( IDONE.LT.IM ) */
425 /* generate the representation tree for the current block and */
426 /* compute the eigenvectors */
429 /* This is a crude protection against infinitely deep trees */
434 /* breadth first processing of the current level of the representation */
435 /* tree: OLDNCL = number of clusters on current level */
437 /* reset NCLUS to count the number of child clusters */
448 /* Process the clusters on the current level */
450 for (i__ = 1; i__ <= i__2; ++i__) {
451 j = oldcls + (i__ << 1);
452 /* OLDFST, OLDLST = first, last index of current cluster. */
453 /* cluster indices start with 1 and are relative */
454 /* to WBEGIN when accessing W, WGAP, WERR, Z */
455 oldfst = iwork[j - 1];
458 /* Retrieve relatively robust representation (RRR) of cluster */
459 /* that has been computed at the previous level */
460 /* The RRR is stored in Z and overwritten once the eigenvectors */
461 /* have been computed or when the cluster is refined */
462 if (*dol == 1 && *dou == *m) {
463 /* Get representation from location of the leftmost evalue */
465 j = wbegin + oldfst - 1;
467 if (wbegin + oldfst - 1 < *dol) {
468 /* Get representation from the left end of Z array */
470 } else if (wbegin + oldfst - 1 > *dou) {
471 /* Get representation from the right end of Z array */
474 j = wbegin + oldfst - 1;
477 dcopy_(&in, &z__[ibegin + j * z_dim1], &c__1, &d__[ibegin]
480 dcopy_(&i__3, &z__[ibegin + (j + 1) * z_dim1], &c__1, &l[
482 sigma = z__[iend + (j + 1) * z_dim1];
483 /* Set the corresponding entries in Z to zero */
484 dlaset_("Full", &in, &c__2, &c_b5, &c_b5, &z__[ibegin + j
487 /* Compute DL and DLL of current RRR */
489 for (j = ibegin; j <= i__3; ++j) {
491 work[indld - 1 + j] = tmp;
492 work[indlld - 1 + j] = tmp * l[j];
496 /* P and Q are index of the first and last eigenvalue to compute */
497 /* within the current block */
498 p = indexw[wbegin - 1 + oldfst];
499 q = indexw[wbegin - 1 + oldlst];
500 /* Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET */
501 /* thru' Q-OFFSET elements of these arrays are to be used. */
502 /* OFFSET = P-OLDFST */
503 offset = indexw[wbegin] - 1;
504 /* perform limited bisection (if necessary) to get approximate */
505 /* eigenvalues to the precision needed. */
506 dlarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &p,
507 &q, rtol1, rtol2, &offset, &work[wbegin], &wgap[
508 wbegin], &werr[wbegin], &work[indwrk], &iwork[
509 iindwk], pivmin, &spdiam, &in, &iinfo);
514 /* We also recompute the extremal gaps. W holds all eigenvalues */
515 /* of the unshifted matrix and must be used for computation */
516 /* of WGAP, the entries of WORK might stem from RRRs with */
517 /* different shifts. The gaps from WBEGIN-1+OLDFST to */
518 /* WBEGIN-1+OLDLST are correctly computed in DLARRB. */
519 /* However, we only allow the gaps to become greater since */
520 /* this is what should happen when we decrease WERR */
523 d__1 = wgap[wbegin + oldfst - 2], d__2 = w[wbegin +
524 oldfst - 1] - werr[wbegin + oldfst - 1] - w[
525 wbegin + oldfst - 2] - werr[wbegin + oldfst -
527 wgap[wbegin + oldfst - 2] = max(d__1,d__2);
529 if (wbegin + oldlst - 1 < wend) {
531 d__1 = wgap[wbegin + oldlst - 1], d__2 = w[wbegin +
532 oldlst] - werr[wbegin + oldlst] - w[wbegin +
533 oldlst - 1] - werr[wbegin + oldlst - 1];
534 wgap[wbegin + oldlst - 1] = max(d__1,d__2);
536 /* Each time the eigenvalues in WORK get refined, we store */
537 /* the newly found approximation with all shifts applied in W */
539 for (j = oldfst; j <= i__3; ++j) {
540 w[wbegin + j - 1] = work[wbegin + j - 1] + sigma;
544 /* Process the current node. */
547 for (j = oldfst; j <= i__3; ++j) {
549 /* we are at the right end of the cluster, this is also the */
550 /* boundary of the child cluster */
552 } else if (wgap[wbegin + j - 1] >= *minrgp * (d__1 = work[
553 wbegin + j - 1], abs(d__1))) {
554 /* the right relative gap is big enough, the child cluster */
555 /* (NEWFST,..,NEWLST) is well separated from the following */
558 /* inside a child cluster, the relative gap is not */
562 /* Compute size of child cluster found */
563 newsiz = newlst - newfst + 1;
564 /* NEWFTT is the place in Z where the new RRR or the computed */
565 /* eigenvector is to be stored */
566 if (*dol == 1 && *dou == *m) {
567 /* Store representation at location of the leftmost evalue */
569 newftt = wbegin + newfst - 1;
571 if (wbegin + newfst - 1 < *dol) {
572 /* Store representation at the left end of Z array */
574 } else if (wbegin + newfst - 1 > *dou) {
575 /* Store representation at the right end of Z array */
578 newftt = wbegin + newfst - 1;
583 /* Current child is not a singleton but a cluster. */
584 /* Compute and store new representation of child. */
587 /* Compute left and right cluster gap. */
589 /* LGAP and RGAP are not computed from WORK because */
590 /* the eigenvalue approximations may stem from RRRs */
591 /* different shifts. However, W hold all eigenvalues */
592 /* of the unshifted matrix. Still, the entries in WGAP */
593 /* have to be computed from WORK since the entries */
594 /* in W might be of the same order so that gaps are not */
595 /* exhibited correctly for very close eigenvalues. */
598 d__1 = 0., d__2 = w[wbegin] - werr[wbegin] - *vl;
599 lgap = max(d__1,d__2);
601 lgap = wgap[wbegin + newfst - 2];
603 rgap = wgap[wbegin + newlst - 1];
605 /* Compute left- and rightmost eigenvalue of child */
606 /* to high precision in order to shift as close */
607 /* as possible and obtain as large relative gaps */
610 for (k = 1; k <= 2; ++k) {
612 p = indexw[wbegin - 1 + newfst];
614 p = indexw[wbegin - 1 + newlst];
616 offset = indexw[wbegin] - 1;
617 dlarrb_(&in, &d__[ibegin], &work[indlld + ibegin
618 - 1], &p, &p, &rqtol, &rqtol, &offset, &
619 work[wbegin], &wgap[wbegin], &werr[wbegin]
620 , &work[indwrk], &iwork[iindwk], pivmin, &
621 spdiam, &in, &iinfo);
625 if (wbegin + newlst - 1 < *dol || wbegin + newfst - 1
627 /* if the cluster contains no desired eigenvalues */
628 /* skip the computation of that branch of the rep. tree */
630 /* We could skip before the refinement of the extremal */
631 /* eigenvalues of the child, but then the representation */
632 /* tree could be different from the one when nothing is */
633 /* skipped. For this reason we skip at this place. */
634 idone = idone + newlst - newfst + 1;
638 /* Compute RRR of child cluster. */
639 /* Note that the new RRR is stored in Z */
641 /* DLARRF needs LWORK = 2*N */
642 dlarrf_(&in, &d__[ibegin], &l[ibegin], &work[indld +
643 ibegin - 1], &newfst, &newlst, &work[wbegin],
644 &wgap[wbegin], &werr[wbegin], &spdiam, &lgap,
645 &rgap, pivmin, &tau, &z__[ibegin + newftt *
646 z_dim1], &z__[ibegin + (newftt + 1) * z_dim1],
647 &work[indwrk], &iinfo);
649 /* a new RRR for the cluster was found by DLARRF */
650 /* update shift and store it */
651 ssigma = sigma + tau;
652 z__[iend + (newftt + 1) * z_dim1] = ssigma;
653 /* WORK() are the midpoints and WERR() the semi-width */
654 /* Note that the entries in W are unchanged. */
656 for (k = newfst; k <= i__4; ++k) {
657 fudge = eps * 3. * (d__1 = work[wbegin + k -
659 work[wbegin + k - 1] -= tau;
660 fudge += eps * 4. * (d__1 = work[wbegin + k -
663 werr[wbegin + k - 1] += fudge;
664 /* Gaps are not fudged. Provided that WERR is small */
665 /* when eigenvalues are close, a zero gap indicates */
666 /* that a new representation is needed for resolving */
667 /* the cluster. A fudge could lead to a wrong decision */
668 /* of judging eigenvalues 'separated' which in */
669 /* reality are not. This could have a negative impact */
670 /* on the orthogonality of the computed eigenvectors. */
674 k = newcls + (nclus << 1);
675 iwork[k - 1] = newfst;
683 /* Compute eigenvector of singleton */
687 tol = log((doublereal) in) * 4. * eps;
690 windex = wbegin + k - 1;
693 windmn = max(i__4,1);
696 windpl = min(i__4,*m);
697 lambda = work[windex];
699 /* Check if eigenvector computation is to be skipped */
700 if (windex < *dol || windex > *dou) {
706 left = work[windex] - werr[windex];
707 right = work[windex] + werr[windex];
708 indeig = indexw[windex];
709 /* Note that since we compute the eigenpairs for a child, */
710 /* all eigenvalue approximations are w.r.t the same shift. */
711 /* In this case, the entries in WORK should be used for */
712 /* computing the gaps since they exhibit even very small */
713 /* differences in the eigenvalues, as opposed to the */
714 /* entries in W which might "look" the same. */
716 /* In the case RANGE='I' and with not much initial */
717 /* accuracy in LAMBDA and VL, the formula */
718 /* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) */
719 /* can lead to an overestimation of the left gap and */
720 /* thus to inadequately early RQI 'convergence'. */
721 /* Prevent this by forcing a small left gap. */
723 d__1 = abs(left), d__2 = abs(right);
724 lgap = eps * max(d__1,d__2);
729 /* In the case RANGE='I' and with not much initial */
730 /* accuracy in LAMBDA and VU, the formula */
731 /* can lead to an overestimation of the right gap and */
732 /* thus to inadequately early RQI 'convergence'. */
733 /* Prevent this by forcing a small right gap. */
735 d__1 = abs(left), d__2 = abs(right);
736 rgap = eps * max(d__1,d__2);
740 gap = min(lgap,rgap);
741 if (k == 1 || k == im) {
742 /* The eigenvector support can become wrong */
743 /* because significant entries could be cut off due to a */
744 /* large GAPTOL parameter in LAR1V. Prevent this. */
751 /* Update WGAP so that it holds the minimum gap */
752 /* to the left or the right. This is crucial in the */
753 /* case where bisection is used to ensure that the */
754 /* eigenvalue is refined up to the required precision. */
755 /* The correct value is restored afterwards. */
756 savgap = wgap[windex];
758 /* We want to use the Rayleigh Quotient Correction */
759 /* as often as possible since it converges quadratically */
760 /* when we are close enough to the desired eigenvalue. */
761 /* However, the Rayleigh Quotient can have the wrong sign */
762 /* and lead us away from the desired eigenvalue. In this */
763 /* case, the best we can do is to use bisection. */
766 /* Bisection is initially turned off unless it is forced */
769 /* Check if bisection should be used to refine eigenvalue */
771 /* Take the bisection as new iterate */
773 itmp1 = iwork[iindr + windex];
774 offset = indexw[wbegin] - 1;
776 dlarrb_(&in, &d__[ibegin], &work[indlld + ibegin
777 - 1], &indeig, &indeig, &c_b5, &d__1, &
778 offset, &work[wbegin], &wgap[wbegin], &
779 werr[wbegin], &work[indwrk], &iwork[
780 iindwk], pivmin, &spdiam, &itmp1, &iinfo);
785 lambda = work[windex];
786 /* Reset twist index from inaccurate LAMBDA to */
787 /* force computation of true MINGMA */
788 iwork[iindr + windex] = 0;
790 /* Given LAMBDA, compute the eigenvector. */
792 dlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &l[
793 ibegin], &work[indld + ibegin - 1], &work[
794 indlld + ibegin - 1], pivmin, &gaptol, &z__[
795 ibegin + windex * z_dim1], &L__1, &negcnt, &
796 ztz, &mingma, &iwork[iindr + windex], &isuppz[
797 (windex << 1) - 1], &nrminv, &resid, &rqcorr,
802 } else if (resid < bstres) {
807 i__4 = isupmn, i__5 = isuppz[(windex << 1) - 1];
808 isupmn = min(i__4,i__5);
810 i__4 = isupmx, i__5 = isuppz[windex * 2];
811 isupmx = max(i__4,i__5);
813 /* sin alpha <= |resid|/gap */
814 /* Note that both the residual and the gap are */
815 /* proportional to the matrix, so ||T|| doesn't play */
816 /* a role in the quotient */
818 /* Convergence test for Rayleigh-Quotient iteration */
819 /* (omitted when Bisection has been used) */
821 if (resid > tol * gap && abs(rqcorr) > rqtol * abs(
822 lambda) && ! usedbs) {
823 /* We need to check that the RQCORR update doesn't */
824 /* move the eigenvalue away from the desired one and */
825 /* towards a neighbor. -> protection with bisection */
826 if (indeig <= negcnt) {
827 /* The wanted eigenvalue lies to the left */
830 /* The wanted eigenvalue lies to the right */
833 /* We only use the RQCORR if it improves the */
834 /* the iterate reasonably. */
835 if (rqcorr * sgndef >= 0. && lambda + rqcorr <=
836 right && lambda + rqcorr >= left) {
838 /* Store new midpoint of bisection interval in WORK */
840 /* The current LAMBDA is on the left of the true */
843 /* We prefer to assume that the error estimate */
844 /* is correct. We could make the interval not */
845 /* as a bracket but to be modified if the RQCORR */
846 /* chooses to. In this case, the RIGHT side should */
847 /* be modified as follows: */
848 /* RIGHT = MAX(RIGHT, LAMBDA + RQCORR) */
850 /* The current LAMBDA is on the right of the true */
853 /* See comment about assuming the error estimate is */
855 /* LEFT = MIN(LEFT, LAMBDA + RQCORR) */
857 work[windex] = (right + left) * .5;
858 /* Take RQCORR since it has the correct sign and */
859 /* improves the iterate reasonably */
861 /* Update width of error interval */
862 werr[windex] = (right - left) * .5;
866 if (right - left < rqtol * abs(lambda)) {
867 /* The eigenvalue is computed to bisection accuracy */
868 /* compute eigenvector and stop */
871 } else if (iter < 10) {
873 } else if (iter == 10) {
882 if (usedrq && usedbs && bstres <= resid) {
887 /* improve error angle by second step */
889 dlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin]
890 , &l[ibegin], &work[indld + ibegin -
891 1], &work[indlld + ibegin - 1],
892 pivmin, &gaptol, &z__[ibegin + windex
893 * z_dim1], &L__1, &negcnt, &ztz, &
894 mingma, &iwork[iindr + windex], &
895 isuppz[(windex << 1) - 1], &nrminv, &
896 resid, &rqcorr, &work[indwrk]);
898 work[windex] = lambda;
901 /* Compute FP-vector support w.r.t. whole matrix */
903 isuppz[(windex << 1) - 1] += oldien;
904 isuppz[windex * 2] += oldien;
905 zfrom = isuppz[(windex << 1) - 1];
906 zto = isuppz[windex * 2];
909 /* Ensure vector is ok if support in the RQI has changed */
910 if (isupmn < zfrom) {
912 for (ii = isupmn; ii <= i__4; ++ii) {
913 z__[ii + windex * z_dim1] = 0.;
919 for (ii = zto + 1; ii <= i__4; ++ii) {
920 z__[ii + windex * z_dim1] = 0.;
924 i__4 = zto - zfrom + 1;
925 dscal_(&i__4, &nrminv, &z__[zfrom + windex * z_dim1],
929 w[windex] = lambda + sigma;
930 /* Recompute the gaps on the left and right */
931 /* But only allow them to become larger and not */
932 /* smaller (which can only happen through "bad" */
933 /* cancellation and doesn't reflect the theory */
934 /* where the initial gaps are underestimated due */
935 /* to WERR being too crude.) */
939 d__1 = wgap[windmn], d__2 = w[windex] - werr[
940 windex] - w[windmn] - werr[windmn];
941 wgap[windmn] = max(d__1,d__2);
945 d__1 = savgap, d__2 = w[windpl] - werr[windpl]
946 - w[windex] - werr[windex];
947 wgap[windex] = max(d__1,d__2);
952 /* here ends the code for the current child */
955 /* Proceed to any remaining child nodes */