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.2) --
8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10: * November 2006
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., th P-OFFSET
443: * thru' Q-OFFSET elements of these arrays are to be used.
444: C OFFSET = P-OLDFST
445: OFFSET = INDEXW( WBEGIN ) - 1
446: * perform limited bisection (if necessary) to get approximate
447: * eigenvalues to the precision needed.
448: CALL 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: C 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlarrv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack