Annotation of rpl/lapack/lapack/dlaebz.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAEBZ
! 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: *> \date November 2011
! 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.9 ! bertrand 322: * -- LAPACK auxiliary routine (version 3.4.0) --
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.9 ! bertrand 325: * November 2011
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>