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.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
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: *DIR$ NOVECTOR
255: DO 30 JI = 1, MINP
256: DO 20 JP = 1, 2
257: TMP1 = D( 1 ) - AB( JI, JP )
258: IF( ABS( TMP1 ).LT.PIVMIN )
259: $ TMP1 = -PIVMIN
260: NAB( JI, JP ) = 0
261: IF( TMP1.LE.ZERO )
262: $ NAB( JI, JP ) = 1
263: *
264: DO 10 J = 2, N
265: TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
266: IF( ABS( TMP1 ).LT.PIVMIN )
267: $ TMP1 = -PIVMIN
268: IF( TMP1.LE.ZERO )
269: $ NAB( JI, JP ) = NAB( JI, JP ) + 1
270: 10 CONTINUE
271: 20 CONTINUE
272: MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
273: 30 CONTINUE
274: RETURN
275: END IF
276: *
277: * Initialize for loop
278: *
279: * KF and KL have the following meaning:
280: * Intervals 1,...,KF-1 have converged.
281: * Intervals KF,...,KL still need to be refined.
282: *
283: KF = 1
284: KL = MINP
285: *
286: * If IJOB=2, initialize C.
287: * If IJOB=3, use the user-supplied starting point.
288: *
289: IF( IJOB.EQ.2 ) THEN
290: DO 40 JI = 1, MINP
291: C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
292: 40 CONTINUE
293: END IF
294: *
295: * Iteration loop
296: *
297: DO 130 JIT = 1, NITMAX
298: *
299: * Loop over intervals
300: *
301: IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
302: *
303: * Begin of Parallel Version of the loop
304: *
305: DO 60 JI = KF, KL
306: *
307: * Compute N(c), the number of eigenvalues less than c
308: *
309: WORK( JI ) = D( 1 ) - C( JI )
310: IWORK( JI ) = 0
311: IF( WORK( JI ).LE.PIVMIN ) THEN
312: IWORK( JI ) = 1
313: WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
314: END IF
315: *
316: DO 50 J = 2, N
317: WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
318: IF( WORK( JI ).LE.PIVMIN ) THEN
319: IWORK( JI ) = IWORK( JI ) + 1
320: WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
321: END IF
322: 50 CONTINUE
323: 60 CONTINUE
324: *
325: IF( IJOB.LE.2 ) THEN
326: *
327: * IJOB=2: Choose all intervals containing eigenvalues.
328: *
329: KLNEW = KL
330: DO 70 JI = KF, KL
331: *
332: * Insure that N(w) is monotone
333: *
334: IWORK( JI ) = MIN( NAB( JI, 2 ),
335: $ MAX( NAB( JI, 1 ), IWORK( JI ) ) )
336: *
337: * Update the Queue -- add intervals if both halves
338: * contain eigenvalues.
339: *
340: IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
341: *
342: * No eigenvalue in the upper interval:
343: * just use the lower interval.
344: *
345: AB( JI, 2 ) = C( JI )
346: *
347: ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
348: *
349: * No eigenvalue in the lower interval:
350: * just use the upper interval.
351: *
352: AB( JI, 1 ) = C( JI )
353: ELSE
354: KLNEW = KLNEW + 1
355: IF( KLNEW.LE.MMAX ) THEN
356: *
357: * Eigenvalue in both intervals -- add upper to
358: * queue.
359: *
360: AB( KLNEW, 2 ) = AB( JI, 2 )
361: NAB( KLNEW, 2 ) = NAB( JI, 2 )
362: AB( KLNEW, 1 ) = C( JI )
363: NAB( KLNEW, 1 ) = IWORK( JI )
364: AB( JI, 2 ) = C( JI )
365: NAB( JI, 2 ) = IWORK( JI )
366: ELSE
367: INFO = MMAX + 1
368: END IF
369: END IF
370: 70 CONTINUE
371: IF( INFO.NE.0 )
372: $ RETURN
373: KL = KLNEW
374: ELSE
375: *
376: * IJOB=3: Binary search. Keep only the interval containing
377: * w s.t. N(w) = NVAL
378: *
379: DO 80 JI = KF, KL
380: IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
381: AB( JI, 1 ) = C( JI )
382: NAB( JI, 1 ) = IWORK( JI )
383: END IF
384: IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
385: AB( JI, 2 ) = C( JI )
386: NAB( JI, 2 ) = IWORK( JI )
387: END IF
388: 80 CONTINUE
389: END IF
390: *
391: ELSE
392: *
393: * End of Parallel Version of the loop
394: *
395: * Begin of Serial Version of the loop
396: *
397: KLNEW = KL
398: DO 100 JI = KF, KL
399: *
400: * Compute N(w), the number of eigenvalues less than w
401: *
402: TMP1 = C( JI )
403: TMP2 = D( 1 ) - TMP1
404: ITMP1 = 0
405: IF( TMP2.LE.PIVMIN ) THEN
406: ITMP1 = 1
407: TMP2 = MIN( TMP2, -PIVMIN )
408: END IF
409: *
410: * A series of compiler directives to defeat vectorization
411: * for the next loop
412: *
413: *$PL$ CMCHAR=' '
414: CDIR$ NEXTSCALAR
415: C$DIR SCALAR
416: CDIR$ NEXT SCALAR
417: CVD$L NOVECTOR
418: CDEC$ NOVECTOR
419: CVD$ NOVECTOR
420: *VDIR NOVECTOR
421: *VOCL LOOP,SCALAR
422: CIBM PREFER SCALAR
423: *$PL$ CMCHAR='*'
424: *
425: DO 90 J = 2, N
426: TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
427: IF( TMP2.LE.PIVMIN ) THEN
428: ITMP1 = ITMP1 + 1
429: TMP2 = MIN( TMP2, -PIVMIN )
430: END IF
431: 90 CONTINUE
432: *
433: IF( IJOB.LE.2 ) THEN
434: *
435: * IJOB=2: Choose all intervals containing eigenvalues.
436: *
437: * Insure that N(w) is monotone
438: *
439: ITMP1 = MIN( NAB( JI, 2 ),
440: $ MAX( NAB( JI, 1 ), ITMP1 ) )
441: *
442: * Update the Queue -- add intervals if both halves
443: * contain eigenvalues.
444: *
445: IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
446: *
447: * No eigenvalue in the upper interval:
448: * just use the lower interval.
449: *
450: AB( JI, 2 ) = TMP1
451: *
452: ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
453: *
454: * No eigenvalue in the lower interval:
455: * just use the upper interval.
456: *
457: AB( JI, 1 ) = TMP1
458: ELSE IF( KLNEW.LT.MMAX ) THEN
459: *
460: * Eigenvalue in both intervals -- add upper to queue.
461: *
462: KLNEW = KLNEW + 1
463: AB( KLNEW, 2 ) = AB( JI, 2 )
464: NAB( KLNEW, 2 ) = NAB( JI, 2 )
465: AB( KLNEW, 1 ) = TMP1
466: NAB( KLNEW, 1 ) = ITMP1
467: AB( JI, 2 ) = TMP1
468: NAB( JI, 2 ) = ITMP1
469: ELSE
470: INFO = MMAX + 1
471: RETURN
472: END IF
473: ELSE
474: *
475: * IJOB=3: Binary search. Keep only the interval
476: * containing w s.t. N(w) = NVAL
477: *
478: IF( ITMP1.LE.NVAL( JI ) ) THEN
479: AB( JI, 1 ) = TMP1
480: NAB( JI, 1 ) = ITMP1
481: END IF
482: IF( ITMP1.GE.NVAL( JI ) ) THEN
483: AB( JI, 2 ) = TMP1
484: NAB( JI, 2 ) = ITMP1
485: END IF
486: END IF
487: 100 CONTINUE
488: KL = KLNEW
489: *
490: * End of Serial Version of the loop
491: *
492: END IF
493: *
494: * Check for convergence
495: *
496: KFNEW = KF
497: DO 110 JI = KF, KL
498: TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
499: TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
500: IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
501: $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
502: *
503: * Converged -- Swap with position KFNEW,
504: * then increment KFNEW
505: *
506: IF( JI.GT.KFNEW ) THEN
507: TMP1 = AB( JI, 1 )
508: TMP2 = AB( JI, 2 )
509: ITMP1 = NAB( JI, 1 )
510: ITMP2 = NAB( JI, 2 )
511: AB( JI, 1 ) = AB( KFNEW, 1 )
512: AB( JI, 2 ) = AB( KFNEW, 2 )
513: NAB( JI, 1 ) = NAB( KFNEW, 1 )
514: NAB( JI, 2 ) = NAB( KFNEW, 2 )
515: AB( KFNEW, 1 ) = TMP1
516: AB( KFNEW, 2 ) = TMP2
517: NAB( KFNEW, 1 ) = ITMP1
518: NAB( KFNEW, 2 ) = ITMP2
519: IF( IJOB.EQ.3 ) THEN
520: ITMP1 = NVAL( JI )
521: NVAL( JI ) = NVAL( KFNEW )
522: NVAL( KFNEW ) = ITMP1
523: END IF
524: END IF
525: KFNEW = KFNEW + 1
526: END IF
527: 110 CONTINUE
528: KF = KFNEW
529: *
530: * Choose Midpoints
531: *
532: DO 120 JI = KF, KL
533: C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
534: 120 CONTINUE
535: *
536: * If no more intervals to refine, quit.
537: *
538: IF( KF.GT.KL )
539: $ GO TO 140
540: 130 CONTINUE
541: *
542: * Converged
543: *
544: 140 CONTINUE
545: INFO = MAX( KL+1-KF, 0 )
546: MOUT = KL
547: *
548: RETURN
549: *
550: * End of DLAEBZ
551: *
552: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaebz.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack