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