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