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: *  -- 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:       REAL              PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
016: *     ..
017: *     .. Array Arguments ..
018:       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
019:      $                   INDEXW( * )
020:       REAL               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, SLARRE 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: *  SSTEMR 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 SLASQ2) to
037: *  conpute all and then discard any unwanted one.
038: *  As an added benefit, SLARRE 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) REAL            
055: *  VU      (input/output) REAL            
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', SLARRE 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) REAL             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) REAL             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) REAL             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) REAL            
089: *  RTOL2   (input) REAL            
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) REAL            
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) REAL             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 ( SLARRE may use the
115: *          remaining N-M elements as workspace).
116: *
117: *  WERR    (output) REAL             array, dimension (N)
118: *          The error bound on the corresponding eigenvalue in W.
119: *
120: *  WGAP    (output) REAL             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) REAL             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) REAL             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 SLARRE.
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 SLARRD.
158: *          = 2:  No base representation could be found in MAXTRY iterations.
159: *                Increasing MAXTRY and recompilation might be a remedy.
160: *          =-3:  Problem in SLARRB when computing the refined root
161: *                representation for SLASQ2.
162: *          =-4:  Problem in SLARRB when preforming bisection on the
163: *                desired part of the spectrum.
164: *          =-5:  Problem in SLASQ2.
165: *          =-6:  Problem in SLASQ2.
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:       REAL               FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
184:      $                   MAXGROWTH, ONE, PERT, TWO, ZERO
185:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
186:      $                     TWO = 2.0E0, FOUR=4.0E0,
187:      $                     HNDRD = 100.0E0,
188:      $                     PERT = 4.0E0,
189:      $                     HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
190:      $                     MAXGROWTH = 64.0E0, FUDGE = 2.0E0 )
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:       REAL               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:       REAL                        SLAMCH
213:       EXTERNAL           SLAMCH, LSAME
214: 
215: *     ..
216: *     .. External Subroutines ..
217:       EXTERNAL           SCOPY, SLARNV, SLARRA, SLARRB, SLARRC, SLARRD,
218:      $                   SLASQ2
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 = SLAMCH( 'S' )
244:       EPS = SLAMCH( 'P' )
245: 
246: *     Set parameters
247:       RTL = HNDRD*EPS
248: *     If one were ever to ask for less initial precision in BSRTOL,
249: *     one should keep in mind that for the subset case, the extremal
250: *     eigenvalues must be at least as accurate as the current setting
251: *     (eigenvalues in the middle need not as much accuracy)
252:       BSRTOL = SQRT(EPS)*(0.5E-3)
253: 
254: *     Treat case of 1x1 matrix for quick return
255:       IF( N.EQ.1 ) THEN
256:          IF( (IRANGE.EQ.ALLRNG).OR.
257:      $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
258:      $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
259:             M = 1
260:             W(1) = D(1)
261: *           The computation error of the eigenvalue is zero
262:             WERR(1) = ZERO
263:             WGAP(1) = ZERO
264:             IBLOCK( 1 ) = 1
265:             INDEXW( 1 ) = 1
266:             GERS(1) = D( 1 )
267:             GERS(2) = D( 1 )
268:          ENDIF
269: *        store the shift for the initial RRR, which is zero in this case
270:          E(1) = ZERO
271:          RETURN
272:       END IF
273: 
274: *     General case: tridiagonal matrix of order > 1
275: *
276: *     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
277: *     Compute maximum off-diagonal entry and pivmin.
278:       GL = D(1)
279:       GU = D(1)
280:       EOLD = ZERO
281:       EMAX = ZERO
282:       E(N) = ZERO
283:       DO 5 I = 1,N
284:          WERR(I) = ZERO
285:          WGAP(I) = ZERO
286:          EABS = ABS( E(I) )
287:          IF( EABS .GE. EMAX ) THEN
288:             EMAX = EABS
289:          END IF
290:          TMP1 = EABS + EOLD
291:          GERS( 2*I-1) = D(I) - TMP1
292:          GL =  MIN( GL, GERS( 2*I - 1))
293:          GERS( 2*I ) = D(I) + TMP1
294:          GU = MAX( GU, GERS(2*I) )
295:          EOLD  = EABS
296:  5    CONTINUE
297: *     The minimum pivot allowed in the Sturm sequence for T
298:       PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
299: *     Compute spectral diameter. The Gerschgorin bounds give an
300: *     estimate that is wrong by at most a factor of SQRT(2)
301:       SPDIAM = GU - GL
302: 
303: *     Compute splitting points
304:       CALL SLARRA( N, D, E, E2, SPLTOL, SPDIAM,
305:      $                    NSPLIT, ISPLIT, IINFO )
306: 
307: *     Can force use of bisection instead of faster DQDS.
308: *     Option left in the code for future multisection work.
309:       FORCEB = .FALSE.
310: 
311: *     Initialize USEDQD, DQDS should be used for ALLRNG unless someone
312: *     explicitly wants bisection.
313:       USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
314: 
315:       IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
316: *        Set interval [VL,VU] that contains all eigenvalues
317:          VL = GL
318:          VU = GU
319:       ELSE
320: *        We call SLARRD to find crude approximations to the eigenvalues
321: *        in the desired range. In case IRANGE = INDRNG, we also obtain the
322: *        interval (VL,VU] that contains all the wanted eigenvalues.
323: *        An interval [LEFT,RIGHT] has converged if
324: *        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
325: *        SLARRD needs a WORK of size 4*N, IWORK of size 3*N
326:          CALL SLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
327:      $                    BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
328:      $                    MM, W, WERR, VL, VU, IBLOCK, INDEXW,
329:      $                    WORK, IWORK, IINFO )
330:          IF( IINFO.NE.0 ) THEN
331:             INFO = -1
332:             RETURN
333:          ENDIF
334: *        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
335:          DO 14 I = MM+1,N
336:             W( I ) = ZERO
337:             WERR( I ) = ZERO
338:             IBLOCK( I ) = 0
339:             INDEXW( I ) = 0
340:  14      CONTINUE
341:       END IF
342: 
343: 
344: ***
345: *     Loop over unreduced blocks
346:       IBEGIN = 1
347:       WBEGIN = 1
348:       DO 170 JBLK = 1, NSPLIT
349:          IEND = ISPLIT( JBLK )
350:          IN = IEND - IBEGIN + 1
351: 
352: *        1 X 1 block
353:          IF( IN.EQ.1 ) THEN
354:             IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
355:      $         ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
356:      $        .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
357:      $        ) THEN
358:                M = M + 1
359:                W( M ) = D( IBEGIN )
360:                WERR(M) = ZERO
361: *              The gap for a single block doesn't matter for the later
362: *              algorithm and is assigned an arbitrary large value
363:                WGAP(M) = ZERO
364:                IBLOCK( M ) = JBLK
365:                INDEXW( M ) = 1
366:                WBEGIN = WBEGIN + 1
367:             ENDIF
368: *           E( IEND ) holds the shift for the initial RRR
369:             E( IEND ) = ZERO
370:             IBEGIN = IEND + 1
371:             GO TO 170
372:          END IF
373: *
374: *        Blocks of size larger than 1x1
375: *
376: *        E( IEND ) will hold the shift for the initial RRR, for now set it =0
377:          E( IEND ) = ZERO
378: *
379: *        Find local outer bounds GL,GU for the block
380:          GL = D(IBEGIN)
381:          GU = D(IBEGIN)
382:          DO 15 I = IBEGIN , IEND
383:             GL = MIN( GERS( 2*I-1 ), GL )
384:             GU = MAX( GERS( 2*I ), GU )
385:  15      CONTINUE
386:          SPDIAM = GU - GL
387: 
388:          IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
389: *           Count the number of eigenvalues in the current block.
390:             MB = 0
391:             DO 20 I = WBEGIN,MM
392:                IF( IBLOCK(I).EQ.JBLK ) THEN
393:                   MB = MB+1
394:                ELSE
395:                   GOTO 21
396:                ENDIF
397:  20         CONTINUE
398:  21         CONTINUE
399: 
400:             IF( MB.EQ.0) THEN
401: *              No eigenvalue in the current block lies in the desired range
402: *              E( IEND ) holds the shift for the initial RRR
403:                E( IEND ) = ZERO
404:                IBEGIN = IEND + 1
405:                GO TO 170
406:             ELSE
407: 
408: *              Decide whether dqds or bisection is more efficient
409:                USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
410:                WEND = WBEGIN + MB - 1
411: *              Calculate gaps for the current block
412: *              In later stages, when representations for individual
413: *              eigenvalues are different, we use SIGMA = E( IEND ).
414:                SIGMA = ZERO
415:                DO 30 I = WBEGIN, WEND - 1
416:                   WGAP( I ) = MAX( ZERO,
417:      $                        W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
418:  30            CONTINUE
419:                WGAP( WEND ) = MAX( ZERO,
420:      $                     VU - SIGMA - (W( WEND )+WERR( WEND )))
421: *              Find local index of the first and last desired evalue.
422:                INDL = INDEXW(WBEGIN)
423:                INDU = INDEXW( WEND )
424:             ENDIF
425:          ENDIF
426:          IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
427: *           Case of DQDS
428: *           Find approximations to the extremal eigenvalues of the block
429:             CALL SLARRK( IN, 1, GL, GU, D(IBEGIN),
430:      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
431:             IF( IINFO.NE.0 ) THEN
432:                INFO = -1
433:                RETURN
434:             ENDIF
435:             ISLEFT = MAX(GL, TMP - TMP1
436:      $               - HNDRD * EPS* ABS(TMP - TMP1))
437: 
438:             CALL SLARRK( IN, IN, GL, GU, D(IBEGIN),
439:      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
440:             IF( IINFO.NE.0 ) THEN
441:                INFO = -1
442:                RETURN
443:             ENDIF
444:             ISRGHT = MIN(GU, TMP + TMP1
445:      $                 + HNDRD * EPS * ABS(TMP + TMP1))
446: *           Improve the estimate of the spectral diameter
447:             SPDIAM = ISRGHT - ISLEFT
448:          ELSE
449: *           Case of bisection
450: *           Find approximations to the wanted extremal eigenvalues
451:             ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
452:      $                  - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
453:             ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
454:      $                  + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
455:          ENDIF
456: 
457: 
458: *        Decide whether the base representation for the current block
459: *        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
460: *        should be on the left or the right end of the current block.
461: *        The strategy is to shift to the end which is "more populated"
462: *        Furthermore, decide whether to use DQDS for the computation of
463: *        the eigenvalue approximations at the end of SLARRE or bisection.
464: *        dqds is chosen if all eigenvalues are desired or the number of
465: *        eigenvalues to be computed is large compared to the blocksize.
466:          IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
467: *           If all the eigenvalues have to be computed, we use dqd
468:             USEDQD = .TRUE.
469: *           INDL is the local index of the first eigenvalue to compute
470:             INDL = 1
471:             INDU = IN
472: *           MB =  number of eigenvalues to compute
473:             MB = IN
474:             WEND = WBEGIN + MB - 1
475: *           Define 1/4 and 3/4 points of the spectrum
476:             S1 = ISLEFT + FOURTH * SPDIAM
477:             S2 = ISRGHT - FOURTH * SPDIAM
478:          ELSE
479: *           SLARRD has computed IBLOCK and INDEXW for each eigenvalue
480: *           approximation.
481: *           choose sigma
482:             IF( USEDQD ) THEN
483:                S1 = ISLEFT + FOURTH * SPDIAM
484:                S2 = ISRGHT - FOURTH * SPDIAM
485:             ELSE
486:                TMP = MIN(ISRGHT,VU) -  MAX(ISLEFT,VL)
487:                S1 =  MAX(ISLEFT,VL) + FOURTH * TMP
488:                S2 =  MIN(ISRGHT,VU) - FOURTH * TMP
489:             ENDIF
490:          ENDIF
491: 
492: *        Compute the negcount at the 1/4 and 3/4 points
493:          IF(MB.GT.1) THEN
494:             CALL SLARRC( 'T', IN, S1, S2, D(IBEGIN),
495:      $                    E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
496:          ENDIF
497: 
498:          IF(MB.EQ.1) THEN
499:             SIGMA = GL
500:             SGNDEF = ONE
501:          ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
502:             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
503:                SIGMA = MAX(ISLEFT,GL)
504:             ELSEIF( USEDQD ) THEN
505: *              use Gerschgorin bound as shift to get pos def matrix
506: *              for dqds
507:                SIGMA = ISLEFT
508:             ELSE
509: *              use approximation of the first desired eigenvalue of the
510: *              block as shift
511:                SIGMA = MAX(ISLEFT,VL)
512:             ENDIF
513:             SGNDEF = ONE
514:          ELSE
515:             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
516:                SIGMA = MIN(ISRGHT,GU)
517:             ELSEIF( USEDQD ) THEN
518: *              use Gerschgorin bound as shift to get neg def matrix
519: *              for dqds
520:                SIGMA = ISRGHT
521:             ELSE
522: *              use approximation of the first desired eigenvalue of the
523: *              block as shift
524:                SIGMA = MIN(ISRGHT,VU)
525:             ENDIF
526:             SGNDEF = -ONE
527:          ENDIF
528: 
529: 
530: *        An initial SIGMA has been chosen that will be used for computing
531: *        T - SIGMA I = L D L^T
532: *        Define the increment TAU of the shift in case the initial shift
533: *        needs to be refined to obtain a factorization with not too much
534: *        element growth.
535:          IF( USEDQD ) THEN
536: *           The initial SIGMA was to the outer end of the spectrum
537: *           the matrix is definite and we need not retreat.
538:             TAU = SPDIAM*EPS*N + TWO*PIVMIN
539:          ELSE
540:             IF(MB.GT.1) THEN
541:                CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
542:                AVGAP = ABS(CLWDTH / REAL(WEND-WBEGIN))
543:                IF( SGNDEF.EQ.ONE ) THEN
544:                   TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
545:                   TAU = MAX(TAU,WERR(WBEGIN))
546:                ELSE
547:                   TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
548:                   TAU = MAX(TAU,WERR(WEND))
549:                ENDIF
550:             ELSE
551:                TAU = WERR(WBEGIN)
552:             ENDIF
553:          ENDIF
554: *
555:          DO 80 IDUM = 1, MAXTRY
556: *           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
557: *           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
558: *           pivots in WORK(2*IN+1:3*IN)
559:             DPIVOT = D( IBEGIN ) - SIGMA
560:             WORK( 1 ) = DPIVOT
561:             DMAX = ABS( WORK(1) )
562:             J = IBEGIN
563:             DO 70 I = 1, IN - 1
564:                WORK( 2*IN+I ) = ONE / WORK( I )
565:                TMP = E( J )*WORK( 2*IN+I )
566:                WORK( IN+I ) = TMP
567:                DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
568:                WORK( I+1 ) = DPIVOT
569:                DMAX = MAX( DMAX, ABS(DPIVOT) )
570:                J = J + 1
571:  70         CONTINUE
572: *           check for element growth
573:             IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
574:                NOREP = .TRUE.
575:             ELSE
576:                NOREP = .FALSE.
577:             ENDIF
578:             IF( USEDQD .AND. .NOT.NOREP ) THEN
579: *              Ensure the definiteness of the representation
580: *              All entries of D (of L D L^T) must have the same sign
581:                DO 71 I = 1, IN
582:                   TMP = SGNDEF*WORK( I )
583:                   IF( TMP.LT.ZERO ) NOREP = .TRUE.
584:  71            CONTINUE
585:             ENDIF
586:             IF(NOREP) THEN
587: *              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
588: *              shift which makes the matrix definite. So we should end up
589: *              here really only in the case of IRANGE = VALRNG or INDRNG.
590:                IF( IDUM.EQ.MAXTRY-1 ) THEN
591:                   IF( SGNDEF.EQ.ONE ) THEN
592: *                    The fudged Gerschgorin shift should succeed
593:                      SIGMA =
594:      $                    GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
595:                   ELSE
596:                      SIGMA =
597:      $                    GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
598:                   END IF
599:                ELSE
600:                   SIGMA = SIGMA - SGNDEF * TAU
601:                   TAU = TWO * TAU
602:                END IF
603:             ELSE
604: *              an initial RRR is found
605:                GO TO 83
606:             END IF
607:  80      CONTINUE
608: *        if the program reaches this point, no base representation could be
609: *        found in MAXTRY iterations.
610:          INFO = 2
611:          RETURN
612: 
613:  83      CONTINUE
614: *        At this point, we have found an initial base representation
615: *        T - SIGMA I = L D L^T with not too much element growth.
616: *        Store the shift.
617:          E( IEND ) = SIGMA
618: *        Store D and L.
619:          CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
620:          CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
621: 
622: 
623:          IF(MB.GT.1 ) THEN
624: *
625: *           Perturb each entry of the base representation by a small
626: *           (but random) relative amount to overcome difficulties with
627: *           glued matrices.
628: *
629:             DO 122 I = 1, 4
630:                ISEED( I ) = 1
631:  122        CONTINUE
632: 
633:             CALL SLARNV(2, ISEED, 2*IN-1, WORK(1))
634:             DO 125 I = 1,IN-1
635:                D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
636:                E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
637:  125        CONTINUE
638:             D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
639: *
640:          ENDIF
641: *
642: *        Don't update the Gerschgorin intervals because keeping track
643: *        of the updates would be too much work in SLARRV.
644: *        We update W instead and use it to locate the proper Gerschgorin
645: *        intervals.
646: 
647: *        Compute the required eigenvalues of L D L' by bisection or dqds
648:          IF ( .NOT.USEDQD ) THEN
649: *           If SLARRD has been used, shift the eigenvalue approximations
650: *           according to their representation. This is necessary for
651: *           a uniform SLARRV since dqds computes eigenvalues of the
652: *           shifted representation. In SLARRV, W will always hold the
653: *           UNshifted eigenvalue approximation.
654:             DO 134 J=WBEGIN,WEND
655:                W(J) = W(J) - SIGMA
656:                WERR(J) = WERR(J) + ABS(W(J)) * EPS
657:  134        CONTINUE
658: *           call SLARRB to reduce eigenvalue error of the approximations
659: *           from SLARRD
660:             DO 135 I = IBEGIN, IEND-1
661:                WORK( I ) = D( I ) * E( I )**2
662:  135        CONTINUE
663: *           use bisection to find EV from INDL to INDU
664:             CALL SLARRB(IN, D(IBEGIN), WORK(IBEGIN),
665:      $                  INDL, INDU, RTOL1, RTOL2, INDL-1,
666:      $                  W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
667:      $                  WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
668:      $                  IN, IINFO )
669:             IF( IINFO .NE. 0 ) THEN
670:                INFO = -4
671:                RETURN
672:             END IF
673: *           SLARRB computes all gaps correctly except for the last one
674: *           Record distance to VU/GU
675:             WGAP( WEND ) = MAX( ZERO,
676:      $           ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
677:             DO 138 I = INDL, INDU
678:                M = M + 1
679:                IBLOCK(M) = JBLK
680:                INDEXW(M) = I
681:  138        CONTINUE
682:          ELSE
683: *           Call dqds to get all eigs (and then possibly delete unwanted
684: *           eigenvalues).
685: *           Note that dqds finds the eigenvalues of the L D L^T representation
686: *           of T to high relative accuracy. High relative accuracy
687: *           might be lost when the shift of the RRR is subtracted to obtain
688: *           the eigenvalues of T. However, T is not guaranteed to define its
689: *           eigenvalues to high relative accuracy anyway.
690: *           Set RTOL to the order of the tolerance used in SLASQ2
691: *           This is an ESTIMATED error, the worst case bound is 4*N*EPS
692: *           which is usually too large and requires unnecessary work to be
693: *           done by bisection when computing the eigenvectors
694:             RTOL = LOG(REAL(IN)) * FOUR * EPS
695:             J = IBEGIN
696:             DO 140 I = 1, IN - 1
697:                WORK( 2*I-1 ) = ABS( D( J ) )
698:                WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
699:                J = J + 1
700:   140       CONTINUE
701:             WORK( 2*IN-1 ) = ABS( D( IEND ) )
702:             WORK( 2*IN ) = ZERO
703:             CALL SLASQ2( IN, WORK, IINFO )
704:             IF( IINFO .NE. 0 ) THEN
705: *              If IINFO = -5 then an index is part of a tight cluster
706: *              and should be changed. The index is in IWORK(1) and the
707: *              gap is in WORK(N+1)
708:                INFO = -5
709:                RETURN
710:             ELSE
711: *              Test that all eigenvalues are positive as expected
712:                DO 149 I = 1, IN
713:                   IF( WORK( I ).LT.ZERO ) THEN
714:                      INFO = -6
715:                      RETURN
716:                   ENDIF
717:  149           CONTINUE
718:             END IF
719:             IF( SGNDEF.GT.ZERO ) THEN
720:                DO 150 I = INDL, INDU
721:                   M = M + 1
722:                   W( M ) = WORK( IN-I+1 )
723:                   IBLOCK( M ) = JBLK
724:                   INDEXW( M ) = I
725:  150           CONTINUE
726:             ELSE
727:                DO 160 I = INDL, INDU
728:                   M = M + 1
729:                   W( M ) = -WORK( I )
730:                   IBLOCK( M ) = JBLK
731:                   INDEXW( M ) = I
732:  160           CONTINUE
733:             END IF
734: 
735:             DO 165 I = M - MB + 1, M
736: *              the value of RTOL below should be the tolerance in SLASQ2
737:                WERR( I ) = RTOL * ABS( W(I) )
738:  165        CONTINUE
739:             DO 166 I = M - MB + 1, M - 1
740: *              compute the right gap between the intervals
741:                WGAP( I ) = MAX( ZERO,
742:      $                          W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
743:  166        CONTINUE
744:             WGAP( M ) = MAX( ZERO,
745:      $           ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
746:          END IF
747: *        proceed with next block
748:          IBEGIN = IEND + 1
749:          WBEGIN = WEND + 1
750:  170  CONTINUE
751: *
752: 
753:       RETURN
754: *
755: *     end of SLARRE
756: *
757:       END
758: