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