Annotation of rpl/lapack/lapack/dstebz.f, revision 1.15
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
1.14 bertrand 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'.
1.9 bertrand 95: *> \endverbatim
96: *>
97: *> \param[in] VU
98: *> \verbatim
99: *> VU is DOUBLE PRECISION
100: *>
1.14 bertrand 101: *> If RANGE='V', the upper bound of the interval to
1.9 bertrand 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
1.14 bertrand 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'.
1.9 bertrand 115: *> \endverbatim
116: *>
117: *> \param[in] IU
118: *> \verbatim
119: *> IU is INTEGER
120: *>
1.14 bertrand 121: *> If RANGE='I', the index of the
122: *> largest eigenvalue to be returned.
1.9 bertrand 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: *
1.14 bertrand 267: *> \date June 2016
1.9 bertrand 268: *
269: *> \ingroup auxOTHERcomputational
270: *
271: * =====================================================================
1.1 bertrand 272: SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
273: $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
274: $ INFO )
275: *
1.14 bertrand 276: * -- LAPACK computational routine (version 3.6.1) --
1.1 bertrand 277: * -- LAPACK is a software package provided by Univ. of Tennessee, --
278: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.14 bertrand 279: * June 2016
1.1 bertrand 280: *
281: * .. Scalar Arguments ..
282: CHARACTER ORDER, RANGE
283: INTEGER IL, INFO, IU, M, N, NSPLIT
284: DOUBLE PRECISION ABSTOL, VL, VU
285: * ..
286: * .. Array Arguments ..
287: INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
288: DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
289: * ..
290: *
291: * =====================================================================
292: *
293: * .. Parameters ..
294: DOUBLE PRECISION ZERO, ONE, TWO, HALF
295: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
296: $ HALF = 1.0D0 / TWO )
297: DOUBLE PRECISION FUDGE, RELFAC
298: PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 )
299: * ..
300: * .. Local Scalars ..
301: LOGICAL NCNVRG, TOOFEW
302: INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
303: $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
304: $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
305: $ NWU
306: DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
307: $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
308: * ..
309: * .. Local Arrays ..
310: INTEGER IDUMMA( 1 )
311: * ..
312: * .. External Functions ..
313: LOGICAL LSAME
314: INTEGER ILAENV
315: DOUBLE PRECISION DLAMCH
316: EXTERNAL LSAME, ILAENV, DLAMCH
317: * ..
318: * .. External Subroutines ..
319: EXTERNAL DLAEBZ, XERBLA
320: * ..
321: * .. Intrinsic Functions ..
322: INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
323: * ..
324: * .. Executable Statements ..
325: *
326: INFO = 0
327: *
328: * Decode RANGE
329: *
330: IF( LSAME( RANGE, 'A' ) ) THEN
331: IRANGE = 1
332: ELSE IF( LSAME( RANGE, 'V' ) ) THEN
333: IRANGE = 2
334: ELSE IF( LSAME( RANGE, 'I' ) ) THEN
335: IRANGE = 3
336: ELSE
337: IRANGE = 0
338: END IF
339: *
340: * Decode ORDER
341: *
342: IF( LSAME( ORDER, 'B' ) ) THEN
343: IORDER = 2
344: ELSE IF( LSAME( ORDER, 'E' ) ) THEN
345: IORDER = 1
346: ELSE
347: IORDER = 0
348: END IF
349: *
350: * Check for Errors
351: *
352: IF( IRANGE.LE.0 ) THEN
353: INFO = -1
354: ELSE IF( IORDER.LE.0 ) THEN
355: INFO = -2
356: ELSE IF( N.LT.0 ) THEN
357: INFO = -3
358: ELSE IF( IRANGE.EQ.2 ) THEN
359: IF( VL.GE.VU )
360: $ INFO = -5
361: ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
362: $ THEN
363: INFO = -6
364: ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
365: $ THEN
366: INFO = -7
367: END IF
368: *
369: IF( INFO.NE.0 ) THEN
370: CALL XERBLA( 'DSTEBZ', -INFO )
371: RETURN
372: END IF
373: *
374: * Initialize error flags
375: *
376: INFO = 0
377: NCNVRG = .FALSE.
378: TOOFEW = .FALSE.
379: *
380: * Quick return if possible
381: *
382: M = 0
383: IF( N.EQ.0 )
384: $ RETURN
385: *
386: * Simplifications:
387: *
388: IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
389: $ IRANGE = 1
390: *
391: * Get machine constants
392: * NB is the minimum vector length for vector bisection, or 0
393: * if only scalar is to be done.
394: *
395: SAFEMN = DLAMCH( 'S' )
396: ULP = DLAMCH( 'P' )
397: RTOLI = ULP*RELFAC
398: NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
399: IF( NB.LE.1 )
400: $ NB = 0
401: *
402: * Special Case when N=1
403: *
404: IF( N.EQ.1 ) THEN
405: NSPLIT = 1
406: ISPLIT( 1 ) = 1
407: IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
408: M = 0
409: ELSE
410: W( 1 ) = D( 1 )
411: IBLOCK( 1 ) = 1
412: M = 1
413: END IF
414: RETURN
415: END IF
416: *
417: * Compute Splitting Points
418: *
419: NSPLIT = 1
420: WORK( N ) = ZERO
421: PIVMIN = ONE
422: *
423: DO 10 J = 2, N
424: TMP1 = E( J-1 )**2
425: IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
426: ISPLIT( NSPLIT ) = J - 1
427: NSPLIT = NSPLIT + 1
428: WORK( J-1 ) = ZERO
429: ELSE
430: WORK( J-1 ) = TMP1
431: PIVMIN = MAX( PIVMIN, TMP1 )
432: END IF
433: 10 CONTINUE
434: ISPLIT( NSPLIT ) = N
435: PIVMIN = PIVMIN*SAFEMN
436: *
437: * Compute Interval and ATOLI
438: *
439: IF( IRANGE.EQ.3 ) THEN
440: *
441: * RANGE='I': Compute the interval containing eigenvalues
442: * IL through IU.
443: *
444: * Compute Gershgorin interval for entire (split) matrix
445: * and use it as the initial interval
446: *
447: GU = D( 1 )
448: GL = D( 1 )
449: TMP1 = ZERO
450: *
451: DO 20 J = 1, N - 1
452: TMP2 = SQRT( WORK( J ) )
453: GU = MAX( GU, D( J )+TMP1+TMP2 )
454: GL = MIN( GL, D( J )-TMP1-TMP2 )
455: TMP1 = TMP2
456: 20 CONTINUE
457: *
458: GU = MAX( GU, D( N )+TMP1 )
459: GL = MIN( GL, D( N )-TMP1 )
460: TNORM = MAX( ABS( GL ), ABS( GU ) )
461: GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
462: GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
463: *
464: * Compute Iteration parameters
465: *
466: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
467: $ LOG( TWO ) ) + 2
468: IF( ABSTOL.LE.ZERO ) THEN
469: ATOLI = ULP*TNORM
470: ELSE
471: ATOLI = ABSTOL
472: END IF
473: *
474: WORK( N+1 ) = GL
475: WORK( N+2 ) = GL
476: WORK( N+3 ) = GU
477: WORK( N+4 ) = GU
478: WORK( N+5 ) = GL
479: WORK( N+6 ) = GU
480: IWORK( 1 ) = -1
481: IWORK( 2 ) = -1
482: IWORK( 3 ) = N + 1
483: IWORK( 4 ) = N + 1
484: IWORK( 5 ) = IL - 1
485: IWORK( 6 ) = IU
486: *
487: CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
488: $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
489: $ IWORK, W, IBLOCK, IINFO )
490: *
491: IF( IWORK( 6 ).EQ.IU ) THEN
492: WL = WORK( N+1 )
493: WLU = WORK( N+3 )
494: NWL = IWORK( 1 )
495: WU = WORK( N+4 )
496: WUL = WORK( N+2 )
497: NWU = IWORK( 4 )
498: ELSE
499: WL = WORK( N+2 )
500: WLU = WORK( N+4 )
501: NWL = IWORK( 2 )
502: WU = WORK( N+3 )
503: WUL = WORK( N+1 )
504: NWU = IWORK( 3 )
505: END IF
506: *
507: IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
508: INFO = 4
509: RETURN
510: END IF
511: ELSE
512: *
513: * RANGE='A' or 'V' -- Set ATOLI
514: *
515: TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
516: $ ABS( D( N ) )+ABS( E( N-1 ) ) )
517: *
518: DO 30 J = 2, N - 1
519: TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
520: $ ABS( E( J ) ) )
521: 30 CONTINUE
522: *
523: IF( ABSTOL.LE.ZERO ) THEN
524: ATOLI = ULP*TNORM
525: ELSE
526: ATOLI = ABSTOL
527: END IF
528: *
529: IF( IRANGE.EQ.2 ) THEN
530: WL = VL
531: WU = VU
532: ELSE
533: WL = ZERO
534: WU = ZERO
535: END IF
536: END IF
537: *
538: * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
539: * NWL accumulates the number of eigenvalues .le. WL,
540: * NWU accumulates the number of eigenvalues .le. WU
541: *
542: M = 0
543: IEND = 0
544: INFO = 0
545: NWL = 0
546: NWU = 0
547: *
548: DO 70 JB = 1, NSPLIT
549: IOFF = IEND
550: IBEGIN = IOFF + 1
551: IEND = ISPLIT( JB )
552: IN = IEND - IOFF
553: *
554: IF( IN.EQ.1 ) THEN
555: *
556: * Special Case -- IN=1
557: *
558: IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
559: $ NWL = NWL + 1
560: IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
561: $ NWU = NWU + 1
562: IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
563: $ D( IBEGIN )-PIVMIN ) ) THEN
564: M = M + 1
565: W( M ) = D( IBEGIN )
566: IBLOCK( M ) = JB
567: END IF
568: ELSE
569: *
570: * General Case -- IN > 1
571: *
572: * Compute Gershgorin Interval
573: * and use it as the initial interval
574: *
575: GU = D( IBEGIN )
576: GL = D( IBEGIN )
577: TMP1 = ZERO
578: *
579: DO 40 J = IBEGIN, IEND - 1
580: TMP2 = ABS( E( J ) )
581: GU = MAX( GU, D( J )+TMP1+TMP2 )
582: GL = MIN( GL, D( J )-TMP1-TMP2 )
583: TMP1 = TMP2
584: 40 CONTINUE
585: *
586: GU = MAX( GU, D( IEND )+TMP1 )
587: GL = MIN( GL, D( IEND )-TMP1 )
588: BNORM = MAX( ABS( GL ), ABS( GU ) )
589: GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
590: GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
591: *
592: * Compute ATOLI for the current submatrix
593: *
594: IF( ABSTOL.LE.ZERO ) THEN
595: ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
596: ELSE
597: ATOLI = ABSTOL
598: END IF
599: *
600: IF( IRANGE.GT.1 ) THEN
601: IF( GU.LT.WL ) THEN
602: NWL = NWL + IN
603: NWU = NWU + IN
604: GO TO 70
605: END IF
606: GL = MAX( GL, WL )
607: GU = MIN( GU, WU )
608: IF( GL.GE.GU )
609: $ GO TO 70
610: END IF
611: *
612: * Set Up Initial Interval
613: *
614: WORK( N+1 ) = GL
615: WORK( N+IN+1 ) = GU
616: CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
617: $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
618: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
619: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
620: *
621: NWL = NWL + IWORK( 1 )
622: NWU = NWU + IWORK( IN+1 )
623: IWOFF = M - IWORK( 1 )
624: *
625: * Compute Eigenvalues
626: *
627: ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
628: $ LOG( TWO ) ) + 2
629: CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
630: $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
631: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
632: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
633: *
634: * Copy Eigenvalues Into W and IBLOCK
635: * Use -JB for block number for unconverged eigenvalues.
636: *
637: DO 60 J = 1, IOUT
638: TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
639: *
640: * Flag non-convergence.
641: *
642: IF( J.GT.IOUT-IINFO ) THEN
643: NCNVRG = .TRUE.
644: IB = -JB
645: ELSE
646: IB = JB
647: END IF
648: DO 50 JE = IWORK( J ) + 1 + IWOFF,
649: $ IWORK( J+IN ) + IWOFF
650: W( JE ) = TMP1
651: IBLOCK( JE ) = IB
652: 50 CONTINUE
653: 60 CONTINUE
654: *
655: M = M + IM
656: END IF
657: 70 CONTINUE
658: *
659: * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
660: * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
661: *
662: IF( IRANGE.EQ.3 ) THEN
663: IM = 0
664: IDISCL = IL - 1 - NWL
665: IDISCU = NWU - IU
666: *
667: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
668: DO 80 JE = 1, M
669: IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
670: IDISCL = IDISCL - 1
671: ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
672: IDISCU = IDISCU - 1
673: ELSE
674: IM = IM + 1
675: W( IM ) = W( JE )
676: IBLOCK( IM ) = IBLOCK( JE )
677: END IF
678: 80 CONTINUE
679: M = IM
680: END IF
681: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
682: *
683: * Code to deal with effects of bad arithmetic:
684: * Some low eigenvalues to be discarded are not in (WL,WLU],
685: * or high eigenvalues to be discarded are not in (WUL,WU]
686: * so just kill off the smallest IDISCL/largest IDISCU
687: * eigenvalues, by simply finding the smallest/largest
688: * eigenvalue(s).
689: *
690: * (If N(w) is monotone non-decreasing, this should never
691: * happen.)
692: *
693: IF( IDISCL.GT.0 ) THEN
694: WKILL = WU
695: DO 100 JDISC = 1, IDISCL
696: IW = 0
697: DO 90 JE = 1, M
698: IF( IBLOCK( JE ).NE.0 .AND.
699: $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
700: IW = JE
701: WKILL = W( JE )
702: END IF
703: 90 CONTINUE
704: IBLOCK( IW ) = 0
705: 100 CONTINUE
706: END IF
707: IF( IDISCU.GT.0 ) THEN
708: *
709: WKILL = WL
710: DO 120 JDISC = 1, IDISCU
711: IW = 0
712: DO 110 JE = 1, M
713: IF( IBLOCK( JE ).NE.0 .AND.
714: $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
715: IW = JE
716: WKILL = W( JE )
717: END IF
718: 110 CONTINUE
719: IBLOCK( IW ) = 0
720: 120 CONTINUE
721: END IF
722: IM = 0
723: DO 130 JE = 1, M
724: IF( IBLOCK( JE ).NE.0 ) THEN
725: IM = IM + 1
726: W( IM ) = W( JE )
727: IBLOCK( IM ) = IBLOCK( JE )
728: END IF
729: 130 CONTINUE
730: M = IM
731: END IF
732: IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
733: TOOFEW = .TRUE.
734: END IF
735: END IF
736: *
737: * If ORDER='B', do nothing -- the eigenvalues are already sorted
738: * by block.
739: * If ORDER='E', sort the eigenvalues from smallest to largest
740: *
741: IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
742: DO 150 JE = 1, M - 1
743: IE = 0
744: TMP1 = W( JE )
745: DO 140 J = JE + 1, M
746: IF( W( J ).LT.TMP1 ) THEN
747: IE = J
748: TMP1 = W( J )
749: END IF
750: 140 CONTINUE
751: *
752: IF( IE.NE.0 ) THEN
753: ITMP1 = IBLOCK( IE )
754: W( IE ) = W( JE )
755: IBLOCK( IE ) = IBLOCK( JE )
756: W( JE ) = TMP1
757: IBLOCK( JE ) = ITMP1
758: END IF
759: 150 CONTINUE
760: END IF
761: *
762: INFO = 0
763: IF( NCNVRG )
764: $ INFO = INFO + 1
765: IF( TOOFEW )
766: $ INFO = INFO + 2
767: RETURN
768: *
769: * End of DSTEBZ
770: *
771: 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