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