Home
Browse Files
  1       SUBROUTINE ZLARRV( N, VL, VU, D, L, PIVMIN,
  2      $                   ISPLIT, M, DOL, DOU, MINRGP,
  3      $                   RTOL1, RTOL2, W, WERR, WGAP,
  4      $                   IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
  5      $                   WORK, IWORK, INFO )
  6 *
  7 *  -- LAPACK auxiliary routine (version 3.3.1) --
  8 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  9 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 10 *  -- April 2011                                                      --
 11 *
 12 *     .. Scalar Arguments ..
 13       INTEGER            DOL, DOU, INFO, LDZ, M, N
 14       DOUBLE PRECISION   MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
 15 *     ..
 16 *     .. Array Arguments ..
 17       INTEGER            IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
 18      $                   ISUPPZ( * ), IWORK( * )
 19       DOUBLE PRECISION   D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
 20      $                   WGAP( * ), WORK( * )
 21       COMPLEX*16        Z( LDZ, * )
 22 *     ..
 23 *
 24 *  Purpose
 25 *  =======
 26 *
 27 *  ZLARRV computes the eigenvectors of the tridiagonal matrix
 28 *  T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
 29 *  The input eigenvalues should have been computed by DLARRE.
 30 *
 31 *  Arguments
 32 *  =========
 33 *
 34 *  N       (input) INTEGER
 35 *          The order of the matrix.  N >= 0.
 36 *
 37 *  VL      (input) DOUBLE PRECISION
 38 *  VU      (input) DOUBLE PRECISION
 39 *          Lower and upper bounds of the interval that contains the desired
 40 *          eigenvalues. VL < VU. Needed to compute gaps on the left or right
 41 *          end of the extremal eigenvalues in the desired RANGE.
 42 *
 43 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
 44 *          On entry, the N diagonal elements of the diagonal matrix D.
 45 *          On exit, D may be overwritten.
 46 *
 47 *  L       (input/output) DOUBLE PRECISION array, dimension (N)
 48 *          On entry, the (N-1) subdiagonal elements of the unit
 49 *          bidiagonal matrix L are in elements 1 to N-1 of L
 50 *          (if the matrix is not splitted.) At the end of each block
 51 *          is stored the corresponding shift as given by DLARRE.
 52 *          On exit, L is overwritten.
 53 *
 54 *  PIVMIN  (in) DOUBLE PRECISION
 55 *          The minimum pivot allowed in the Sturm sequence.
 56 *
 57 *  ISPLIT  (input) INTEGER array, dimension (N)
 58 *          The splitting points, at which T breaks up into blocks.
 59 *          The first block consists of rows/columns 1 to
 60 *          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
 61 *          through ISPLIT( 2 ), etc.
 62 *
 63 *  M       (input) INTEGER
 64 *          The total number of input eigenvalues.  0 <= M <= N.
 65 *
 66 *  DOL     (input) INTEGER
 67 *  DOU     (input) INTEGER
 68 *          If the user wants to compute only selected eigenvectors from all
 69 *          the eigenvalues supplied, he can specify an index range DOL:DOU.
 70 *          Or else the setting DOL=1, DOU=M should be applied.
 71 *          Note that DOL and DOU refer to the order in which the eigenvalues
 72 *          are stored in W.
 73 *          If the user wants to compute only selected eigenpairs, then
 74 *          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
 75 *          computed eigenvectors. All other columns of Z are set to zero.
 76 *
 77 *  MINRGP  (input) DOUBLE PRECISION
 78 *
 79 *  RTOL1   (input) DOUBLE PRECISION
 80 *  RTOL2   (input) DOUBLE PRECISION
 81 *           Parameters for bisection.
 82 *           An interval [LEFT,RIGHT] has converged if
 83 *           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
 84 *
 85 *  W       (input/output) DOUBLE PRECISION array, dimension (N)
 86 *          The first M elements of W contain the APPROXIMATE eigenvalues for
 87 *          which eigenvectors are to be computed.  The eigenvalues
 88 *          should be grouped by split-off block and ordered from
 89 *          smallest to largest within the block ( The output array
 90 *          W from DLARRE is expected here ). Furthermore, they are with
 91 *          respect to the shift of the corresponding root representation
 92 *          for their block. On exit, W holds the eigenvalues of the
 93 *          UNshifted matrix.
 94 *
 95 *  WERR    (input/output) DOUBLE PRECISION array, dimension (N)
 96 *          The first M elements contain the semiwidth of the uncertainty
 97 *          interval of the corresponding eigenvalue in W
 98 *
 99 *  WGAP    (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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*16       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) DOUBLE PRECISION 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 ZLARRV.
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 DLARRB when refining a child's eigenvalues.
149 *          =-2:  Problem in DLARRF 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 DLARRB 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*16         CZERO
177       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
178       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, HALF
179       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
180      $                     TWO = 2.0D0, THREE = 3.0D0,
181      $                     FOUR = 4.0D0, HALF = 0.5D0)
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       DOUBLE PRECISION   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       DOUBLE PRECISION   DLAMCH
202       EXTERNAL           DLAMCH
203 *     ..
204 *     .. External Subroutines ..
205       EXTERNAL           DCOPY, DLARRB, DLARRF, ZDSCAL, ZLAR1V,
206      $                   ZLASET
207 *     ..
208 *     .. Intrinsic Functions ..
209       INTRINSIC ABS, DBLE, MAX, MIN
210       INTRINSIC DCMPLX
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 ZLASET( 'Full', N, ZUSEDW, CZERO, CZERO,
256      $                    Z(1,ZUSEDL), LDZ )
257 
258       EPS = DLAMCH( '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 ) = DCMPLX( 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 DCOPY( 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 ) = DBLE( Z( IBEGIN+K-1,
418      $                                 J ) )
419                      L( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1,
420      $                                 J+1 ) )
421    45             CONTINUE
422                   D( IEND ) = DBLE( Z( IEND, J ) )
423                   SIGMA = DBLE( Z( IEND, J+1 ) )
424 
425 *                 Set the corresponding entries in Z to zero
426                   CALL ZLASET( '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., the P-OFFSET
443 *                 through the Q-OFFSET elements of these arrays are to be used.
444 *                  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 DLARRB( 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 DLARRB.
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 DLARRB( 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 *                    DLARRF needs LWORK = 2*N
584                      CALL DLARRF( 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, DLARRF 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      $                     DCMPLX( WORK( INDIN1+K-1 ), ZERO )
597                         Z( IBEGIN+K-1, NEWFTT+1 ) =
598      $                     DCMPLX( WORK( INDIN2+K-1 ), ZERO )
599    56                CONTINUE
600                      Z( IEND, NEWFTT ) =
601      $                  DCMPLX( WORK( INDIN1+IN-1 ), ZERO )
602                      IF( IINFO.EQ.0 ) THEN
603 *                       a new RRR for the cluster was found by DLARRF
604 *                       update shift and store it
605                         SSIGMA = SIGMA + TAU
606                         Z( IEND, NEWFTT+1 ) = DCMPLX( 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(DBLE(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 DLARRB( 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 ZLAR1V( 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 ZLAR1V( 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 ZDSCAL( 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 ZLARRV
916 *
917       END
Home