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.
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: *
276: *> \ingroup OTHERauxiliary
277: *
278: *> \par Further Details:
279: * =====================
280: *>
281: *> \verbatim
282: *>
283: *> This routine is intended to be called only by other LAPACK
284: *> routines, thus the interface is less user-friendly. It is intended
285: *> for two purposes:
286: *>
287: *> (a) finding eigenvalues. In this case, DLAEBZ should have one or
288: *> more initial intervals set up in AB, and DLAEBZ should be called
289: *> with IJOB=1. This sets up NAB, and also counts the eigenvalues.
290: *> Intervals with no eigenvalues would usually be thrown out at
291: *> this point. Also, if not all the eigenvalues in an interval i
292: *> are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
293: *> For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
294: *> eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX
295: *> no smaller than the value of MOUT returned by the call with
296: *> IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
297: *> through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
298: *> tolerance specified by ABSTOL and RELTOL.
299: *>
300: *> (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
301: *> In this case, start with a Gershgorin interval (a,b). Set up
302: *> AB to contain 2 search intervals, both initially (a,b). One
303: *> NVAL element should contain f-1 and the other should contain l
304: *> , while C should contain a and b, resp. NAB(i,1) should be -1
305: *> and NAB(i,2) should be N+1, to flag an error if the desired
306: *> interval does not lie in (a,b). DLAEBZ is then called with
307: *> IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
308: *> j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
309: *> if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
310: *> >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
311: *> N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
312: *> w(l-r)=...=w(l+k) are handled similarly.
313: *> \endverbatim
314: *>
315: * =====================================================================
316: SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
317: $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
318: $ NAB, WORK, IWORK, INFO )
319: *
320: * -- LAPACK auxiliary routine --
321: * -- LAPACK is a software package provided by Univ. of Tennessee, --
322: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
323: *
324: * .. Scalar Arguments ..
325: INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
326: DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
327: * ..
328: * .. Array Arguments ..
329: INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
330: DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
331: $ WORK( * )
332: * ..
333: *
334: * =====================================================================
335: *
336: * .. Parameters ..
337: DOUBLE PRECISION ZERO, TWO, HALF
338: PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
339: $ HALF = 1.0D0 / TWO )
340: * ..
341: * .. Local Scalars ..
342: INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
343: $ KLNEW
344: DOUBLE PRECISION TMP1, TMP2
345: * ..
346: * .. Intrinsic Functions ..
347: INTRINSIC ABS, MAX, MIN
348: * ..
349: * .. Executable Statements ..
350: *
351: * Check for Errors
352: *
353: INFO = 0
354: IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
355: INFO = -1
356: RETURN
357: END IF
358: *
359: * Initialize NAB
360: *
361: IF( IJOB.EQ.1 ) THEN
362: *
363: * Compute the number of eigenvalues in the initial intervals.
364: *
365: MOUT = 0
366: DO 30 JI = 1, MINP
367: DO 20 JP = 1, 2
368: TMP1 = D( 1 ) - AB( JI, JP )
369: IF( ABS( TMP1 ).LT.PIVMIN )
370: $ TMP1 = -PIVMIN
371: NAB( JI, JP ) = 0
372: IF( TMP1.LE.ZERO )
373: $ NAB( JI, JP ) = 1
374: *
375: DO 10 J = 2, N
376: TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
377: IF( ABS( TMP1 ).LT.PIVMIN )
378: $ TMP1 = -PIVMIN
379: IF( TMP1.LE.ZERO )
380: $ NAB( JI, JP ) = NAB( JI, JP ) + 1
381: 10 CONTINUE
382: 20 CONTINUE
383: MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
384: 30 CONTINUE
385: RETURN
386: END IF
387: *
388: * Initialize for loop
389: *
390: * KF and KL have the following meaning:
391: * Intervals 1,...,KF-1 have converged.
392: * Intervals KF,...,KL still need to be refined.
393: *
394: KF = 1
395: KL = MINP
396: *
397: * If IJOB=2, initialize C.
398: * If IJOB=3, use the user-supplied starting point.
399: *
400: IF( IJOB.EQ.2 ) THEN
401: DO 40 JI = 1, MINP
402: C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
403: 40 CONTINUE
404: END IF
405: *
406: * Iteration loop
407: *
408: DO 130 JIT = 1, NITMAX
409: *
410: * Loop over intervals
411: *
412: IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
413: *
414: * Begin of Parallel Version of the loop
415: *
416: DO 60 JI = KF, KL
417: *
418: * Compute N(c), the number of eigenvalues less than c
419: *
420: WORK( JI ) = D( 1 ) - C( JI )
421: IWORK( JI ) = 0
422: IF( WORK( JI ).LE.PIVMIN ) THEN
423: IWORK( JI ) = 1
424: WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
425: END IF
426: *
427: DO 50 J = 2, N
428: WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
429: IF( WORK( JI ).LE.PIVMIN ) THEN
430: IWORK( JI ) = IWORK( JI ) + 1
431: WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
432: END IF
433: 50 CONTINUE
434: 60 CONTINUE
435: *
436: IF( IJOB.LE.2 ) THEN
437: *
438: * IJOB=2: Choose all intervals containing eigenvalues.
439: *
440: KLNEW = KL
441: DO 70 JI = KF, KL
442: *
443: * Insure that N(w) is monotone
444: *
445: IWORK( JI ) = MIN( NAB( JI, 2 ),
446: $ MAX( NAB( JI, 1 ), IWORK( JI ) ) )
447: *
448: * Update the Queue -- add intervals if both halves
449: * contain eigenvalues.
450: *
451: IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
452: *
453: * No eigenvalue in the upper interval:
454: * just use the lower interval.
455: *
456: AB( JI, 2 ) = C( JI )
457: *
458: ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
459: *
460: * No eigenvalue in the lower interval:
461: * just use the upper interval.
462: *
463: AB( JI, 1 ) = C( JI )
464: ELSE
465: KLNEW = KLNEW + 1
466: IF( KLNEW.LE.MMAX ) THEN
467: *
468: * Eigenvalue in both intervals -- add upper to
469: * queue.
470: *
471: AB( KLNEW, 2 ) = AB( JI, 2 )
472: NAB( KLNEW, 2 ) = NAB( JI, 2 )
473: AB( KLNEW, 1 ) = C( JI )
474: NAB( KLNEW, 1 ) = IWORK( JI )
475: AB( JI, 2 ) = C( JI )
476: NAB( JI, 2 ) = IWORK( JI )
477: ELSE
478: INFO = MMAX + 1
479: END IF
480: END IF
481: 70 CONTINUE
482: IF( INFO.NE.0 )
483: $ RETURN
484: KL = KLNEW
485: ELSE
486: *
487: * IJOB=3: Binary search. Keep only the interval containing
488: * w s.t. N(w) = NVAL
489: *
490: DO 80 JI = KF, KL
491: IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
492: AB( JI, 1 ) = C( JI )
493: NAB( JI, 1 ) = IWORK( JI )
494: END IF
495: IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
496: AB( JI, 2 ) = C( JI )
497: NAB( JI, 2 ) = IWORK( JI )
498: END IF
499: 80 CONTINUE
500: END IF
501: *
502: ELSE
503: *
504: * End of Parallel Version of the loop
505: *
506: * Begin of Serial Version of the loop
507: *
508: KLNEW = KL
509: DO 100 JI = KF, KL
510: *
511: * Compute N(w), the number of eigenvalues less than w
512: *
513: TMP1 = C( JI )
514: TMP2 = D( 1 ) - TMP1
515: ITMP1 = 0
516: IF( TMP2.LE.PIVMIN ) THEN
517: ITMP1 = 1
518: TMP2 = MIN( TMP2, -PIVMIN )
519: END IF
520: *
521: DO 90 J = 2, N
522: TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
523: IF( TMP2.LE.PIVMIN ) THEN
524: ITMP1 = ITMP1 + 1
525: TMP2 = MIN( TMP2, -PIVMIN )
526: END IF
527: 90 CONTINUE
528: *
529: IF( IJOB.LE.2 ) THEN
530: *
531: * IJOB=2: Choose all intervals containing eigenvalues.
532: *
533: * Insure that N(w) is monotone
534: *
535: ITMP1 = MIN( NAB( JI, 2 ),
536: $ MAX( NAB( JI, 1 ), ITMP1 ) )
537: *
538: * Update the Queue -- add intervals if both halves
539: * contain eigenvalues.
540: *
541: IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
542: *
543: * No eigenvalue in the upper interval:
544: * just use the lower interval.
545: *
546: AB( JI, 2 ) = TMP1
547: *
548: ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
549: *
550: * No eigenvalue in the lower interval:
551: * just use the upper interval.
552: *
553: AB( JI, 1 ) = TMP1
554: ELSE IF( KLNEW.LT.MMAX ) THEN
555: *
556: * Eigenvalue in both intervals -- add upper to queue.
557: *
558: KLNEW = KLNEW + 1
559: AB( KLNEW, 2 ) = AB( JI, 2 )
560: NAB( KLNEW, 2 ) = NAB( JI, 2 )
561: AB( KLNEW, 1 ) = TMP1
562: NAB( KLNEW, 1 ) = ITMP1
563: AB( JI, 2 ) = TMP1
564: NAB( JI, 2 ) = ITMP1
565: ELSE
566: INFO = MMAX + 1
567: RETURN
568: END IF
569: ELSE
570: *
571: * IJOB=3: Binary search. Keep only the interval
572: * containing w s.t. N(w) = NVAL
573: *
574: IF( ITMP1.LE.NVAL( JI ) ) THEN
575: AB( JI, 1 ) = TMP1
576: NAB( JI, 1 ) = ITMP1
577: END IF
578: IF( ITMP1.GE.NVAL( JI ) ) THEN
579: AB( JI, 2 ) = TMP1
580: NAB( JI, 2 ) = ITMP1
581: END IF
582: END IF
583: 100 CONTINUE
584: KL = KLNEW
585: *
586: END IF
587: *
588: * Check for convergence
589: *
590: KFNEW = KF
591: DO 110 JI = KF, KL
592: TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
593: TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
594: IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
595: $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
596: *
597: * Converged -- Swap with position KFNEW,
598: * then increment KFNEW
599: *
600: IF( JI.GT.KFNEW ) THEN
601: TMP1 = AB( JI, 1 )
602: TMP2 = AB( JI, 2 )
603: ITMP1 = NAB( JI, 1 )
604: ITMP2 = NAB( JI, 2 )
605: AB( JI, 1 ) = AB( KFNEW, 1 )
606: AB( JI, 2 ) = AB( KFNEW, 2 )
607: NAB( JI, 1 ) = NAB( KFNEW, 1 )
608: NAB( JI, 2 ) = NAB( KFNEW, 2 )
609: AB( KFNEW, 1 ) = TMP1
610: AB( KFNEW, 2 ) = TMP2
611: NAB( KFNEW, 1 ) = ITMP1
612: NAB( KFNEW, 2 ) = ITMP2
613: IF( IJOB.EQ.3 ) THEN
614: ITMP1 = NVAL( JI )
615: NVAL( JI ) = NVAL( KFNEW )
616: NVAL( KFNEW ) = ITMP1
617: END IF
618: END IF
619: KFNEW = KFNEW + 1
620: END IF
621: 110 CONTINUE
622: KF = KFNEW
623: *
624: * Choose Midpoints
625: *
626: DO 120 JI = KF, KL
627: C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
628: 120 CONTINUE
629: *
630: * If no more intervals to refine, quit.
631: *
632: IF( KF.GT.KL )
633: $ GO TO 140
634: 130 CONTINUE
635: *
636: * Converged
637: *
638: 140 CONTINUE
639: INFO = MAX( KL+1-KF, 0 )
640: MOUT = KL
641: *
642: RETURN
643: *
644: * End of DLAEBZ
645: *
646: 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