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: * =====================================================================
262: SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
263: $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
264: $ INFO )
265: *
266: * -- LAPACK computational routine (version 3.4.0) --
267: * -- LAPACK is a software package provided by Univ. of Tennessee, --
268: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269: * November 2011
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