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