Annotation of rpl/lapack/lapack/dstebz.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DSTEBZ
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DSTEBZ + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
! 22: * M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
! 23: * INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER ORDER, RANGE
! 27: * INTEGER IL, INFO, IU, M, N, NSPLIT
! 28: * DOUBLE PRECISION ABSTOL, VL, VU
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
! 32: * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
! 33: * ..
! 34: *
! 35: *
! 36: *> \par Purpose:
! 37: * =============
! 38: *>
! 39: *> \verbatim
! 40: *>
! 41: *> DSTEBZ computes the eigenvalues of a symmetric tridiagonal
! 42: *> matrix T. The user may ask for all eigenvalues, all eigenvalues
! 43: *> in the half-open interval (VL, VU], or the IL-th through IU-th
! 44: *> eigenvalues.
! 45: *>
! 46: *> To avoid overflow, the matrix must be scaled so that its
! 47: *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
! 48: *> accuracy, it should not be much smaller than that.
! 49: *>
! 50: *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
! 51: *> Matrix", Report CS41, Computer Science Dept., Stanford
! 52: *> University, July 21, 1966.
! 53: *> \endverbatim
! 54: *
! 55: * Arguments:
! 56: * ==========
! 57: *
! 58: *> \param[in] RANGE
! 59: *> \verbatim
! 60: *> RANGE is CHARACTER*1
! 61: *> = 'A': ("All") all eigenvalues will be found.
! 62: *> = 'V': ("Value") all eigenvalues in the half-open interval
! 63: *> (VL, VU] will be found.
! 64: *> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
! 65: *> entire matrix) will be found.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[in] ORDER
! 69: *> \verbatim
! 70: *> ORDER is CHARACTER*1
! 71: *> = 'B': ("By Block") the eigenvalues will be grouped by
! 72: *> split-off block (see IBLOCK, ISPLIT) and
! 73: *> ordered from smallest to largest within
! 74: *> the block.
! 75: *> = 'E': ("Entire matrix")
! 76: *> the eigenvalues for the entire matrix
! 77: *> will be ordered from smallest to
! 78: *> largest.
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] N
! 82: *> \verbatim
! 83: *> N is INTEGER
! 84: *> The order of the tridiagonal matrix T. N >= 0.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in] VL
! 88: *> \verbatim
! 89: *> VL is DOUBLE PRECISION
! 90: *> \endverbatim
! 91: *>
! 92: *> \param[in] VU
! 93: *> \verbatim
! 94: *> VU is DOUBLE PRECISION
! 95: *>
! 96: *> If RANGE='V', the lower and upper bounds of the interval to
! 97: *> be searched for eigenvalues. Eigenvalues less than or equal
! 98: *> to VL, or greater than VU, will not be returned. VL < VU.
! 99: *> Not referenced if RANGE = 'A' or 'I'.
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in] IL
! 103: *> \verbatim
! 104: *> IL is INTEGER
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[in] IU
! 108: *> \verbatim
! 109: *> IU is INTEGER
! 110: *>
! 111: *> If RANGE='I', the indices (in ascending order) of the
! 112: *> smallest and largest eigenvalues to be returned.
! 113: *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
! 114: *> Not referenced if RANGE = 'A' or 'V'.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[in] ABSTOL
! 118: *> \verbatim
! 119: *> ABSTOL is DOUBLE PRECISION
! 120: *> The absolute tolerance for the eigenvalues. An eigenvalue
! 121: *> (or cluster) is considered to be located if it has been
! 122: *> determined to lie in an interval whose width is ABSTOL or
! 123: *> less. If ABSTOL is less than or equal to zero, then ULP*|T|
! 124: *> will be used, where |T| means the 1-norm of T.
! 125: *>
! 126: *> Eigenvalues will be computed most accurately when ABSTOL is
! 127: *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
! 128: *> \endverbatim
! 129: *>
! 130: *> \param[in] D
! 131: *> \verbatim
! 132: *> D is DOUBLE PRECISION array, dimension (N)
! 133: *> The n diagonal elements of the tridiagonal matrix T.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[in] E
! 137: *> \verbatim
! 138: *> E is DOUBLE PRECISION array, dimension (N-1)
! 139: *> The (n-1) off-diagonal elements of the tridiagonal matrix T.
! 140: *> \endverbatim
! 141: *>
! 142: *> \param[out] M
! 143: *> \verbatim
! 144: *> M is INTEGER
! 145: *> The actual number of eigenvalues found. 0 <= M <= N.
! 146: *> (See also the description of INFO=2,3.)
! 147: *> \endverbatim
! 148: *>
! 149: *> \param[out] NSPLIT
! 150: *> \verbatim
! 151: *> NSPLIT is INTEGER
! 152: *> The number of diagonal blocks in the matrix T.
! 153: *> 1 <= NSPLIT <= N.
! 154: *> \endverbatim
! 155: *>
! 156: *> \param[out] W
! 157: *> \verbatim
! 158: *> W is DOUBLE PRECISION array, dimension (N)
! 159: *> On exit, the first M elements of W will contain the
! 160: *> eigenvalues. (DSTEBZ may use the remaining N-M elements as
! 161: *> workspace.)
! 162: *> \endverbatim
! 163: *>
! 164: *> \param[out] IBLOCK
! 165: *> \verbatim
! 166: *> IBLOCK is INTEGER array, dimension (N)
! 167: *> At each row/column j where E(j) is zero or small, the
! 168: *> matrix T is considered to split into a block diagonal
! 169: *> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
! 170: *> block (from 1 to the number of blocks) the eigenvalue W(i)
! 171: *> belongs. (DSTEBZ may use the remaining N-M elements as
! 172: *> workspace.)
! 173: *> \endverbatim
! 174: *>
! 175: *> \param[out] ISPLIT
! 176: *> \verbatim
! 177: *> ISPLIT is INTEGER array, dimension (N)
! 178: *> The splitting points, at which T breaks up into submatrices.
! 179: *> The first submatrix consists of rows/columns 1 to ISPLIT(1),
! 180: *> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
! 181: *> etc., and the NSPLIT-th consists of rows/columns
! 182: *> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
! 183: *> (Only the first NSPLIT elements will actually be used, but
! 184: *> since the user cannot know a priori what value NSPLIT will
! 185: *> have, N words must be reserved for ISPLIT.)
! 186: *> \endverbatim
! 187: *>
! 188: *> \param[out] WORK
! 189: *> \verbatim
! 190: *> WORK is DOUBLE PRECISION array, dimension (4*N)
! 191: *> \endverbatim
! 192: *>
! 193: *> \param[out] IWORK
! 194: *> \verbatim
! 195: *> IWORK is INTEGER array, dimension (3*N)
! 196: *> \endverbatim
! 197: *>
! 198: *> \param[out] INFO
! 199: *> \verbatim
! 200: *> INFO is INTEGER
! 201: *> = 0: successful exit
! 202: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 203: *> > 0: some or all of the eigenvalues failed to converge or
! 204: *> were not computed:
! 205: *> =1 or 3: Bisection failed to converge for some
! 206: *> eigenvalues; these eigenvalues are flagged by a
! 207: *> negative block number. The effect is that the
! 208: *> eigenvalues may not be as accurate as the
! 209: *> absolute and relative tolerances. This is
! 210: *> generally caused by unexpectedly inaccurate
! 211: *> arithmetic.
! 212: *> =2 or 3: RANGE='I' only: Not all of the eigenvalues
! 213: *> IL:IU were found.
! 214: *> Effect: M < IU+1-IL
! 215: *> Cause: non-monotonic arithmetic, causing the
! 216: *> Sturm sequence to be non-monotonic.
! 217: *> Cure: recalculate, using RANGE='A', and pick
! 218: *> out eigenvalues IL:IU. In some cases,
! 219: *> increasing the PARAMETER "FUDGE" may
! 220: *> make things work.
! 221: *> = 4: RANGE='I', and the Gershgorin interval
! 222: *> initially used was too small. No eigenvalues
! 223: *> were computed.
! 224: *> Probable cause: your machine has sloppy
! 225: *> floating-point arithmetic.
! 226: *> Cure: Increase the PARAMETER "FUDGE",
! 227: *> recompile, and try again.
! 228: *> \endverbatim
! 229: *
! 230: *> \par Internal Parameters:
! 231: * =========================
! 232: *>
! 233: *> \verbatim
! 234: *> RELFAC DOUBLE PRECISION, default = 2.0e0
! 235: *> The relative tolerance. An interval (a,b] lies within
! 236: *> "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
! 237: *> where "ulp" is the machine precision (distance from 1 to
! 238: *> the next larger floating point number.)
! 239: *>
! 240: *> FUDGE DOUBLE PRECISION, default = 2
! 241: *> A "fudge factor" to widen the Gershgorin intervals. Ideally,
! 242: *> a value of 1 should work, but on machines with sloppy
! 243: *> arithmetic, this needs to be larger. The default for
! 244: *> publicly released versions should be large enough to handle
! 245: *> the worst machine around. Note that this has no effect
! 246: *> on accuracy of the solution.
! 247: *> \endverbatim
! 248: *
! 249: * Authors:
! 250: * ========
! 251: *
! 252: *> \author Univ. of Tennessee
! 253: *> \author Univ. of California Berkeley
! 254: *> \author Univ. of Colorado Denver
! 255: *> \author NAG Ltd.
! 256: *
! 257: *> \date November 2011
! 258: *
! 259: *> \ingroup auxOTHERcomputational
! 260: *
! 261: * =====================================================================
1.1 bertrand 262: SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
263: $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
264: $ INFO )
265: *
1.9 ! bertrand 266: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 267: * -- LAPACK is a software package provided by Univ. of Tennessee, --
268: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 269: * November 2011
1.1 bertrand 270: *
271: * .. Scalar Arguments ..
272: CHARACTER ORDER, RANGE
273: INTEGER IL, INFO, IU, M, N, NSPLIT
274: DOUBLE PRECISION ABSTOL, VL, VU
275: * ..
276: * .. Array Arguments ..
277: INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
278: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
279: * ..
280: *
281: * =====================================================================
282: *
283: * .. Parameters ..
284: DOUBLE PRECISION ZERO, ONE, TWO, HALF
285: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
286: $ HALF = 1.0D0 / TWO )
287: DOUBLE PRECISION FUDGE, RELFAC
288: PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 )
289: * ..
290: * .. Local Scalars ..
291: LOGICAL NCNVRG, TOOFEW
292: INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
293: $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
294: $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
295: $ NWU
296: DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
297: $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
298: * ..
299: * .. Local Arrays ..
300: INTEGER IDUMMA( 1 )
301: * ..
302: * .. External Functions ..
303: LOGICAL LSAME
304: INTEGER ILAENV
305: DOUBLE PRECISION DLAMCH
306: EXTERNAL LSAME, ILAENV, DLAMCH
307: * ..
308: * .. External Subroutines ..
309: EXTERNAL DLAEBZ, XERBLA
310: * ..
311: * .. Intrinsic Functions ..
312: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
313: * ..
314: * .. Executable Statements ..
315: *
316: INFO = 0
317: *
318: * Decode RANGE
319: *
320: IF( LSAME( RANGE, 'A' ) ) THEN
321: IRANGE = 1
322: ELSE IF( LSAME( RANGE, 'V' ) ) THEN
323: IRANGE = 2
324: ELSE IF( LSAME( RANGE, 'I' ) ) THEN
325: IRANGE = 3
326: ELSE
327: IRANGE = 0
328: END IF
329: *
330: * Decode ORDER
331: *
332: IF( LSAME( ORDER, 'B' ) ) THEN
333: IORDER = 2
334: ELSE IF( LSAME( ORDER, 'E' ) ) THEN
335: IORDER = 1
336: ELSE
337: IORDER = 0
338: END IF
339: *
340: * Check for Errors
341: *
342: IF( IRANGE.LE.0 ) THEN
343: INFO = -1
344: ELSE IF( IORDER.LE.0 ) THEN
345: INFO = -2
346: ELSE IF( N.LT.0 ) THEN
347: INFO = -3
348: ELSE IF( IRANGE.EQ.2 ) THEN
349: IF( VL.GE.VU )
350: $ INFO = -5
351: ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
352: $ THEN
353: INFO = -6
354: ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
355: $ THEN
356: INFO = -7
357: END IF
358: *
359: IF( INFO.NE.0 ) THEN
360: CALL XERBLA( 'DSTEBZ', -INFO )
361: RETURN
362: END IF
363: *
364: * Initialize error flags
365: *
366: INFO = 0
367: NCNVRG = .FALSE.
368: TOOFEW = .FALSE.
369: *
370: * Quick return if possible
371: *
372: M = 0
373: IF( N.EQ.0 )
374: $ RETURN
375: *
376: * Simplifications:
377: *
378: IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
379: $ IRANGE = 1
380: *
381: * Get machine constants
382: * NB is the minimum vector length for vector bisection, or 0
383: * if only scalar is to be done.
384: *
385: SAFEMN = DLAMCH( 'S' )
386: ULP = DLAMCH( 'P' )
387: RTOLI = ULP*RELFAC
388: NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
389: IF( NB.LE.1 )
390: $ NB = 0
391: *
392: * Special Case when N=1
393: *
394: IF( N.EQ.1 ) THEN
395: NSPLIT = 1
396: ISPLIT( 1 ) = 1
397: IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
398: M = 0
399: ELSE
400: W( 1 ) = D( 1 )
401: IBLOCK( 1 ) = 1
402: M = 1
403: END IF
404: RETURN
405: END IF
406: *
407: * Compute Splitting Points
408: *
409: NSPLIT = 1
410: WORK( N ) = ZERO
411: PIVMIN = ONE
412: *
413: DO 10 J = 2, N
414: TMP1 = E( J-1 )**2
415: IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
416: ISPLIT( NSPLIT ) = J - 1
417: NSPLIT = NSPLIT + 1
418: WORK( J-1 ) = ZERO
419: ELSE
420: WORK( J-1 ) = TMP1
421: PIVMIN = MAX( PIVMIN, TMP1 )
422: END IF
423: 10 CONTINUE
424: ISPLIT( NSPLIT ) = N
425: PIVMIN = PIVMIN*SAFEMN
426: *
427: * Compute Interval and ATOLI
428: *
429: IF( IRANGE.EQ.3 ) THEN
430: *
431: * RANGE='I': Compute the interval containing eigenvalues
432: * IL through IU.
433: *
434: * Compute Gershgorin interval for entire (split) matrix
435: * and use it as the initial interval
436: *
437: GU = D( 1 )
438: GL = D( 1 )
439: TMP1 = ZERO
440: *
441: DO 20 J = 1, N - 1
442: TMP2 = SQRT( WORK( J ) )
443: GU = MAX( GU, D( J )+TMP1+TMP2 )
444: GL = MIN( GL, D( J )-TMP1-TMP2 )
445: TMP1 = TMP2
446: 20 CONTINUE
447: *
448: GU = MAX( GU, D( N )+TMP1 )
449: GL = MIN( GL, D( N )-TMP1 )
450: TNORM = MAX( ABS( GL ), ABS( GU ) )
451: GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
452: GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
453: *
454: * Compute Iteration parameters
455: *
456: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
457: $ LOG( TWO ) ) + 2
458: IF( ABSTOL.LE.ZERO ) THEN
459: ATOLI = ULP*TNORM
460: ELSE
461: ATOLI = ABSTOL
462: END IF
463: *
464: WORK( N+1 ) = GL
465: WORK( N+2 ) = GL
466: WORK( N+3 ) = GU
467: WORK( N+4 ) = GU
468: WORK( N+5 ) = GL
469: WORK( N+6 ) = GU
470: IWORK( 1 ) = -1
471: IWORK( 2 ) = -1
472: IWORK( 3 ) = N + 1
473: IWORK( 4 ) = N + 1
474: IWORK( 5 ) = IL - 1
475: IWORK( 6 ) = IU
476: *
477: CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
478: $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
479: $ IWORK, W, IBLOCK, IINFO )
480: *
481: IF( IWORK( 6 ).EQ.IU ) THEN
482: WL = WORK( N+1 )
483: WLU = WORK( N+3 )
484: NWL = IWORK( 1 )
485: WU = WORK( N+4 )
486: WUL = WORK( N+2 )
487: NWU = IWORK( 4 )
488: ELSE
489: WL = WORK( N+2 )
490: WLU = WORK( N+4 )
491: NWL = IWORK( 2 )
492: WU = WORK( N+3 )
493: WUL = WORK( N+1 )
494: NWU = IWORK( 3 )
495: END IF
496: *
497: IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
498: INFO = 4
499: RETURN
500: END IF
501: ELSE
502: *
503: * RANGE='A' or 'V' -- Set ATOLI
504: *
505: TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
506: $ ABS( D( N ) )+ABS( E( N-1 ) ) )
507: *
508: DO 30 J = 2, N - 1
509: TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
510: $ ABS( E( J ) ) )
511: 30 CONTINUE
512: *
513: IF( ABSTOL.LE.ZERO ) THEN
514: ATOLI = ULP*TNORM
515: ELSE
516: ATOLI = ABSTOL
517: END IF
518: *
519: IF( IRANGE.EQ.2 ) THEN
520: WL = VL
521: WU = VU
522: ELSE
523: WL = ZERO
524: WU = ZERO
525: END IF
526: END IF
527: *
528: * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
529: * NWL accumulates the number of eigenvalues .le. WL,
530: * NWU accumulates the number of eigenvalues .le. WU
531: *
532: M = 0
533: IEND = 0
534: INFO = 0
535: NWL = 0
536: NWU = 0
537: *
538: DO 70 JB = 1, NSPLIT
539: IOFF = IEND
540: IBEGIN = IOFF + 1
541: IEND = ISPLIT( JB )
542: IN = IEND - IOFF
543: *
544: IF( IN.EQ.1 ) THEN
545: *
546: * Special Case -- IN=1
547: *
548: IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
549: $ NWL = NWL + 1
550: IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
551: $ NWU = NWU + 1
552: IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
553: $ D( IBEGIN )-PIVMIN ) ) THEN
554: M = M + 1
555: W( M ) = D( IBEGIN )
556: IBLOCK( M ) = JB
557: END IF
558: ELSE
559: *
560: * General Case -- IN > 1
561: *
562: * Compute Gershgorin Interval
563: * and use it as the initial interval
564: *
565: GU = D( IBEGIN )
566: GL = D( IBEGIN )
567: TMP1 = ZERO
568: *
569: DO 40 J = IBEGIN, IEND - 1
570: TMP2 = ABS( E( J ) )
571: GU = MAX( GU, D( J )+TMP1+TMP2 )
572: GL = MIN( GL, D( J )-TMP1-TMP2 )
573: TMP1 = TMP2
574: 40 CONTINUE
575: *
576: GU = MAX( GU, D( IEND )+TMP1 )
577: GL = MIN( GL, D( IEND )-TMP1 )
578: BNORM = MAX( ABS( GL ), ABS( GU ) )
579: GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
580: GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
581: *
582: * Compute ATOLI for the current submatrix
583: *
584: IF( ABSTOL.LE.ZERO ) THEN
585: ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
586: ELSE
587: ATOLI = ABSTOL
588: END IF
589: *
590: IF( IRANGE.GT.1 ) THEN
591: IF( GU.LT.WL ) THEN
592: NWL = NWL + IN
593: NWU = NWU + IN
594: GO TO 70
595: END IF
596: GL = MAX( GL, WL )
597: GU = MIN( GU, WU )
598: IF( GL.GE.GU )
599: $ GO TO 70
600: END IF
601: *
602: * Set Up Initial Interval
603: *
604: WORK( N+1 ) = GL
605: WORK( N+IN+1 ) = GU
606: CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
607: $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
608: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
609: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
610: *
611: NWL = NWL + IWORK( 1 )
612: NWU = NWU + IWORK( IN+1 )
613: IWOFF = M - IWORK( 1 )
614: *
615: * Compute Eigenvalues
616: *
617: ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
618: $ LOG( TWO ) ) + 2
619: CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
620: $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
621: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
622: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
623: *
624: * Copy Eigenvalues Into W and IBLOCK
625: * Use -JB for block number for unconverged eigenvalues.
626: *
627: DO 60 J = 1, IOUT
628: TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
629: *
630: * Flag non-convergence.
631: *
632: IF( J.GT.IOUT-IINFO ) THEN
633: NCNVRG = .TRUE.
634: IB = -JB
635: ELSE
636: IB = JB
637: END IF
638: DO 50 JE = IWORK( J ) + 1 + IWOFF,
639: $ IWORK( J+IN ) + IWOFF
640: W( JE ) = TMP1
641: IBLOCK( JE ) = IB
642: 50 CONTINUE
643: 60 CONTINUE
644: *
645: M = M + IM
646: END IF
647: 70 CONTINUE
648: *
649: * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
650: * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
651: *
652: IF( IRANGE.EQ.3 ) THEN
653: IM = 0
654: IDISCL = IL - 1 - NWL
655: IDISCU = NWU - IU
656: *
657: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
658: DO 80 JE = 1, M
659: IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
660: IDISCL = IDISCL - 1
661: ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
662: IDISCU = IDISCU - 1
663: ELSE
664: IM = IM + 1
665: W( IM ) = W( JE )
666: IBLOCK( IM ) = IBLOCK( JE )
667: END IF
668: 80 CONTINUE
669: M = IM
670: END IF
671: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
672: *
673: * Code to deal with effects of bad arithmetic:
674: * Some low eigenvalues to be discarded are not in (WL,WLU],
675: * or high eigenvalues to be discarded are not in (WUL,WU]
676: * so just kill off the smallest IDISCL/largest IDISCU
677: * eigenvalues, by simply finding the smallest/largest
678: * eigenvalue(s).
679: *
680: * (If N(w) is monotone non-decreasing, this should never
681: * happen.)
682: *
683: IF( IDISCL.GT.0 ) THEN
684: WKILL = WU
685: DO 100 JDISC = 1, IDISCL
686: IW = 0
687: DO 90 JE = 1, M
688: IF( IBLOCK( JE ).NE.0 .AND.
689: $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
690: IW = JE
691: WKILL = W( JE )
692: END IF
693: 90 CONTINUE
694: IBLOCK( IW ) = 0
695: 100 CONTINUE
696: END IF
697: IF( IDISCU.GT.0 ) THEN
698: *
699: WKILL = WL
700: DO 120 JDISC = 1, IDISCU
701: IW = 0
702: DO 110 JE = 1, M
703: IF( IBLOCK( JE ).NE.0 .AND.
704: $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
705: IW = JE
706: WKILL = W( JE )
707: END IF
708: 110 CONTINUE
709: IBLOCK( IW ) = 0
710: 120 CONTINUE
711: END IF
712: IM = 0
713: DO 130 JE = 1, M
714: IF( IBLOCK( JE ).NE.0 ) THEN
715: IM = IM + 1
716: W( IM ) = W( JE )
717: IBLOCK( IM ) = IBLOCK( JE )
718: END IF
719: 130 CONTINUE
720: M = IM
721: END IF
722: IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
723: TOOFEW = .TRUE.
724: END IF
725: END IF
726: *
727: * If ORDER='B', do nothing -- the eigenvalues are already sorted
728: * by block.
729: * If ORDER='E', sort the eigenvalues from smallest to largest
730: *
731: IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
732: DO 150 JE = 1, M - 1
733: IE = 0
734: TMP1 = W( JE )
735: DO 140 J = JE + 1, M
736: IF( W( J ).LT.TMP1 ) THEN
737: IE = J
738: TMP1 = W( J )
739: END IF
740: 140 CONTINUE
741: *
742: IF( IE.NE.0 ) THEN
743: ITMP1 = IBLOCK( IE )
744: W( IE ) = W( JE )
745: IBLOCK( IE ) = IBLOCK( JE )
746: W( JE ) = TMP1
747: IBLOCK( JE ) = ITMP1
748: END IF
749: 150 CONTINUE
750: END IF
751: *
752: INFO = 0
753: IF( NCNVRG )
754: $ INFO = INFO + 1
755: IF( TOOFEW )
756: $ INFO = INFO + 2
757: RETURN
758: *
759: * End of DSTEBZ
760: *
761: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dstebz.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack