1 SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
2 $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
3 $ NAB, WORK, IWORK, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
12 DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
16 DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
17 $ WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DLAEBZ contains the iteration loops which compute and use the
24 * function N(w), which is the count of eigenvalues of a symmetric
25 * tridiagonal matrix T less than or equal to its argument w. It
26 * performs a choice of two types of loops:
27 *
28 * IJOB=1, followed by
29 * IJOB=2: It takes as input a list of intervals and returns a list of
30 * sufficiently small intervals whose union contains the same
31 * eigenvalues as the union of the original intervals.
32 * The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
33 * The output interval (AB(j,1),AB(j,2)] will contain
34 * eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
35 *
36 * IJOB=3: It performs a binary search in each input interval
37 * (AB(j,1),AB(j,2)] for a point w(j) such that
38 * N(w(j))=NVAL(j), and uses C(j) as the starting point of
39 * the search. If such a w(j) is found, then on output
40 * AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output
41 * (AB(j,1),AB(j,2)] will be a small interval containing the
42 * point where N(w) jumps through NVAL(j), unless that point
43 * lies outside the initial interval.
44 *
45 * Note that the intervals are in all cases half-open intervals,
46 * i.e., of the form (a,b] , which includes b but not a .
47 *
48 * To avoid underflow, the matrix should be scaled so that its largest
49 * element is no greater than overflow**(1/2) * underflow**(1/4)
50 * in absolute value. To assure the most accurate computation
51 * of small eigenvalues, the matrix should be scaled to be
52 * not much smaller than that, either.
53 *
54 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
55 * Matrix", Report CS41, Computer Science Dept., Stanford
56 * University, July 21, 1966
57 *
58 * Note: the arguments are, in general, *not* checked for unreasonable
59 * values.
60 *
61 * Arguments
62 * =========
63 *
64 * IJOB (input) INTEGER
65 * Specifies what is to be done:
66 * = 1: Compute NAB for the initial intervals.
67 * = 2: Perform bisection iteration to find eigenvalues of T.
68 * = 3: Perform bisection iteration to invert N(w), i.e.,
69 * to find a point which has a specified number of
70 * eigenvalues of T to its left.
71 * Other values will cause DLAEBZ to return with INFO=-1.
72 *
73 * NITMAX (input) INTEGER
74 * The maximum number of "levels" of bisection to be
75 * performed, i.e., an interval of width W will not be made
76 * smaller than 2^(-NITMAX) * W. If not all intervals
77 * have converged after NITMAX iterations, then INFO is set
78 * to the number of non-converged intervals.
79 *
80 * N (input) INTEGER
81 * The dimension n of the tridiagonal matrix T. It must be at
82 * least 1.
83 *
84 * MMAX (input) INTEGER
85 * The maximum number of intervals. If more than MMAX intervals
86 * are generated, then DLAEBZ will quit with INFO=MMAX+1.
87 *
88 * MINP (input) INTEGER
89 * The initial number of intervals. It may not be greater than
90 * MMAX.
91 *
92 * NBMIN (input) INTEGER
93 * The smallest number of intervals that should be processed
94 * using a vector loop. If zero, then only the scalar loop
95 * will be used.
96 *
97 * ABSTOL (input) DOUBLE PRECISION
98 * The minimum (absolute) width of an interval. When an
99 * interval is narrower than ABSTOL, or than RELTOL times the
100 * larger (in magnitude) endpoint, then it is considered to be
101 * sufficiently small, i.e., converged. This must be at least
102 * zero.
103 *
104 * RELTOL (input) DOUBLE PRECISION
105 * The minimum relative width of an interval. When an interval
106 * is narrower than ABSTOL, or than RELTOL times the larger (in
107 * magnitude) endpoint, then it is considered to be
108 * sufficiently small, i.e., converged. Note: this should
109 * always be at least radix*machine epsilon.
110 *
111 * PIVMIN (input) DOUBLE PRECISION
112 * The minimum absolute value of a "pivot" in the Sturm
113 * sequence loop. This *must* be at least max |e(j)**2| *
114 * safe_min and at least safe_min, where safe_min is at least
115 * the smallest number that can divide one without overflow.
116 *
117 * D (input) DOUBLE PRECISION array, dimension (N)
118 * The diagonal elements of the tridiagonal matrix T.
119 *
120 * E (input) DOUBLE PRECISION array, dimension (N)
121 * The offdiagonal elements of the tridiagonal matrix T in
122 * positions 1 through N-1. E(N) is arbitrary.
123 *
124 * E2 (input) DOUBLE PRECISION array, dimension (N)
125 * The squares of the offdiagonal elements of the tridiagonal
126 * matrix T. E2(N) is ignored.
127 *
128 * NVAL (input/output) INTEGER array, dimension (MINP)
129 * If IJOB=1 or 2, not referenced.
130 * If IJOB=3, the desired values of N(w). The elements of NVAL
131 * will be reordered to correspond with the intervals in AB.
132 * Thus, NVAL(j) on output will not, in general be the same as
133 * NVAL(j) on input, but it will correspond with the interval
134 * (AB(j,1),AB(j,2)] on output.
135 *
136 * AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
137 * The endpoints of the intervals. AB(j,1) is a(j), the left
138 * endpoint of the j-th interval, and AB(j,2) is b(j), the
139 * right endpoint of the j-th interval. The input intervals
140 * will, in general, be modified, split, and reordered by the
141 * calculation.
142 *
143 * C (input/output) DOUBLE PRECISION array, dimension (MMAX)
144 * If IJOB=1, ignored.
145 * If IJOB=2, workspace.
146 * If IJOB=3, then on input C(j) should be initialized to the
147 * first search point in the binary search.
148 *
149 * MOUT (output) INTEGER
150 * If IJOB=1, the number of eigenvalues in the intervals.
151 * If IJOB=2 or 3, the number of intervals output.
152 * If IJOB=3, MOUT will equal MINP.
153 *
154 * NAB (input/output) INTEGER array, dimension (MMAX,2)
155 * If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
156 * If IJOB=2, then on input, NAB(i,j) should be set. It must
157 * satisfy the condition:
158 * N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
159 * which means that in interval i only eigenvalues
160 * NAB(i,1)+1,...,NAB(i,2) will be considered. Usually,
161 * NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
162 * IJOB=1.
163 * On output, NAB(i,j) will contain
164 * max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
165 * the input interval that the output interval
166 * (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
167 * the input values of NAB(k,1) and NAB(k,2).
168 * If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
169 * unless N(w) > NVAL(i) for all search points w , in which
170 * case NAB(i,1) will not be modified, i.e., the output
171 * value will be the same as the input value (modulo
172 * reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
173 * for all search points w , in which case NAB(i,2) will
174 * not be modified. Normally, NAB should be set to some
175 * distinctive value(s) before DLAEBZ is called.
176 *
177 * WORK (workspace) DOUBLE PRECISION array, dimension (MMAX)
178 * Workspace.
179 *
180 * IWORK (workspace) INTEGER array, dimension (MMAX)
181 * Workspace.
182 *
183 * INFO (output) INTEGER
184 * = 0: All intervals converged.
185 * = 1--MMAX: The last INFO intervals did not converge.
186 * = MMAX+1: More than MMAX intervals were generated.
187 *
188 * Further Details
189 * ===============
190 *
191 * This routine is intended to be called only by other LAPACK
192 * routines, thus the interface is less user-friendly. It is intended
193 * for two purposes:
194 *
195 * (a) finding eigenvalues. In this case, DLAEBZ should have one or
196 * more initial intervals set up in AB, and DLAEBZ should be called
197 * with IJOB=1. This sets up NAB, and also counts the eigenvalues.
198 * Intervals with no eigenvalues would usually be thrown out at
199 * this point. Also, if not all the eigenvalues in an interval i
200 * are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
201 * For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
202 * eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX
203 * no smaller than the value of MOUT returned by the call with
204 * IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
205 * through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
206 * tolerance specified by ABSTOL and RELTOL.
207 *
208 * (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
209 * In this case, start with a Gershgorin interval (a,b). Set up
210 * AB to contain 2 search intervals, both initially (a,b). One
211 * NVAL element should contain f-1 and the other should contain l
212 * , while C should contain a and b, resp. NAB(i,1) should be -1
213 * and NAB(i,2) should be N+1, to flag an error if the desired
214 * interval does not lie in (a,b). DLAEBZ is then called with
215 * IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
216 * j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
217 * if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
218 * >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
219 * N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
220 * w(l-r)=...=w(l+k) are handled similarly.
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225 DOUBLE PRECISION ZERO, TWO, HALF
226 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
227 $ HALF = 1.0D0 / TWO )
228 * ..
229 * .. Local Scalars ..
230 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
231 $ KLNEW
232 DOUBLE PRECISION TMP1, TMP2
233 * ..
234 * .. Intrinsic Functions ..
235 INTRINSIC ABS, MAX, MIN
236 * ..
237 * .. Executable Statements ..
238 *
239 * Check for Errors
240 *
241 INFO = 0
242 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
243 INFO = -1
244 RETURN
245 END IF
246 *
247 * Initialize NAB
248 *
249 IF( IJOB.EQ.1 ) THEN
250 *
251 * Compute the number of eigenvalues in the initial intervals.
252 *
253 MOUT = 0
254 DO 30 JI = 1, MINP
255 DO 20 JP = 1, 2
256 TMP1 = D( 1 ) - AB( JI, JP )
257 IF( ABS( TMP1 ).LT.PIVMIN )
258 $ TMP1 = -PIVMIN
259 NAB( JI, JP ) = 0
260 IF( TMP1.LE.ZERO )
261 $ NAB( JI, JP ) = 1
262 *
263 DO 10 J = 2, N
264 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
265 IF( ABS( TMP1 ).LT.PIVMIN )
266 $ TMP1 = -PIVMIN
267 IF( TMP1.LE.ZERO )
268 $ NAB( JI, JP ) = NAB( JI, JP ) + 1
269 10 CONTINUE
270 20 CONTINUE
271 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
272 30 CONTINUE
273 RETURN
274 END IF
275 *
276 * Initialize for loop
277 *
278 * KF and KL have the following meaning:
279 * Intervals 1,...,KF-1 have converged.
280 * Intervals KF,...,KL still need to be refined.
281 *
282 KF = 1
283 KL = MINP
284 *
285 * If IJOB=2, initialize C.
286 * If IJOB=3, use the user-supplied starting point.
287 *
288 IF( IJOB.EQ.2 ) THEN
289 DO 40 JI = 1, MINP
290 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
291 40 CONTINUE
292 END IF
293 *
294 * Iteration loop
295 *
296 DO 130 JIT = 1, NITMAX
297 *
298 * Loop over intervals
299 *
300 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
301 *
302 * Begin of Parallel Version of the loop
303 *
304 DO 60 JI = KF, KL
305 *
306 * Compute N(c), the number of eigenvalues less than c
307 *
308 WORK( JI ) = D( 1 ) - C( JI )
309 IWORK( JI ) = 0
310 IF( WORK( JI ).LE.PIVMIN ) THEN
311 IWORK( JI ) = 1
312 WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
313 END IF
314 *
315 DO 50 J = 2, N
316 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
317 IF( WORK( JI ).LE.PIVMIN ) THEN
318 IWORK( JI ) = IWORK( JI ) + 1
319 WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
320 END IF
321 50 CONTINUE
322 60 CONTINUE
323 *
324 IF( IJOB.LE.2 ) THEN
325 *
326 * IJOB=2: Choose all intervals containing eigenvalues.
327 *
328 KLNEW = KL
329 DO 70 JI = KF, KL
330 *
331 * Insure that N(w) is monotone
332 *
333 IWORK( JI ) = MIN( NAB( JI, 2 ),
334 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) )
335 *
336 * Update the Queue -- add intervals if both halves
337 * contain eigenvalues.
338 *
339 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
340 *
341 * No eigenvalue in the upper interval:
342 * just use the lower interval.
343 *
344 AB( JI, 2 ) = C( JI )
345 *
346 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
347 *
348 * No eigenvalue in the lower interval:
349 * just use the upper interval.
350 *
351 AB( JI, 1 ) = C( JI )
352 ELSE
353 KLNEW = KLNEW + 1
354 IF( KLNEW.LE.MMAX ) THEN
355 *
356 * Eigenvalue in both intervals -- add upper to
357 * queue.
358 *
359 AB( KLNEW, 2 ) = AB( JI, 2 )
360 NAB( KLNEW, 2 ) = NAB( JI, 2 )
361 AB( KLNEW, 1 ) = C( JI )
362 NAB( KLNEW, 1 ) = IWORK( JI )
363 AB( JI, 2 ) = C( JI )
364 NAB( JI, 2 ) = IWORK( JI )
365 ELSE
366 INFO = MMAX + 1
367 END IF
368 END IF
369 70 CONTINUE
370 IF( INFO.NE.0 )
371 $ RETURN
372 KL = KLNEW
373 ELSE
374 *
375 * IJOB=3: Binary search. Keep only the interval containing
376 * w s.t. N(w) = NVAL
377 *
378 DO 80 JI = KF, KL
379 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
380 AB( JI, 1 ) = C( JI )
381 NAB( JI, 1 ) = IWORK( JI )
382 END IF
383 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
384 AB( JI, 2 ) = C( JI )
385 NAB( JI, 2 ) = IWORK( JI )
386 END IF
387 80 CONTINUE
388 END IF
389 *
390 ELSE
391 *
392 * End of Parallel Version of the loop
393 *
394 * Begin of Serial Version of the loop
395 *
396 KLNEW = KL
397 DO 100 JI = KF, KL
398 *
399 * Compute N(w), the number of eigenvalues less than w
400 *
401 TMP1 = C( JI )
402 TMP2 = D( 1 ) - TMP1
403 ITMP1 = 0
404 IF( TMP2.LE.PIVMIN ) THEN
405 ITMP1 = 1
406 TMP2 = MIN( TMP2, -PIVMIN )
407 END IF
408 *
409 DO 90 J = 2, N
410 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
411 IF( TMP2.LE.PIVMIN ) THEN
412 ITMP1 = ITMP1 + 1
413 TMP2 = MIN( TMP2, -PIVMIN )
414 END IF
415 90 CONTINUE
416 *
417 IF( IJOB.LE.2 ) THEN
418 *
419 * IJOB=2: Choose all intervals containing eigenvalues.
420 *
421 * Insure that N(w) is monotone
422 *
423 ITMP1 = MIN( NAB( JI, 2 ),
424 $ MAX( NAB( JI, 1 ), ITMP1 ) )
425 *
426 * Update the Queue -- add intervals if both halves
427 * contain eigenvalues.
428 *
429 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
430 *
431 * No eigenvalue in the upper interval:
432 * just use the lower interval.
433 *
434 AB( JI, 2 ) = TMP1
435 *
436 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
437 *
438 * No eigenvalue in the lower interval:
439 * just use the upper interval.
440 *
441 AB( JI, 1 ) = TMP1
442 ELSE IF( KLNEW.LT.MMAX ) THEN
443 *
444 * Eigenvalue in both intervals -- add upper to queue.
445 *
446 KLNEW = KLNEW + 1
447 AB( KLNEW, 2 ) = AB( JI, 2 )
448 NAB( KLNEW, 2 ) = NAB( JI, 2 )
449 AB( KLNEW, 1 ) = TMP1
450 NAB( KLNEW, 1 ) = ITMP1
451 AB( JI, 2 ) = TMP1
452 NAB( JI, 2 ) = ITMP1
453 ELSE
454 INFO = MMAX + 1
455 RETURN
456 END IF
457 ELSE
458 *
459 * IJOB=3: Binary search. Keep only the interval
460 * containing w s.t. N(w) = NVAL
461 *
462 IF( ITMP1.LE.NVAL( JI ) ) THEN
463 AB( JI, 1 ) = TMP1
464 NAB( JI, 1 ) = ITMP1
465 END IF
466 IF( ITMP1.GE.NVAL( JI ) ) THEN
467 AB( JI, 2 ) = TMP1
468 NAB( JI, 2 ) = ITMP1
469 END IF
470 END IF
471 100 CONTINUE
472 KL = KLNEW
473 *
474 END IF
475 *
476 * Check for convergence
477 *
478 KFNEW = KF
479 DO 110 JI = KF, KL
480 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
481 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
482 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
483 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
484 *
485 * Converged -- Swap with position KFNEW,
486 * then increment KFNEW
487 *
488 IF( JI.GT.KFNEW ) THEN
489 TMP1 = AB( JI, 1 )
490 TMP2 = AB( JI, 2 )
491 ITMP1 = NAB( JI, 1 )
492 ITMP2 = NAB( JI, 2 )
493 AB( JI, 1 ) = AB( KFNEW, 1 )
494 AB( JI, 2 ) = AB( KFNEW, 2 )
495 NAB( JI, 1 ) = NAB( KFNEW, 1 )
496 NAB( JI, 2 ) = NAB( KFNEW, 2 )
497 AB( KFNEW, 1 ) = TMP1
498 AB( KFNEW, 2 ) = TMP2
499 NAB( KFNEW, 1 ) = ITMP1
500 NAB( KFNEW, 2 ) = ITMP2
501 IF( IJOB.EQ.3 ) THEN
502 ITMP1 = NVAL( JI )
503 NVAL( JI ) = NVAL( KFNEW )
504 NVAL( KFNEW ) = ITMP1
505 END IF
506 END IF
507 KFNEW = KFNEW + 1
508 END IF
509 110 CONTINUE
510 KF = KFNEW
511 *
512 * Choose Midpoints
513 *
514 DO 120 JI = KF, KL
515 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
516 120 CONTINUE
517 *
518 * If no more intervals to refine, quit.
519 *
520 IF( KF.GT.KL )
521 $ GO TO 140
522 130 CONTINUE
523 *
524 * Converged
525 *
526 140 CONTINUE
527 INFO = MAX( KL+1-KF, 0 )
528 MOUT = KL
529 *
530 RETURN
531 *
532 * End of DLAEBZ
533 *
534 END
2 $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
3 $ NAB, WORK, IWORK, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
12 DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
16 DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
17 $ WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DLAEBZ contains the iteration loops which compute and use the
24 * function N(w), which is the count of eigenvalues of a symmetric
25 * tridiagonal matrix T less than or equal to its argument w. It
26 * performs a choice of two types of loops:
27 *
28 * IJOB=1, followed by
29 * IJOB=2: It takes as input a list of intervals and returns a list of
30 * sufficiently small intervals whose union contains the same
31 * eigenvalues as the union of the original intervals.
32 * The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
33 * The output interval (AB(j,1),AB(j,2)] will contain
34 * eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
35 *
36 * IJOB=3: It performs a binary search in each input interval
37 * (AB(j,1),AB(j,2)] for a point w(j) such that
38 * N(w(j))=NVAL(j), and uses C(j) as the starting point of
39 * the search. If such a w(j) is found, then on output
40 * AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output
41 * (AB(j,1),AB(j,2)] will be a small interval containing the
42 * point where N(w) jumps through NVAL(j), unless that point
43 * lies outside the initial interval.
44 *
45 * Note that the intervals are in all cases half-open intervals,
46 * i.e., of the form (a,b] , which includes b but not a .
47 *
48 * To avoid underflow, the matrix should be scaled so that its largest
49 * element is no greater than overflow**(1/2) * underflow**(1/4)
50 * in absolute value. To assure the most accurate computation
51 * of small eigenvalues, the matrix should be scaled to be
52 * not much smaller than that, either.
53 *
54 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
55 * Matrix", Report CS41, Computer Science Dept., Stanford
56 * University, July 21, 1966
57 *
58 * Note: the arguments are, in general, *not* checked for unreasonable
59 * values.
60 *
61 * Arguments
62 * =========
63 *
64 * IJOB (input) INTEGER
65 * Specifies what is to be done:
66 * = 1: Compute NAB for the initial intervals.
67 * = 2: Perform bisection iteration to find eigenvalues of T.
68 * = 3: Perform bisection iteration to invert N(w), i.e.,
69 * to find a point which has a specified number of
70 * eigenvalues of T to its left.
71 * Other values will cause DLAEBZ to return with INFO=-1.
72 *
73 * NITMAX (input) INTEGER
74 * The maximum number of "levels" of bisection to be
75 * performed, i.e., an interval of width W will not be made
76 * smaller than 2^(-NITMAX) * W. If not all intervals
77 * have converged after NITMAX iterations, then INFO is set
78 * to the number of non-converged intervals.
79 *
80 * N (input) INTEGER
81 * The dimension n of the tridiagonal matrix T. It must be at
82 * least 1.
83 *
84 * MMAX (input) INTEGER
85 * The maximum number of intervals. If more than MMAX intervals
86 * are generated, then DLAEBZ will quit with INFO=MMAX+1.
87 *
88 * MINP (input) INTEGER
89 * The initial number of intervals. It may not be greater than
90 * MMAX.
91 *
92 * NBMIN (input) INTEGER
93 * The smallest number of intervals that should be processed
94 * using a vector loop. If zero, then only the scalar loop
95 * will be used.
96 *
97 * ABSTOL (input) DOUBLE PRECISION
98 * The minimum (absolute) width of an interval. When an
99 * interval is narrower than ABSTOL, or than RELTOL times the
100 * larger (in magnitude) endpoint, then it is considered to be
101 * sufficiently small, i.e., converged. This must be at least
102 * zero.
103 *
104 * RELTOL (input) DOUBLE PRECISION
105 * The minimum relative width of an interval. When an interval
106 * is narrower than ABSTOL, or than RELTOL times the larger (in
107 * magnitude) endpoint, then it is considered to be
108 * sufficiently small, i.e., converged. Note: this should
109 * always be at least radix*machine epsilon.
110 *
111 * PIVMIN (input) DOUBLE PRECISION
112 * The minimum absolute value of a "pivot" in the Sturm
113 * sequence loop. This *must* be at least max |e(j)**2| *
114 * safe_min and at least safe_min, where safe_min is at least
115 * the smallest number that can divide one without overflow.
116 *
117 * D (input) DOUBLE PRECISION array, dimension (N)
118 * The diagonal elements of the tridiagonal matrix T.
119 *
120 * E (input) DOUBLE PRECISION array, dimension (N)
121 * The offdiagonal elements of the tridiagonal matrix T in
122 * positions 1 through N-1. E(N) is arbitrary.
123 *
124 * E2 (input) DOUBLE PRECISION array, dimension (N)
125 * The squares of the offdiagonal elements of the tridiagonal
126 * matrix T. E2(N) is ignored.
127 *
128 * NVAL (input/output) INTEGER array, dimension (MINP)
129 * If IJOB=1 or 2, not referenced.
130 * If IJOB=3, the desired values of N(w). The elements of NVAL
131 * will be reordered to correspond with the intervals in AB.
132 * Thus, NVAL(j) on output will not, in general be the same as
133 * NVAL(j) on input, but it will correspond with the interval
134 * (AB(j,1),AB(j,2)] on output.
135 *
136 * AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
137 * The endpoints of the intervals. AB(j,1) is a(j), the left
138 * endpoint of the j-th interval, and AB(j,2) is b(j), the
139 * right endpoint of the j-th interval. The input intervals
140 * will, in general, be modified, split, and reordered by the
141 * calculation.
142 *
143 * C (input/output) DOUBLE PRECISION array, dimension (MMAX)
144 * If IJOB=1, ignored.
145 * If IJOB=2, workspace.
146 * If IJOB=3, then on input C(j) should be initialized to the
147 * first search point in the binary search.
148 *
149 * MOUT (output) INTEGER
150 * If IJOB=1, the number of eigenvalues in the intervals.
151 * If IJOB=2 or 3, the number of intervals output.
152 * If IJOB=3, MOUT will equal MINP.
153 *
154 * NAB (input/output) INTEGER array, dimension (MMAX,2)
155 * If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
156 * If IJOB=2, then on input, NAB(i,j) should be set. It must
157 * satisfy the condition:
158 * N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
159 * which means that in interval i only eigenvalues
160 * NAB(i,1)+1,...,NAB(i,2) will be considered. Usually,
161 * NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
162 * IJOB=1.
163 * On output, NAB(i,j) will contain
164 * max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
165 * the input interval that the output interval
166 * (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
167 * the input values of NAB(k,1) and NAB(k,2).
168 * If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
169 * unless N(w) > NVAL(i) for all search points w , in which
170 * case NAB(i,1) will not be modified, i.e., the output
171 * value will be the same as the input value (modulo
172 * reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
173 * for all search points w , in which case NAB(i,2) will
174 * not be modified. Normally, NAB should be set to some
175 * distinctive value(s) before DLAEBZ is called.
176 *
177 * WORK (workspace) DOUBLE PRECISION array, dimension (MMAX)
178 * Workspace.
179 *
180 * IWORK (workspace) INTEGER array, dimension (MMAX)
181 * Workspace.
182 *
183 * INFO (output) INTEGER
184 * = 0: All intervals converged.
185 * = 1--MMAX: The last INFO intervals did not converge.
186 * = MMAX+1: More than MMAX intervals were generated.
187 *
188 * Further Details
189 * ===============
190 *
191 * This routine is intended to be called only by other LAPACK
192 * routines, thus the interface is less user-friendly. It is intended
193 * for two purposes:
194 *
195 * (a) finding eigenvalues. In this case, DLAEBZ should have one or
196 * more initial intervals set up in AB, and DLAEBZ should be called
197 * with IJOB=1. This sets up NAB, and also counts the eigenvalues.
198 * Intervals with no eigenvalues would usually be thrown out at
199 * this point. Also, if not all the eigenvalues in an interval i
200 * are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
201 * For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
202 * eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX
203 * no smaller than the value of MOUT returned by the call with
204 * IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
205 * through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
206 * tolerance specified by ABSTOL and RELTOL.
207 *
208 * (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
209 * In this case, start with a Gershgorin interval (a,b). Set up
210 * AB to contain 2 search intervals, both initially (a,b). One
211 * NVAL element should contain f-1 and the other should contain l
212 * , while C should contain a and b, resp. NAB(i,1) should be -1
213 * and NAB(i,2) should be N+1, to flag an error if the desired
214 * interval does not lie in (a,b). DLAEBZ is then called with
215 * IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
216 * j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
217 * if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
218 * >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
219 * N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
220 * w(l-r)=...=w(l+k) are handled similarly.
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225 DOUBLE PRECISION ZERO, TWO, HALF
226 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
227 $ HALF = 1.0D0 / TWO )
228 * ..
229 * .. Local Scalars ..
230 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
231 $ KLNEW
232 DOUBLE PRECISION TMP1, TMP2
233 * ..
234 * .. Intrinsic Functions ..
235 INTRINSIC ABS, MAX, MIN
236 * ..
237 * .. Executable Statements ..
238 *
239 * Check for Errors
240 *
241 INFO = 0
242 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
243 INFO = -1
244 RETURN
245 END IF
246 *
247 * Initialize NAB
248 *
249 IF( IJOB.EQ.1 ) THEN
250 *
251 * Compute the number of eigenvalues in the initial intervals.
252 *
253 MOUT = 0
254 DO 30 JI = 1, MINP
255 DO 20 JP = 1, 2
256 TMP1 = D( 1 ) - AB( JI, JP )
257 IF( ABS( TMP1 ).LT.PIVMIN )
258 $ TMP1 = -PIVMIN
259 NAB( JI, JP ) = 0
260 IF( TMP1.LE.ZERO )
261 $ NAB( JI, JP ) = 1
262 *
263 DO 10 J = 2, N
264 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
265 IF( ABS( TMP1 ).LT.PIVMIN )
266 $ TMP1 = -PIVMIN
267 IF( TMP1.LE.ZERO )
268 $ NAB( JI, JP ) = NAB( JI, JP ) + 1
269 10 CONTINUE
270 20 CONTINUE
271 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
272 30 CONTINUE
273 RETURN
274 END IF
275 *
276 * Initialize for loop
277 *
278 * KF and KL have the following meaning:
279 * Intervals 1,...,KF-1 have converged.
280 * Intervals KF,...,KL still need to be refined.
281 *
282 KF = 1
283 KL = MINP
284 *
285 * If IJOB=2, initialize C.
286 * If IJOB=3, use the user-supplied starting point.
287 *
288 IF( IJOB.EQ.2 ) THEN
289 DO 40 JI = 1, MINP
290 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
291 40 CONTINUE
292 END IF
293 *
294 * Iteration loop
295 *
296 DO 130 JIT = 1, NITMAX
297 *
298 * Loop over intervals
299 *
300 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
301 *
302 * Begin of Parallel Version of the loop
303 *
304 DO 60 JI = KF, KL
305 *
306 * Compute N(c), the number of eigenvalues less than c
307 *
308 WORK( JI ) = D( 1 ) - C( JI )
309 IWORK( JI ) = 0
310 IF( WORK( JI ).LE.PIVMIN ) THEN
311 IWORK( JI ) = 1
312 WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
313 END IF
314 *
315 DO 50 J = 2, N
316 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
317 IF( WORK( JI ).LE.PIVMIN ) THEN
318 IWORK( JI ) = IWORK( JI ) + 1
319 WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
320 END IF
321 50 CONTINUE
322 60 CONTINUE
323 *
324 IF( IJOB.LE.2 ) THEN
325 *
326 * IJOB=2: Choose all intervals containing eigenvalues.
327 *
328 KLNEW = KL
329 DO 70 JI = KF, KL
330 *
331 * Insure that N(w) is monotone
332 *
333 IWORK( JI ) = MIN( NAB( JI, 2 ),
334 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) )
335 *
336 * Update the Queue -- add intervals if both halves
337 * contain eigenvalues.
338 *
339 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
340 *
341 * No eigenvalue in the upper interval:
342 * just use the lower interval.
343 *
344 AB( JI, 2 ) = C( JI )
345 *
346 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
347 *
348 * No eigenvalue in the lower interval:
349 * just use the upper interval.
350 *
351 AB( JI, 1 ) = C( JI )
352 ELSE
353 KLNEW = KLNEW + 1
354 IF( KLNEW.LE.MMAX ) THEN
355 *
356 * Eigenvalue in both intervals -- add upper to
357 * queue.
358 *
359 AB( KLNEW, 2 ) = AB( JI, 2 )
360 NAB( KLNEW, 2 ) = NAB( JI, 2 )
361 AB( KLNEW, 1 ) = C( JI )
362 NAB( KLNEW, 1 ) = IWORK( JI )
363 AB( JI, 2 ) = C( JI )
364 NAB( JI, 2 ) = IWORK( JI )
365 ELSE
366 INFO = MMAX + 1
367 END IF
368 END IF
369 70 CONTINUE
370 IF( INFO.NE.0 )
371 $ RETURN
372 KL = KLNEW
373 ELSE
374 *
375 * IJOB=3: Binary search. Keep only the interval containing
376 * w s.t. N(w) = NVAL
377 *
378 DO 80 JI = KF, KL
379 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
380 AB( JI, 1 ) = C( JI )
381 NAB( JI, 1 ) = IWORK( JI )
382 END IF
383 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
384 AB( JI, 2 ) = C( JI )
385 NAB( JI, 2 ) = IWORK( JI )
386 END IF
387 80 CONTINUE
388 END IF
389 *
390 ELSE
391 *
392 * End of Parallel Version of the loop
393 *
394 * Begin of Serial Version of the loop
395 *
396 KLNEW = KL
397 DO 100 JI = KF, KL
398 *
399 * Compute N(w), the number of eigenvalues less than w
400 *
401 TMP1 = C( JI )
402 TMP2 = D( 1 ) - TMP1
403 ITMP1 = 0
404 IF( TMP2.LE.PIVMIN ) THEN
405 ITMP1 = 1
406 TMP2 = MIN( TMP2, -PIVMIN )
407 END IF
408 *
409 DO 90 J = 2, N
410 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
411 IF( TMP2.LE.PIVMIN ) THEN
412 ITMP1 = ITMP1 + 1
413 TMP2 = MIN( TMP2, -PIVMIN )
414 END IF
415 90 CONTINUE
416 *
417 IF( IJOB.LE.2 ) THEN
418 *
419 * IJOB=2: Choose all intervals containing eigenvalues.
420 *
421 * Insure that N(w) is monotone
422 *
423 ITMP1 = MIN( NAB( JI, 2 ),
424 $ MAX( NAB( JI, 1 ), ITMP1 ) )
425 *
426 * Update the Queue -- add intervals if both halves
427 * contain eigenvalues.
428 *
429 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
430 *
431 * No eigenvalue in the upper interval:
432 * just use the lower interval.
433 *
434 AB( JI, 2 ) = TMP1
435 *
436 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
437 *
438 * No eigenvalue in the lower interval:
439 * just use the upper interval.
440 *
441 AB( JI, 1 ) = TMP1
442 ELSE IF( KLNEW.LT.MMAX ) THEN
443 *
444 * Eigenvalue in both intervals -- add upper to queue.
445 *
446 KLNEW = KLNEW + 1
447 AB( KLNEW, 2 ) = AB( JI, 2 )
448 NAB( KLNEW, 2 ) = NAB( JI, 2 )
449 AB( KLNEW, 1 ) = TMP1
450 NAB( KLNEW, 1 ) = ITMP1
451 AB( JI, 2 ) = TMP1
452 NAB( JI, 2 ) = ITMP1
453 ELSE
454 INFO = MMAX + 1
455 RETURN
456 END IF
457 ELSE
458 *
459 * IJOB=3: Binary search. Keep only the interval
460 * containing w s.t. N(w) = NVAL
461 *
462 IF( ITMP1.LE.NVAL( JI ) ) THEN
463 AB( JI, 1 ) = TMP1
464 NAB( JI, 1 ) = ITMP1
465 END IF
466 IF( ITMP1.GE.NVAL( JI ) ) THEN
467 AB( JI, 2 ) = TMP1
468 NAB( JI, 2 ) = ITMP1
469 END IF
470 END IF
471 100 CONTINUE
472 KL = KLNEW
473 *
474 END IF
475 *
476 * Check for convergence
477 *
478 KFNEW = KF
479 DO 110 JI = KF, KL
480 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
481 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
482 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
483 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
484 *
485 * Converged -- Swap with position KFNEW,
486 * then increment KFNEW
487 *
488 IF( JI.GT.KFNEW ) THEN
489 TMP1 = AB( JI, 1 )
490 TMP2 = AB( JI, 2 )
491 ITMP1 = NAB( JI, 1 )
492 ITMP2 = NAB( JI, 2 )
493 AB( JI, 1 ) = AB( KFNEW, 1 )
494 AB( JI, 2 ) = AB( KFNEW, 2 )
495 NAB( JI, 1 ) = NAB( KFNEW, 1 )
496 NAB( JI, 2 ) = NAB( KFNEW, 2 )
497 AB( KFNEW, 1 ) = TMP1
498 AB( KFNEW, 2 ) = TMP2
499 NAB( KFNEW, 1 ) = ITMP1
500 NAB( KFNEW, 2 ) = ITMP2
501 IF( IJOB.EQ.3 ) THEN
502 ITMP1 = NVAL( JI )
503 NVAL( JI ) = NVAL( KFNEW )
504 NVAL( KFNEW ) = ITMP1
505 END IF
506 END IF
507 KFNEW = KFNEW + 1
508 END IF
509 110 CONTINUE
510 KF = KFNEW
511 *
512 * Choose Midpoints
513 *
514 DO 120 JI = KF, KL
515 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
516 120 CONTINUE
517 *
518 * If no more intervals to refine, quit.
519 *
520 IF( KF.GT.KL )
521 $ GO TO 140
522 130 CONTINUE
523 *
524 * Converged
525 *
526 140 CONTINUE
527 INFO = MAX( KL+1-KF, 0 )
528 MOUT = KL
529 *
530 RETURN
531 *
532 * End of DLAEBZ
533 *
534 END