LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clarrv.f
Go to the documentation of this file.
1*> \brief \b CLARRV 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 CLARRV + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarrv.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarrv.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarrv.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLARRV( 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* REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
28* ..
29* .. Array Arguments ..
30* INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
31* $ ISUPPZ( * ), IWORK( * )
32* REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
33* $ WGAP( * ), WORK( * )
34* COMPLEX Z( LDZ, * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> CLARRV 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 SLARRE.
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 REAL
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 REAL
68*> Upper bound of the interval that contains the desired
69*> eigenvalues. VL < VU. Needed to compute gaps on the left or right
70*> end of the extremal eigenvalues in the desired RANGE.
71*> \endverbatim
72*>
73*> \param[in,out] D
74*> \verbatim
75*> D is REAL array, dimension (N)
76*> On entry, the N diagonal elements of the diagonal matrix D.
77*> On exit, D may be overwritten.
78*> \endverbatim
79*>
80*> \param[in,out] L
81*> \verbatim
82*> L is REAL array, dimension (N)
83*> On entry, the (N-1) subdiagonal elements of the unit
84*> bidiagonal matrix L are in elements 1 to N-1 of L
85*> (if the matrix is not split.) At the end of each block
86*> is stored the corresponding shift as given by SLARRE.
87*> On exit, L is overwritten.
88*> \endverbatim
89*>
90*> \param[in] PIVMIN
91*> \verbatim
92*> PIVMIN is REAL
93*> The minimum pivot allowed in the Sturm sequence.
94*> \endverbatim
95*>
96*> \param[in] ISPLIT
97*> \verbatim
98*> ISPLIT is INTEGER array, dimension (N)
99*> The splitting points, at which T breaks up into blocks.
100*> The first block consists of rows/columns 1 to
101*> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
102*> through ISPLIT( 2 ), etc.
103*> \endverbatim
104*>
105*> \param[in] M
106*> \verbatim
107*> M is INTEGER
108*> The total number of input eigenvalues. 0 <= M <= N.
109*> \endverbatim
110*>
111*> \param[in] DOL
112*> \verbatim
113*> DOL is INTEGER
114*> \endverbatim
115*>
116*> \param[in] DOU
117*> \verbatim
118*> DOU is INTEGER
119*> If the user wants to compute only selected eigenvectors from all
120*> the eigenvalues supplied, he can specify an index range DOL:DOU.
121*> Or else the setting DOL=1, DOU=M should be applied.
122*> Note that DOL and DOU refer to the order in which the eigenvalues
123*> are stored in W.
124*> If the user wants to compute only selected eigenpairs, then
125*> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
126*> computed eigenvectors. All other columns of Z are set to zero.
127*> \endverbatim
128*>
129*> \param[in] MINRGP
130*> \verbatim
131*> MINRGP is REAL
132*> \endverbatim
133*>
134*> \param[in] RTOL1
135*> \verbatim
136*> RTOL1 is REAL
137*> \endverbatim
138*>
139*> \param[in] RTOL2
140*> \verbatim
141*> RTOL2 is REAL
142*> Parameters for bisection.
143*> An interval [LEFT,RIGHT] has converged if
144*> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
145*> \endverbatim
146*>
147*> \param[in,out] W
148*> \verbatim
149*> W is REAL array, dimension (N)
150*> The first M elements of W contain the APPROXIMATE eigenvalues for
151*> which eigenvectors are to be computed. The eigenvalues
152*> should be grouped by split-off block and ordered from
153*> smallest to largest within the block ( The output array
154*> W from SLARRE is expected here ). Furthermore, they are with
155*> respect to the shift of the corresponding root representation
156*> for their block. On exit, W holds the eigenvalues of the
157*> UNshifted matrix.
158*> \endverbatim
159*>
160*> \param[in,out] WERR
161*> \verbatim
162*> WERR is REAL array, dimension (N)
163*> The first M elements contain the semiwidth of the uncertainty
164*> interval of the corresponding eigenvalue in W
165*> \endverbatim
166*>
167*> \param[in,out] WGAP
168*> \verbatim
169*> WGAP is REAL array, dimension (N)
170*> The separation from the right neighbor eigenvalue in W.
171*> \endverbatim
172*>
173*> \param[in] IBLOCK
174*> \verbatim
175*> IBLOCK is INTEGER array, dimension (N)
176*> The indices of the blocks (submatrices) associated with the
177*> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
178*> W(i) belongs to the first block from the top, =2 if W(i)
179*> belongs to the second block, etc.
180*> \endverbatim
181*>
182*> \param[in] INDEXW
183*> \verbatim
184*> INDEXW is INTEGER array, dimension (N)
185*> The indices of the eigenvalues within each block (submatrix);
186*> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
187*> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
188*> \endverbatim
189*>
190*> \param[in] GERS
191*> \verbatim
192*> GERS is REAL array, dimension (2*N)
193*> The N Gerschgorin intervals (the i-th Gerschgorin interval
194*> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
195*> be computed from the original UNshifted matrix.
196*> \endverbatim
197*>
198*> \param[out] Z
199*> \verbatim
200*> Z is COMPLEX array, dimension (LDZ, max(1,M) )
201*> If INFO = 0, the first M columns of Z contain the
202*> orthonormal eigenvectors of the matrix T
203*> corresponding to the input eigenvalues, with the i-th
204*> column of Z holding the eigenvector associated with W(i).
205*> Note: the user must ensure that at least max(1,M) columns are
206*> supplied in the array Z.
207*> \endverbatim
208*>
209*> \param[in] LDZ
210*> \verbatim
211*> LDZ is INTEGER
212*> The leading dimension of the array Z. LDZ >= 1, and if
213*> JOBZ = 'V', LDZ >= max(1,N).
214*> \endverbatim
215*>
216*> \param[out] ISUPPZ
217*> \verbatim
218*> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
219*> The support of the eigenvectors in Z, i.e., the indices
220*> indicating the nonzero elements in Z. The I-th eigenvector
221*> is nonzero only in elements ISUPPZ( 2*I-1 ) through
222*> ISUPPZ( 2*I ).
223*> \endverbatim
224*>
225*> \param[out] WORK
226*> \verbatim
227*> WORK is REAL array, dimension (12*N)
228*> \endverbatim
229*>
230*> \param[out] IWORK
231*> \verbatim
232*> IWORK is INTEGER array, dimension (7*N)
233*> \endverbatim
234*>
235*> \param[out] INFO
236*> \verbatim
237*> INFO is INTEGER
238*> = 0: successful exit
239*>
240*> > 0: A problem occurred in CLARRV.
241*> < 0: One of the called subroutines signaled an internal problem.
242*> Needs inspection of the corresponding parameter IINFO
243*> for further information.
244*>
245*> =-1: Problem in SLARRB when refining a child's eigenvalues.
246*> =-2: Problem in SLARRF when computing the RRR of a child.
247*> When a child is inside a tight cluster, it can be difficult
248*> to find an RRR. A partial remedy from the user's point of
249*> view is to make the parameter MINRGP smaller and recompile.
250*> However, as the orthogonality of the computed vectors is
251*> proportional to 1/MINRGP, the user should be aware that
252*> he might be trading in precision when he decreases MINRGP.
253*> =-3: Problem in SLARRB when refining a single eigenvalue
254*> after the Rayleigh correction was rejected.
255*> = 5: The Rayleigh Quotient Iteration failed to converge to
256*> full accuracy in MAXITR steps.
257*> \endverbatim
258*
259* Authors:
260* ========
261*
262*> \author Univ. of Tennessee
263*> \author Univ. of California Berkeley
264*> \author Univ. of Colorado Denver
265*> \author NAG Ltd.
266*
267*> \ingroup larrv
268*
269*> \par Contributors:
270* ==================
271*>
272*> Beresford Parlett, University of California, Berkeley, USA \n
273*> Jim Demmel, University of California, Berkeley, USA \n
274*> Inderjit Dhillon, University of Texas, Austin, USA \n
275*> Osni Marques, LBNL/NERSC, USA \n
276*> Christof Voemel, University of California, Berkeley, USA
277*
278* =====================================================================
279 SUBROUTINE clarrv( N, VL, VU, D, L, PIVMIN,
280 $ ISPLIT, M, DOL, DOU, MINRGP,
281 $ RTOL1, RTOL2, W, WERR, WGAP,
282 $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ,
283 $ WORK, IWORK, INFO )
284*
285* -- LAPACK auxiliary routine --
286* -- LAPACK is a software package provided by Univ. of Tennessee, --
287* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288*
289* .. Scalar Arguments ..
290 INTEGER DOL, DOU, INFO, LDZ, M, N
291 REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
292* ..
293* .. Array Arguments ..
294 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
295 $ isuppz( * ), iwork( * )
296 REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
297 $ WGAP( * ), WORK( * )
298 COMPLEX Z( LDZ, * )
299* ..
300*
301* =====================================================================
302*
303* .. Parameters ..
304 INTEGER MAXITR
305 PARAMETER ( MAXITR = 10 )
306 complex czero
307 parameter( czero = ( 0.0e0, 0.0e0 ) )
308 REAL ZERO, ONE, TWO, THREE, FOUR, HALF
309 parameter( zero = 0.0e0, one = 1.0e0,
310 $ two = 2.0e0, three = 3.0e0,
311 $ four = 4.0e0, half = 0.5e0)
312* ..
313* .. Local Scalars ..
314 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
315 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
316 $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG,
317 $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER,
318 $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
319 $ ndepth, negcnt, newcls, newfst, newftt, newlst,
320 $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
321 $ oldncl, p, parity, q, wbegin, wend, windex,
322 $ windmn, windpl, zfrom, zto, zusedl, zusedu,
323 $ zusedw
324 INTEGER INDIN1, INDIN2
325 REAL BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
326 $ LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID,
327 $ RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF,
328 $ sigma, spdiam, ssigma, tau, tmp, tol, ztz
329* ..
330* .. External Functions ..
331 REAL SLAMCH
332 EXTERNAL SLAMCH
333* ..
334* .. External Subroutines ..
335 EXTERNAL clar1v, claset, csscal, scopy,
336 $ slarrb,
337 $ slarrf
338* ..
339* .. Intrinsic Functions ..
340 INTRINSIC abs, real, max, min
341 INTRINSIC cmplx
342* ..
343* .. Executable Statements ..
344* ..
345
346 info = 0
347*
348* Quick return if possible
349*
350 IF( (n.LE.0).OR.(m.LE.0) ) THEN
351 RETURN
352 END IF
353*
354* The first N entries of WORK are reserved for the eigenvalues
355 indld = n+1
356 indlld= 2*n+1
357 indin1 = 3*n + 1
358 indin2 = 4*n + 1
359 indwrk = 5*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 claset( 'Full', n, zusedw, czero, czero,
395 $ z(1,zusedl), ldz )
396
397 eps = slamch( '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 ) = cmplx( one, zero )
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 scopy( 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 DO 45 k = 1, in - 1
556 d( ibegin+k-1 ) = real( z( ibegin+k-1,
557 $ j ) )
558 l( ibegin+k-1 ) = real( z( ibegin+k-1,
559 $ j+1 ) )
560 45 CONTINUE
561 d( iend ) = real( z( iend, j ) )
562 sigma = real( z( iend, j+1 ) )
563
564* Set the corresponding entries in Z to zero
565 CALL claset( 'Full', in, 2, czero, czero,
566 $ z( ibegin, j), ldz )
567 END IF
568
569* Compute DL and DLL of current RRR
570 DO 50 j = ibegin, iend-1
571 tmp = d( j )*l( j )
572 work( indld-1+j ) = tmp
573 work( indlld-1+j ) = tmp*l( j )
574 50 CONTINUE
575
576 IF( ndepth.GT.0 ) THEN
577* P and Q are index of the first and last eigenvalue to compute
578* within the current block
579 p = indexw( wbegin-1+oldfst )
580 q = indexw( wbegin-1+oldlst )
581* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET
582* through the Q-OFFSET elements of these arrays are to be used.
583* OFFSET = P-OLDFST
584 offset = indexw( wbegin ) - 1
585* perform limited bisection (if necessary) to get approximate
586* eigenvalues to the precision needed.
587 CALL slarrb( in, d( ibegin ),
588 $ work(indlld+ibegin-1),
589 $ p, q, rtol1, rtol2, offset,
590 $ work(wbegin),wgap(wbegin),werr(wbegin),
591 $ work( indwrk ), iwork( iindwk ),
592 $ pivmin, spdiam, in, iinfo )
593 IF( iinfo.NE.0 ) THEN
594 info = -1
595 RETURN
596 ENDIF
597* We also recompute the extremal gaps. W holds all eigenvalues
598* of the unshifted matrix and must be used for computation
599* of WGAP, the entries of WORK might stem from RRRs with
600* different shifts. The gaps from WBEGIN-1+OLDFST to
601* WBEGIN-1+OLDLST are correctly computed in SLARRB.
602* However, we only allow the gaps to become greater since
603* this is what should happen when we decrease WERR
604 IF( oldfst.GT.1) THEN
605 wgap( wbegin+oldfst-2 ) =
606 $ max(wgap(wbegin+oldfst-2),
607 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
608 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
609 ENDIF
610 IF( wbegin + oldlst -1 .LT. wend ) THEN
611 wgap( wbegin+oldlst-1 ) =
612 $ max(wgap(wbegin+oldlst-1),
613 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
614 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
615 ENDIF
616* Each time the eigenvalues in WORK get refined, we store
617* the newly found approximation with all shifts applied in W
618 DO 53 j=oldfst,oldlst
619 w(wbegin+j-1) = work(wbegin+j-1)+sigma
620 53 CONTINUE
621 END IF
622
623* Process the current node.
624 newfst = oldfst
625 DO 140 j = oldfst, oldlst
626 IF( j.EQ.oldlst ) THEN
627* we are at the right end of the cluster, this is also the
628* boundary of the child cluster
629 newlst = j
630 ELSE IF ( wgap( wbegin + j -1).GE.
631 $ minrgp* abs( work(wbegin + j -1) ) ) THEN
632* the right relative gap is big enough, the child cluster
633* (NEWFST,..,NEWLST) is well separated from the following
634 newlst = j
635 ELSE
636* inside a child cluster, the relative gap is not
637* big enough.
638 GOTO 140
639 END IF
640
641* Compute size of child cluster found
642 newsiz = newlst - newfst + 1
643
644* NEWFTT is the place in Z where the new RRR or the computed
645* eigenvector is to be stored
646 IF((dol.EQ.1).AND.(dou.EQ.m)) THEN
647* Store representation at location of the leftmost evalue
648* of the cluster
649 newftt = wbegin + newfst - 1
650 ELSE
651 IF(wbegin+newfst-1.LT.dol) THEN
652* Store representation at the left end of Z array
653 newftt = dol - 1
654 ELSEIF(wbegin+newfst-1.GT.dou) THEN
655* Store representation at the right end of Z array
656 newftt = dou
657 ELSE
658 newftt = wbegin + newfst - 1
659 ENDIF
660 ENDIF
661
662 IF( newsiz.GT.1) THEN
663*
664* Current child is not a singleton but a cluster.
665* Compute and store new representation of child.
666*
667*
668* Compute left and right cluster gap.
669*
670* LGAP and RGAP are not computed from WORK because
671* the eigenvalue approximations may stem from RRRs
672* different shifts. However, W hold all eigenvalues
673* of the unshifted matrix. Still, the entries in WGAP
674* have to be computed from WORK since the entries
675* in W might be of the same order so that gaps are not
676* exhibited correctly for very close eigenvalues.
677 IF( newfst.EQ.1 ) THEN
678 lgap = max( zero,
679 $ w(wbegin)-werr(wbegin) - vl )
680 ELSE
681 lgap = wgap( wbegin+newfst-2 )
682 ENDIF
683 rgap = wgap( wbegin+newlst-1 )
684*
685* Compute left- and rightmost eigenvalue of child
686* to high precision in order to shift as close
687* as possible and obtain as large relative gaps
688* as possible
689*
690 DO 55 k =1,2
691 IF(k.EQ.1) THEN
692 p = indexw( wbegin-1+newfst )
693 ELSE
694 p = indexw( wbegin-1+newlst )
695 ENDIF
696 offset = indexw( wbegin ) - 1
697 CALL slarrb( in, d(ibegin),
698 $ work( indlld+ibegin-1 ),p,p,
699 $ rqtol, rqtol, offset,
700 $ work(wbegin),wgap(wbegin),
701 $ werr(wbegin),work( indwrk ),
702 $ iwork( iindwk ), pivmin, spdiam,
703 $ in, iinfo )
704 55 CONTINUE
705*
706 IF((wbegin+newlst-1.LT.dol).OR.
707 $ (wbegin+newfst-1.GT.dou)) THEN
708* if the cluster contains no desired eigenvalues
709* skip the computation of that branch of the rep. tree
710*
711* We could skip before the refinement of the extremal
712* eigenvalues of the child, but then the representation
713* tree could be different from the one when nothing is
714* skipped. For this reason we skip at this place.
715 idone = idone + newlst - newfst + 1
716 GOTO 139
717 ENDIF
718*
719* Compute RRR of child cluster.
720* Note that the new RRR is stored in Z
721*
722* SLARRF needs LWORK = 2*N
723 CALL slarrf( in, d( ibegin ), l( ibegin ),
724 $ work(indld+ibegin-1),
725 $ newfst, newlst, work(wbegin),
726 $ wgap(wbegin), werr(wbegin),
727 $ spdiam, lgap, rgap, pivmin, tau,
728 $ work( indin1 ), work( indin2 ),
729 $ work( indwrk ), iinfo )
730* In the complex case, SLARRF cannot write
731* the new RRR directly into Z and needs an intermediate
732* workspace
733 DO 56 k = 1, in-1
734 z( ibegin+k-1, newftt ) =
735 $ cmplx( work( indin1+k-1 ), zero )
736 z( ibegin+k-1, newftt+1 ) =
737 $ cmplx( work( indin2+k-1 ), zero )
738 56 CONTINUE
739 z( iend, newftt ) =
740 $ cmplx( work( indin1+in-1 ), zero )
741 IF( iinfo.EQ.0 ) THEN
742* a new RRR for the cluster was found by SLARRF
743* update shift and store it
744 ssigma = sigma + tau
745 z( iend, newftt+1 ) = cmplx( ssigma, zero )
746* WORK() are the midpoints and WERR() the semi-width
747* Note that the entries in W are unchanged.
748 DO 116 k = newfst, newlst
749 fudge =
750 $ three*eps*abs(work(wbegin+k-1))
751 work( wbegin + k - 1 ) =
752 $ work( wbegin + k - 1) - tau
753 fudge = fudge +
754 $ four*eps*abs(work(wbegin+k-1))
755* Fudge errors
756 werr( wbegin + k - 1 ) =
757 $ werr( wbegin + k - 1 ) + fudge
758* Gaps are not fudged. Provided that WERR is small
759* when eigenvalues are close, a zero gap indicates
760* that a new representation is needed for resolving
761* the cluster. A fudge could lead to a wrong decision
762* of judging eigenvalues 'separated' which in
763* reality are not. This could have a negative impact
764* on the orthogonality of the computed eigenvectors.
765 116 CONTINUE
766
767 nclus = nclus + 1
768 k = newcls + 2*nclus
769 iwork( k-1 ) = newfst
770 iwork( k ) = newlst
771 ELSE
772 info = -2
773 RETURN
774 ENDIF
775 ELSE
776*
777* Compute eigenvector of singleton
778*
779 iter = 0
780*
781 tol = four * log(real(in)) * eps
782*
783 k = newfst
784 windex = wbegin + k - 1
785 windmn = max(windex - 1,1)
786 windpl = min(windex + 1,m)
787 lambda = work( windex )
788 done = done + 1
789* Check if eigenvector computation is to be skipped
790 IF((windex.LT.dol).OR.
791 $ (windex.GT.dou)) THEN
792 eskip = .true.
793 GOTO 125
794 ELSE
795 eskip = .false.
796 ENDIF
797 left = work( windex ) - werr( windex )
798 right = work( windex ) + werr( windex )
799 indeig = indexw( windex )
800* Note that since we compute the eigenpairs for a child,
801* all eigenvalue approximations are w.r.t the same shift.
802* In this case, the entries in WORK should be used for
803* computing the gaps since they exhibit even very small
804* differences in the eigenvalues, as opposed to the
805* entries in W which might "look" the same.
806
807 IF( k .EQ. 1) THEN
808* In the case RANGE='I' and with not much initial
809* accuracy in LAMBDA and VL, the formula
810* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA )
811* can lead to an overestimation of the left gap and
812* thus to inadequately early RQI 'convergence'.
813* Prevent this by forcing a small left gap.
814 lgap = eps*max(abs(left),abs(right))
815 ELSE
816 lgap = wgap(windmn)
817 ENDIF
818 IF( k .EQ. im) THEN
819* In the case RANGE='I' and with not much initial
820* accuracy in LAMBDA and VU, the formula
821* can lead to an overestimation of the right gap and
822* thus to inadequately early RQI 'convergence'.
823* Prevent this by forcing a small right gap.
824 rgap = eps*max(abs(left),abs(right))
825 ELSE
826 rgap = wgap(windex)
827 ENDIF
828 gap = min( lgap, rgap )
829 IF(( k .EQ. 1).OR.(k .EQ. im)) THEN
830* The eigenvector support can become wrong
831* because significant entries could be cut off due to a
832* large GAPTOL parameter in LAR1V. Prevent this.
833 gaptol = zero
834 ELSE
835 gaptol = gap * eps
836 ENDIF
837 isupmn = in
838 isupmx = 1
839* Update WGAP so that it holds the minimum gap
840* to the left or the right. This is crucial in the
841* case where bisection is used to ensure that the
842* eigenvalue is refined up to the required precision.
843* The correct value is restored afterwards.
844 savgap = wgap(windex)
845 wgap(windex) = gap
846* We want to use the Rayleigh Quotient Correction
847* as often as possible since it converges quadratically
848* when we are close enough to the desired eigenvalue.
849* However, the Rayleigh Quotient can have the wrong sign
850* and lead us away from the desired eigenvalue. In this
851* case, the best we can do is to use bisection.
852 usedbs = .false.
853 usedrq = .false.
854* Bisection is initially turned off unless it is forced
855 needbs = .NOT.tryrqc
856 120 CONTINUE
857* Check if bisection should be used to refine eigenvalue
858 IF(needbs) THEN
859* Take the bisection as new iterate
860 usedbs = .true.
861 itmp1 = iwork( iindr+windex )
862 offset = indexw( wbegin ) - 1
863 CALL slarrb( in, d(ibegin),
864 $ work(indlld+ibegin-1),indeig,indeig,
865 $ zero, two*eps, offset,
866 $ work(wbegin),wgap(wbegin),
867 $ werr(wbegin),work( indwrk ),
868 $ iwork( iindwk ), pivmin, spdiam,
869 $ itmp1, iinfo )
870 IF( iinfo.NE.0 ) THEN
871 info = -3
872 RETURN
873 ENDIF
874 lambda = work( windex )
875* Reset twist index from inaccurate LAMBDA to
876* force computation of true MINGMA
877 iwork( iindr+windex ) = 0
878 ENDIF
879* Given LAMBDA, compute the eigenvector.
880 CALL clar1v( in, 1, in, lambda, d( ibegin ),
881 $ l( ibegin ), work(indld+ibegin-1),
882 $ work(indlld+ibegin-1),
883 $ pivmin, gaptol, z( ibegin, windex ),
884 $ .NOT.usedbs, negcnt, ztz, mingma,
885 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
886 $ nrminv, resid, rqcorr, work( indwrk ) )
887 IF(iter .EQ. 0) THEN
888 bstres = resid
889 bstw = lambda
890 ELSEIF(resid.LT.bstres) THEN
891 bstres = resid
892 bstw = lambda
893 ENDIF
894 isupmn = min(isupmn,isuppz( 2*windex-1 ))
895 isupmx = max(isupmx,isuppz( 2*windex ))
896 iter = iter + 1
897
898* sin alpha <= |resid|/gap
899* Note that both the residual and the gap are
900* proportional to the matrix, so ||T|| doesn't play
901* a role in the quotient
902
903*
904* Convergence test for Rayleigh-Quotient iteration
905* (omitted when Bisection has been used)
906*
907 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
908 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
909 $ THEN
910* We need to check that the RQCORR update doesn't
911* move the eigenvalue away from the desired one and
912* towards a neighbor. -> protection with bisection
913 IF(indeig.LE.negcnt) THEN
914* The wanted eigenvalue lies to the left
915 sgndef = -one
916 ELSE
917* The wanted eigenvalue lies to the right
918 sgndef = one
919 ENDIF
920* We only use the RQCORR if it improves the
921* the iterate reasonably.
922 IF( ( rqcorr*sgndef.GE.zero )
923 $ .AND.( lambda + rqcorr.LE. right)
924 $ .AND.( lambda + rqcorr.GE. left)
925 $ ) THEN
926 usedrq = .true.
927* Store new midpoint of bisection interval in WORK
928 IF(sgndef.EQ.one) THEN
929* The current LAMBDA is on the left of the true
930* eigenvalue
931 left = lambda
932* We prefer to assume that the error estimate
933* is correct. We could make the interval not
934* as a bracket but to be modified if the RQCORR
935* chooses to. In this case, the RIGHT side should
936* be modified as follows:
937* RIGHT = MAX(RIGHT, LAMBDA + RQCORR)
938 ELSE
939* The current LAMBDA is on the right of the true
940* eigenvalue
941 right = lambda
942* See comment about assuming the error estimate is
943* correct above.
944* LEFT = MIN(LEFT, LAMBDA + RQCORR)
945 ENDIF
946 work( windex ) =
947 $ half * (right + left)
948* Take RQCORR since it has the correct sign and
949* improves the iterate reasonably
950 lambda = lambda + rqcorr
951* Update width of error interval
952 werr( windex ) =
953 $ half * (right-left)
954 ELSE
955 needbs = .true.
956 ENDIF
957 IF(right-left.LT.rqtol*abs(lambda)) THEN
958* The eigenvalue is computed to bisection accuracy
959* compute eigenvector and stop
960 usedbs = .true.
961 GOTO 120
962 ELSEIF( iter.LT.maxitr ) THEN
963 GOTO 120
964 ELSEIF( iter.EQ.maxitr ) THEN
965 needbs = .true.
966 GOTO 120
967 ELSE
968 info = 5
969 RETURN
970 END IF
971 ELSE
972 stp2ii = .false.
973 IF(usedrq .AND. usedbs .AND.
974 $ bstres.LE.resid) THEN
975 lambda = bstw
976 stp2ii = .true.
977 ENDIF
978 IF (stp2ii) THEN
979* improve error angle by second step
980 CALL clar1v( in, 1, in, lambda,
981 $ d( ibegin ), l( ibegin ),
982 $ work(indld+ibegin-1),
983 $ work(indlld+ibegin-1),
984 $ pivmin, gaptol, z( ibegin, windex ),
985 $ .NOT.usedbs, negcnt, ztz, mingma,
986 $ iwork( iindr+windex ),
987 $ isuppz( 2*windex-1 ),
988 $ nrminv, resid, rqcorr, work( indwrk ) )
989 ENDIF
990 work( windex ) = lambda
991 END IF
992*
993* Compute FP-vector support w.r.t. whole matrix
994*
995 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
996 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
997 zfrom = isuppz( 2*windex-1 )
998 zto = isuppz( 2*windex )
999 isupmn = isupmn + oldien
1000 isupmx = isupmx + oldien
1001* Ensure vector is ok if support in the RQI has changed
1002 IF(isupmn.LT.zfrom) THEN
1003 DO 122 ii = isupmn,zfrom-1
1004 z( ii, windex ) = zero
1005 122 CONTINUE
1006 ENDIF
1007 IF(isupmx.GT.zto) THEN
1008 DO 123 ii = zto+1,isupmx
1009 z( ii, windex ) = zero
1010 123 CONTINUE
1011 ENDIF
1012 CALL csscal( zto-zfrom+1, nrminv,
1013 $ z( zfrom, windex ), 1 )
1014 125 CONTINUE
1015* Update W
1016 w( windex ) = lambda+sigma
1017* Recompute the gaps on the left and right
1018* But only allow them to become larger and not
1019* smaller (which can only happen through "bad"
1020* cancellation and doesn't reflect the theory
1021* where the initial gaps are underestimated due
1022* to WERR being too crude.)
1023 IF(.NOT.eskip) THEN
1024 IF( k.GT.1) THEN
1025 wgap( windmn ) = max( wgap(windmn),
1026 $ w(windex)-werr(windex)
1027 $ - w(windmn)-werr(windmn) )
1028 ENDIF
1029 IF( windex.LT.wend ) THEN
1030 wgap( windex ) = max( savgap,
1031 $ w( windpl )-werr( windpl )
1032 $ - w( windex )-werr( windex) )
1033 ENDIF
1034 ENDIF
1035 idone = idone + 1
1036 ENDIF
1037* here ends the code for the current child
1038*
1039 139 CONTINUE
1040* Proceed to any remaining child nodes
1041 newfst = j + 1
1042 140 CONTINUE
1043 150 CONTINUE
1044 ndepth = ndepth + 1
1045 GO TO 40
1046 END IF
1047 ibegin = iend + 1
1048 wbegin = wend + 1
1049 170 CONTINUE
1050*
1051
1052 RETURN
1053*
1054* End of CLARRV
1055*
1056 END
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
clar1v
subroutine clar1v(n, b1, bn, lambda, d, l, ld, lld, pivmin, gaptol, z, wantnc, negcnt, ztz, mingma, r, isuppz, nrminv, resid, rqcorr, work)
CLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
Definition clar1v.f:228
slarrb
subroutine slarrb(n, d, lld, ifirst, ilast, rtol1, rtol2, offset, w, wgap, werr, work, iwork, pivmin, spdiam, twist, info)
SLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition slarrb.f:194
slarrf
subroutine slarrf(n, d, l, ld, clstrt, clend, w, wgap, werr, spdiam, clgapl, clgapr, pivmin, sigma, dplus, lplus, work, info)
SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is rela...
Definition slarrf.f:191
clarrv
subroutine clarrv(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)
CLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition clarrv.f:284
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:104
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78