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