1: *> \brief \b DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLARRV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrv.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrv.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrv.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
22: * ISPLIT, M, DOL, DOU, MINRGP,
23: * RTOL1, RTOL2, W, WERR, WGAP,
24: * IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
25: * WORK, IWORK, INFO )
26: *
27: * .. Scalar Arguments ..
28: * INTEGER DOL, DOU, INFO, LDZ, M, N
29: * DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
30: * ..
31: * .. Array Arguments ..
32: * INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
33: * $ ISUPPZ( * ), IWORK( * )
34: * DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
35: * $ WGAP( * ), WORK( * )
36: * DOUBLE PRECISION Z( LDZ, * )
37: * ..
38: *
39: *
40: *> \par Purpose:
41: * =============
42: *>
43: *> \verbatim
44: *>
45: *> DLARRV computes the eigenvectors of the tridiagonal matrix
46: *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
47: *> The input eigenvalues should have been computed by DLARRE.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] N
54: *> \verbatim
55: *> N is INTEGER
56: *> The order of the matrix. N >= 0.
57: *> \endverbatim
58: *>
59: *> \param[in] VL
60: *> \verbatim
61: *> VL is DOUBLE PRECISION
62: *> Lower bound of the interval that contains the desired
63: *> eigenvalues. VL < VU. Needed to compute gaps on the left or right
64: *> end of the extremal eigenvalues in the desired RANGE.
65: *> \endverbatim
66: *>
67: *> \param[in] VU
68: *> \verbatim
69: *> VU is DOUBLE PRECISION
70: *> Upper bound of the interval that contains the desired
71: *> eigenvalues. VL < VU.
72: *> Note: VU is currently not used by this implementation of DLARRV, VU is
73: *> passed to DLARRV because it could be used compute gaps on the right end
74: *> of the extremal eigenvalues. However, with not much initial accuracy in
75: *> LAMBDA and VU, the formula can lead to an overestimation of the right gap
76: *> and thus to inadequately early RQI 'convergence'. This is currently
77: *> prevented this by forcing a small right gap. And so it turns out that VU
78: *> is currently not used by this implementation of DLARRV.
79: *> \endverbatim
80: *>
81: *> \param[in,out] D
82: *> \verbatim
83: *> D is DOUBLE PRECISION array, dimension (N)
84: *> On entry, the N diagonal elements of the diagonal matrix D.
85: *> On exit, D may be overwritten.
86: *> \endverbatim
87: *>
88: *> \param[in,out] L
89: *> \verbatim
90: *> L is DOUBLE PRECISION array, dimension (N)
91: *> On entry, the (N-1) subdiagonal elements of the unit
92: *> bidiagonal matrix L are in elements 1 to N-1 of L
93: *> (if the matrix is not split.) At the end of each block
94: *> is stored the corresponding shift as given by DLARRE.
95: *> On exit, L is overwritten.
96: *> \endverbatim
97: *>
98: *> \param[in] PIVMIN
99: *> \verbatim
100: *> PIVMIN is DOUBLE PRECISION
101: *> The minimum pivot allowed in the Sturm sequence.
102: *> \endverbatim
103: *>
104: *> \param[in] ISPLIT
105: *> \verbatim
106: *> ISPLIT is INTEGER array, dimension (N)
107: *> The splitting points, at which T breaks up into blocks.
108: *> The first block consists of rows/columns 1 to
109: *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
110: *> through ISPLIT( 2 ), etc.
111: *> \endverbatim
112: *>
113: *> \param[in] M
114: *> \verbatim
115: *> M is INTEGER
116: *> The total number of input eigenvalues. 0 <= M <= N.
117: *> \endverbatim
118: *>
119: *> \param[in] DOL
120: *> \verbatim
121: *> DOL is INTEGER
122: *> \endverbatim
123: *>
124: *> \param[in] DOU
125: *> \verbatim
126: *> DOU is INTEGER
127: *> If the user wants to compute only selected eigenvectors from all
128: *> the eigenvalues supplied, he can specify an index range DOL:DOU.
129: *> Or else the setting DOL=1, DOU=M should be applied.
130: *> Note that DOL and DOU refer to the order in which the eigenvalues
131: *> are stored in W.
132: *> If the user wants to compute only selected eigenpairs, then
133: *> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
134: *> computed eigenvectors. All other columns of Z are set to zero.
135: *> \endverbatim
136: *>
137: *> \param[in] MINRGP
138: *> \verbatim
139: *> MINRGP is DOUBLE PRECISION
140: *> \endverbatim
141: *>
142: *> \param[in] RTOL1
143: *> \verbatim
144: *> RTOL1 is DOUBLE PRECISION
145: *> \endverbatim
146: *>
147: *> \param[in] RTOL2
148: *> \verbatim
149: *> RTOL2 is DOUBLE PRECISION
150: *> Parameters for bisection.
151: *> An interval [LEFT,RIGHT] has converged if
152: *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
153: *> \endverbatim
154: *>
155: *> \param[in,out] W
156: *> \verbatim
157: *> W is DOUBLE PRECISION array, dimension (N)
158: *> The first M elements of W contain the APPROXIMATE eigenvalues for
159: *> which eigenvectors are to be computed. The eigenvalues
160: *> should be grouped by split-off block and ordered from
161: *> smallest to largest within the block ( The output array
162: *> W from DLARRE is expected here ). Furthermore, they are with
163: *> respect to the shift of the corresponding root representation
164: *> for their block. On exit, W holds the eigenvalues of the
165: *> UNshifted matrix.
166: *> \endverbatim
167: *>
168: *> \param[in,out] WERR
169: *> \verbatim
170: *> WERR is DOUBLE PRECISION array, dimension (N)
171: *> The first M elements contain the semiwidth of the uncertainty
172: *> interval of the corresponding eigenvalue in W
173: *> \endverbatim
174: *>
175: *> \param[in,out] WGAP
176: *> \verbatim
177: *> WGAP is DOUBLE PRECISION array, dimension (N)
178: *> The separation from the right neighbor eigenvalue in W.
179: *> \endverbatim
180: *>
181: *> \param[in] IBLOCK
182: *> \verbatim
183: *> IBLOCK is INTEGER array, dimension (N)
184: *> The indices of the blocks (submatrices) associated with the
185: *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
186: *> W(i) belongs to the first block from the top, =2 if W(i)
187: *> belongs to the second block, etc.
188: *> \endverbatim
189: *>
190: *> \param[in] INDEXW
191: *> \verbatim
192: *> INDEXW is INTEGER array, dimension (N)
193: *> The indices of the eigenvalues within each block (submatrix);
194: *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
195: *> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
196: *> \endverbatim
197: *>
198: *> \param[in] GERS
199: *> \verbatim
200: *> GERS is DOUBLE PRECISION array, dimension (2*N)
201: *> The N Gerschgorin intervals (the i-th Gerschgorin interval
202: *> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
203: *> be computed from the original UNshifted matrix.
204: *> \endverbatim
205: *>
206: *> \param[out] Z
207: *> \verbatim
208: *> Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
209: *> If INFO = 0, the first M columns of Z contain the
210: *> orthonormal eigenvectors of the matrix T
211: *> corresponding to the input eigenvalues, with the i-th
212: *> column of Z holding the eigenvector associated with W(i).
213: *> Note: the user must ensure that at least max(1,M) columns are
214: *> supplied in the array Z.
215: *> \endverbatim
216: *>
217: *> \param[in] LDZ
218: *> \verbatim
219: *> LDZ is INTEGER
220: *> The leading dimension of the array Z. LDZ >= 1, and if
221: *> JOBZ = 'V', LDZ >= max(1,N).
222: *> \endverbatim
223: *>
224: *> \param[out] ISUPPZ
225: *> \verbatim
226: *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
227: *> The support of the eigenvectors in Z, i.e., the indices
228: *> indicating the nonzero elements in Z. The I-th eigenvector
229: *> is nonzero only in elements ISUPPZ( 2*I-1 ) through
230: *> ISUPPZ( 2*I ).
231: *> \endverbatim
232: *>
233: *> \param[out] WORK
234: *> \verbatim
235: *> WORK is DOUBLE PRECISION array, dimension (12*N)
236: *> \endverbatim
237: *>
238: *> \param[out] IWORK
239: *> \verbatim
240: *> IWORK is INTEGER array, dimension (7*N)
241: *> \endverbatim
242: *>
243: *> \param[out] INFO
244: *> \verbatim
245: *> INFO is INTEGER
246: *> = 0: successful exit
247: *>
248: *> > 0: A problem occurred in DLARRV.
249: *> < 0: One of the called subroutines signaled an internal problem.
250: *> Needs inspection of the corresponding parameter IINFO
251: *> for further information.
252: *>
253: *> =-1: Problem in DLARRB when refining a child's eigenvalues.
254: *> =-2: Problem in DLARRF when computing the RRR of a child.
255: *> When a child is inside a tight cluster, it can be difficult
256: *> to find an RRR. A partial remedy from the user's point of
257: *> view is to make the parameter MINRGP smaller and recompile.
258: *> However, as the orthogonality of the computed vectors is
259: *> proportional to 1/MINRGP, the user should be aware that
260: *> he might be trading in precision when he decreases MINRGP.
261: *> =-3: Problem in DLARRB when refining a single eigenvalue
262: *> after the Rayleigh correction was rejected.
263: *> = 5: The Rayleigh Quotient Iteration failed to converge to
264: *> full accuracy in MAXITR steps.
265: *> \endverbatim
266: *
267: * Authors:
268: * ========
269: *
270: *> \author Univ. of Tennessee
271: *> \author Univ. of California Berkeley
272: *> \author Univ. of Colorado Denver
273: *> \author NAG Ltd.
274: *
275: *> \date June 2016
276: *
277: *> \ingroup doubleOTHERauxiliary
278: *
279: *> \par Contributors:
280: * ==================
281: *>
282: *> Beresford Parlett, University of California, Berkeley, USA \n
283: *> Jim Demmel, University of California, Berkeley, USA \n
284: *> Inderjit Dhillon, University of Texas, Austin, USA \n
285: *> Osni Marques, LBNL/NERSC, USA \n
286: *> Christof Voemel, University of California, Berkeley, USA
287: *
288: * =====================================================================
289: SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN,
290: $ ISPLIT, M, DOL, DOU, MINRGP,
291: $ RTOL1, RTOL2, W, WERR, WGAP,
292: $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
293: $ WORK, IWORK, INFO )
294: *
295: * -- LAPACK auxiliary routine (version 3.8.0) --
296: * -- LAPACK is a software package provided by Univ. of Tennessee, --
297: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
298: * June 2016
299: *
300: * .. Scalar Arguments ..
301: INTEGER DOL, DOU, INFO, LDZ, M, N
302: DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
303: * ..
304: * .. Array Arguments ..
305: INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
306: $ ISUPPZ( * ), IWORK( * )
307: DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
308: $ WGAP( * ), WORK( * )
309: DOUBLE PRECISION Z( LDZ, * )
310: * ..
311: *
312: * =====================================================================
313: *
314: * .. Parameters ..
315: INTEGER MAXITR
316: PARAMETER ( MAXITR = 10 )
317: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF
318: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0,
319: $ TWO = 2.0D0, THREE = 3.0D0,
320: $ FOUR = 4.0D0, HALF = 0.5D0)
321: * ..
322: * .. Local Scalars ..
323: LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
324: INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
325: $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
326: $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
327: $ ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS,
328: $ NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST,
329: $ NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST,
330: $ OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX,
331: $ WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU,
332: $ ZUSEDW
333: DOUBLE PRECISION BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
334: $ LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
335: $ RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
336: $ SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
337: * ..
338: * .. External Functions ..
339: DOUBLE PRECISION DLAMCH
340: EXTERNAL DLAMCH
341: * ..
342: * .. External Subroutines ..
343: EXTERNAL DCOPY, DLAR1V, DLARRB, DLARRF, DLASET,
344: $ DSCAL
345: * ..
346: * .. Intrinsic Functions ..
347: INTRINSIC ABS, DBLE, MAX, MIN
348: * ..
349: * .. Executable Statements ..
350: * ..
351:
352: INFO = 0
353: *
354: * Quick return if possible
355: *
356: IF( N.LE.0 ) THEN
357: RETURN
358: END IF
359: *
360: * The first N entries of WORK are reserved for the eigenvalues
361: INDLD = N+1
362: INDLLD= 2*N+1
363: INDWRK= 3*N+1
364: MINWSIZE = 12 * N
365:
366: DO 5 I= 1,MINWSIZE
367: WORK( I ) = ZERO
368: 5 CONTINUE
369:
370: * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
371: * factorization used to compute the FP vector
372: IINDR = 0
373: * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
374: * layer and the one above.
375: IINDC1 = N
376: IINDC2 = 2*N
377: IINDWK = 3*N + 1
378:
379: MINIWSIZE = 7 * N
380: DO 10 I= 1,MINIWSIZE
381: IWORK( I ) = 0
382: 10 CONTINUE
383:
384: ZUSEDL = 1
385: IF(DOL.GT.1) THEN
386: * Set lower bound for use of Z
387: ZUSEDL = DOL-1
388: ENDIF
389: ZUSEDU = M
390: IF(DOU.LT.M) THEN
391: * Set lower bound for use of Z
392: ZUSEDU = DOU+1
393: ENDIF
394: * The width of the part of Z that is used
395: ZUSEDW = ZUSEDU - ZUSEDL + 1
396:
397:
398: CALL DLASET( 'Full', N, ZUSEDW, ZERO, ZERO,
399: $ Z(1,ZUSEDL), LDZ )
400:
401: EPS = DLAMCH( 'Precision' )
402: RQTOL = TWO * EPS
403: *
404: * Set expert flags for standard code.
405: TRYRQC = .TRUE.
406:
407: IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
408: ELSE
409: * Only selected eigenpairs are computed. Since the other evalues
410: * are not refined by RQ iteration, bisection has to compute to full
411: * accuracy.
412: RTOL1 = FOUR * EPS
413: RTOL2 = FOUR * EPS
414: ENDIF
415:
416: * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
417: * desired eigenvalues. The support of the nonzero eigenvector
418: * entries is contained in the interval IBEGIN:IEND.
419: * Remark that if k eigenpairs are desired, then the eigenvectors
420: * are stored in k contiguous columns of Z.
421:
422: * DONE is the number of eigenvectors already computed
423: DONE = 0
424: IBEGIN = 1
425: WBEGIN = 1
426: DO 170 JBLK = 1, IBLOCK( M )
427: IEND = ISPLIT( JBLK )
428: SIGMA = L( IEND )
429: * Find the eigenvectors of the submatrix indexed IBEGIN
430: * through IEND.
431: WEND = WBEGIN - 1
432: 15 CONTINUE
433: IF( WEND.LT.M ) THEN
434: IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN
435: WEND = WEND + 1
436: GO TO 15
437: END IF
438: END IF
439: IF( WEND.LT.WBEGIN ) THEN
440: IBEGIN = IEND + 1
441: GO TO 170
442: ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN
443: IBEGIN = IEND + 1
444: WBEGIN = WEND + 1
445: GO TO 170
446: END IF
447:
448: * Find local spectral diameter of the block
449: GL = GERS( 2*IBEGIN-1 )
450: GU = GERS( 2*IBEGIN )
451: DO 20 I = IBEGIN+1 , IEND
452: GL = MIN( GERS( 2*I-1 ), GL )
453: GU = MAX( GERS( 2*I ), GU )
454: 20 CONTINUE
455: SPDIAM = GU - GL
456:
457: * OLDIEN is the last index of the previous block
458: OLDIEN = IBEGIN - 1
459: * Calculate the size of the current block
460: IN = IEND - IBEGIN + 1
461: * The number of eigenvalues in the current block
462: IM = WEND - WBEGIN + 1
463:
464: * This is for a 1x1 block
465: IF( IBEGIN.EQ.IEND ) THEN
466: DONE = DONE+1
467: Z( IBEGIN, WBEGIN ) = ONE
468: ISUPPZ( 2*WBEGIN-1 ) = IBEGIN
469: ISUPPZ( 2*WBEGIN ) = IBEGIN
470: W( WBEGIN ) = W( WBEGIN ) + SIGMA
471: WORK( WBEGIN ) = W( WBEGIN )
472: IBEGIN = IEND + 1
473: WBEGIN = WBEGIN + 1
474: GO TO 170
475: END IF
476:
477: * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
478: * Note that these can be approximations, in this case, the corresp.
479: * entries of WERR give the size of the uncertainty interval.
480: * The eigenvalue approximations will be refined when necessary as
481: * high relative accuracy is required for the computation of the
482: * corresponding eigenvectors.
483: CALL DCOPY( IM, W( WBEGIN ), 1,
484: $ WORK( WBEGIN ), 1 )
485:
486: * We store in W the eigenvalue approximations w.r.t. the original
487: * matrix T.
488: DO 30 I=1,IM
489: W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA
490: 30 CONTINUE
491:
492:
493: * NDEPTH is the current depth of the representation tree
494: NDEPTH = 0
495: * PARITY is either 1 or 0
496: PARITY = 1
497: * NCLUS is the number of clusters for the next level of the
498: * representation tree, we start with NCLUS = 1 for the root
499: NCLUS = 1
500: IWORK( IINDC1+1 ) = 1
501: IWORK( IINDC1+2 ) = IM
502:
503: * IDONE is the number of eigenvectors already computed in the current
504: * block
505: IDONE = 0
506: * loop while( IDONE.LT.IM )
507: * generate the representation tree for the current block and
508: * compute the eigenvectors
509: 40 CONTINUE
510: IF( IDONE.LT.IM ) THEN
511: * This is a crude protection against infinitely deep trees
512: IF( NDEPTH.GT.M ) THEN
513: INFO = -2
514: RETURN
515: ENDIF
516: * breadth first processing of the current level of the representation
517: * tree: OLDNCL = number of clusters on current level
518: OLDNCL = NCLUS
519: * reset NCLUS to count the number of child clusters
520: NCLUS = 0
521: *
522: PARITY = 1 - PARITY
523: IF( PARITY.EQ.0 ) THEN
524: OLDCLS = IINDC1
525: NEWCLS = IINDC2
526: ELSE
527: OLDCLS = IINDC2
528: NEWCLS = IINDC1
529: END IF
530: * Process the clusters on the current level
531: DO 150 I = 1, OLDNCL
532: J = OLDCLS + 2*I
533: * OLDFST, OLDLST = first, last index of current cluster.
534: * cluster indices start with 1 and are relative
535: * to WBEGIN when accessing W, WGAP, WERR, Z
536: OLDFST = IWORK( J-1 )
537: OLDLST = IWORK( J )
538: IF( NDEPTH.GT.0 ) THEN
539: * Retrieve relatively robust representation (RRR) of cluster
540: * that has been computed at the previous level
541: * The RRR is stored in Z and overwritten once the eigenvectors
542: * have been computed or when the cluster is refined
543:
544: IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
545: * Get representation from location of the leftmost evalue
546: * of the cluster
547: J = WBEGIN + OLDFST - 1
548: ELSE
549: IF(WBEGIN+OLDFST-1.LT.DOL) THEN
550: * Get representation from the left end of Z array
551: J = DOL - 1
552: ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN
553: * Get representation from the right end of Z array
554: J = DOU
555: ELSE
556: J = WBEGIN + OLDFST - 1
557: ENDIF
558: ENDIF
559: CALL DCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 )
560: CALL DCOPY( IN-1, Z( IBEGIN, J+1 ), 1, L( IBEGIN ),
561: $ 1 )
562: SIGMA = Z( IEND, J+1 )
563:
564: * Set the corresponding entries in Z to zero
565: CALL DLASET( 'Full', IN, 2, ZERO, ZERO,
566: $ Z( IBEGIN, J), LDZ )
567: END IF
568:
569: * Compute DL and DLL of current RRR
570: DO 50 J = IBEGIN, IEND-1
571: TMP = D( J )*L( J )
572: WORK( INDLD-1+J ) = TMP
573: WORK( INDLLD-1+J ) = TMP*L( J )
574: 50 CONTINUE
575:
576: IF( NDEPTH.GT.0 ) THEN
577: * P and Q are index of the first and last eigenvalue to compute
578: * within the current block
579: P = INDEXW( WBEGIN-1+OLDFST )
580: Q = INDEXW( WBEGIN-1+OLDLST )
581: * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
582: * through the Q-OFFSET elements of these arrays are to be used.
583: * OFFSET = P-OLDFST
584: OFFSET = INDEXW( WBEGIN ) - 1
585: * perform limited bisection (if necessary) to get approximate
586: * eigenvalues to the precision needed.
587: CALL DLARRB( IN, D( IBEGIN ),
588: $ WORK(INDLLD+IBEGIN-1),
589: $ P, Q, RTOL1, RTOL2, OFFSET,
590: $ WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN),
591: $ WORK( INDWRK ), IWORK( IINDWK ),
592: $ PIVMIN, SPDIAM, IN, IINFO )
593: IF( IINFO.NE.0 ) THEN
594: INFO = -1
595: RETURN
596: ENDIF
597: * We also recompute the extremal gaps. W holds all eigenvalues
598: * of the unshifted matrix and must be used for computation
599: * of WGAP, the entries of WORK might stem from RRRs with
600: * different shifts. The gaps from WBEGIN-1+OLDFST to
601: * WBEGIN-1+OLDLST are correctly computed in DLARRB.
602: * However, we only allow the gaps to become greater since
603: * this is what should happen when we decrease WERR
604: IF( OLDFST.GT.1) THEN
605: WGAP( WBEGIN+OLDFST-2 ) =
606: $ MAX(WGAP(WBEGIN+OLDFST-2),
607: $ W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1)
608: $ - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) )
609: ENDIF
610: IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN
611: WGAP( WBEGIN+OLDLST-1 ) =
612: $ MAX(WGAP(WBEGIN+OLDLST-1),
613: $ W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST)
614: $ - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) )
615: ENDIF
616: * Each time the eigenvalues in WORK get refined, we store
617: * the newly found approximation with all shifts applied in W
618: DO 53 J=OLDFST,OLDLST
619: W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA
620: 53 CONTINUE
621: END IF
622:
623: * Process the current node.
624: NEWFST = OLDFST
625: DO 140 J = OLDFST, OLDLST
626: IF( J.EQ.OLDLST ) THEN
627: * we are at the right end of the cluster, this is also the
628: * boundary of the child cluster
629: NEWLST = J
630: ELSE IF ( WGAP( WBEGIN + J -1).GE.
631: $ MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN
632: * the right relative gap is big enough, the child cluster
633: * (NEWFST,..,NEWLST) is well separated from the following
634: NEWLST = J
635: ELSE
636: * inside a child cluster, the relative gap is not
637: * big enough.
638: GOTO 140
639: END IF
640:
641: * Compute size of child cluster found
642: NEWSIZ = NEWLST - NEWFST + 1
643:
644: * NEWFTT is the place in Z where the new RRR or the computed
645: * eigenvector is to be stored
646: IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN
647: * Store representation at location of the leftmost evalue
648: * of the cluster
649: NEWFTT = WBEGIN + NEWFST - 1
650: ELSE
651: IF(WBEGIN+NEWFST-1.LT.DOL) THEN
652: * Store representation at the left end of Z array
653: NEWFTT = DOL - 1
654: ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN
655: * Store representation at the right end of Z array
656: NEWFTT = DOU
657: ELSE
658: NEWFTT = WBEGIN + NEWFST - 1
659: ENDIF
660: ENDIF
661:
662: IF( NEWSIZ.GT.1) THEN
663: *
664: * Current child is not a singleton but a cluster.
665: * Compute and store new representation of child.
666: *
667: *
668: * Compute left and right cluster gap.
669: *
670: * LGAP and RGAP are not computed from WORK because
671: * the eigenvalue approximations may stem from RRRs
672: * different shifts. However, W hold all eigenvalues
673: * of the unshifted matrix. Still, the entries in WGAP
674: * have to be computed from WORK since the entries
675: * in W might be of the same order so that gaps are not
676: * exhibited correctly for very close eigenvalues.
677: IF( NEWFST.EQ.1 ) THEN
678: LGAP = MAX( ZERO,
679: $ W(WBEGIN)-WERR(WBEGIN) - VL )
680: ELSE
681: LGAP = WGAP( WBEGIN+NEWFST-2 )
682: ENDIF
683: RGAP = WGAP( WBEGIN+NEWLST-1 )
684: *
685: * Compute left- and rightmost eigenvalue of child
686: * to high precision in order to shift as close
687: * as possible and obtain as large relative gaps
688: * as possible
689: *
690: DO 55 K =1,2
691: IF(K.EQ.1) THEN
692: P = INDEXW( WBEGIN-1+NEWFST )
693: ELSE
694: P = INDEXW( WBEGIN-1+NEWLST )
695: ENDIF
696: OFFSET = INDEXW( WBEGIN ) - 1
697: CALL DLARRB( IN, D(IBEGIN),
698: $ WORK( INDLLD+IBEGIN-1 ),P,P,
699: $ RQTOL, RQTOL, OFFSET,
700: $ WORK(WBEGIN),WGAP(WBEGIN),
701: $ WERR(WBEGIN),WORK( INDWRK ),
702: $ IWORK( IINDWK ), PIVMIN, SPDIAM,
703: $ IN, IINFO )
704: 55 CONTINUE
705: *
706: IF((WBEGIN+NEWLST-1.LT.DOL).OR.
707: $ (WBEGIN+NEWFST-1.GT.DOU)) THEN
708: * if the cluster contains no desired eigenvalues
709: * skip the computation of that branch of the rep. tree
710: *
711: * We could skip before the refinement of the extremal
712: * eigenvalues of the child, but then the representation
713: * tree could be different from the one when nothing is
714: * skipped. For this reason we skip at this place.
715: IDONE = IDONE + NEWLST - NEWFST + 1
716: GOTO 139
717: ENDIF
718: *
719: * Compute RRR of child cluster.
720: * Note that the new RRR is stored in Z
721: *
722: * DLARRF needs LWORK = 2*N
723: CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ),
724: $ WORK(INDLD+IBEGIN-1),
725: $ NEWFST, NEWLST, WORK(WBEGIN),
726: $ WGAP(WBEGIN), WERR(WBEGIN),
727: $ SPDIAM, LGAP, RGAP, PIVMIN, TAU,
728: $ Z(IBEGIN, NEWFTT),Z(IBEGIN, NEWFTT+1),
729: $ WORK( INDWRK ), IINFO )
730: IF( IINFO.EQ.0 ) THEN
731: * a new RRR for the cluster was found by DLARRF
732: * update shift and store it
733: SSIGMA = SIGMA + TAU
734: Z( IEND, NEWFTT+1 ) = SSIGMA
735: * WORK() are the midpoints and WERR() the semi-width
736: * Note that the entries in W are unchanged.
737: DO 116 K = NEWFST, NEWLST
738: FUDGE =
739: $ THREE*EPS*ABS(WORK(WBEGIN+K-1))
740: WORK( WBEGIN + K - 1 ) =
741: $ WORK( WBEGIN + K - 1) - TAU
742: FUDGE = FUDGE +
743: $ FOUR*EPS*ABS(WORK(WBEGIN+K-1))
744: * Fudge errors
745: WERR( WBEGIN + K - 1 ) =
746: $ WERR( WBEGIN + K - 1 ) + FUDGE
747: * Gaps are not fudged. Provided that WERR is small
748: * when eigenvalues are close, a zero gap indicates
749: * that a new representation is needed for resolving
750: * the cluster. A fudge could lead to a wrong decision
751: * of judging eigenvalues 'separated' which in
752: * reality are not. This could have a negative impact
753: * on the orthogonality of the computed eigenvectors.
754: 116 CONTINUE
755:
756: NCLUS = NCLUS + 1
757: K = NEWCLS + 2*NCLUS
758: IWORK( K-1 ) = NEWFST
759: IWORK( K ) = NEWLST
760: ELSE
761: INFO = -2
762: RETURN
763: ENDIF
764: ELSE
765: *
766: * Compute eigenvector of singleton
767: *
768: ITER = 0
769: *
770: TOL = FOUR * LOG(DBLE(IN)) * EPS
771: *
772: K = NEWFST
773: WINDEX = WBEGIN + K - 1
774: WINDMN = MAX(WINDEX - 1,1)
775: WINDPL = MIN(WINDEX + 1,M)
776: LAMBDA = WORK( WINDEX )
777: DONE = DONE + 1
778: * Check if eigenvector computation is to be skipped
779: IF((WINDEX.LT.DOL).OR.
780: $ (WINDEX.GT.DOU)) THEN
781: ESKIP = .TRUE.
782: GOTO 125
783: ELSE
784: ESKIP = .FALSE.
785: ENDIF
786: LEFT = WORK( WINDEX ) - WERR( WINDEX )
787: RIGHT = WORK( WINDEX ) + WERR( WINDEX )
788: INDEIG = INDEXW( WINDEX )
789: * Note that since we compute the eigenpairs for a child,
790: * all eigenvalue approximations are w.r.t the same shift.
791: * In this case, the entries in WORK should be used for
792: * computing the gaps since they exhibit even very small
793: * differences in the eigenvalues, as opposed to the
794: * entries in W which might "look" the same.
795:
796: IF( K .EQ. 1) THEN
797: * In the case RANGE='I' and with not much initial
798: * accuracy in LAMBDA and VL, the formula
799: * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
800: * can lead to an overestimation of the left gap and
801: * thus to inadequately early RQI 'convergence'.
802: * Prevent this by forcing a small left gap.
803: LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
804: ELSE
805: LGAP = WGAP(WINDMN)
806: ENDIF
807: IF( K .EQ. IM) THEN
808: * In the case RANGE='I' and with not much initial
809: * accuracy in LAMBDA and VU, the formula
810: * can lead to an overestimation of the right gap and
811: * thus to inadequately early RQI 'convergence'.
812: * Prevent this by forcing a small right gap.
813: RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT))
814: ELSE
815: RGAP = WGAP(WINDEX)
816: ENDIF
817: GAP = MIN( LGAP, RGAP )
818: IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN
819: * The eigenvector support can become wrong
820: * because significant entries could be cut off due to a
821: * large GAPTOL parameter in LAR1V. Prevent this.
822: GAPTOL = ZERO
823: ELSE
824: GAPTOL = GAP * EPS
825: ENDIF
826: ISUPMN = IN
827: ISUPMX = 1
828: * Update WGAP so that it holds the minimum gap
829: * to the left or the right. This is crucial in the
830: * case where bisection is used to ensure that the
831: * eigenvalue is refined up to the required precision.
832: * The correct value is restored afterwards.
833: SAVGAP = WGAP(WINDEX)
834: WGAP(WINDEX) = GAP
835: * We want to use the Rayleigh Quotient Correction
836: * as often as possible since it converges quadratically
837: * when we are close enough to the desired eigenvalue.
838: * However, the Rayleigh Quotient can have the wrong sign
839: * and lead us away from the desired eigenvalue. In this
840: * case, the best we can do is to use bisection.
841: USEDBS = .FALSE.
842: USEDRQ = .FALSE.
843: * Bisection is initially turned off unless it is forced
844: NEEDBS = .NOT.TRYRQC
845: 120 CONTINUE
846: * Check if bisection should be used to refine eigenvalue
847: IF(NEEDBS) THEN
848: * Take the bisection as new iterate
849: USEDBS = .TRUE.
850: ITMP1 = IWORK( IINDR+WINDEX )
851: OFFSET = INDEXW( WBEGIN ) - 1
852: CALL DLARRB( IN, D(IBEGIN),
853: $ WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG,
854: $ ZERO, TWO*EPS, OFFSET,
855: $ WORK(WBEGIN),WGAP(WBEGIN),
856: $ WERR(WBEGIN),WORK( INDWRK ),
857: $ IWORK( IINDWK ), PIVMIN, SPDIAM,
858: $ ITMP1, IINFO )
859: IF( IINFO.NE.0 ) THEN
860: INFO = -3
861: RETURN
862: ENDIF
863: LAMBDA = WORK( WINDEX )
864: * Reset twist index from inaccurate LAMBDA to
865: * force computation of true MINGMA
866: IWORK( IINDR+WINDEX ) = 0
867: ENDIF
868: * Given LAMBDA, compute the eigenvector.
869: CALL DLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ),
870: $ L( IBEGIN ), WORK(INDLD+IBEGIN-1),
871: $ WORK(INDLLD+IBEGIN-1),
872: $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
873: $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
874: $ IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ),
875: $ NRMINV, RESID, RQCORR, WORK( INDWRK ) )
876: IF(ITER .EQ. 0) THEN
877: BSTRES = RESID
878: BSTW = LAMBDA
879: ELSEIF(RESID.LT.BSTRES) THEN
880: BSTRES = RESID
881: BSTW = LAMBDA
882: ENDIF
883: ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 ))
884: ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX ))
885: ITER = ITER + 1
886:
887: * sin alpha <= |resid|/gap
888: * Note that both the residual and the gap are
889: * proportional to the matrix, so ||T|| doesn't play
890: * a role in the quotient
891:
892: *
893: * Convergence test for Rayleigh-Quotient iteration
894: * (omitted when Bisection has been used)
895: *
896: IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT.
897: $ RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS)
898: $ THEN
899: * We need to check that the RQCORR update doesn't
900: * move the eigenvalue away from the desired one and
901: * towards a neighbor. -> protection with bisection
902: IF(INDEIG.LE.NEGCNT) THEN
903: * The wanted eigenvalue lies to the left
904: SGNDEF = -ONE
905: ELSE
906: * The wanted eigenvalue lies to the right
907: SGNDEF = ONE
908: ENDIF
909: * We only use the RQCORR if it improves the
910: * the iterate reasonably.
911: IF( ( RQCORR*SGNDEF.GE.ZERO )
912: $ .AND.( LAMBDA + RQCORR.LE. RIGHT)
913: $ .AND.( LAMBDA + RQCORR.GE. LEFT)
914: $ ) THEN
915: USEDRQ = .TRUE.
916: * Store new midpoint of bisection interval in WORK
917: IF(SGNDEF.EQ.ONE) THEN
918: * The current LAMBDA is on the left of the true
919: * eigenvalue
920: LEFT = LAMBDA
921: * We prefer to assume that the error estimate
922: * is correct. We could make the interval not
923: * as a bracket but to be modified if the RQCORR
924: * chooses to. In this case, the RIGHT side should
925: * be modified as follows:
926: * RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
927: ELSE
928: * The current LAMBDA is on the right of the true
929: * eigenvalue
930: RIGHT = LAMBDA
931: * See comment about assuming the error estimate is
932: * correct above.
933: * LEFT = MIN(LEFT, LAMBDA + RQCORR)
934: ENDIF
935: WORK( WINDEX ) =
936: $ HALF * (RIGHT + LEFT)
937: * Take RQCORR since it has the correct sign and
938: * improves the iterate reasonably
939: LAMBDA = LAMBDA + RQCORR
940: * Update width of error interval
941: WERR( WINDEX ) =
942: $ HALF * (RIGHT-LEFT)
943: ELSE
944: NEEDBS = .TRUE.
945: ENDIF
946: IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN
947: * The eigenvalue is computed to bisection accuracy
948: * compute eigenvector and stop
949: USEDBS = .TRUE.
950: GOTO 120
951: ELSEIF( ITER.LT.MAXITR ) THEN
952: GOTO 120
953: ELSEIF( ITER.EQ.MAXITR ) THEN
954: NEEDBS = .TRUE.
955: GOTO 120
956: ELSE
957: INFO = 5
958: RETURN
959: END IF
960: ELSE
961: STP2II = .FALSE.
962: IF(USEDRQ .AND. USEDBS .AND.
963: $ BSTRES.LE.RESID) THEN
964: LAMBDA = BSTW
965: STP2II = .TRUE.
966: ENDIF
967: IF (STP2II) THEN
968: * improve error angle by second step
969: CALL DLAR1V( IN, 1, IN, LAMBDA,
970: $ D( IBEGIN ), L( IBEGIN ),
971: $ WORK(INDLD+IBEGIN-1),
972: $ WORK(INDLLD+IBEGIN-1),
973: $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ),
974: $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA,
975: $ IWORK( IINDR+WINDEX ),
976: $ ISUPPZ( 2*WINDEX-1 ),
977: $ NRMINV, RESID, RQCORR, WORK( INDWRK ) )
978: ENDIF
979: WORK( WINDEX ) = LAMBDA
980: END IF
981: *
982: * Compute FP-vector support w.r.t. whole matrix
983: *
984: ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN
985: ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN
986: ZFROM = ISUPPZ( 2*WINDEX-1 )
987: ZTO = ISUPPZ( 2*WINDEX )
988: ISUPMN = ISUPMN + OLDIEN
989: ISUPMX = ISUPMX + OLDIEN
990: * Ensure vector is ok if support in the RQI has changed
991: IF(ISUPMN.LT.ZFROM) THEN
992: DO 122 II = ISUPMN,ZFROM-1
993: Z( II, WINDEX ) = ZERO
994: 122 CONTINUE
995: ENDIF
996: IF(ISUPMX.GT.ZTO) THEN
997: DO 123 II = ZTO+1,ISUPMX
998: Z( II, WINDEX ) = ZERO
999: 123 CONTINUE
1000: ENDIF
1001: CALL DSCAL( ZTO-ZFROM+1, NRMINV,
1002: $ Z( ZFROM, WINDEX ), 1 )
1003: 125 CONTINUE
1004: * Update W
1005: W( WINDEX ) = LAMBDA+SIGMA
1006: * Recompute the gaps on the left and right
1007: * But only allow them to become larger and not
1008: * smaller (which can only happen through "bad"
1009: * cancellation and doesn't reflect the theory
1010: * where the initial gaps are underestimated due
1011: * to WERR being too crude.)
1012: IF(.NOT.ESKIP) THEN
1013: IF( K.GT.1) THEN
1014: WGAP( WINDMN ) = MAX( WGAP(WINDMN),
1015: $ W(WINDEX)-WERR(WINDEX)
1016: $ - W(WINDMN)-WERR(WINDMN) )
1017: ENDIF
1018: IF( WINDEX.LT.WEND ) THEN
1019: WGAP( WINDEX ) = MAX( SAVGAP,
1020: $ W( WINDPL )-WERR( WINDPL )
1021: $ - W( WINDEX )-WERR( WINDEX) )
1022: ENDIF
1023: ENDIF
1024: IDONE = IDONE + 1
1025: ENDIF
1026: * here ends the code for the current child
1027: *
1028: 139 CONTINUE
1029: * Proceed to any remaining child nodes
1030: NEWFST = J + 1
1031: 140 CONTINUE
1032: 150 CONTINUE
1033: NDEPTH = NDEPTH + 1
1034: GO TO 40
1035: END IF
1036: IBEGIN = IEND + 1
1037: WBEGIN = WEND + 1
1038: 170 CONTINUE
1039: *
1040:
1041: RETURN
1042: *
1043: * End of DLARRV
1044: *
1045: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack