3 /* Subroutine */ int slarrb_(integer *n, real *d__, real *lld, integer *
4 ifirst, integer *ilast, real *rtol1, real *rtol2, integer *offset,
5 real *w, real *wgap, real *werr, real *work, integer *iwork, real *
6 pivmin, real *spdiam, integer *twist, integer *info)
8 /* System generated locals */
12 /* Builtin functions */
13 double log(doublereal);
16 integer i__, k, r__, i1, ii, ip;
17 real gap, mid, tmp, back, lgap, rgap, left;
18 integer iter, nint, prev, next;
19 real cvrgd, right, width;
20 extern integer slaneg_(integer *, real *, real *, real *, real *, integer
24 integer olnint, maxitr;
27 /* -- LAPACK auxiliary routine (version 3.1) -- */
28 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
31 /* .. Scalar Arguments .. */
33 /* .. Array Arguments .. */
39 /* Given the relatively robust representation(RRR) L D L^T, SLARRB */
40 /* does "limited" bisection to refine the eigenvalues of L D L^T, */
41 /* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial */
42 /* guesses for these eigenvalues are input in W, the corresponding estimate */
43 /* of the error in these guesses and their gaps are input in WERR */
44 /* and WGAP, respectively. During bisection, intervals */
45 /* [left, right] are maintained by storing their mid-points and */
46 /* semi-widths in the arrays W and WERR respectively. */
51 /* N (input) INTEGER */
52 /* The order of the matrix. */
54 /* D (input) REAL array, dimension (N) */
55 /* The N diagonal elements of the diagonal matrix D. */
57 /* LLD (input) REAL array, dimension (N-1) */
58 /* The (N-1) elements L(i)*L(i)*D(i). */
60 /* IFIRST (input) INTEGER */
61 /* The index of the first eigenvalue to be computed. */
63 /* ILAST (input) INTEGER */
64 /* The index of the last eigenvalue to be computed. */
66 /* RTOL1 (input) REAL */
67 /* RTOL2 (input) REAL */
68 /* Tolerance for the convergence of the bisection intervals. */
69 /* An interval [LEFT,RIGHT] has converged if */
70 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
71 /* where GAP is the (estimated) distance to the nearest */
74 /* OFFSET (input) INTEGER */
75 /* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET */
76 /* through ILAST-OFFSET elements of these arrays are to be used. */
78 /* W (input/output) REAL array, dimension (N) */
79 /* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are */
80 /* estimates of the eigenvalues of L D L^T indexed IFIRST throug */
82 /* On output, these estimates are refined. */
84 /* WGAP (input/output) REAL array, dimension (N-1) */
85 /* On input, the (estimated) gaps between consecutive */
86 /* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between */
87 /* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST */
88 /* then WGAP(IFIRST-OFFSET) must be set to ZERO. */
89 /* On output, these gaps are refined. */
91 /* WERR (input/output) REAL array, dimension (N) */
92 /* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are */
93 /* the errors in the estimates of the corresponding elements in W. */
94 /* On output, these errors are refined. */
96 /* WORK (workspace) REAL array, dimension (2*N) */
99 /* IWORK (workspace) INTEGER array, dimension (2*N) */
102 /* PIVMIN (input) DOUBLE PRECISION */
103 /* The minimum pivot in the Sturm sequence. */
105 /* SPDIAM (input) DOUBLE PRECISION */
106 /* The spectral diameter of the matrix. */
108 /* TWIST (input) INTEGER */
109 /* The twist index for the twisted factorization that is used */
110 /* for the negcount. */
111 /* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T */
112 /* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T */
113 /* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) */
115 /* INFO (output) INTEGER */
118 /* Further Details */
119 /* =============== */
121 /* Based on contributions by */
122 /* Beresford Parlett, University of California, Berkeley, USA */
123 /* Jim Demmel, University of California, Berkeley, USA */
124 /* Inderjit Dhillon, University of Texas, Austin, USA */
125 /* Osni Marques, LBNL/NERSC, USA */
126 /* Christof Voemel, University of California, Berkeley, USA */
128 /* ===================================================================== */
130 /* .. Parameters .. */
132 /* .. Local Scalars .. */
134 /* .. External Functions .. */
137 /* .. Intrinsic Functions .. */
139 /* .. Executable Statements .. */
141 /* Parameter adjustments */
153 maxitr = (integer) ((log(*spdiam + *pivmin) - log(*pivmin)) / log(2.f)) +
155 mnwdth = *pivmin * 2.f;
158 if (r__ < 1 || r__ > *n) {
162 /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. */
163 /* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while */
164 /* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) */
165 /* for an unconverged interval is set to the index of the next unconverged */
166 /* interval, and is -1 or 0 for a converged interval. Thus a linked */
167 /* list of unconverged intervals is set up. */
170 /* The number of unconverged intervals */
172 /* The last unconverged interval found */
174 rgap = wgap[i1 - *offset];
176 for (i__ = i1; i__ <= i__1; ++i__) {
179 left = w[ii] - werr[ii];
180 right = w[ii] + werr[ii];
183 gap = dmin(lgap,rgap);
184 /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue */
185 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT */
187 /* Do while( NEGCNT(LEFT).GT.I-1 ) */
191 negcnt = slaneg_(n, &d__[1], &lld[1], &left, pivmin, &r__);
192 if (negcnt > i__ - 1) {
198 /* Do while( NEGCNT(RIGHT).LT.I ) */
199 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT */
203 negcnt = slaneg_(n, &d__[1], &lld[1], &right, pivmin, &r__);
209 width = (r__1 = left - right, dabs(r__1)) * .5f;
211 r__1 = dabs(left), r__2 = dabs(right);
212 tmp = dmax(r__1,r__2);
214 r__1 = *rtol1 * gap, r__2 = *rtol2 * tmp;
215 cvrgd = dmax(r__1,r__2);
216 if (width <= cvrgd || width <= mnwdth) {
217 /* This interval has already converged and does not need refinement. */
218 /* (Note that the gaps might change through refining the */
219 /* eigenvalues, however, they can only get bigger.) */
220 /* Remove it from the list. */
222 /* Make sure that I1 always points to the first unconverged interval */
223 if (i__ == i1 && i__ < *ilast) {
226 if (prev >= i1 && i__ <= *ilast) {
227 iwork[(prev << 1) - 1] = i__ + 1;
230 /* unconverged interval found */
233 iwork[k - 1] = i__ + 1;
241 /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals */
242 /* and while (ITER.LT.MAXITR) */
250 for (ip = 1; ip <= i__1; ++ip) {
258 gap = dmin(lgap,rgap);
262 mid = (left + right) * .5f;
263 /* semiwidth of interval */
266 r__1 = dabs(left), r__2 = dabs(right);
267 tmp = dmax(r__1,r__2);
269 r__1 = *rtol1 * gap, r__2 = *rtol2 * tmp;
270 cvrgd = dmax(r__1,r__2);
271 if (width <= cvrgd || width <= mnwdth || iter == maxitr) {
272 /* reduce number of unconverged intervals */
274 /* Mark interval as converged. */
279 /* Prev holds the last unconverged interval previously examined */
281 iwork[(prev << 1) - 1] = next;
289 /* Perform one bisection step */
291 negcnt = slaneg_(n, &d__[1], &lld[1], &mid, pivmin, &r__);
292 if (negcnt <= i__ - 1) {
302 /* do another loop if there are still unconverged intervals */
303 /* However, in the last iteration, all intervals are accepted */
304 /* since this is the best we can do. */
305 if (nint > 0 && iter <= maxitr) {
310 /* At this point, all the intervals have converged */
312 for (i__ = *ifirst; i__ <= i__1; ++i__) {
315 /* All intervals marked by '0' have been refined. */
316 if (iwork[k - 1] == 0) {
317 w[ii] = (work[k - 1] + work[k]) * .5f;
318 werr[ii] = work[k] - w[ii];
324 for (i__ = *ifirst + 1; i__ <= i__1; ++i__) {
328 r__1 = 0.f, r__2 = w[ii] - werr[ii] - w[ii - 1] - werr[ii - 1];
329 wgap[ii - 1] = dmax(r__1,r__2);