3 /* Subroutine */ int slaebz_(integer *ijob, integer *nitmax, integer *n,
4 integer *mmax, integer *minp, integer *nbmin, real *abstol, real *
5 reltol, real *pivmin, real *d__, real *e, real *e2, integer *nval,
6 real *ab, real *c__, integer *mout, integer *nab, real *work, integer
9 /* System generated locals */
10 integer nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4,
12 real r__1, r__2, r__3, r__4;
15 integer j, kf, ji, kl, jp, jit;
17 integer itmp1, itmp2, kfnew, klnew;
20 /* -- LAPACK auxiliary routine (version 3.1) -- */
21 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
24 /* .. Scalar Arguments .. */
26 /* .. Array Arguments .. */
32 /* SLAEBZ contains the iteration loops which compute and use the */
33 /* function N(w), which is the count of eigenvalues of a symmetric */
34 /* tridiagonal matrix T less than or equal to its argument w. It */
35 /* performs a choice of two types of loops: */
37 /* IJOB=1, followed by */
38 /* IJOB=2: It takes as input a list of intervals and returns a list of */
39 /* sufficiently small intervals whose union contains the same */
40 /* eigenvalues as the union of the original intervals. */
41 /* The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. */
42 /* The output interval (AB(j,1),AB(j,2)] will contain */
43 /* eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. */
45 /* IJOB=3: It performs a binary search in each input interval */
46 /* (AB(j,1),AB(j,2)] for a point w(j) such that */
47 /* N(w(j))=NVAL(j), and uses C(j) as the starting point of */
48 /* the search. If such a w(j) is found, then on output */
49 /* AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output */
50 /* (AB(j,1),AB(j,2)] will be a small interval containing the */
51 /* point where N(w) jumps through NVAL(j), unless that point */
52 /* lies outside the initial interval. */
54 /* Note that the intervals are in all cases half-open intervals, */
55 /* i.e., of the form (a,b] , which includes b but not a . */
57 /* To avoid underflow, the matrix should be scaled so that its largest */
58 /* element is no greater than overflow**(1/2) * underflow**(1/4) */
59 /* in absolute value. To assure the most accurate computation */
60 /* of small eigenvalues, the matrix should be scaled to be */
61 /* not much smaller than that, either. */
63 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
64 /* Matrix", Report CS41, Computer Science Dept., Stanford */
65 /* University, July 21, 1966 */
67 /* Note: the arguments are, in general, *not* checked for unreasonable */
73 /* IJOB (input) INTEGER */
74 /* Specifies what is to be done: */
75 /* = 1: Compute NAB for the initial intervals. */
76 /* = 2: Perform bisection iteration to find eigenvalues of T. */
77 /* = 3: Perform bisection iteration to invert N(w), i.e., */
78 /* to find a point which has a specified number of */
79 /* eigenvalues of T to its left. */
80 /* Other values will cause SLAEBZ to return with INFO=-1. */
82 /* NITMAX (input) INTEGER */
83 /* The maximum number of "levels" of bisection to be */
84 /* performed, i.e., an interval of width W will not be made */
85 /* smaller than 2^(-NITMAX) * W. If not all intervals */
86 /* have converged after NITMAX iterations, then INFO is set */
87 /* to the number of non-converged intervals. */
89 /* N (input) INTEGER */
90 /* The dimension n of the tridiagonal matrix T. It must be at */
93 /* MMAX (input) INTEGER */
94 /* The maximum number of intervals. If more than MMAX intervals */
95 /* are generated, then SLAEBZ will quit with INFO=MMAX+1. */
97 /* MINP (input) INTEGER */
98 /* The initial number of intervals. It may not be greater than */
101 /* NBMIN (input) INTEGER */
102 /* The smallest number of intervals that should be processed */
103 /* using a vector loop. If zero, then only the scalar loop */
106 /* ABSTOL (input) REAL */
107 /* The minimum (absolute) width of an interval. When an */
108 /* interval is narrower than ABSTOL, or than RELTOL times the */
109 /* larger (in magnitude) endpoint, then it is considered to be */
110 /* sufficiently small, i.e., converged. This must be at least */
113 /* RELTOL (input) REAL */
114 /* The minimum relative width of an interval. When an interval */
115 /* is narrower than ABSTOL, or than RELTOL times the larger (in */
116 /* magnitude) endpoint, then it is considered to be */
117 /* sufficiently small, i.e., converged. Note: this should */
118 /* always be at least radix*machine epsilon. */
120 /* PIVMIN (input) REAL */
121 /* The minimum absolute value of a "pivot" in the Sturm */
122 /* sequence loop. This *must* be at least max |e(j)**2| * */
123 /* safe_min and at least safe_min, where safe_min is at least */
124 /* the smallest number that can divide one without overflow. */
126 /* D (input) REAL array, dimension (N) */
127 /* The diagonal elements of the tridiagonal matrix T. */
129 /* E (input) REAL array, dimension (N) */
130 /* The offdiagonal elements of the tridiagonal matrix T in */
131 /* positions 1 through N-1. E(N) is arbitrary. */
133 /* E2 (input) REAL array, dimension (N) */
134 /* The squares of the offdiagonal elements of the tridiagonal */
135 /* matrix T. E2(N) is ignored. */
137 /* NVAL (input/output) INTEGER array, dimension (MINP) */
138 /* If IJOB=1 or 2, not referenced. */
139 /* If IJOB=3, the desired values of N(w). The elements of NVAL */
140 /* will be reordered to correspond with the intervals in AB. */
141 /* Thus, NVAL(j) on output will not, in general be the same as */
142 /* NVAL(j) on input, but it will correspond with the interval */
143 /* (AB(j,1),AB(j,2)] on output. */
145 /* AB (input/output) REAL array, dimension (MMAX,2) */
146 /* The endpoints of the intervals. AB(j,1) is a(j), the left */
147 /* endpoint of the j-th interval, and AB(j,2) is b(j), the */
148 /* right endpoint of the j-th interval. The input intervals */
149 /* will, in general, be modified, split, and reordered by the */
152 /* C (input/output) REAL array, dimension (MMAX) */
153 /* If IJOB=1, ignored. */
154 /* If IJOB=2, workspace. */
155 /* If IJOB=3, then on input C(j) should be initialized to the */
156 /* first search point in the binary search. */
158 /* MOUT (output) INTEGER */
159 /* If IJOB=1, the number of eigenvalues in the intervals. */
160 /* If IJOB=2 or 3, the number of intervals output. */
161 /* If IJOB=3, MOUT will equal MINP. */
163 /* NAB (input/output) INTEGER array, dimension (MMAX,2) */
164 /* If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). */
165 /* If IJOB=2, then on input, NAB(i,j) should be set. It must */
166 /* satisfy the condition: */
167 /* N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), */
168 /* which means that in interval i only eigenvalues */
169 /* NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, */
170 /* NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with */
172 /* On output, NAB(i,j) will contain */
173 /* max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of */
174 /* the input interval that the output interval */
175 /* (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the */
176 /* the input values of NAB(k,1) and NAB(k,2). */
177 /* If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), */
178 /* unless N(w) > NVAL(i) for all search points w , in which */
179 /* case NAB(i,1) will not be modified, i.e., the output */
180 /* value will be the same as the input value (modulo */
181 /* reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) */
182 /* for all search points w , in which case NAB(i,2) will */
183 /* not be modified. Normally, NAB should be set to some */
184 /* distinctive value(s) before SLAEBZ is called. */
186 /* WORK (workspace) REAL array, dimension (MMAX) */
189 /* IWORK (workspace) INTEGER array, dimension (MMAX) */
192 /* INFO (output) INTEGER */
193 /* = 0: All intervals converged. */
194 /* = 1--MMAX: The last INFO intervals did not converge. */
195 /* = MMAX+1: More than MMAX intervals were generated. */
197 /* Further Details */
198 /* =============== */
200 /* This routine is intended to be called only by other LAPACK */
201 /* routines, thus the interface is less user-friendly. It is intended */
202 /* for two purposes: */
204 /* (a) finding eigenvalues. In this case, SLAEBZ should have one or */
205 /* more initial intervals set up in AB, and SLAEBZ should be called */
206 /* with IJOB=1. This sets up NAB, and also counts the eigenvalues. */
207 /* Intervals with no eigenvalues would usually be thrown out at */
208 /* this point. Also, if not all the eigenvalues in an interval i */
209 /* are desired, NAB(i,1) can be increased or NAB(i,2) decreased. */
210 /* For example, set NAB(i,1)=NAB(i,2)-1 to get the largest */
211 /* eigenvalue. SLAEBZ is then called with IJOB=2 and MMAX */
212 /* no smaller than the value of MOUT returned by the call with */
213 /* IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 */
214 /* through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the */
215 /* tolerance specified by ABSTOL and RELTOL. */
217 /* (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). */
218 /* In this case, start with a Gershgorin interval (a,b). Set up */
219 /* AB to contain 2 search intervals, both initially (a,b). One */
220 /* NVAL element should contain f-1 and the other should contain l */
221 /* , while C should contain a and b, resp. NAB(i,1) should be -1 */
222 /* and NAB(i,2) should be N+1, to flag an error if the desired */
223 /* interval does not lie in (a,b). SLAEBZ is then called with */
224 /* IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- */
225 /* j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while */
226 /* if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r */
227 /* >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and */
228 /* N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and */
229 /* w(l-r)=...=w(l+k) are handled similarly. */
231 /* ===================================================================== */
233 /* .. Parameters .. */
235 /* .. Local Scalars .. */
237 /* .. Intrinsic Functions .. */
239 /* .. Executable Statements .. */
241 /* Check for Errors */
243 /* Parameter adjustments */
245 nab_offset = 1 + nab_dim1;
248 ab_offset = 1 + ab_dim1;
260 if (*ijob < 1 || *ijob > 3) {
269 /* Compute the number of eigenvalues in the initial intervals. */
274 for (ji = 1; ji <= i__1; ++ji) {
275 for (jp = 1; jp <= 2; ++jp) {
276 tmp1 = d__[1] - ab[ji + jp * ab_dim1];
277 if (dabs(tmp1) < *pivmin) {
280 nab[ji + jp * nab_dim1] = 0;
282 nab[ji + jp * nab_dim1] = 1;
286 for (j = 2; j <= i__2; ++j) {
287 tmp1 = d__[j] - e2[j - 1] / tmp1 - ab[ji + jp * ab_dim1];
288 if (dabs(tmp1) < *pivmin) {
292 ++nab[ji + jp * nab_dim1];
298 *mout = *mout + nab[ji + (nab_dim1 << 1)] - nab[ji + nab_dim1];
304 /* Initialize for loop */
306 /* KF and KL have the following meaning: */
307 /* Intervals 1,...,KF-1 have converged. */
308 /* Intervals KF,...,KL still need to be refined. */
313 /* If IJOB=2, initialize C. */
314 /* If IJOB=3, use the user-supplied starting point. */
318 for (ji = 1; ji <= i__1; ++ji) {
319 c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5f;
327 for (jit = 1; jit <= i__1; ++jit) {
329 /* Loop over intervals */
331 if (kl - kf + 1 >= *nbmin && *nbmin > 0) {
333 /* Begin of Parallel Version of the loop */
336 for (ji = kf; ji <= i__2; ++ji) {
338 /* Compute N(c), the number of eigenvalues less than c */
340 work[ji] = d__[1] - c__[ji];
342 if (work[ji] <= *pivmin) {
345 r__1 = work[ji], r__2 = -(*pivmin);
346 work[ji] = dmin(r__1,r__2);
350 for (j = 2; j <= i__3; ++j) {
351 work[ji] = d__[j] - e2[j - 1] / work[ji] - c__[ji];
352 if (work[ji] <= *pivmin) {
355 r__1 = work[ji], r__2 = -(*pivmin);
356 work[ji] = dmin(r__1,r__2);
365 /* IJOB=2: Choose all intervals containing eigenvalues. */
369 for (ji = kf; ji <= i__2; ++ji) {
371 /* Insure that N(w) is monotone */
375 i__5 = nab[ji + nab_dim1], i__6 = iwork[ji];
376 i__3 = nab[ji + (nab_dim1 << 1)], i__4 = max(i__5,i__6);
377 iwork[ji] = min(i__3,i__4);
379 /* Update the Queue -- add intervals if both halves */
380 /* contain eigenvalues. */
382 if (iwork[ji] == nab[ji + (nab_dim1 << 1)]) {
384 /* No eigenvalue in the upper interval: */
385 /* just use the lower interval. */
387 ab[ji + (ab_dim1 << 1)] = c__[ji];
389 } else if (iwork[ji] == nab[ji + nab_dim1]) {
391 /* No eigenvalue in the lower interval: */
392 /* just use the upper interval. */
394 ab[ji + ab_dim1] = c__[ji];
397 if (klnew <= *mmax) {
399 /* Eigenvalue in both intervals -- add upper to */
402 ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 <<
404 nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1
406 ab[klnew + ab_dim1] = c__[ji];
407 nab[klnew + nab_dim1] = iwork[ji];
408 ab[ji + (ab_dim1 << 1)] = c__[ji];
409 nab[ji + (nab_dim1 << 1)] = iwork[ji];
422 /* IJOB=3: Binary search. Keep only the interval containing */
423 /* w s.t. N(w) = NVAL */
426 for (ji = kf; ji <= i__2; ++ji) {
427 if (iwork[ji] <= nval[ji]) {
428 ab[ji + ab_dim1] = c__[ji];
429 nab[ji + nab_dim1] = iwork[ji];
431 if (iwork[ji] >= nval[ji]) {
432 ab[ji + (ab_dim1 << 1)] = c__[ji];
433 nab[ji + (nab_dim1 << 1)] = iwork[ji];
441 /* End of Parallel Version of the loop */
443 /* Begin of Serial Version of the loop */
447 for (ji = kf; ji <= i__2; ++ji) {
449 /* Compute N(w), the number of eigenvalues less than w */
452 tmp2 = d__[1] - tmp1;
454 if (tmp2 <= *pivmin) {
457 r__1 = tmp2, r__2 = -(*pivmin);
458 tmp2 = dmin(r__1,r__2);
461 /* A series of compiler directives to defeat vectorization */
462 /* for the next loop */
464 /* $PL$ CMCHAR=' ' */
465 /* DIR$ NEXTSCALAR */
467 /* DIR$ NEXT SCALAR */
472 /* VOCL LOOP,SCALAR */
473 /* IBM PREFER SCALAR */
474 /* $PL$ CMCHAR='*' */
477 for (j = 2; j <= i__3; ++j) {
478 tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1;
479 if (tmp2 <= *pivmin) {
482 r__1 = tmp2, r__2 = -(*pivmin);
483 tmp2 = dmin(r__1,r__2);
490 /* IJOB=2: Choose all intervals containing eigenvalues. */
492 /* Insure that N(w) is monotone */
496 i__5 = nab[ji + nab_dim1];
497 i__3 = nab[ji + (nab_dim1 << 1)], i__4 = max(i__5,itmp1);
498 itmp1 = min(i__3,i__4);
500 /* Update the Queue -- add intervals if both halves */
501 /* contain eigenvalues. */
503 if (itmp1 == nab[ji + (nab_dim1 << 1)]) {
505 /* No eigenvalue in the upper interval: */
506 /* just use the lower interval. */
508 ab[ji + (ab_dim1 << 1)] = tmp1;
510 } else if (itmp1 == nab[ji + nab_dim1]) {
512 /* No eigenvalue in the lower interval: */
513 /* just use the upper interval. */
515 ab[ji + ab_dim1] = tmp1;
516 } else if (klnew < *mmax) {
518 /* Eigenvalue in both intervals -- add upper to queue. */
521 ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << 1)];
522 nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 <<
524 ab[klnew + ab_dim1] = tmp1;
525 nab[klnew + nab_dim1] = itmp1;
526 ab[ji + (ab_dim1 << 1)] = tmp1;
527 nab[ji + (nab_dim1 << 1)] = itmp1;
534 /* IJOB=3: Binary search. Keep only the interval */
535 /* containing w s.t. N(w) = NVAL */
537 if (itmp1 <= nval[ji]) {
538 ab[ji + ab_dim1] = tmp1;
539 nab[ji + nab_dim1] = itmp1;
541 if (itmp1 >= nval[ji]) {
542 ab[ji + (ab_dim1 << 1)] = tmp1;
543 nab[ji + (nab_dim1 << 1)] = itmp1;
550 /* End of Serial Version of the loop */
554 /* Check for convergence */
558 for (ji = kf; ji <= i__2; ++ji) {
559 tmp1 = (r__1 = ab[ji + (ab_dim1 << 1)] - ab[ji + ab_dim1], dabs(
562 r__3 = (r__1 = ab[ji + (ab_dim1 << 1)], dabs(r__1)), r__4 = (r__2
563 = ab[ji + ab_dim1], dabs(r__2));
564 tmp2 = dmax(r__3,r__4);
566 r__1 = max(*abstol,*pivmin), r__2 = *reltol * tmp2;
567 if (tmp1 < dmax(r__1,r__2) || nab[ji + nab_dim1] >= nab[ji + (
570 /* Converged -- Swap with position KFNEW, */
571 /* then increment KFNEW */
574 tmp1 = ab[ji + ab_dim1];
575 tmp2 = ab[ji + (ab_dim1 << 1)];
576 itmp1 = nab[ji + nab_dim1];
577 itmp2 = nab[ji + (nab_dim1 << 1)];
578 ab[ji + ab_dim1] = ab[kfnew + ab_dim1];
579 ab[ji + (ab_dim1 << 1)] = ab[kfnew + (ab_dim1 << 1)];
580 nab[ji + nab_dim1] = nab[kfnew + nab_dim1];
581 nab[ji + (nab_dim1 << 1)] = nab[kfnew + (nab_dim1 << 1)];
582 ab[kfnew + ab_dim1] = tmp1;
583 ab[kfnew + (ab_dim1 << 1)] = tmp2;
584 nab[kfnew + nab_dim1] = itmp1;
585 nab[kfnew + (nab_dim1 << 1)] = itmp2;
588 nval[ji] = nval[kfnew];
598 /* Choose Midpoints */
601 for (ji = kf; ji <= i__2; ++ji) {
602 c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5f;
606 /* If no more intervals to refine, quit. */