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