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