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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *     8-18-00:  Increase FUDGE factor for T3E (eca)
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          ORDER, RANGE
013:       INTEGER            IL, INFO, IU, M, N, NSPLIT
014:       REAL               ABSTOL, VL, VU
015: *     ..
016: *     .. Array Arguments ..
017:       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
018:       REAL               D( * ), E( * ), W( * ), WORK( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  SSTEBZ computes the eigenvalues of a symmetric tridiagonal
025: *  matrix T.  The user may ask for all eigenvalues, all eigenvalues
026: *  in the half-open interval (VL, VU], or the IL-th through IU-th
027: *  eigenvalues.
028: *
029: *  To avoid overflow, the matrix must be scaled so that its
030: *  largest element is no greater than overflow**(1/2) *
031: *  underflow**(1/4) in absolute value, and for greatest
032: *  accuracy, it should not be much smaller than that.
033: *
034: *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
035: *  Matrix", Report CS41, Computer Science Dept., Stanford
036: *  University, July 21, 1966.
037: *
038: *  Arguments
039: *  =========
040: *
041: *  RANGE   (input) CHARACTER*1
042: *          = 'A': ("All")   all eigenvalues will be found.
043: *          = 'V': ("Value") all eigenvalues in the half-open interval
044: *                           (VL, VU] will be found.
045: *          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
046: *                           entire matrix) will be found.
047: *
048: *  ORDER   (input) CHARACTER*1
049: *          = 'B': ("By Block") the eigenvalues will be grouped by
050: *                              split-off block (see IBLOCK, ISPLIT) and
051: *                              ordered from smallest to largest within
052: *                              the block.
053: *          = 'E': ("Entire matrix")
054: *                              the eigenvalues for the entire matrix
055: *                              will be ordered from smallest to
056: *                              largest.
057: *
058: *  N       (input) INTEGER
059: *          The order of the tridiagonal matrix T.  N >= 0.
060: *
061: *  VL      (input) REAL
062: *  VU      (input) REAL
063: *          If RANGE='V', the lower and upper bounds of the interval to
064: *          be searched for eigenvalues.  Eigenvalues less than or equal
065: *          to VL, or greater than VU, will not be returned.  VL < VU.
066: *          Not referenced if RANGE = 'A' or 'I'.
067: *
068: *  IL      (input) INTEGER
069: *  IU      (input) INTEGER
070: *          If RANGE='I', the indices (in ascending order) of the
071: *          smallest and largest eigenvalues to be returned.
072: *          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
073: *          Not referenced if RANGE = 'A' or 'V'.
074: *
075: *  ABSTOL  (input) REAL
076: *          The absolute tolerance for the eigenvalues.  An eigenvalue
077: *          (or cluster) is considered to be located if it has been
078: *          determined to lie in an interval whose width is ABSTOL or
079: *          less.  If ABSTOL is less than or equal to zero, then ULP*|T|
080: *          will be used, where |T| means the 1-norm of T.
081: *
082: *          Eigenvalues will be computed most accurately when ABSTOL is
083: *          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
084: *
085: *  D       (input) REAL array, dimension (N)
086: *          The n diagonal elements of the tridiagonal matrix T.
087: *
088: *  E       (input) REAL array, dimension (N-1)
089: *          The (n-1) off-diagonal elements of the tridiagonal matrix T.
090: *
091: *  M       (output) INTEGER
092: *          The actual number of eigenvalues found. 0 <= M <= N.
093: *          (See also the description of INFO=2,3.)
094: *
095: *  NSPLIT  (output) INTEGER
096: *          The number of diagonal blocks in the matrix T.
097: *          1 <= NSPLIT <= N.
098: *
099: *  W       (output) REAL array, dimension (N)
100: *          On exit, the first M elements of W will contain the
101: *          eigenvalues.  (SSTEBZ 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.  (SSTEBZ 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) REAL 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  REAL, 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   REAL, 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:       REAL               ZERO, ONE, TWO, HALF
176:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
177:      $                   HALF = 1.0E0 / TWO )
178:       REAL               FUDGE, RELFAC
179:       PARAMETER          ( FUDGE = 2.1E0, RELFAC = 2.0E0 )
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:       REAL               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:       REAL               SLAMCH
197:       EXTERNAL           LSAME, ILAENV, SLAMCH
198: *     ..
199: *     .. External Subroutines ..
200:       EXTERNAL           SLAEBZ, 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 ) INFO = -5
241:       ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
242:      $          THEN
243:          INFO = -6
244:       ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
245:      $          THEN
246:          INFO = -7
247:       END IF
248: *
249:       IF( INFO.NE.0 ) THEN
250:          CALL XERBLA( 'SSTEBZ', -INFO )
251:          RETURN
252:       END IF
253: *
254: *     Initialize error flags
255: *
256:       INFO = 0
257:       NCNVRG = .FALSE.
258:       TOOFEW = .FALSE.
259: *
260: *     Quick return if possible
261: *
262:       M = 0
263:       IF( N.EQ.0 )
264:      $   RETURN
265: *
266: *     Simplifications:
267: *
268:       IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
269:      $   IRANGE = 1
270: *
271: *     Get machine constants
272: *     NB is the minimum vector length for vector bisection, or 0
273: *     if only scalar is to be done.
274: *
275:       SAFEMN = SLAMCH( 'S' )
276:       ULP = SLAMCH( 'P' )
277:       RTOLI = ULP*RELFAC
278:       NB = ILAENV( 1, 'SSTEBZ', ' ', N, -1, -1, -1 )
279:       IF( NB.LE.1 )
280:      $   NB = 0
281: *
282: *     Special Case when N=1
283: *
284:       IF( N.EQ.1 ) THEN
285:          NSPLIT = 1
286:          ISPLIT( 1 ) = 1
287:          IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
288:             M = 0
289:          ELSE
290:             W( 1 ) = D( 1 )
291:             IBLOCK( 1 ) = 1
292:             M = 1
293:          END IF
294:          RETURN
295:       END IF
296: *
297: *     Compute Splitting Points
298: *
299:       NSPLIT = 1
300:       WORK( N ) = ZERO
301:       PIVMIN = ONE
302: *
303: CDIR$ NOVECTOR
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 SLAEBZ( 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 SLAEBZ( 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 SLAEBZ( 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 SSTEBZ
651: *
652:       END
653: