LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlarrv.f
Go to the documentation of this file.
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 < 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*> \ingroup larrv
276*
277*> \par Contributors:
278* ==================
279*>
280*> Beresford Parlett, University of California, Berkeley, USA \n
281*> Jim Demmel, University of California, Berkeley, USA \n
282*> Inderjit Dhillon, University of Texas, Austin, USA \n
283*> Osni Marques, LBNL/NERSC, USA \n
284*> Christof Voemel, University of California, Berkeley, USA
285*
286* =====================================================================
287 SUBROUTINE dlarrv( N, VL, VU, D, L, PIVMIN,
288 $ ISPLIT, M, DOL, DOU, MINRGP,
289 $ RTOL1, RTOL2, W, WERR, WGAP,
290 $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
291 $ WORK, IWORK, INFO )
292*
293* -- LAPACK auxiliary routine --
294* -- LAPACK is a software package provided by Univ. of Tennessee, --
295* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
296*
297* .. Scalar Arguments ..
298 INTEGER DOL, DOU, INFO, LDZ, M, N
299 DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
300* ..
301* .. Array Arguments ..
302 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
303 $ isuppz( * ), iwork( * )
304 DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
305 $ WGAP( * ), WORK( * )
306 DOUBLE PRECISION Z( LDZ, * )
307* ..
308*
309* =====================================================================
310*
311* .. Parameters ..
312 INTEGER MAXITR
313 PARAMETER ( MAXITR = 10 )
314 double precision zero, one, two, three, four, half
315 parameter( zero = 0.0d0, one = 1.0d0,
316 $ two = 2.0d0, three = 3.0d0,
317 $ four = 4.0d0, half = 0.5d0)
318* ..
319* .. Local Scalars ..
320 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
321 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
322 $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
323 $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
324 $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
325 $ ndepth, negcnt, newcls, newfst, newftt, newlst,
326 $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
327 $ oldncl, p, parity, q, wbegin, wend, windex,
328 $ windmn, windpl, zfrom, zto, zusedl, zusedu,
329 $ zusedw
330 DOUBLE PRECISION BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
331 $ LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
332 $ RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
333 $ SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ
334* ..
335* .. External Functions ..
336 DOUBLE PRECISION DLAMCH
337 EXTERNAL DLAMCH
338* ..
339* .. External Subroutines ..
340 EXTERNAL dcopy, dlar1v, dlarrb, dlarrf, dlaset,
341 $ dscal
342* ..
343* .. Intrinsic Functions ..
344 INTRINSIC abs, dble, max, min
345* ..
346* .. Executable Statements ..
347* ..
348
349 info = 0
350*
351* Quick return if possible
352*
353 IF( (n.LE.0).OR.(m.LE.0) ) THEN
354 RETURN
355 END IF
356*
357* The first N entries of WORK are reserved for the eigenvalues
358 indld = n+1
359 indlld= 2*n+1
360 indwrk= 3*n+1
361 minwsize = 12 * n
362
363 DO 5 i= 1,minwsize
364 work( i ) = zero
365 5 CONTINUE
366
367* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the
368* factorization used to compute the FP vector
369 iindr = 0
370* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current
371* layer and the one above.
372 iindc1 = n
373 iindc2 = 2*n
374 iindwk = 3*n + 1
375
376 miniwsize = 7 * n
377 DO 10 i= 1,miniwsize
378 iwork( i ) = 0
379 10 CONTINUE
380
381 zusedl = 1
382 IF(dol.GT.1) THEN
383* Set lower bound for use of Z
384 zusedl = dol-1
385 ENDIF
386 zusedu = m
387 IF(dou.LT.m) THEN
388* Set lower bound for use of Z
389 zusedu = dou+1
390 ENDIF
391* The width of the part of Z that is used
392 zusedw = zusedu - zusedl + 1
393
394
395 CALL dlaset( 'Full', n, zusedw, zero, zero,
396 $ z(1,zusedl), ldz )
397
398 eps = dlamch( 'Precision' )
399 rqtol = two * eps
400*
401* Set expert flags for standard code.
402 tryrqc = .true.
403
404 IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
405 ELSE
406* Only selected eigenpairs are computed. Since the other evalues
407* are not refined by RQ iteration, bisection has to compute to full
408* accuracy.
409 rtol1 = four * eps
410 rtol2 = four * eps
411 ENDIF
412
413* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the
414* desired eigenvalues. The support of the nonzero eigenvector
415* entries is contained in the interval IBEGIN:IEND.
416* Remark that if k eigenpairs are desired, then the eigenvectors
417* are stored in k contiguous columns of Z.
418
419* DONE is the number of eigenvectors already computed
420 done = 0
421 ibegin = 1
422 wbegin = 1
423 DO 170 jblk = 1, iblock( m )
424 iend = isplit( jblk )
425 sigma = l( iend )
426* Find the eigenvectors of the submatrix indexed IBEGIN
427* through IEND.
428 wend = wbegin - 1
429 15 CONTINUE
430 IF( wend.LT.m ) THEN
431 IF( iblock( wend+1 ).EQ.jblk ) THEN
432 wend = wend + 1
433 GO TO 15
434 END IF
435 END IF
436 IF( wend.LT.wbegin ) THEN
437 ibegin = iend + 1
438 GO TO 170
439 ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) ) THEN
440 ibegin = iend + 1
441 wbegin = wend + 1
442 GO TO 170
443 END IF
444
445* Find local spectral diameter of the block
446 gl = gers( 2*ibegin-1 )
447 gu = gers( 2*ibegin )
448 DO 20 i = ibegin+1 , iend
449 gl = min( gers( 2*i-1 ), gl )
450 gu = max( gers( 2*i ), gu )
451 20 CONTINUE
452 spdiam = gu - gl
453
454* OLDIEN is the last index of the previous block
455 oldien = ibegin - 1
456* Calculate the size of the current block
457 in = iend - ibegin + 1
458* The number of eigenvalues in the current block
459 im = wend - wbegin + 1
460
461* This is for a 1x1 block
462 IF( ibegin.EQ.iend ) THEN
463 done = done+1
464 z( ibegin, wbegin ) = one
465 isuppz( 2*wbegin-1 ) = ibegin
466 isuppz( 2*wbegin ) = ibegin
467 w( wbegin ) = w( wbegin ) + sigma
468 work( wbegin ) = w( wbegin )
469 ibegin = iend + 1
470 wbegin = wbegin + 1
471 GO TO 170
472 END IF
473
474* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND)
475* Note that these can be approximations, in this case, the corresp.
476* entries of WERR give the size of the uncertainty interval.
477* The eigenvalue approximations will be refined when necessary as
478* high relative accuracy is required for the computation of the
479* corresponding eigenvectors.
480 CALL dcopy( im, w( wbegin ), 1,
481 $ work( wbegin ), 1 )
482
483* We store in W the eigenvalue approximations w.r.t. the original
484* matrix T.
485 DO 30 i=1,im
486 w(wbegin+i-1) = w(wbegin+i-1)+sigma
487 30 CONTINUE
488
489
490* NDEPTH is the current depth of the representation tree
491 ndepth = 0
492* PARITY is either 1 or 0
493 parity = 1
494* NCLUS is the number of clusters for the next level of the
495* representation tree, we start with NCLUS = 1 for the root
496 nclus = 1
497 iwork( iindc1+1 ) = 1
498 iwork( iindc1+2 ) = im
499
500* IDONE is the number of eigenvectors already computed in the current
501* block
502 idone = 0
503* loop while( IDONE.LT.IM )
504* generate the representation tree for the current block and
505* compute the eigenvectors
506 40 CONTINUE
507 IF( idone.LT.im ) THEN
508* This is a crude protection against infinitely deep trees
509 IF( ndepth.GT.m ) THEN
510 info = -2
511 RETURN
512 ENDIF
513* breadth first processing of the current level of the representation
514* tree: OLDNCL = number of clusters on current level
515 oldncl = nclus
516* reset NCLUS to count the number of child clusters
517 nclus = 0
518*
519 parity = 1 - parity
520 IF( parity.EQ.0 ) THEN
521 oldcls = iindc1
522 newcls = iindc2
523 ELSE
524 oldcls = iindc2
525 newcls = iindc1
526 END IF
527* Process the clusters on the current level
528 DO 150 i = 1, oldncl
529 j = oldcls + 2*i
530* OLDFST, OLDLST = first, last index of current cluster.
531* cluster indices start with 1 and are relative
532* to WBEGIN when accessing W, WGAP, WERR, Z
533 oldfst = iwork( j-1 )
534 oldlst = iwork( j )
535 IF( ndepth.GT.0 ) THEN
536* Retrieve relatively robust representation (RRR) of cluster
537* that has been computed at the previous level
538* The RRR is stored in Z and overwritten once the eigenvectors
539* have been computed or when the cluster is refined
540
541 IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
542* Get representation from location of the leftmost evalue
543* of the cluster
544 j = wbegin + oldfst - 1
545 ELSE
546 IF(wbegin+oldfst-1.LT.dol) THEN
547* Get representation from the left end of Z array
548 j = dol - 1
549 ELSEIF(wbegin+oldfst-1.GT.dou) THEN
550* Get representation from the right end of Z array
551 j = dou
552 ELSE
553 j = wbegin + oldfst - 1
554 ENDIF
555 ENDIF
556 CALL dcopy( in, z( ibegin, j ), 1, d( ibegin ), 1 )
557 CALL dcopy( in-1, z( ibegin, j+1 ), 1, l( ibegin ),
558 $ 1 )
559 sigma = z( iend, j+1 )
560
561* Set the corresponding entries in Z to zero
562 CALL dlaset( 'Full', in, 2, zero, zero,
563 $ z( ibegin, j), ldz )
564 END IF
565
566* Compute DL and DLL of current RRR
567 DO 50 j = ibegin, iend-1
568 tmp = d( j )*l( j )
569 work( indld-1+j ) = tmp
570 work( indlld-1+j ) = tmp*l( j )
571 50 CONTINUE
572
573 IF( ndepth.GT.0 ) THEN
574* P and Q are index of the first and last eigenvalue to compute
575* within the current block
576 p = indexw( wbegin-1+oldfst )
577 q = indexw( wbegin-1+oldlst )
578* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
579* through the Q-OFFSET elements of these arrays are to be used.
580* OFFSET = P-OLDFST
581 offset = indexw( wbegin ) - 1
582* perform limited bisection (if necessary) to get approximate
583* eigenvalues to the precision needed.
584 CALL dlarrb( in, d( ibegin ),
585 $ work(indlld+ibegin-1),
586 $ p, q, rtol1, rtol2, offset,
587 $ work(wbegin),wgap(wbegin),werr(wbegin),
588 $ work( indwrk ), iwork( iindwk ),
589 $ pivmin, spdiam, in, iinfo )
590 IF( iinfo.NE.0 ) THEN
591 info = -1
592 RETURN
593 ENDIF
594* We also recompute the extremal gaps. W holds all eigenvalues
595* of the unshifted matrix and must be used for computation
596* of WGAP, the entries of WORK might stem from RRRs with
597* different shifts. The gaps from WBEGIN-1+OLDFST to
598* WBEGIN-1+OLDLST are correctly computed in DLARRB.
599* However, we only allow the gaps to become greater since
600* this is what should happen when we decrease WERR
601 IF( oldfst.GT.1) THEN
602 wgap( wbegin+oldfst-2 ) =
603 $ max(wgap(wbegin+oldfst-2),
604 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
605 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
606 ENDIF
607 IF( wbegin + oldlst -1 .LT. wend ) THEN
608 wgap( wbegin+oldlst-1 ) =
609 $ max(wgap(wbegin+oldlst-1),
610 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
611 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
612 ENDIF
613* Each time the eigenvalues in WORK get refined, we store
614* the newly found approximation with all shifts applied in W
615 DO 53 j=oldfst,oldlst
616 w(wbegin+j-1) = work(wbegin+j-1)+sigma
617 53 CONTINUE
618 END IF
619
620* Process the current node.
621 newfst = oldfst
622 DO 140 j = oldfst, oldlst
623 IF( j.EQ.oldlst ) THEN
624* we are at the right end of the cluster, this is also the
625* boundary of the child cluster
626 newlst = j
627 ELSE IF ( wgap( wbegin + j -1).GE.
628 $ minrgp* abs( work(wbegin + j -1) ) ) THEN
629* the right relative gap is big enough, the child cluster
630* (NEWFST,..,NEWLST) is well separated from the following
631 newlst = j
632 ELSE
633* inside a child cluster, the relative gap is not
634* big enough.
635 GOTO 140
636 END IF
637
638* Compute size of child cluster found
639 newsiz = newlst - newfst + 1
640
641* NEWFTT is the place in Z where the new RRR or the computed
642* eigenvector is to be stored
643 IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
644* Store representation at location of the leftmost evalue
645* of the cluster
646 newftt = wbegin + newfst - 1
647 ELSE
648 IF(wbegin+newfst-1.LT.dol) THEN
649* Store representation at the left end of Z array
650 newftt = dol - 1
651 ELSEIF(wbegin+newfst-1.GT.dou) THEN
652* Store representation at the right end of Z array
653 newftt = dou
654 ELSE
655 newftt = wbegin + newfst - 1
656 ENDIF
657 ENDIF
658
659 IF( newsiz.GT.1) THEN
660*
661* Current child is not a singleton but a cluster.
662* Compute and store new representation of child.
663*
664*
665* Compute left and right cluster gap.
666*
667* LGAP and RGAP are not computed from WORK because
668* the eigenvalue approximations may stem from RRRs
669* different shifts. However, W hold all eigenvalues
670* of the unshifted matrix. Still, the entries in WGAP
671* have to be computed from WORK since the entries
672* in W might be of the same order so that gaps are not
673* exhibited correctly for very close eigenvalues.
674 IF( newfst.EQ.1 ) THEN
675 lgap = max( zero,
676 $ w(wbegin)-werr(wbegin) - vl )
677 ELSE
678 lgap = wgap( wbegin+newfst-2 )
679 ENDIF
680 rgap = wgap( wbegin+newlst-1 )
681*
682* Compute left- and rightmost eigenvalue of child
683* to high precision in order to shift as close
684* as possible and obtain as large relative gaps
685* as possible
686*
687 DO 55 k =1,2
688 IF(k.EQ.1) THEN
689 p = indexw( wbegin-1+newfst )
690 ELSE
691 p = indexw( wbegin-1+newlst )
692 ENDIF
693 offset = indexw( wbegin ) - 1
694 CALL dlarrb( in, d(ibegin),
695 $ work( indlld+ibegin-1 ),p,p,
696 $ rqtol, rqtol, offset,
697 $ work(wbegin),wgap(wbegin),
698 $ werr(wbegin),work( indwrk ),
699 $ iwork( iindwk ), pivmin, spdiam,
700 $ in, iinfo )
701 55 CONTINUE
702*
703 IF((wbegin+newlst-1.LT.dol).OR.
704 $ (wbegin+newfst-1.GT.dou)) THEN
705* if the cluster contains no desired eigenvalues
706* skip the computation of that branch of the rep. tree
707*
708* We could skip before the refinement of the extremal
709* eigenvalues of the child, but then the representation
710* tree could be different from the one when nothing is
711* skipped. For this reason we skip at this place.
712 idone = idone + newlst - newfst + 1
713 GOTO 139
714 ENDIF
715*
716* Compute RRR of child cluster.
717* Note that the new RRR is stored in Z
718*
719* DLARRF needs LWORK = 2*N
720 CALL dlarrf( in, d( ibegin ), l( ibegin ),
721 $ work(indld+ibegin-1),
722 $ newfst, newlst, work(wbegin),
723 $ wgap(wbegin), werr(wbegin),
724 $ spdiam, lgap, rgap, pivmin, tau,
725 $ z(ibegin, newftt),z(ibegin, newftt+1),
726 $ work( indwrk ), iinfo )
727 IF( iinfo.EQ.0 ) THEN
728* a new RRR for the cluster was found by DLARRF
729* update shift and store it
730 ssigma = sigma + tau
731 z( iend, newftt+1 ) = ssigma
732* WORK() are the midpoints and WERR() the semi-width
733* Note that the entries in W are unchanged.
734 DO 116 k = newfst, newlst
735 fudge =
736 $ three*eps*abs(work(wbegin+k-1))
737 work( wbegin + k - 1 ) =
738 $ work( wbegin + k - 1) - tau
739 fudge = fudge +
740 $ four*eps*abs(work(wbegin+k-1))
741* Fudge errors
742 werr( wbegin + k - 1 ) =
743 $ werr( wbegin + k - 1 ) + fudge
744* Gaps are not fudged. Provided that WERR is small
745* when eigenvalues are close, a zero gap indicates
746* that a new representation is needed for resolving
747* the cluster. A fudge could lead to a wrong decision
748* of judging eigenvalues 'separated' which in
749* reality are not. This could have a negative impact
750* on the orthogonality of the computed eigenvectors.
751 116 CONTINUE
752
753 nclus = nclus + 1
754 k = newcls + 2*nclus
755 iwork( k-1 ) = newfst
756 iwork( k ) = newlst
757 ELSE
758 info = -2
759 RETURN
760 ENDIF
761 ELSE
762*
763* Compute eigenvector of singleton
764*
765 iter = 0
766*
767 tol = four * log(dble(in)) * eps
768*
769 k = newfst
770 windex = wbegin + k - 1
771 windmn = max(windex - 1,1)
772 windpl = min(windex + 1,m)
773 lambda = work( windex )
774 done = done + 1
775* Check if eigenvector computation is to be skipped
776 IF((windex.LT.dol).OR.
777 $ (windex.GT.dou)) THEN
778 eskip = .true.
779 GOTO 125
780 ELSE
781 eskip = .false.
782 ENDIF
783 left = work( windex ) - werr( windex )
784 right = work( windex ) + werr( windex )
785 indeig = indexw( windex )
786* Note that since we compute the eigenpairs for a child,
787* all eigenvalue approximations are w.r.t the same shift.
788* In this case, the entries in WORK should be used for
789* computing the gaps since they exhibit even very small
790* differences in the eigenvalues, as opposed to the
791* entries in W which might "look" the same.
792
793 IF( k .EQ. 1) THEN
794* In the case RANGE='I' and with not much initial
795* accuracy in LAMBDA and VL, the formula
796* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
797* can lead to an overestimation of the left gap and
798* thus to inadequately early RQI 'convergence'.
799* Prevent this by forcing a small left gap.
800 lgap = eps*max(abs(left),abs(right))
801 ELSE
802 lgap = wgap(windmn)
803 ENDIF
804 IF( k .EQ. im) THEN
805* In the case RANGE='I' and with not much initial
806* accuracy in LAMBDA and VU, the formula
807* can lead to an overestimation of the right gap and
808* thus to inadequately early RQI 'convergence'.
809* Prevent this by forcing a small right gap.
810 rgap = eps*max(abs(left),abs(right))
811 ELSE
812 rgap = wgap(windex)
813 ENDIF
814 gap = min( lgap, rgap )
815 IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
816* The eigenvector support can become wrong
817* because significant entries could be cut off due to a
818* large GAPTOL parameter in LAR1V. Prevent this.
819 gaptol = zero
820 ELSE
821 gaptol = gap * eps
822 ENDIF
823 isupmn = in
824 isupmx = 1
825* Update WGAP so that it holds the minimum gap
826* to the left or the right. This is crucial in the
827* case where bisection is used to ensure that the
828* eigenvalue is refined up to the required precision.
829* The correct value is restored afterwards.
830 savgap = wgap(windex)
831 wgap(windex) = gap
832* We want to use the Rayleigh Quotient Correction
833* as often as possible since it converges quadratically
834* when we are close enough to the desired eigenvalue.
835* However, the Rayleigh Quotient can have the wrong sign
836* and lead us away from the desired eigenvalue. In this
837* case, the best we can do is to use bisection.
838 usedbs = .false.
839 usedrq = .false.
840* Bisection is initially turned off unless it is forced
841 needbs = .NOT.tryrqc
842 120 CONTINUE
843* Check if bisection should be used to refine eigenvalue
844 IF(needbs) THEN
845* Take the bisection as new iterate
846 usedbs = .true.
847 itmp1 = iwork( iindr+windex )
848 offset = indexw( wbegin ) - 1
849 CALL dlarrb( in, d(ibegin),
850 $ work(indlld+ibegin-1),indeig,indeig,
851 $ zero, two*eps, offset,
852 $ work(wbegin),wgap(wbegin),
853 $ werr(wbegin),work( indwrk ),
854 $ iwork( iindwk ), pivmin, spdiam,
855 $ itmp1, iinfo )
856 IF( iinfo.NE.0 ) THEN
857 info = -3
858 RETURN
859 ENDIF
860 lambda = work( windex )
861* Reset twist index from inaccurate LAMBDA to
862* force computation of true MINGMA
863 iwork( iindr+windex ) = 0
864 ENDIF
865* Given LAMBDA, compute the eigenvector.
866 CALL dlar1v( in, 1, in, lambda, d( ibegin ),
867 $ l( ibegin ), work(indld+ibegin-1),
868 $ work(indlld+ibegin-1),
869 $ pivmin, gaptol, z( ibegin, windex ),
870 $ .NOT.usedbs, negcnt, ztz, mingma,
871 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
872 $ nrminv, resid, rqcorr, work( indwrk ) )
873 IF(iter .EQ. 0) THEN
874 bstres = resid
875 bstw = lambda
876 ELSEIF(resid.LT.bstres) THEN
877 bstres = resid
878 bstw = lambda
879 ENDIF
880 isupmn = min(isupmn,isuppz( 2*windex-1 ))
881 isupmx = max(isupmx,isuppz( 2*windex ))
882 iter = iter + 1
883
884* sin alpha <= |resid|/gap
885* Note that both the residual and the gap are
886* proportional to the matrix, so ||T|| doesn't play
887* a role in the quotient
888
889*
890* Convergence test for Rayleigh-Quotient iteration
891* (omitted when Bisection has been used)
892*
893 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
894 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
895 $ THEN
896* We need to check that the RQCORR update doesn't
897* move the eigenvalue away from the desired one and
898* towards a neighbor. -> protection with bisection
899 IF(indeig.LE.negcnt) THEN
900* The wanted eigenvalue lies to the left
901 sgndef = -one
902 ELSE
903* The wanted eigenvalue lies to the right
904 sgndef = one
905 ENDIF
906* We only use the RQCORR if it improves the
907* the iterate reasonably.
908 IF( ( rqcorr*sgndef.GE.zero )
909 $ .AND.( lambda + rqcorr.LE. right)
910 $ .AND.( lambda + rqcorr.GE. left)
911 $ ) THEN
912 usedrq = .true.
913* Store new midpoint of bisection interval in WORK
914 IF(sgndef.EQ.one) THEN
915* The current LAMBDA is on the left of the true
916* eigenvalue
917 left = lambda
918* We prefer to assume that the error estimate
919* is correct. We could make the interval not
920* as a bracket but to be modified if the RQCORR
921* chooses to. In this case, the RIGHT side should
922* be modified as follows:
923* RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
924 ELSE
925* The current LAMBDA is on the right of the true
926* eigenvalue
927 right = lambda
928* See comment about assuming the error estimate is
929* correct above.
930* LEFT = MIN(LEFT, LAMBDA + RQCORR)
931 ENDIF
932 work( windex ) =
933 $ half * (right + left)
934* Take RQCORR since it has the correct sign and
935* improves the iterate reasonably
936 lambda = lambda + rqcorr
937* Update width of error interval
938 werr( windex ) =
939 $ half * (right-left)
940 ELSE
941 needbs = .true.
942 ENDIF
943 IF(right-left.LT.rqtol*abs(lambda)) THEN
944* The eigenvalue is computed to bisection accuracy
945* compute eigenvector and stop
946 usedbs = .true.
947 GOTO 120
948 ELSEIF( iter.LT.maxitr ) THEN
949 GOTO 120
950 ELSEIF( iter.EQ.maxitr ) THEN
951 needbs = .true.
952 GOTO 120
953 ELSE
954 info = 5
955 RETURN
956 END IF
957 ELSE
958 stp2ii = .false.
959 IF(usedrq .AND. usedbs .AND.
960 $ bstres.LE.resid) THEN
961 lambda = bstw
962 stp2ii = .true.
963 ENDIF
964 IF (stp2ii) THEN
965* improve error angle by second step
966 CALL dlar1v( in, 1, in, lambda,
967 $ d( ibegin ), l( ibegin ),
968 $ work(indld+ibegin-1),
969 $ work(indlld+ibegin-1),
970 $ pivmin, gaptol, z( ibegin, windex ),
971 $ .NOT.usedbs, negcnt, ztz, mingma,
972 $ iwork( iindr+windex ),
973 $ isuppz( 2*windex-1 ),
974 $ nrminv, resid, rqcorr, work( indwrk ) )
975 ENDIF
976 work( windex ) = lambda
977 END IF
978*
979* Compute FP-vector support w.r.t. whole matrix
980*
981 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
982 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
983 zfrom = isuppz( 2*windex-1 )
984 zto = isuppz( 2*windex )
985 isupmn = isupmn + oldien
986 isupmx = isupmx + oldien
987* Ensure vector is ok if support in the RQI has changed
988 IF(isupmn.LT.zfrom) THEN
989 DO 122 ii = isupmn,zfrom-1
990 z( ii, windex ) = zero
991 122 CONTINUE
992 ENDIF
993 IF(isupmx.GT.zto) THEN
994 DO 123 ii = zto+1,isupmx
995 z( ii, windex ) = zero
996 123 CONTINUE
997 ENDIF
998 CALL dscal( zto-zfrom+1, nrminv,
999 $ z( zfrom, windex ), 1 )
1000 125 CONTINUE
1001* Update W
1002 w( windex ) = lambda+sigma
1003* Recompute the gaps on the left and right
1004* But only allow them to become larger and not
1005* smaller (which can only happen through "bad"
1006* cancellation and doesn't reflect the theory
1007* where the initial gaps are underestimated due
1008* to WERR being too crude.)
1009 IF(.NOT.eskip) THEN
1010 IF( k.GT.1) THEN
1011 wgap( windmn ) = max( wgap(windmn),
1012 $ w(windex)-werr(windex)
1013 $ - w(windmn)-werr(windmn) )
1014 ENDIF
1015 IF( windex.LT.wend ) THEN
1016 wgap( windex ) = max( savgap,
1017 $ w( windpl )-werr( windpl )
1018 $ - w( windex )-werr( windex) )
1019 ENDIF
1020 ENDIF
1021 idone = idone + 1
1022 ENDIF
1023* here ends the code for the current child
1024*
1025 139 CONTINUE
1026* Proceed to any remaining child nodes
1027 newfst = j + 1
1028 140 CONTINUE
1029 150 CONTINUE
1030 ndepth = ndepth + 1
1031 GO TO 40
1032 END IF
1033 ibegin = iend + 1
1034 wbegin = wend + 1
1035 170 CONTINUE
1036*
1037
1038 RETURN
1039*
1040* End of DLARRV
1041*
1042 END
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlar1v
subroutine dlar1v(n, b1, bn, lambda, d, l, ld, lld, pivmin, gaptol, z, wantnc, negcnt, ztz, mingma, r, isuppz, nrminv, resid, rqcorr, work)
DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
Definition dlar1v.f:230
dlarrb
subroutine dlarrb(n, d, lld, ifirst, ilast, rtol1, rtol2, offset, w, wgap, werr, work, iwork, pivmin, spdiam, twist, info)
DLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition dlarrb.f:196
dlarrf
subroutine dlarrf(n, d, l, ld, clstrt, clend, w, wgap, werr, spdiam, clgapl, clgapr, pivmin, sigma, dplus, lplus, work, info)
DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is rela...
Definition dlarrf.f:193
dlarrv
subroutine dlarrv(n, vl, vu, d, l, pivmin, isplit, m, dol, dou, minrgp, rtol1, rtol2, w, werr, wgap, iblock, indexw, gers, z, ldz, isuppz, work, iwork, info)
DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition dlarrv.f:292
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79