3 /* Subroutine */ int dlarrb_(integer *n, doublereal *d__, doublereal *lld,
4 integer *ifirst, integer *ilast, doublereal *rtol1, doublereal *rtol2,
5 integer *offset, doublereal *w, doublereal *wgap, doublereal *werr,
6 doublereal *work, integer *iwork, doublereal *pivmin, doublereal *
7 spdiam, integer *twist, integer *info)
9 /* System generated locals */
11 doublereal d__1, d__2;
13 /* Builtin functions */
14 double log(doublereal);
17 integer i__, k, r__, i1, ii, ip;
18 doublereal gap, mid, tmp, back, lgap, rgap, left;
19 integer iter, nint, prev, next;
20 doublereal cvrgd, right, width;
21 extern integer dlaneg_(integer *, doublereal *, doublereal *, doublereal *
22 , doublereal *, integer *);
25 integer olnint, maxitr;
28 /* -- LAPACK auxiliary routine (version 3.1) -- */
29 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
32 /* .. Scalar Arguments .. */
34 /* .. Array Arguments .. */
40 /* Given the relatively robust representation(RRR) L D L^T, DLARRB */
41 /* does "limited" bisection to refine the eigenvalues of L D L^T, */
42 /* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial */
43 /* guesses for these eigenvalues are input in W, the corresponding estimate */
44 /* of the error in these guesses and their gaps are input in WERR */
45 /* and WGAP, respectively. During bisection, intervals */
46 /* [left, right] are maintained by storing their mid-points and */
47 /* semi-widths in the arrays W and WERR respectively. */
52 /* N (input) INTEGER */
53 /* The order of the matrix. */
55 /* D (input) DOUBLE PRECISION array, dimension (N) */
56 /* The N diagonal elements of the diagonal matrix D. */
58 /* LLD (input) DOUBLE PRECISION array, dimension (N-1) */
59 /* The (N-1) elements L(i)*L(i)*D(i). */
61 /* IFIRST (input) INTEGER */
62 /* The index of the first eigenvalue to be computed. */
64 /* ILAST (input) INTEGER */
65 /* The index of the last eigenvalue to be computed. */
67 /* RTOL1 (input) DOUBLE PRECISION */
68 /* RTOL2 (input) DOUBLE PRECISION */
69 /* Tolerance for the convergence of the bisection intervals. */
70 /* An interval [LEFT,RIGHT] has converged if */
71 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
72 /* where GAP is the (estimated) distance to the nearest */
75 /* OFFSET (input) INTEGER */
76 /* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET */
77 /* through ILAST-OFFSET elements of these arrays are to be used. */
79 /* W (input/output) DOUBLE PRECISION array, dimension (N) */
80 /* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are */
81 /* estimates of the eigenvalues of L D L^T indexed IFIRST throug */
83 /* On output, these estimates are refined. */
85 /* WGAP (input/output) DOUBLE PRECISION array, dimension (N-1) */
86 /* On input, the (estimated) gaps between consecutive */
87 /* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between */
88 /* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST */
89 /* then WGAP(IFIRST-OFFSET) must be set to ZERO. */
90 /* On output, these gaps are refined. */
92 /* WERR (input/output) DOUBLE PRECISION array, dimension (N) */
93 /* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are */
94 /* the errors in the estimates of the corresponding elements in W. */
95 /* On output, these errors are refined. */
97 /* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */
100 /* IWORK (workspace) INTEGER array, dimension (2*N) */
103 /* PIVMIN (input) DOUBLE PRECISION */
104 /* The minimum pivot in the Sturm sequence. */
106 /* SPDIAM (input) DOUBLE PRECISION */
107 /* The spectral diameter of the matrix. */
109 /* TWIST (input) INTEGER */
110 /* The twist index for the twisted factorization that is used */
111 /* for the negcount. */
112 /* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T */
113 /* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T */
114 /* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) */
116 /* INFO (output) INTEGER */
119 /* Further Details */
120 /* =============== */
122 /* Based on contributions by */
123 /* Beresford Parlett, University of California, Berkeley, USA */
124 /* Jim Demmel, University of California, Berkeley, USA */
125 /* Inderjit Dhillon, University of Texas, Austin, USA */
126 /* Osni Marques, LBNL/NERSC, USA */
127 /* Christof Voemel, University of California, Berkeley, USA */
129 /* ===================================================================== */
131 /* .. Parameters .. */
133 /* .. Local Scalars .. */
135 /* .. External Functions .. */
138 /* .. Intrinsic Functions .. */
140 /* .. Executable Statements .. */
142 /* Parameter adjustments */
154 maxitr = (integer) ((log(*spdiam + *pivmin) - log(*pivmin)) / log(2.)) +
156 mnwdth = *pivmin * 2.;
159 if (r__ < 1 || r__ > *n) {
163 /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. */
164 /* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while */
165 /* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) */
166 /* for an unconverged interval is set to the index of the next unconverged */
167 /* interval, and is -1 or 0 for a converged interval. Thus a linked */
168 /* list of unconverged intervals is set up. */
171 /* The number of unconverged intervals */
173 /* The last unconverged interval found */
175 rgap = wgap[i1 - *offset];
177 for (i__ = i1; i__ <= i__1; ++i__) {
180 left = w[ii] - werr[ii];
181 right = w[ii] + werr[ii];
184 gap = min(lgap,rgap);
185 /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue */
186 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT */
188 /* Do while( NEGCNT(LEFT).GT.I-1 ) */
192 negcnt = dlaneg_(n, &d__[1], &lld[1], &left, pivmin, &r__);
193 if (negcnt > i__ - 1) {
199 /* Do while( NEGCNT(RIGHT).LT.I ) */
200 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT */
204 negcnt = dlaneg_(n, &d__[1], &lld[1], &right, pivmin, &r__);
210 width = (d__1 = left - right, abs(d__1)) * .5;
212 d__1 = abs(left), d__2 = abs(right);
213 tmp = max(d__1,d__2);
215 d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
216 cvrgd = max(d__1,d__2);
217 if (width <= cvrgd || width <= mnwdth) {
218 /* This interval has already converged and does not need refinement. */
219 /* (Note that the gaps might change through refining the */
220 /* eigenvalues, however, they can only get bigger.) */
221 /* Remove it from the list. */
223 /* Make sure that I1 always points to the first unconverged interval */
224 if (i__ == i1 && i__ < *ilast) {
227 if (prev >= i1 && i__ <= *ilast) {
228 iwork[(prev << 1) - 1] = i__ + 1;
231 /* unconverged interval found */
234 iwork[k - 1] = i__ + 1;
242 /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals */
243 /* and while (ITER.LT.MAXITR) */
251 for (ip = 1; ip <= i__1; ++ip) {
259 gap = min(lgap,rgap);
263 mid = (left + right) * .5;
264 /* semiwidth of interval */
267 d__1 = abs(left), d__2 = abs(right);
268 tmp = max(d__1,d__2);
270 d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
271 cvrgd = max(d__1,d__2);
272 if (width <= cvrgd || width <= mnwdth || iter == maxitr) {
273 /* reduce number of unconverged intervals */
275 /* Mark interval as converged. */
280 /* Prev holds the last unconverged interval previously examined */
282 iwork[(prev << 1) - 1] = next;
290 /* Perform one bisection step */
292 negcnt = dlaneg_(n, &d__[1], &lld[1], &mid, pivmin, &r__);
293 if (negcnt <= i__ - 1) {
303 /* do another loop if there are still unconverged intervals */
304 /* However, in the last iteration, all intervals are accepted */
305 /* since this is the best we can do. */
306 if (nint > 0 && iter <= maxitr) {
311 /* At this point, all the intervals have converged */
313 for (i__ = *ifirst; i__ <= i__1; ++i__) {
316 /* All intervals marked by '0' have been refined. */
317 if (iwork[k - 1] == 0) {
318 w[ii] = (work[k - 1] + work[k]) * .5;
319 werr[ii] = work[k] - w[ii];
325 for (i__ = *ifirst + 1; i__ <= i__1; ++i__) {
329 d__1 = 0., d__2 = w[ii] - werr[ii] - w[ii - 1] - werr[ii - 1];
330 wgap[ii - 1] = max(d__1,d__2);