001:       SUBROUTINE SLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
002:      $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
003:      $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
004:      $                    WORK, IWORK, INFO )
005:       IMPLICIT NONE
006: *
007: *  -- LAPACK auxiliary routine (version 3.2) --
008: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
009: *     November 2006
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          RANGE
013:       INTEGER            IL, INFO, IU, M, N, NSPLIT
014:       REAL              PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
015: *     ..
016: *     .. Array Arguments ..
017:       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
018:      $                   INDEXW( * )
019:       REAL               D( * ), E( * ), E2( * ), GERS( * ),
020:      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  To find the desired eigenvalues of a given real symmetric
027: *  tridiagonal matrix T, SLARRE sets any "small" off-diagonal
028: *  elements to zero, and for each unreduced block T_i, it finds
029: *  (a) a suitable shift at one end of the block's spectrum,
030: *  (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
031: *  (c) eigenvalues of each L_i D_i L_i^T.
032: *  The representations and eigenvalues found are then used by
033: *  SSTEMR to compute the eigenvectors of T.
034: *  The accuracy varies depending on whether bisection is used to
035: *  find a few eigenvalues or the dqds algorithm (subroutine SLASQ2) to
036: *  conpute all and then discard any unwanted one.
037: *  As an added benefit, SLARRE also outputs the n
038: *  Gerschgorin intervals for the matrices L_i D_i L_i^T.
039: *
040: *  Arguments
041: *  =========
042: *
043: *  RANGE   (input) CHARACTER
044: *          = 'A': ("All")   all eigenvalues will be found.
045: *          = 'V': ("Value") all eigenvalues in the half-open interval
046: *                           (VL, VU] will be found.
047: *          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
048: *                           entire matrix) will be found.
049: *
050: *  N       (input) INTEGER
051: *          The order of the matrix. N > 0.
052: *
053: *  VL      (input/output) REAL            
054: *  VU      (input/output) REAL            
055: *          If RANGE='V', the lower and upper bounds for the eigenvalues.
056: *          Eigenvalues less than or equal to VL, or greater than VU,
057: *          will not be returned.  VL < VU.
058: *          If RANGE='I' or ='A', SLARRE computes bounds on the desired
059: *          part of the spectrum.
060: *
061: *  IL      (input) INTEGER
062: *  IU      (input) INTEGER
063: *          If RANGE='I', the indices (in ascending order) of the
064: *          smallest and largest eigenvalues to be returned.
065: *          1 <= IL <= IU <= N.
066: *
067: *  D       (input/output) REAL             array, dimension (N)
068: *          On entry, the N diagonal elements of the tridiagonal
069: *          matrix T.
070: *          On exit, the N diagonal elements of the diagonal
071: *          matrices D_i.
072: *
073: *  E       (input/output) REAL             array, dimension (N)
074: *          On entry, the first (N-1) entries contain the subdiagonal
075: *          elements of the tridiagonal matrix T; E(N) need not be set.
076: *          On exit, E contains the subdiagonal elements of the unit
077: *          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
078: *          1 <= I <= NSPLIT, contain the base points sigma_i on output.
079: *
080: *  E2      (input/output) REAL             array, dimension (N)
081: *          On entry, the first (N-1) entries contain the SQUARES of the
082: *          subdiagonal elements of the tridiagonal matrix T;
083: *          E2(N) need not be set.
084: *          On exit, the entries E2( ISPLIT( I ) ),
085: *          1 <= I <= NSPLIT, have been set to zero
086: *
087: *  RTOL1   (input) REAL            
088: *  RTOL2   (input) REAL            
089: *           Parameters for bisection.
090: *           An interval [LEFT,RIGHT] has converged if
091: *           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
092: *
093: *  SPLTOL (input) REAL            
094: *          The threshold for splitting.
095: *
096: *  NSPLIT  (output) INTEGER
097: *          The number of blocks T splits into. 1 <= NSPLIT <= N.
098: *
099: *  ISPLIT  (output) INTEGER array, dimension (N)
100: *          The splitting points, at which T breaks up into blocks.
101: *          The first block consists of rows/columns 1 to ISPLIT(1),
102: *          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
103: *          etc., and the NSPLIT-th consists of rows/columns
104: *          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
105: *
106: *  M       (output) INTEGER
107: *          The total number of eigenvalues (of all L_i D_i L_i^T)
108: *          found.
109: *
110: *  W       (output) REAL             array, dimension (N)
111: *          The first M elements contain the eigenvalues. The
112: *          eigenvalues of each of the blocks, L_i D_i L_i^T, are
113: *          sorted in ascending order ( SLARRE may use the
114: *          remaining N-M elements as workspace).
115: *
116: *  WERR    (output) REAL             array, dimension (N)
117: *          The error bound on the corresponding eigenvalue in W.
118: *
119: *  WGAP    (output) REAL             array, dimension (N)
120: *          The separation from the right neighbor eigenvalue in W.
121: *          The gap is only with respect to the eigenvalues of the same block
122: *          as each block has its own representation tree.
123: *          Exception: at the right end of a block we store the left gap
124: *
125: *  IBLOCK  (output) INTEGER array, dimension (N)
126: *          The indices of the blocks (submatrices) associated with the
127: *          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
128: *          W(i) belongs to the first block from the top, =2 if W(i)
129: *          belongs to the second block, etc.
130: *
131: *  INDEXW  (output) INTEGER array, dimension (N)
132: *          The indices of the eigenvalues within each block (submatrix);
133: *          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
134: *          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
135: *
136: *  GERS    (output) REAL             array, dimension (2*N)
137: *          The N Gerschgorin intervals (the i-th Gerschgorin interval
138: *          is (GERS(2*i-1), GERS(2*i)).
139: *
140: *  PIVMIN  (output) DOUBLE PRECISION
141: *          The minimum pivot in the Sturm sequence for T.
142: *
143: *  WORK    (workspace) REAL             array, dimension (6*N)
144: *          Workspace.
145: *
146: *  IWORK   (workspace) INTEGER array, dimension (5*N)
147: *          Workspace.
148: *
149: *  INFO    (output) INTEGER
150: *          = 0:  successful exit
151: *          > 0:  A problem occured in SLARRE.
152: *          < 0:  One of the called subroutines signaled an internal problem.
153: *                Needs inspection of the corresponding parameter IINFO
154: *                for further information.
155: *
156: *          =-1:  Problem in SLARRD.
157: *          = 2:  No base representation could be found in MAXTRY iterations.
158: *                Increasing MAXTRY and recompilation might be a remedy.
159: *          =-3:  Problem in SLARRB when computing the refined root
160: *                representation for SLASQ2.
161: *          =-4:  Problem in SLARRB when preforming bisection on the
162: *                desired part of the spectrum.
163: *          =-5:  Problem in SLASQ2.
164: *          =-6:  Problem in SLASQ2.
165: *
166: *  Further Details
167: *  The base representations are required to suffer very little
168: *  element growth and consequently define all their eigenvalues to
169: *  high relative accuracy.
170: *  ===============
171: *
172: *  Based on contributions by
173: *     Beresford Parlett, University of California, Berkeley, USA
174: *     Jim Demmel, University of California, Berkeley, USA
175: *     Inderjit Dhillon, University of Texas, Austin, USA
176: *     Osni Marques, LBNL/NERSC, USA
177: *     Christof Voemel, University of California, Berkeley, USA
178: *
179: *  =====================================================================
180: *
181: *     .. Parameters ..
182:       REAL               FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
183:      $                   MAXGROWTH, ONE, PERT, TWO, ZERO
184:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
185:      $                     TWO = 2.0E0, FOUR=4.0E0,
186:      $                     HNDRD = 100.0E0,
187:      $                     PERT = 4.0E0,
188:      $                     HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
189:      $                     MAXGROWTH = 64.0E0, FUDGE = 2.0E0 )
190:       INTEGER            MAXTRY, ALLRNG, INDRNG, VALRNG
191:       PARAMETER          ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
192:      $                     VALRNG = 3 )
193: *     ..
194: *     .. Local Scalars ..
195:       LOGICAL            FORCEB, NOREP, USEDQD
196:       INTEGER            CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
197:      $                   IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
198:      $                   WBEGIN, WEND
199:       REAL               AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
200:      $                   EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
201:      $                   RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
202:      $                   TAU, TMP, TMP1
203: 
204: 
205: *     ..
206: *     .. Local Arrays ..
207:       INTEGER            ISEED( 4 )
208: *     ..
209: *     .. External Functions ..
210:       LOGICAL            LSAME
211:       REAL                        SLAMCH
212:       EXTERNAL           SLAMCH, LSAME
213: 
214: *     ..
215: *     .. External Subroutines ..
216:       EXTERNAL           SCOPY, SLARNV, SLARRA, SLARRB, SLARRC, SLARRD,
217:      $                   SLASQ2
218: *     ..
219: *     .. Intrinsic Functions ..
220:       INTRINSIC          ABS, MAX, MIN
221: 
222: *     ..
223: *     .. Executable Statements ..
224: *
225: 
226:       INFO = 0
227: 
228: *
229: *     Decode RANGE
230: *
231:       IF( LSAME( RANGE, 'A' ) ) THEN
232:          IRANGE = ALLRNG
233:       ELSE IF( LSAME( RANGE, 'V' ) ) THEN
234:          IRANGE = VALRNG
235:       ELSE IF( LSAME( RANGE, 'I' ) ) THEN
236:          IRANGE = INDRNG
237:       END IF
238: 
239:       M = 0
240: 
241: *     Get machine constants
242:       SAFMIN = SLAMCH( 'S' )
243:       EPS = SLAMCH( 'P' )
244: 
245: *     Set parameters
246:       RTL = HNDRD*EPS
247: *     If one were ever to ask for less initial precision in BSRTOL,
248: *     one should keep in mind that for the subset case, the extremal
249: *     eigenvalues must be at least as accurate as the current setting
250: *     (eigenvalues in the middle need not as much accuracy)
251:       BSRTOL = SQRT(EPS)*(0.5E-3)
252: 
253: *     Treat case of 1x1 matrix for quick return
254:       IF( N.EQ.1 ) THEN
255:          IF( (IRANGE.EQ.ALLRNG).OR.
256:      $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
257:      $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
258:             M = 1
259:             W(1) = D(1)
260: *           The computation error of the eigenvalue is zero
261:             WERR(1) = ZERO
262:             WGAP(1) = ZERO
263:             IBLOCK( 1 ) = 1
264:             INDEXW( 1 ) = 1
265:             GERS(1) = D( 1 )
266:             GERS(2) = D( 1 )
267:          ENDIF
268: *        store the shift for the initial RRR, which is zero in this case
269:          E(1) = ZERO
270:          RETURN
271:       END IF
272: 
273: *     General case: tridiagonal matrix of order > 1
274: *
275: *     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
276: *     Compute maximum off-diagonal entry and pivmin.
277:       GL = D(1)
278:       GU = D(1)
279:       EOLD = ZERO
280:       EMAX = ZERO
281:       E(N) = ZERO
282:       DO 5 I = 1,N
283:          WERR(I) = ZERO
284:          WGAP(I) = ZERO
285:          EABS = ABS( E(I) )
286:          IF( EABS .GE. EMAX ) THEN
287:             EMAX = EABS
288:          END IF
289:          TMP1 = EABS + EOLD
290:          GERS( 2*I-1) = D(I) - TMP1
291:          GL =  MIN( GL, GERS( 2*I - 1))
292:          GERS( 2*I ) = D(I) + TMP1
293:          GU = MAX( GU, GERS(2*I) )
294:          EOLD  = EABS
295:  5    CONTINUE
296: *     The minimum pivot allowed in the Sturm sequence for T
297:       PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
298: *     Compute spectral diameter. The Gerschgorin bounds give an
299: *     estimate that is wrong by at most a factor of SQRT(2)
300:       SPDIAM = GU - GL
301: 
302: *     Compute splitting points
303:       CALL SLARRA( N, D, E, E2, SPLTOL, SPDIAM,
304:      $                    NSPLIT, ISPLIT, IINFO )
305: 
306: *     Can force use of bisection instead of faster DQDS.
307: *     Option left in the code for future multisection work.
308:       FORCEB = .FALSE.
309: 
310: *     Initialize USEDQD, DQDS should be used for ALLRNG unless someone
311: *     explicitly wants bisection.
312:       USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
313: 
314:       IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
315: *        Set interval [VL,VU] that contains all eigenvalues
316:          VL = GL
317:          VU = GU
318:       ELSE
319: *        We call SLARRD to find crude approximations to the eigenvalues
320: *        in the desired range. In case IRANGE = INDRNG, we also obtain the
321: *        interval (VL,VU] that contains all the wanted eigenvalues.
322: *        An interval [LEFT,RIGHT] has converged if
323: *        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
324: *        SLARRD needs a WORK of size 4*N, IWORK of size 3*N
325:          CALL SLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
326:      $                    BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
327:      $                    MM, W, WERR, VL, VU, IBLOCK, INDEXW,
328:      $                    WORK, IWORK, IINFO )
329:          IF( IINFO.NE.0 ) THEN
330:             INFO = -1
331:             RETURN
332:          ENDIF
333: *        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
334:          DO 14 I = MM+1,N
335:             W( I ) = ZERO
336:             WERR( I ) = ZERO
337:             IBLOCK( I ) = 0
338:             INDEXW( I ) = 0
339:  14      CONTINUE
340:       END IF
341: 
342: 
343: ***
344: *     Loop over unreduced blocks
345:       IBEGIN = 1
346:       WBEGIN = 1
347:       DO 170 JBLK = 1, NSPLIT
348:          IEND = ISPLIT( JBLK )
349:          IN = IEND - IBEGIN + 1
350: 
351: *        1 X 1 block
352:          IF( IN.EQ.1 ) THEN
353:             IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
354:      $         ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
355:      $        .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
356:      $        ) THEN
357:                M = M + 1
358:                W( M ) = D( IBEGIN )
359:                WERR(M) = ZERO
360: *              The gap for a single block doesn't matter for the later
361: *              algorithm and is assigned an arbitrary large value
362:                WGAP(M) = ZERO
363:                IBLOCK( M ) = JBLK
364:                INDEXW( M ) = 1
365:                WBEGIN = WBEGIN + 1
366:             ENDIF
367: *           E( IEND ) holds the shift for the initial RRR
368:             E( IEND ) = ZERO
369:             IBEGIN = IEND + 1
370:             GO TO 170
371:          END IF
372: *
373: *        Blocks of size larger than 1x1
374: *
375: *        E( IEND ) will hold the shift for the initial RRR, for now set it =0
376:          E( IEND ) = ZERO
377: *
378: *        Find local outer bounds GL,GU for the block
379:          GL = D(IBEGIN)
380:          GU = D(IBEGIN)
381:          DO 15 I = IBEGIN , IEND
382:             GL = MIN( GERS( 2*I-1 ), GL )
383:             GU = MAX( GERS( 2*I ), GU )
384:  15      CONTINUE
385:          SPDIAM = GU - GL
386: 
387:          IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
388: *           Count the number of eigenvalues in the current block.
389:             MB = 0
390:             DO 20 I = WBEGIN,MM
391:                IF( IBLOCK(I).EQ.JBLK ) THEN
392:                   MB = MB+1
393:                ELSE
394:                   GOTO 21
395:                ENDIF
396:  20         CONTINUE
397:  21         CONTINUE
398: 
399:             IF( MB.EQ.0) THEN
400: *              No eigenvalue in the current block lies in the desired range
401: *              E( IEND ) holds the shift for the initial RRR
402:                E( IEND ) = ZERO
403:                IBEGIN = IEND + 1
404:                GO TO 170
405:             ELSE
406: 
407: *              Decide whether dqds or bisection is more efficient
408:                USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
409:                WEND = WBEGIN + MB - 1
410: *              Calculate gaps for the current block
411: *              In later stages, when representations for individual
412: *              eigenvalues are different, we use SIGMA = E( IEND ).
413:                SIGMA = ZERO
414:                DO 30 I = WBEGIN, WEND - 1
415:                   WGAP( I ) = MAX( ZERO,
416:      $                        W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
417:  30            CONTINUE
418:                WGAP( WEND ) = MAX( ZERO,
419:      $                     VU - SIGMA - (W( WEND )+WERR( WEND )))
420: *              Find local index of the first and last desired evalue.
421:                INDL = INDEXW(WBEGIN)
422:                INDU = INDEXW( WEND )
423:             ENDIF
424:          ENDIF
425:          IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
426: *           Case of DQDS
427: *           Find approximations to the extremal eigenvalues of the block
428:             CALL SLARRK( IN, 1, GL, GU, D(IBEGIN),
429:      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
430:             IF( IINFO.NE.0 ) THEN
431:                INFO = -1
432:                RETURN
433:             ENDIF
434:             ISLEFT = MAX(GL, TMP - TMP1
435:      $               - HNDRD * EPS* ABS(TMP - TMP1))
436: 
437:             CALL SLARRK( IN, IN, GL, GU, D(IBEGIN),
438:      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
439:             IF( IINFO.NE.0 ) THEN
440:                INFO = -1
441:                RETURN
442:             ENDIF
443:             ISRGHT = MIN(GU, TMP + TMP1
444:      $                 + HNDRD * EPS * ABS(TMP + TMP1))
445: *           Improve the estimate of the spectral diameter
446:             SPDIAM = ISRGHT - ISLEFT
447:          ELSE
448: *           Case of bisection
449: *           Find approximations to the wanted extremal eigenvalues
450:             ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
451:      $                  - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
452:             ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
453:      $                  + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
454:          ENDIF
455: 
456: 
457: *        Decide whether the base representation for the current block
458: *        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
459: *        should be on the left or the right end of the current block.
460: *        The strategy is to shift to the end which is "more populated"
461: *        Furthermore, decide whether to use DQDS for the computation of
462: *        the eigenvalue approximations at the end of SLARRE or bisection.
463: *        dqds is chosen if all eigenvalues are desired or the number of
464: *        eigenvalues to be computed is large compared to the blocksize.
465:          IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
466: *           If all the eigenvalues have to be computed, we use dqd
467:             USEDQD = .TRUE.
468: *           INDL is the local index of the first eigenvalue to compute
469:             INDL = 1
470:             INDU = IN
471: *           MB =  number of eigenvalues to compute
472:             MB = IN
473:             WEND = WBEGIN + MB - 1
474: *           Define 1/4 and 3/4 points of the spectrum
475:             S1 = ISLEFT + FOURTH * SPDIAM
476:             S2 = ISRGHT - FOURTH * SPDIAM
477:          ELSE
478: *           SLARRD has computed IBLOCK and INDEXW for each eigenvalue
479: *           approximation.
480: *           choose sigma
481:             IF( USEDQD ) THEN
482:                S1 = ISLEFT + FOURTH * SPDIAM
483:                S2 = ISRGHT - FOURTH * SPDIAM
484:             ELSE
485:                TMP = MIN(ISRGHT,VU) -  MAX(ISLEFT,VL)
486:                S1 =  MAX(ISLEFT,VL) + FOURTH * TMP
487:                S2 =  MIN(ISRGHT,VU) - FOURTH * TMP
488:             ENDIF
489:          ENDIF
490: 
491: *        Compute the negcount at the 1/4 and 3/4 points
492:          IF(MB.GT.1) THEN
493:             CALL SLARRC( 'T', IN, S1, S2, D(IBEGIN),
494:      $                    E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
495:          ENDIF
496: 
497:          IF(MB.EQ.1) THEN
498:             SIGMA = GL
499:             SGNDEF = ONE
500:          ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
501:             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
502:                SIGMA = MAX(ISLEFT,GL)
503:             ELSEIF( USEDQD ) THEN
504: *              use Gerschgorin bound as shift to get pos def matrix
505: *              for dqds
506:                SIGMA = ISLEFT
507:             ELSE
508: *              use approximation of the first desired eigenvalue of the
509: *              block as shift
510:                SIGMA = MAX(ISLEFT,VL)
511:             ENDIF
512:             SGNDEF = ONE
513:          ELSE
514:             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
515:                SIGMA = MIN(ISRGHT,GU)
516:             ELSEIF( USEDQD ) THEN
517: *              use Gerschgorin bound as shift to get neg def matrix
518: *              for dqds
519:                SIGMA = ISRGHT
520:             ELSE
521: *              use approximation of the first desired eigenvalue of the
522: *              block as shift
523:                SIGMA = MIN(ISRGHT,VU)
524:             ENDIF
525:             SGNDEF = -ONE
526:          ENDIF
527: 
528: 
529: *        An initial SIGMA has been chosen that will be used for computing
530: *        T - SIGMA I = L D L^T
531: *        Define the increment TAU of the shift in case the initial shift
532: *        needs to be refined to obtain a factorization with not too much
533: *        element growth.
534:          IF( USEDQD ) THEN
535: *           The initial SIGMA was to the outer end of the spectrum
536: *           the matrix is definite and we need not retreat.
537:             TAU = SPDIAM*EPS*N + TWO*PIVMIN
538:          ELSE
539:             IF(MB.GT.1) THEN
540:                CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
541:                AVGAP = ABS(CLWDTH / REAL(WEND-WBEGIN))
542:                IF( SGNDEF.EQ.ONE ) THEN
543:                   TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
544:                   TAU = MAX(TAU,WERR(WBEGIN))
545:                ELSE
546:                   TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
547:                   TAU = MAX(TAU,WERR(WEND))
548:                ENDIF
549:             ELSE
550:                TAU = WERR(WBEGIN)
551:             ENDIF
552:          ENDIF
553: *
554:          DO 80 IDUM = 1, MAXTRY
555: *           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
556: *           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
557: *           pivots in WORK(2*IN+1:3*IN)
558:             DPIVOT = D( IBEGIN ) - SIGMA
559:             WORK( 1 ) = DPIVOT
560:             DMAX = ABS( WORK(1) )
561:             J = IBEGIN
562:             DO 70 I = 1, IN - 1
563:                WORK( 2*IN+I ) = ONE / WORK( I )
564:                TMP = E( J )*WORK( 2*IN+I )
565:                WORK( IN+I ) = TMP
566:                DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
567:                WORK( I+1 ) = DPIVOT
568:                DMAX = MAX( DMAX, ABS(DPIVOT) )
569:                J = J + 1
570:  70         CONTINUE
571: *           check for element growth
572:             IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
573:                NOREP = .TRUE.
574:             ELSE
575:                NOREP = .FALSE.
576:             ENDIF
577:             IF( USEDQD .AND. .NOT.NOREP ) THEN
578: *              Ensure the definiteness of the representation
579: *              All entries of D (of L D L^T) must have the same sign
580:                DO 71 I = 1, IN
581:                   TMP = SGNDEF*WORK( I )
582:                   IF( TMP.LT.ZERO ) NOREP = .TRUE.
583:  71            CONTINUE
584:             ENDIF
585:             IF(NOREP) THEN
586: *              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
587: *              shift which makes the matrix definite. So we should end up
588: *              here really only in the case of IRANGE = VALRNG or INDRNG.
589:                IF( IDUM.EQ.MAXTRY-1 ) THEN
590:                   IF( SGNDEF.EQ.ONE ) THEN
591: *                    The fudged Gerschgorin shift should succeed
592:                      SIGMA =
593:      $                    GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
594:                   ELSE
595:                      SIGMA =
596:      $                    GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
597:                   END IF
598:                ELSE
599:                   SIGMA = SIGMA - SGNDEF * TAU
600:                   TAU = TWO * TAU
601:                END IF
602:             ELSE
603: *              an initial RRR is found
604:                GO TO 83
605:             END IF
606:  80      CONTINUE
607: *        if the program reaches this point, no base representation could be
608: *        found in MAXTRY iterations.
609:          INFO = 2
610:          RETURN
611: 
612:  83      CONTINUE
613: *        At this point, we have found an initial base representation
614: *        T - SIGMA I = L D L^T with not too much element growth.
615: *        Store the shift.
616:          E( IEND ) = SIGMA
617: *        Store D and L.
618:          CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
619:          CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
620: 
621: 
622:          IF(MB.GT.1 ) THEN
623: *
624: *           Perturb each entry of the base representation by a small
625: *           (but random) relative amount to overcome difficulties with
626: *           glued matrices.
627: *
628:             DO 122 I = 1, 4
629:                ISEED( I ) = 1
630:  122        CONTINUE
631: 
632:             CALL SLARNV(2, ISEED, 2*IN-1, WORK(1))
633:             DO 125 I = 1,IN-1
634:                D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
635:                E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
636:  125        CONTINUE
637:             D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
638: *
639:          ENDIF
640: *
641: *        Don't update the Gerschgorin intervals because keeping track
642: *        of the updates would be too much work in SLARRV.
643: *        We update W instead and use it to locate the proper Gerschgorin
644: *        intervals.
645: 
646: *        Compute the required eigenvalues of L D L' by bisection or dqds
647:          IF ( .NOT.USEDQD ) THEN
648: *           If SLARRD has been used, shift the eigenvalue approximations
649: *           according to their representation. This is necessary for
650: *           a uniform SLARRV since dqds computes eigenvalues of the
651: *           shifted representation. In SLARRV, W will always hold the
652: *           UNshifted eigenvalue approximation.
653:             DO 134 J=WBEGIN,WEND
654:                W(J) = W(J) - SIGMA
655:                WERR(J) = WERR(J) + ABS(W(J)) * EPS
656:  134        CONTINUE
657: *           call SLARRB to reduce eigenvalue error of the approximations
658: *           from SLARRD
659:             DO 135 I = IBEGIN, IEND-1
660:                WORK( I ) = D( I ) * E( I )**2
661:  135        CONTINUE
662: *           use bisection to find EV from INDL to INDU
663:             CALL SLARRB(IN, D(IBEGIN), WORK(IBEGIN),
664:      $                  INDL, INDU, RTOL1, RTOL2, INDL-1,
665:      $                  W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
666:      $                  WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
667:      $                  IN, IINFO )
668:             IF( IINFO .NE. 0 ) THEN
669:                INFO = -4
670:                RETURN
671:             END IF
672: *           SLARRB computes all gaps correctly except for the last one
673: *           Record distance to VU/GU
674:             WGAP( WEND ) = MAX( ZERO,
675:      $           ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
676:             DO 138 I = INDL, INDU
677:                M = M + 1
678:                IBLOCK(M) = JBLK
679:                INDEXW(M) = I
680:  138        CONTINUE
681:          ELSE
682: *           Call dqds to get all eigs (and then possibly delete unwanted
683: *           eigenvalues).
684: *           Note that dqds finds the eigenvalues of the L D L^T representation
685: *           of T to high relative accuracy. High relative accuracy
686: *           might be lost when the shift of the RRR is subtracted to obtain
687: *           the eigenvalues of T. However, T is not guaranteed to define its
688: *           eigenvalues to high relative accuracy anyway.
689: *           Set RTOL to the order of the tolerance used in SLASQ2
690: *           This is an ESTIMATED error, the worst case bound is 4*N*EPS
691: *           which is usually too large and requires unnecessary work to be
692: *           done by bisection when computing the eigenvectors
693:             RTOL = LOG(REAL(IN)) * FOUR * EPS
694:             J = IBEGIN
695:             DO 140 I = 1, IN - 1
696:                WORK( 2*I-1 ) = ABS( D( J ) )
697:                WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
698:                J = J + 1
699:   140       CONTINUE
700:             WORK( 2*IN-1 ) = ABS( D( IEND ) )
701:             WORK( 2*IN ) = ZERO
702:             CALL SLASQ2( IN, WORK, IINFO )
703:             IF( IINFO .NE. 0 ) THEN
704: *              If IINFO = -5 then an index is part of a tight cluster
705: *              and should be changed. The index is in IWORK(1) and the
706: *              gap is in WORK(N+1)
707:                INFO = -5
708:                RETURN
709:             ELSE
710: *              Test that all eigenvalues are positive as expected
711:                DO 149 I = 1, IN
712:                   IF( WORK( I ).LT.ZERO ) THEN
713:                      INFO = -6
714:                      RETURN
715:                   ENDIF
716:  149           CONTINUE
717:             END IF
718:             IF( SGNDEF.GT.ZERO ) THEN
719:                DO 150 I = INDL, INDU
720:                   M = M + 1
721:                   W( M ) = WORK( IN-I+1 )
722:                   IBLOCK( M ) = JBLK
723:                   INDEXW( M ) = I
724:  150           CONTINUE
725:             ELSE
726:                DO 160 I = INDL, INDU
727:                   M = M + 1
728:                   W( M ) = -WORK( I )
729:                   IBLOCK( M ) = JBLK
730:                   INDEXW( M ) = I
731:  160           CONTINUE
732:             END IF
733: 
734:             DO 165 I = M - MB + 1, M
735: *              the value of RTOL below should be the tolerance in SLASQ2
736:                WERR( I ) = RTOL * ABS( W(I) )
737:  165        CONTINUE
738:             DO 166 I = M - MB + 1, M - 1
739: *              compute the right gap between the intervals
740:                WGAP( I ) = MAX( ZERO,
741:      $                          W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
742:  166        CONTINUE
743:             WGAP( M ) = MAX( ZERO,
744:      $           ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
745:          END IF
746: *        proceed with next block
747:          IBEGIN = IEND + 1
748:          WBEGIN = WEND + 1
749:  170  CONTINUE
750: *
751: 
752:       RETURN
753: *
754: *     end of SLARRE
755: *
756:       END
757: