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