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