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