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