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.3.1) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * -- April 2011 --
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: DO 10 J = 2, N
305: TMP1 = E( J-1 )**2
306: IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
307: ISPLIT( NSPLIT ) = J - 1
308: NSPLIT = NSPLIT + 1
309: WORK( J-1 ) = ZERO
310: ELSE
311: WORK( J-1 ) = TMP1
312: PIVMIN = MAX( PIVMIN, TMP1 )
313: END IF
314: 10 CONTINUE
315: ISPLIT( NSPLIT ) = N
316: PIVMIN = PIVMIN*SAFEMN
317: *
318: * Compute Interval and ATOLI
319: *
320: IF( IRANGE.EQ.3 ) THEN
321: *
322: * RANGE='I': Compute the interval containing eigenvalues
323: * IL through IU.
324: *
325: * Compute Gershgorin interval for entire (split) matrix
326: * and use it as the initial interval
327: *
328: GU = D( 1 )
329: GL = D( 1 )
330: TMP1 = ZERO
331: *
332: DO 20 J = 1, N - 1
333: TMP2 = SQRT( WORK( J ) )
334: GU = MAX( GU, D( J )+TMP1+TMP2 )
335: GL = MIN( GL, D( J )-TMP1-TMP2 )
336: TMP1 = TMP2
337: 20 CONTINUE
338: *
339: GU = MAX( GU, D( N )+TMP1 )
340: GL = MIN( GL, D( N )-TMP1 )
341: TNORM = MAX( ABS( GL ), ABS( GU ) )
342: GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
343: GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
344: *
345: * Compute Iteration parameters
346: *
347: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
348: $ LOG( TWO ) ) + 2
349: IF( ABSTOL.LE.ZERO ) THEN
350: ATOLI = ULP*TNORM
351: ELSE
352: ATOLI = ABSTOL
353: END IF
354: *
355: WORK( N+1 ) = GL
356: WORK( N+2 ) = GL
357: WORK( N+3 ) = GU
358: WORK( N+4 ) = GU
359: WORK( N+5 ) = GL
360: WORK( N+6 ) = GU
361: IWORK( 1 ) = -1
362: IWORK( 2 ) = -1
363: IWORK( 3 ) = N + 1
364: IWORK( 4 ) = N + 1
365: IWORK( 5 ) = IL - 1
366: IWORK( 6 ) = IU
367: *
368: CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
369: $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
370: $ IWORK, W, IBLOCK, IINFO )
371: *
372: IF( IWORK( 6 ).EQ.IU ) THEN
373: WL = WORK( N+1 )
374: WLU = WORK( N+3 )
375: NWL = IWORK( 1 )
376: WU = WORK( N+4 )
377: WUL = WORK( N+2 )
378: NWU = IWORK( 4 )
379: ELSE
380: WL = WORK( N+2 )
381: WLU = WORK( N+4 )
382: NWL = IWORK( 2 )
383: WU = WORK( N+3 )
384: WUL = WORK( N+1 )
385: NWU = IWORK( 3 )
386: END IF
387: *
388: IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
389: INFO = 4
390: RETURN
391: END IF
392: ELSE
393: *
394: * RANGE='A' or 'V' -- Set ATOLI
395: *
396: TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
397: $ ABS( D( N ) )+ABS( E( N-1 ) ) )
398: *
399: DO 30 J = 2, N - 1
400: TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
401: $ ABS( E( J ) ) )
402: 30 CONTINUE
403: *
404: IF( ABSTOL.LE.ZERO ) THEN
405: ATOLI = ULP*TNORM
406: ELSE
407: ATOLI = ABSTOL
408: END IF
409: *
410: IF( IRANGE.EQ.2 ) THEN
411: WL = VL
412: WU = VU
413: ELSE
414: WL = ZERO
415: WU = ZERO
416: END IF
417: END IF
418: *
419: * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
420: * NWL accumulates the number of eigenvalues .le. WL,
421: * NWU accumulates the number of eigenvalues .le. WU
422: *
423: M = 0
424: IEND = 0
425: INFO = 0
426: NWL = 0
427: NWU = 0
428: *
429: DO 70 JB = 1, NSPLIT
430: IOFF = IEND
431: IBEGIN = IOFF + 1
432: IEND = ISPLIT( JB )
433: IN = IEND - IOFF
434: *
435: IF( IN.EQ.1 ) THEN
436: *
437: * Special Case -- IN=1
438: *
439: IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
440: $ NWL = NWL + 1
441: IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
442: $ NWU = NWU + 1
443: IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
444: $ D( IBEGIN )-PIVMIN ) ) THEN
445: M = M + 1
446: W( M ) = D( IBEGIN )
447: IBLOCK( M ) = JB
448: END IF
449: ELSE
450: *
451: * General Case -- IN > 1
452: *
453: * Compute Gershgorin Interval
454: * and use it as the initial interval
455: *
456: GU = D( IBEGIN )
457: GL = D( IBEGIN )
458: TMP1 = ZERO
459: *
460: DO 40 J = IBEGIN, IEND - 1
461: TMP2 = ABS( E( J ) )
462: GU = MAX( GU, D( J )+TMP1+TMP2 )
463: GL = MIN( GL, D( J )-TMP1-TMP2 )
464: TMP1 = TMP2
465: 40 CONTINUE
466: *
467: GU = MAX( GU, D( IEND )+TMP1 )
468: GL = MIN( GL, D( IEND )-TMP1 )
469: BNORM = MAX( ABS( GL ), ABS( GU ) )
470: GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
471: GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
472: *
473: * Compute ATOLI for the current submatrix
474: *
475: IF( ABSTOL.LE.ZERO ) THEN
476: ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
477: ELSE
478: ATOLI = ABSTOL
479: END IF
480: *
481: IF( IRANGE.GT.1 ) THEN
482: IF( GU.LT.WL ) THEN
483: NWL = NWL + IN
484: NWU = NWU + IN
485: GO TO 70
486: END IF
487: GL = MAX( GL, WL )
488: GU = MIN( GU, WU )
489: IF( GL.GE.GU )
490: $ GO TO 70
491: END IF
492: *
493: * Set Up Initial Interval
494: *
495: WORK( N+1 ) = GL
496: WORK( N+IN+1 ) = GU
497: CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
498: $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
499: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
500: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
501: *
502: NWL = NWL + IWORK( 1 )
503: NWU = NWU + IWORK( IN+1 )
504: IWOFF = M - IWORK( 1 )
505: *
506: * Compute Eigenvalues
507: *
508: ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
509: $ LOG( TWO ) ) + 2
510: CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
511: $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
512: $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
513: $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
514: *
515: * Copy Eigenvalues Into W and IBLOCK
516: * Use -JB for block number for unconverged eigenvalues.
517: *
518: DO 60 J = 1, IOUT
519: TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
520: *
521: * Flag non-convergence.
522: *
523: IF( J.GT.IOUT-IINFO ) THEN
524: NCNVRG = .TRUE.
525: IB = -JB
526: ELSE
527: IB = JB
528: END IF
529: DO 50 JE = IWORK( J ) + 1 + IWOFF,
530: $ IWORK( J+IN ) + IWOFF
531: W( JE ) = TMP1
532: IBLOCK( JE ) = IB
533: 50 CONTINUE
534: 60 CONTINUE
535: *
536: M = M + IM
537: END IF
538: 70 CONTINUE
539: *
540: * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
541: * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
542: *
543: IF( IRANGE.EQ.3 ) THEN
544: IM = 0
545: IDISCL = IL - 1 - NWL
546: IDISCU = NWU - IU
547: *
548: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
549: DO 80 JE = 1, M
550: IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
551: IDISCL = IDISCL - 1
552: ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
553: IDISCU = IDISCU - 1
554: ELSE
555: IM = IM + 1
556: W( IM ) = W( JE )
557: IBLOCK( IM ) = IBLOCK( JE )
558: END IF
559: 80 CONTINUE
560: M = IM
561: END IF
562: IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
563: *
564: * Code to deal with effects of bad arithmetic:
565: * Some low eigenvalues to be discarded are not in (WL,WLU],
566: * or high eigenvalues to be discarded are not in (WUL,WU]
567: * so just kill off the smallest IDISCL/largest IDISCU
568: * eigenvalues, by simply finding the smallest/largest
569: * eigenvalue(s).
570: *
571: * (If N(w) is monotone non-decreasing, this should never
572: * happen.)
573: *
574: IF( IDISCL.GT.0 ) THEN
575: WKILL = WU
576: DO 100 JDISC = 1, IDISCL
577: IW = 0
578: DO 90 JE = 1, M
579: IF( IBLOCK( JE ).NE.0 .AND.
580: $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
581: IW = JE
582: WKILL = W( JE )
583: END IF
584: 90 CONTINUE
585: IBLOCK( IW ) = 0
586: 100 CONTINUE
587: END IF
588: IF( IDISCU.GT.0 ) THEN
589: *
590: WKILL = WL
591: DO 120 JDISC = 1, IDISCU
592: IW = 0
593: DO 110 JE = 1, M
594: IF( IBLOCK( JE ).NE.0 .AND.
595: $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
596: IW = JE
597: WKILL = W( JE )
598: END IF
599: 110 CONTINUE
600: IBLOCK( IW ) = 0
601: 120 CONTINUE
602: END IF
603: IM = 0
604: DO 130 JE = 1, M
605: IF( IBLOCK( JE ).NE.0 ) THEN
606: IM = IM + 1
607: W( IM ) = W( JE )
608: IBLOCK( IM ) = IBLOCK( JE )
609: END IF
610: 130 CONTINUE
611: M = IM
612: END IF
613: IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
614: TOOFEW = .TRUE.
615: END IF
616: END IF
617: *
618: * If ORDER='B', do nothing -- the eigenvalues are already sorted
619: * by block.
620: * If ORDER='E', sort the eigenvalues from smallest to largest
621: *
622: IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
623: DO 150 JE = 1, M - 1
624: IE = 0
625: TMP1 = W( JE )
626: DO 140 J = JE + 1, M
627: IF( W( J ).LT.TMP1 ) THEN
628: IE = J
629: TMP1 = W( J )
630: END IF
631: 140 CONTINUE
632: *
633: IF( IE.NE.0 ) THEN
634: ITMP1 = IBLOCK( IE )
635: W( IE ) = W( JE )
636: IBLOCK( IE ) = IBLOCK( JE )
637: W( JE ) = TMP1
638: IBLOCK( JE ) = ITMP1
639: END IF
640: 150 CONTINUE
641: END IF
642: *
643: INFO = 0
644: IF( NCNVRG )
645: $ INFO = INFO + 1
646: IF( TOOFEW )
647: $ INFO = INFO + 2
648: RETURN
649: *
650: * End of DSTEBZ
651: *
652: END
CVSweb interface <joel.bertrand@systella.fr>