001:       SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
002:      $                   ISPLIT, M, DOL, DOU, MINRGP,
003:      $                   RTOL1, RTOL2, W, WERR, WGAP,
004:      $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
005:      $                   WORK, IWORK, INFO )
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:       INTEGER            DOL, DOU, INFO, LDZ, M, N
013:       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
017:      $                   ISUPPZ( * ), IWORK( * )
018:       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
019:      $                   WGAP( * ), WORK( * )
020:       DOUBLE PRECISION  Z( LDZ, * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  DLARRV computes the eigenvectors of the tridiagonal matrix
027: *  T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T.
028: *  The input eigenvalues should have been computed by DLARRE.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  N       (input) INTEGER
034: *          The order of the matrix.  N >= 0.
035: *
036: *  VL      (input) DOUBLE PRECISION
037: *  VU      (input) DOUBLE PRECISION
038: *          Lower and upper bounds of the interval that contains the desired
039: *          eigenvalues. VL < VU. Needed to compute gaps on the left or right
040: *          end of the extremal eigenvalues in the desired RANGE.
041: *
042: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
043: *          On entry, the N diagonal elements of the diagonal matrix D.
044: *          On exit, D may be overwritten.
045: *
046: *  L       (input/output) DOUBLE PRECISION array, dimension (N)
047: *          On entry, the (N-1) subdiagonal elements of the unit
048: *          bidiagonal matrix L are in elements 1 to N-1 of L
049: *          (if the matrix is not splitted.) At the end of each block
050: *          is stored the corresponding shift as given by DLARRE.
051: *          On exit, L is overwritten.
052: *
053: *  PIVMIN  (in) DOUBLE PRECISION
054: *          The minimum pivot allowed in the Sturm sequence.
055: *
056: *  ISPLIT  (input) INTEGER array, dimension (N)
057: *          The splitting points, at which T breaks up into blocks.
058: *          The first block consists of rows/columns 1 to
059: *          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
060: *          through ISPLIT( 2 ), etc.
061: *
062: *  M       (input) INTEGER
063: *          The total number of input eigenvalues.  0 <= M <= N.
064: *
065: *  DOL     (input) INTEGER
066: *  DOU     (input) INTEGER
067: *          If the user wants to compute only selected eigenvectors from all
068: *          the eigenvalues supplied, he can specify an index range DOL:DOU.
069: *          Or else the setting DOL=1, DOU=M should be applied.
070: *          Note that DOL and DOU refer to the order in which the eigenvalues
071: *          are stored in W.
072: *          If the user wants to compute only selected eigenpairs, then
073: *          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
074: *          computed eigenvectors. All other columns of Z are set to zero.
075: *
076: *  MINRGP  (input) DOUBLE PRECISION
077: *
078: *  RTOL1   (input) DOUBLE PRECISION
079: *  RTOL2   (input) DOUBLE PRECISION
080: *           Parameters for bisection.
081: *           An interval [LEFT,RIGHT] has converged if
082: *           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
083: *
084: *  W       (input/output) DOUBLE PRECISION array, dimension (N)
085: *          The first M elements of W contain the APPROXIMATE eigenvalues for
086: *          which eigenvectors are to be computed.  The eigenvalues
087: *          should be grouped by split-off block and ordered from
088: *          smallest to largest within the block ( The output array
089: *          W from DLARRE is expected here ). Furthermore, they are with
090: *          respect to the shift of the corresponding root representation
091: *          for their block. On exit, W holds the eigenvalues of the
092: *          UNshifted matrix.
093: *
094: *  WERR    (input/output) DOUBLE PRECISION array, dimension (N)
095: *          The first M elements contain the semiwidth of the uncertainty
096: *          interval of the corresponding eigenvalue in W
097: *
098: *  WGAP    (input/output) DOUBLE PRECISION array, dimension (N)
099: *          The separation from the right neighbor eigenvalue in W.
100: *
101: *  IBLOCK  (input) INTEGER array, dimension (N)
102: *          The indices of the blocks (submatrices) associated with the
103: *          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
104: *          W(i) belongs to the first block from the top, =2 if W(i)
105: *          belongs to the second block, etc.
106: *
107: *  INDEXW  (input) INTEGER array, dimension (N)
108: *          The indices of the eigenvalues within each block (submatrix);
109: *          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
110: *          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
111: *
112: *  GERS    (input) DOUBLE PRECISION array, dimension (2*N)
113: *          The N Gerschgorin intervals (the i-th Gerschgorin interval
114: *          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
115: *          be computed from the original UNshifted matrix.
116: *
117: *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
118: *          If INFO = 0, the first M columns of Z contain the
119: *          orthonormal eigenvectors of the matrix T
120: *          corresponding to the input eigenvalues, with the i-th
121: *          column of Z holding the eigenvector associated with W(i).
122: *          Note: the user must ensure that at least max(1,M) columns are
123: *          supplied in the array Z.
124: *
125: *  LDZ     (input) INTEGER
126: *          The leading dimension of the array Z.  LDZ >= 1, and if
127: *          JOBZ = 'V', LDZ >= max(1,N).
128: *
129: *  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) )
130: *          The support of the eigenvectors in Z, i.e., the indices
131: *          indicating the nonzero elements in Z. The I-th eigenvector
132: *          is nonzero only in elements ISUPPZ( 2*I-1 ) through
133: *          ISUPPZ( 2*I ).
134: *
135: *  WORK    (workspace) DOUBLE PRECISION array, dimension (12*N)
136: *
137: *  IWORK   (workspace) INTEGER array, dimension (7*N)
138: *
139: *  INFO    (output) INTEGER
140: *          = 0:  successful exit
141: *
142: *          > 0:  A problem occured in DLARRV.
143: *          < 0:  One of the called subroutines signaled an internal problem.
144: *                Needs inspection of the corresponding parameter IINFO
145: *                for further information.
146: *
147: *          =-1:  Problem in DLARRB when refining a child's eigenvalues.
148: *          =-2:  Problem in DLARRF when computing the RRR of a child.
149: *                When a child is inside a tight cluster, it can be difficult
150: *                to find an RRR. A partial remedy from the user's point of
151: *                view is to make the parameter MINRGP smaller and recompile.
152: *                However, as the orthogonality of the computed vectors is
153: *                proportional to 1/MINRGP, the user should be aware that
154: *                he might be trading in precision when he decreases MINRGP.
155: *          =-3:  Problem in DLARRB when refining a single eigenvalue
156: *                after the Rayleigh correction was rejected.
157: *          = 5:  The Rayleigh Quotient Iteration failed to converge to
158: *                full accuracy in MAXITR steps.
159: *
160: *  Further Details
161: *  ===============
162: *
163: *  Based on contributions by
164: *     Beresford Parlett, University of California, Berkeley, USA
165: *     Jim Demmel, University of California, Berkeley, USA
166: *     Inderjit Dhillon, University of Texas, Austin, USA
167: *     Osni Marques, LBNL/NERSC, USA
168: *     Christof Voemel, University of California, Berkeley, USA
169: *
170: *  =====================================================================
171: *
172: *     .. Parameters ..
173:       INTEGER            MAXITR
174:       PARAMETER          ( MAXITR = 10 )
175:       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
176:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
177:      $                     TWO = 2.0D0, THREE = 3.0D0,
178:      $                     FOUR = 4.0D0, HALF = 0.5D0)
179: *     ..
180: *     .. Local Scalars ..
181:       LOGICAL            ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
182:       INTEGER            DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
183:      $                   IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
184:      $                   INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
185:      $                   ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
186:      $                   NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
187:      $                   NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
188:      $                   OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
189:      $                   WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
190:      $                   ZUSEDW
191:       DOUBLE PRECISION   BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
192:      $                   LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
193:      $                   RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
194:      $                   SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
195: *     ..
196: *     .. External Functions ..
197:       DOUBLE PRECISION   DLAMCH
198:       EXTERNAL           DLAMCH
199: *     ..
200: *     .. External Subroutines ..
201:       EXTERNAL           DCOPY, DLAR1V, DLARRB, DLARRF, DLASET,
202:      $                   DSCAL
203: *     ..
204: *     .. Intrinsic Functions ..
205:       INTRINSIC ABS, DBLE, MAX, MIN
206: *     ..
207: *     .. Executable Statements ..
208: *     ..
209: 
210: *     The first N entries of WORK are reserved for the eigenvalues
211:       INDLD = N+1
212:       INDLLD= 2*N+1
213:       INDWRK= 3*N+1
214:       MINWSIZE = 12 * N
215: 
216:       DO 5 I= 1,MINWSIZE
217:          WORK( I ) = ZERO
218:  5    CONTINUE
219: 
220: *     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
221: *     factorization used to compute the FP vector
222:       IINDR = 0
223: *     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
224: *     layer and the one above.
225:       IINDC1 = N
226:       IINDC2 = 2*N
227:       IINDWK = 3*N + 1
228: 
229:       MINIWSIZE = 7 * N
230:       DO 10 I= 1,MINIWSIZE
231:          IWORK( I ) = 0
232:  10   CONTINUE
233: 
234:       ZUSEDL = 1
235:       IF(DOL.GT.1) THEN
236: *        Set lower bound for use of Z
237:          ZUSEDL = DOL-1
238:       ENDIF
239:       ZUSEDU = M
240:       IF(DOU.LT.M) THEN
241: *        Set lower bound for use of Z
242:          ZUSEDU = DOU+1
243:       ENDIF
244: *     The width of the part of Z that is used
245:       ZUSEDW = ZUSEDU - ZUSEDL + 1
246: 
247: 
248:       CALL DLASET( 'Full', N, ZUSEDW, ZERO, ZERO,
249:      $                    Z(1,ZUSEDL), LDZ )
250: 
251:       EPS = DLAMCH( 'Precision' )
252:       RQTOL = TWO * EPS
253: *
254: *     Set expert flags for standard code.
255:       TRYRQC = .TRUE.
256: 
257:       IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
258:       ELSE
259: *        Only selected eigenpairs are computed. Since the other evalues
260: *        are not refined by RQ iteration, bisection has to compute to full
261: *        accuracy.
262:          RTOL1 = FOUR * EPS
263:          RTOL2 = FOUR * EPS
264:       ENDIF
265: 
266: *     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
267: *     desired eigenvalues. The support of the nonzero eigenvector
268: *     entries is contained in the interval IBEGIN:IEND.
269: *     Remark that if k eigenpairs are desired, then the eigenvectors
270: *     are stored in k contiguous columns of Z.
271: 
272: *     DONE is the number of eigenvectors already computed
273:       DONE = 0
274:       IBEGIN = 1
275:       WBEGIN = 1
276:       DO 170 JBLK = 1, IBLOCK( M )
277:          IEND = ISPLIT( JBLK )
278:          SIGMA = L( IEND )
279: *        Find the eigenvectors of the submatrix indexed IBEGIN
280: *        through IEND.
281:          WEND = WBEGIN - 1
282:  15      CONTINUE
283:          IF( WEND.LT.M ) THEN
284:             IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
285:                WEND = WEND + 1
286:                GO TO 15
287:             END IF
288:          END IF
289:          IF( WEND.LT.WBEGIN ) THEN
290:             IBEGIN = IEND + 1
291:             GO TO 170
292:          ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
293:             IBEGIN = IEND + 1
294:             WBEGIN = WEND + 1
295:             GO TO 170
296:          END IF
297: 
298: *        Find local spectral diameter of the block
299:          GL = GERS( 2*IBEGIN-1 )
300:          GU = GERS( 2*IBEGIN )
301:          DO 20 I = IBEGIN+1 , IEND
302:             GL = MIN( GERS( 2*I-1 ), GL )
303:             GU = MAX( GERS( 2*I ), GU )
304:  20      CONTINUE
305:          SPDIAM = GU - GL
306: 
307: *        OLDIEN is the last index of the previous block
308:          OLDIEN = IBEGIN - 1
309: *        Calculate the size of the current block
310:          IN = IEND - IBEGIN + 1
311: *        The number of eigenvalues in the current block
312:          IM = WEND - WBEGIN + 1
313: 
314: *        This is for a 1x1 block
315:          IF( IBEGIN.EQ.IEND ) THEN
316:             DONE = DONE+1
317:             Z( IBEGIN, WBEGIN ) = ONE
318:             ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
319:             ISUPPZ( 2*WBEGIN ) = IBEGIN
320:             W( WBEGIN ) = W( WBEGIN ) + SIGMA
321:             WORK( WBEGIN ) = W( WBEGIN )
322:             IBEGIN = IEND + 1
323:             WBEGIN = WBEGIN + 1
324:             GO TO 170
325:          END IF
326: 
327: *        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
328: *        Note that these can be approximations, in this case, the corresp.
329: *        entries of WERR give the size of the uncertainty interval.
330: *        The eigenvalue approximations will be refined when necessary as
331: *        high relative accuracy is required for the computation of the
332: *        corresponding eigenvectors.
333:          CALL DCOPY( IM, W( WBEGIN ), 1,
334:      &                   WORK( WBEGIN ), 1 )
335: 
336: *        We store in W the eigenvalue approximations w.r.t. the original
337: *        matrix T.
338:          DO 30 I=1,IM
339:             W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
340:  30      CONTINUE
341: 
342: 
343: *        NDEPTH is the current depth of the representation tree
344:          NDEPTH = 0
345: *        PARITY is either 1 or 0
346:          PARITY = 1
347: *        NCLUS is the number of clusters for the next level of the
348: *        representation tree, we start with NCLUS = 1 for the root
349:          NCLUS = 1
350:          IWORK( IINDC1+1 ) = 1
351:          IWORK( IINDC1+2 ) = IM
352: 
353: *        IDONE is the number of eigenvectors already computed in the current
354: *        block
355:          IDONE = 0
356: *        loop while( IDONE.LT.IM )
357: *        generate the representation tree for the current block and
358: *        compute the eigenvectors
359:    40    CONTINUE
360:          IF( IDONE.LT.IM ) THEN
361: *           This is a crude protection against infinitely deep trees
362:             IF( NDEPTH.GT.M ) THEN
363:                INFO = -2
364:                RETURN
365:             ENDIF
366: *           breadth first processing of the current level of the representation
367: *           tree: OLDNCL = number of clusters on current level
368:             OLDNCL = NCLUS
369: *           reset NCLUS to count the number of child clusters
370:             NCLUS = 0
371: *
372:             PARITY = 1 - PARITY
373:             IF( PARITY.EQ.0 ) THEN
374:                OLDCLS = IINDC1
375:                NEWCLS = IINDC2
376:             ELSE
377:                OLDCLS = IINDC2
378:                NEWCLS = IINDC1
379:             END IF
380: *           Process the clusters on the current level
381:             DO 150 I = 1, OLDNCL
382:                J = OLDCLS + 2*I
383: *              OLDFST, OLDLST = first, last index of current cluster.
384: *                               cluster indices start with 1 and are relative
385: *                               to WBEGIN when accessing W, WGAP, WERR, Z
386:                OLDFST = IWORK( J-1 )
387:                OLDLST = IWORK( J )
388:                IF( NDEPTH.GT.0 ) THEN
389: *                 Retrieve relatively robust representation (RRR) of cluster
390: *                 that has been computed at the previous level
391: *                 The RRR is stored in Z and overwritten once the eigenvectors
392: *                 have been computed or when the cluster is refined
393: 
394:                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
395: *                    Get representation from location of the leftmost evalue
396: *                    of the cluster
397:                      J = WBEGIN + OLDFST - 1
398:                   ELSE
399:                      IF(WBEGIN+OLDFST-1.LT.DOL) THEN
400: *                       Get representation from the left end of Z array
401:                         J = DOL - 1
402:                      ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
403: *                       Get representation from the right end of Z array
404:                         J = DOU
405:                      ELSE
406:                         J = WBEGIN + OLDFST - 1
407:                      ENDIF
408:                   ENDIF
409:                   CALL DCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 )
410:                   CALL DCOPY( IN-1, Z( IBEGIN, J+1 ), 1, L( IBEGIN ),
411:      $               1 )
412:                   SIGMA = Z( IEND, J+1 )
413: 
414: *                 Set the corresponding entries in Z to zero
415:                   CALL DLASET( 'Full', IN, 2, ZERO, ZERO,
416:      $                         Z( IBEGIN, J), LDZ )
417:                END IF
418: 
419: *              Compute DL and DLL of current RRR
420:                DO 50 J = IBEGIN, IEND-1
421:                   TMP = D( J )*L( J )
422:                   WORK( INDLD-1+J ) = TMP
423:                   WORK( INDLLD-1+J ) = TMP*L( J )
424:    50          CONTINUE
425: 
426:                IF( NDEPTH.GT.0 ) THEN
427: *                 P and Q are index of the first and last eigenvalue to compute
428: *                 within the current block
429:                   P = INDEXW( WBEGIN-1+OLDFST )
430:                   Q = INDEXW( WBEGIN-1+OLDLST )
431: *                 Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET
432: *                 thru' Q-OFFSET elements of these arrays are to be used.
433: C                  OFFSET = P-OLDFST
434:                   OFFSET = INDEXW( WBEGIN ) - 1
435: *                 perform limited bisection (if necessary) to get approximate
436: *                 eigenvalues to the precision needed.
437:                   CALL DLARRB( IN, D( IBEGIN ),
438:      $                         WORK(INDLLD+IBEGIN-1),
439:      $                         P, Q, RTOL1, RTOL2, OFFSET,
440:      $                         WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
441:      $                         WORK( INDWRK ), IWORK( IINDWK ),
442:      $                         PIVMIN, SPDIAM, IN, IINFO )
443:                   IF( IINFO.NE.0 ) THEN
444:                      INFO = -1
445:                      RETURN
446:                   ENDIF
447: *                 We also recompute the extremal gaps. W holds all eigenvalues
448: *                 of the unshifted matrix and must be used for computation
449: *                 of WGAP, the entries of WORK might stem from RRRs with
450: *                 different shifts. The gaps from WBEGIN-1+OLDFST to
451: *                 WBEGIN-1+OLDLST are correctly computed in DLARRB.
452: *                 However, we only allow the gaps to become greater since
453: *                 this is what should happen when we decrease WERR
454:                   IF( OLDFST.GT.1) THEN
455:                      WGAP( WBEGIN+OLDFST-2 ) =
456:      $             MAX(WGAP(WBEGIN+OLDFST-2),
457:      $                 W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
458:      $                 - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
459:                   ENDIF
460:                   IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
461:                      WGAP( WBEGIN+OLDLST-1 ) =
462:      $               MAX(WGAP(WBEGIN+OLDLST-1),
463:      $                   W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
464:      $                   - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
465:                   ENDIF
466: *                 Each time the eigenvalues in WORK get refined, we store
467: *                 the newly found approximation with all shifts applied in W
468:                   DO 53 J=OLDFST,OLDLST
469:                      W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
470:  53               CONTINUE
471:                END IF
472: 
473: *              Process the current node.
474:                NEWFST = OLDFST
475:                DO 140 J = OLDFST, OLDLST
476:                   IF( J.EQ.OLDLST ) THEN
477: *                    we are at the right end of the cluster, this is also the
478: *                    boundary of the child cluster
479:                      NEWLST = J
480:                   ELSE IF ( WGAP( WBEGIN + J -1).GE.
481:      $                    MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
482: *                    the right relative gap is big enough, the child cluster
483: *                    (NEWFST,..,NEWLST) is well separated from the following
484:                      NEWLST = J
485:                    ELSE
486: *                    inside a child cluster, the relative gap is not
487: *                    big enough.
488:                      GOTO 140
489:                   END IF
490: 
491: *                 Compute size of child cluster found
492:                   NEWSIZ = NEWLST - NEWFST + 1
493: 
494: *                 NEWFTT is the place in Z where the new RRR or the computed
495: *                 eigenvector is to be stored
496:                   IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
497: *                    Store representation at location of the leftmost evalue
498: *                    of the cluster
499:                      NEWFTT = WBEGIN + NEWFST - 1
500:                   ELSE
501:                      IF(WBEGIN+NEWFST-1.LT.DOL) THEN
502: *                       Store representation at the left end of Z array
503:                         NEWFTT = DOL - 1
504:                      ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
505: *                       Store representation at the right end of Z array
506:                         NEWFTT = DOU
507:                      ELSE
508:                         NEWFTT = WBEGIN + NEWFST - 1
509:                      ENDIF
510:                   ENDIF
511: 
512:                   IF( NEWSIZ.GT.1) THEN
513: *
514: *                    Current child is not a singleton but a cluster.
515: *                    Compute and store new representation of child.
516: *
517: *
518: *                    Compute left and right cluster gap.
519: *
520: *                    LGAP and RGAP are not computed from WORK because
521: *                    the eigenvalue approximations may stem from RRRs
522: *                    different shifts. However, W hold all eigenvalues
523: *                    of the unshifted matrix. Still, the entries in WGAP
524: *                    have to be computed from WORK since the entries
525: *                    in W might be of the same order so that gaps are not
526: *                    exhibited correctly for very close eigenvalues.
527:                      IF( NEWFST.EQ.1 ) THEN
528:                         LGAP = MAX( ZERO,
529:      $                       W(WBEGIN)-WERR(WBEGIN) - VL )
530:                     ELSE
531:                         LGAP = WGAP( WBEGIN+NEWFST-2 )
532:                      ENDIF
533:                      RGAP = WGAP( WBEGIN+NEWLST-1 )
534: *
535: *                    Compute left- and rightmost eigenvalue of child
536: *                    to high precision in order to shift as close
537: *                    as possible and obtain as large relative gaps
538: *                    as possible
539: *
540:                      DO 55 K =1,2
541:                         IF(K.EQ.1) THEN
542:                            P = INDEXW( WBEGIN-1+NEWFST )
543:                         ELSE
544:                            P = INDEXW( WBEGIN-1+NEWLST )
545:                         ENDIF
546:                         OFFSET = INDEXW( WBEGIN ) - 1
547:                         CALL DLARRB( IN, D(IBEGIN),
548:      $                       WORK( INDLLD+IBEGIN-1 ),P,P,
549:      $                       RQTOL, RQTOL, OFFSET,
550:      $                       WORK(WBEGIN),WGAP(WBEGIN),
551:      $                       WERR(WBEGIN),WORK( INDWRK ),
552:      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
553:      $                       IN, IINFO )
554:  55                  CONTINUE
555: *
556:                      IF((WBEGIN+NEWLST-1.LT.DOL).OR.
557:      $                  (WBEGIN+NEWFST-1.GT.DOU)) THEN
558: *                       if the cluster contains no desired eigenvalues
559: *                       skip the computation of that branch of the rep. tree
560: *
561: *                       We could skip before the refinement of the extremal
562: *                       eigenvalues of the child, but then the representation
563: *                       tree could be different from the one when nothing is
564: *                       skipped. For this reason we skip at this place.
565:                         IDONE = IDONE + NEWLST - NEWFST + 1
566:                         GOTO 139
567:                      ENDIF
568: *
569: *                    Compute RRR of child cluster.
570: *                    Note that the new RRR is stored in Z
571: *
572: C                    DLARRF needs LWORK = 2*N
573:                      CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),
574:      $                         WORK(INDLD+IBEGIN-1),
575:      $                         NEWFST, NEWLST, WORK(WBEGIN),
576:      $                         WGAP(WBEGIN), WERR(WBEGIN),
577:      $                         SPDIAM, LGAP, RGAP, PIVMIN, TAU,
578:      $                         Z(IBEGIN, NEWFTT),Z(IBEGIN, NEWFTT+1),
579:      $                         WORK( INDWRK ), IINFO )
580:                      IF( IINFO.EQ.0 ) THEN
581: *                       a new RRR for the cluster was found by DLARRF
582: *                       update shift and store it
583:                         SSIGMA = SIGMA + TAU
584:                         Z( IEND, NEWFTT+1 ) = SSIGMA
585: *                       WORK() are the midpoints and WERR() the semi-width
586: *                       Note that the entries in W are unchanged.
587:                         DO 116 K = NEWFST, NEWLST
588:                            FUDGE =
589:      $                          THREE*EPS*ABS(WORK(WBEGIN+K-1))
590:                            WORK( WBEGIN + K - 1 ) =
591:      $                          WORK( WBEGIN + K - 1) - TAU
592:                            FUDGE = FUDGE +
593:      $                          FOUR*EPS*ABS(WORK(WBEGIN+K-1))
594: *                          Fudge errors
595:                            WERR( WBEGIN + K - 1 ) =
596:      $                          WERR( WBEGIN + K - 1 ) + FUDGE
597: *                          Gaps are not fudged. Provided that WERR is small
598: *                          when eigenvalues are close, a zero gap indicates
599: *                          that a new representation is needed for resolving
600: *                          the cluster. A fudge could lead to a wrong decision
601: *                          of judging eigenvalues 'separated' which in
602: *                          reality are not. This could have a negative impact
603: *                          on the orthogonality of the computed eigenvectors.
604:  116                    CONTINUE
605: 
606:                         NCLUS = NCLUS + 1
607:                         K = NEWCLS + 2*NCLUS
608:                         IWORK( K-1 ) = NEWFST
609:                         IWORK( K ) = NEWLST
610:                      ELSE
611:                         INFO = -2
612:                         RETURN
613:                      ENDIF
614:                   ELSE
615: *
616: *                    Compute eigenvector of singleton
617: *
618:                      ITER = 0
619: *
620:                      TOL = FOUR * LOG(DBLE(IN)) * EPS
621: *
622:                      K = NEWFST
623:                      WINDEX = WBEGIN + K - 1
624:                      WINDMN = MAX(WINDEX - 1,1)
625:                      WINDPL = MIN(WINDEX + 1,M)
626:                      LAMBDA = WORK( WINDEX )
627:                      DONE = DONE + 1
628: *                    Check if eigenvector computation is to be skipped
629:                      IF((WINDEX.LT.DOL).OR.
630:      $                  (WINDEX.GT.DOU)) THEN
631:                         ESKIP = .TRUE.
632:                         GOTO 125
633:                      ELSE
634:                         ESKIP = .FALSE.
635:                      ENDIF
636:                      LEFT = WORK( WINDEX ) - WERR( WINDEX )
637:                      RIGHT = WORK( WINDEX ) + WERR( WINDEX )
638:                      INDEIG = INDEXW( WINDEX )
639: *                    Note that since we compute the eigenpairs for a child,
640: *                    all eigenvalue approximations are w.r.t the same shift.
641: *                    In this case, the entries in WORK should be used for
642: *                    computing the gaps since they exhibit even very small
643: *                    differences in the eigenvalues, as opposed to the
644: *                    entries in W which might "look" the same.
645: 
646:                      IF( K .EQ. 1) THEN
647: *                       In the case RANGE='I' and with not much initial
648: *                       accuracy in LAMBDA and VL, the formula
649: *                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
650: *                       can lead to an overestimation of the left gap and
651: *                       thus to inadequately early RQI 'convergence'.
652: *                       Prevent this by forcing a small left gap.
653:                         LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
654:                      ELSE
655:                         LGAP = WGAP(WINDMN)
656:                      ENDIF
657:                      IF( K .EQ. IM) THEN
658: *                       In the case RANGE='I' and with not much initial
659: *                       accuracy in LAMBDA and VU, the formula
660: *                       can lead to an overestimation of the right gap and
661: *                       thus to inadequately early RQI 'convergence'.
662: *                       Prevent this by forcing a small right gap.
663:                         RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
664:                      ELSE
665:                         RGAP = WGAP(WINDEX)
666:                      ENDIF
667:                      GAP = MIN( LGAP, RGAP )
668:                      IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
669: *                       The eigenvector support can become wrong
670: *                       because significant entries could be cut off due to a
671: *                       large GAPTOL parameter in LAR1V. Prevent this.
672:                         GAPTOL = ZERO
673:                      ELSE
674:                         GAPTOL = GAP * EPS
675:                      ENDIF
676:                      ISUPMN = IN
677:                      ISUPMX = 1
678: *                    Update WGAP so that it holds the minimum gap
679: *                    to the left or the right. This is crucial in the
680: *                    case where bisection is used to ensure that the
681: *                    eigenvalue is refined up to the required precision.
682: *                    The correct value is restored afterwards.
683:                      SAVGAP = WGAP(WINDEX)
684:                      WGAP(WINDEX) = GAP
685: *                    We want to use the Rayleigh Quotient Correction
686: *                    as often as possible since it converges quadratically
687: *                    when we are close enough to the desired eigenvalue.
688: *                    However, the Rayleigh Quotient can have the wrong sign
689: *                    and lead us away from the desired eigenvalue. In this
690: *                    case, the best we can do is to use bisection.
691:                      USEDBS = .FALSE.
692:                      USEDRQ = .FALSE.
693: *                    Bisection is initially turned off unless it is forced
694:                      NEEDBS =  .NOT.TRYRQC
695:  120                 CONTINUE
696: *                    Check if bisection should be used to refine eigenvalue
697:                      IF(NEEDBS) THEN
698: *                       Take the bisection as new iterate
699:                         USEDBS = .TRUE.
700:                         ITMP1 = IWORK( IINDR+WINDEX )
701:                         OFFSET = INDEXW( WBEGIN ) - 1
702:                         CALL DLARRB( IN, D(IBEGIN),
703:      $                       WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
704:      $                       ZERO, TWO*EPS, OFFSET,
705:      $                       WORK(WBEGIN),WGAP(WBEGIN),
706:      $                       WERR(WBEGIN),WORK( INDWRK ),
707:      $                       IWORK( IINDWK ), PIVMIN, SPDIAM,
708:      $                       ITMP1, IINFO )
709:                         IF( IINFO.NE.0 ) THEN
710:                            INFO = -3
711:                            RETURN
712:                         ENDIF
713:                         LAMBDA = WORK( WINDEX )
714: *                       Reset twist index from inaccurate LAMBDA to
715: *                       force computation of true MINGMA
716:                         IWORK( IINDR+WINDEX ) = 0
717:                      ENDIF
718: *                    Given LAMBDA, compute the eigenvector.
719:                      CALL DLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
720:      $                    L( IBEGIN ), WORK(INDLD+IBEGIN-1),
721:      $                    WORK(INDLLD+IBEGIN-1),
722:      $                    PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
723:      $                    .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
724:      $                    IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
725:      $                    NRMINV, RESID, RQCORR, WORK( INDWRK ) )
726:                      IF(ITER .EQ. 0) THEN
727:                         BSTRES = RESID
728:                         BSTW = LAMBDA
729:                      ELSEIF(RESID.LT.BSTRES) THEN
730:                         BSTRES = RESID
731:                         BSTW = LAMBDA
732:                      ENDIF
733:                      ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
734:                      ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
735:                      ITER = ITER + 1
736: 
737: *                    sin alpha <= |resid|/gap
738: *                    Note that both the residual and the gap are
739: *                    proportional to the matrix, so ||T|| doesn't play
740: *                    a role in the quotient
741: 
742: *
743: *                    Convergence test for Rayleigh-Quotient iteration
744: *                    (omitted when Bisection has been used)
745: *
746:                      IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
747:      $                    RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
748:      $                    THEN
749: *                       We need to check that the RQCORR update doesn't
750: *                       move the eigenvalue away from the desired one and
751: *                       towards a neighbor. -> protection with bisection
752:                         IF(INDEIG.LE.NEGCNT) THEN
753: *                          The wanted eigenvalue lies to the left
754:                            SGNDEF = -ONE
755:                         ELSE
756: *                          The wanted eigenvalue lies to the right
757:                            SGNDEF = ONE
758:                         ENDIF
759: *                       We only use the RQCORR if it improves the
760: *                       the iterate reasonably.
761:                         IF( ( RQCORR*SGNDEF.GE.ZERO )
762:      $                       .AND.( LAMBDA + RQCORR.LE. RIGHT)
763:      $                       .AND.( LAMBDA + RQCORR.GE. LEFT)
764:      $                       ) THEN
765:                            USEDRQ = .TRUE.
766: *                          Store new midpoint of bisection interval in WORK
767:                            IF(SGNDEF.EQ.ONE) THEN
768: *                             The current LAMBDA is on the left of the true
769: *                             eigenvalue
770:                               LEFT = LAMBDA
771: *                             We prefer to assume that the error estimate
772: *                             is correct. We could make the interval not
773: *                             as a bracket but to be modified if the RQCORR
774: *                             chooses to. In this case, the RIGHT side should
775: *                             be modified as follows:
776: *                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
777:                            ELSE
778: *                             The current LAMBDA is on the right of the true
779: *                             eigenvalue
780:                               RIGHT = LAMBDA
781: *                             See comment about assuming the error estimate is
782: *                             correct above.
783: *                              LEFT = MIN(LEFT, LAMBDA + RQCORR)
784:                            ENDIF
785:                            WORK( WINDEX ) =
786:      $                       HALF * (RIGHT + LEFT)
787: *                          Take RQCORR since it has the correct sign and
788: *                          improves the iterate reasonably
789:                            LAMBDA = LAMBDA + RQCORR
790: *                          Update width of error interval
791:                            WERR( WINDEX ) =
792:      $                             HALF * (RIGHT-LEFT)
793:                         ELSE
794:                            NEEDBS = .TRUE.
795:                         ENDIF
796:                         IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
797: *                             The eigenvalue is computed to bisection accuracy
798: *                             compute eigenvector and stop
799:                            USEDBS = .TRUE.
800:                            GOTO 120
801:                         ELSEIF( ITER.LT.MAXITR ) THEN
802:                            GOTO 120
803:                         ELSEIF( ITER.EQ.MAXITR ) THEN
804:                            NEEDBS = .TRUE.
805:                            GOTO 120
806:                         ELSE
807:                            INFO = 5
808:                            RETURN
809:                         END IF
810:                      ELSE
811:                         STP2II = .FALSE.
812:         IF(USEDRQ .AND. USEDBS .AND.
813:      $                     BSTRES.LE.RESID) THEN
814:                            LAMBDA = BSTW
815:                            STP2II = .TRUE.
816:                         ENDIF
817:                         IF (STP2II) THEN
818: *                          improve error angle by second step
819:                            CALL DLAR1V( IN, 1, IN, LAMBDA,
820:      $                          D( IBEGIN ), L( IBEGIN ),
821:      $                          WORK(INDLD+IBEGIN-1),
822:      $                          WORK(INDLLD+IBEGIN-1),
823:      $                          PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
824:      $                          .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
825:      $                          IWORK( IINDR+WINDEX ),
826:      $                          ISUPPZ( 2*WINDEX-1 ),
827:      $                          NRMINV, RESID, RQCORR, WORK( INDWRK ) )
828:                         ENDIF
829:                         WORK( WINDEX ) = LAMBDA
830:                      END IF
831: *
832: *                    Compute FP-vector support w.r.t. whole matrix
833: *
834:                      ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
835:                      ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
836:                      ZFROM = ISUPPZ( 2*WINDEX-1 )
837:                      ZTO = ISUPPZ( 2*WINDEX )
838:                      ISUPMN = ISUPMN + OLDIEN
839:                      ISUPMX = ISUPMX + OLDIEN
840: *                    Ensure vector is ok if support in the RQI has changed
841:                      IF(ISUPMN.LT.ZFROM) THEN
842:                         DO 122 II = ISUPMN,ZFROM-1
843:                            Z( II, WINDEX ) = ZERO
844:  122                    CONTINUE
845:                      ENDIF
846:                      IF(ISUPMX.GT.ZTO) THEN
847:                         DO 123 II = ZTO+1,ISUPMX
848:                            Z( II, WINDEX ) = ZERO
849:  123                    CONTINUE
850:                      ENDIF
851:                      CALL DSCAL( ZTO-ZFROM+1, NRMINV,
852:      $                       Z( ZFROM, WINDEX ), 1 )
853:  125                 CONTINUE
854: *                    Update W
855:                      W( WINDEX ) = LAMBDA+SIGMA
856: *                    Recompute the gaps on the left and right
857: *                    But only allow them to become larger and not
858: *                    smaller (which can only happen through "bad"
859: *                    cancellation and doesn't reflect the theory
860: *                    where the initial gaps are underestimated due
861: *                    to WERR being too crude.)
862:                      IF(.NOT.ESKIP) THEN
863:                         IF( K.GT.1) THEN
864:                            WGAP( WINDMN ) = MAX( WGAP(WINDMN),
865:      $                          W(WINDEX)-WERR(WINDEX)
866:      $                          - W(WINDMN)-WERR(WINDMN) )
867:                         ENDIF
868:                         IF( WINDEX.LT.WEND ) THEN
869:                            WGAP( WINDEX ) = MAX( SAVGAP,
870:      $                          W( WINDPL )-WERR( WINDPL )
871:      $                          - W( WINDEX )-WERR( WINDEX) )
872:                         ENDIF
873:                      ENDIF
874:                      IDONE = IDONE + 1
875:                   ENDIF
876: *                 here ends the code for the current child
877: *
878:  139              CONTINUE
879: *                 Proceed to any remaining child nodes
880:                   NEWFST = J + 1
881:  140           CONTINUE
882:  150        CONTINUE
883:             NDEPTH = NDEPTH + 1
884:             GO TO 40
885:          END IF
886:          IBEGIN = IEND + 1
887:          WBEGIN = WEND + 1
888:  170  CONTINUE
889: *
890: 
891:       RETURN
892: *
893: *     End of DLARRV
894: *
895:       END
896: