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