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