LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlarrv()

subroutine dlarrv ( integer  n,
double precision  vl,
double precision  vu,
double precision, dimension( * )  d,
double precision, dimension( * )  l,
double precision  pivmin,
integer, dimension( * )  isplit,
integer  m,
integer  dol,
integer  dou,
double precision  minrgp,
double precision  rtol1,
double precision  rtol2,
double precision, dimension( * )  w,
double precision, dimension( * )  werr,
double precision, dimension( * )  wgap,
integer, dimension( * )  iblock,
integer, dimension( * )  indexw,
double precision, dimension( * )  gers,
double precision, dimension( ldz, * )  z,
integer  ldz,
integer, dimension( * )  isuppz,
double precision, dimension( * )  work,
integer, dimension( * )  iwork,
integer  info 
)

DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.

Download DLARRV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLARRV computes the eigenvectors of the tridiagonal matrix
 T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.
 The input eigenvalues should have been computed by DLARRE.
Parameters
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in]VL
          VL is DOUBLE PRECISION
          Lower bound of the interval that contains the desired
          eigenvalues. VL < VU. Needed to compute gaps on the left or right
          end of the extremal eigenvalues in the desired RANGE.
[in]VU
          VU is DOUBLE PRECISION
          Upper bound of the interval that contains the desired
          eigenvalues. VL < VU. 
          Note: VU is currently not used by this implementation of DLARRV, VU is
          passed to DLARRV because it could be used compute gaps on the right end
          of the extremal eigenvalues. However, with not much initial accuracy in
          LAMBDA and VU, the formula can lead to an overestimation of the right gap
          and thus to inadequately early RQI 'convergence'. This is currently
          prevented this by forcing a small right gap. And so it turns out that VU
          is currently not used by this implementation of DLARRV.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the N diagonal elements of the diagonal matrix D.
          On exit, D may be overwritten.
[in,out]L
          L is DOUBLE PRECISION array, dimension (N)
          On entry, the (N-1) subdiagonal elements of the unit
          bidiagonal matrix L are in elements 1 to N-1 of L
          (if the matrix is not split.) At the end of each block
          is stored the corresponding shift as given by DLARRE.
          On exit, L is overwritten.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot allowed in the Sturm sequence.
[in]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into blocks.
          The first block consists of rows/columns 1 to
          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
          through ISPLIT( 2 ), etc.
[in]M
          M is INTEGER
          The total number of input eigenvalues.  0 <= M <= N.
[in]DOL
          DOL is INTEGER
[in]DOU
          DOU is INTEGER
          If the user wants to compute only selected eigenvectors from all
          the eigenvalues supplied, he can specify an index range DOL:DOU.
          Or else the setting DOL=1, DOU=M should be applied.
          Note that DOL and DOU refer to the order in which the eigenvalues
          are stored in W.
          If the user wants to compute only selected eigenpairs, then
          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the
          computed eigenvectors. All other columns of Z are set to zero.
[in]MINRGP
          MINRGP is DOUBLE PRECISION
[in]RTOL1
          RTOL1 is DOUBLE PRECISION
[in]RTOL2
          RTOL2 is DOUBLE PRECISION
           Parameters for bisection.
           An interval [LEFT,RIGHT] has converged if
           RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
[in,out]W
          W is DOUBLE PRECISION array, dimension (N)
          The first M elements of W contain the APPROXIMATE eigenvalues for
          which eigenvectors are to be computed.  The eigenvalues
          should be grouped by split-off block and ordered from
          smallest to largest within the block ( The output array
          W from DLARRE is expected here ). Furthermore, they are with
          respect to the shift of the corresponding root representation
          for their block. On exit, W holds the eigenvalues of the
          UNshifted matrix.
[in,out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          The first M elements contain the semiwidth of the uncertainty
          interval of the corresponding eigenvalue in W
[in,out]WGAP
          WGAP is DOUBLE PRECISION array, dimension (N)
          The separation from the right neighbor eigenvalue in W.
[in]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          The indices of the blocks (submatrices) associated with the
          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
          W(i) belongs to the first block from the top, =2 if W(i)
          belongs to the second block, etc.
[in]INDEXW
          INDEXW is INTEGER array, dimension (N)
          The indices of the eigenvalues within each block (submatrix);
          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.
[in]GERS
          GERS is DOUBLE PRECISION array, dimension (2*N)
          The N Gerschgorin intervals (the i-th Gerschgorin interval
          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should
          be computed from the original UNshifted matrix.
[out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
          If INFO = 0, the first M columns of Z contain the
          orthonormal eigenvectors of the matrix T
          corresponding to the input eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]ISUPPZ
          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
          The support of the eigenvectors in Z, i.e., the indices
          indicating the nonzero elements in Z. The I-th eigenvector
          is nonzero only in elements ISUPPZ( 2*I-1 ) through
          ISUPPZ( 2*I ).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (12*N)
[out]IWORK
          IWORK is INTEGER array, dimension (7*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit

          > 0:  A problem occurred in DLARRV.
          < 0:  One of the called subroutines signaled an internal problem.
                Needs inspection of the corresponding parameter IINFO
                for further information.

          =-1:  Problem in DLARRB when refining a child's eigenvalues.
          =-2:  Problem in DLARRF when computing the RRR of a child.
                When a child is inside a tight cluster, it can be difficult
                to find an RRR. A partial remedy from the user's point of
                view is to make the parameter MINRGP smaller and recompile.
                However, as the orthogonality of the computed vectors is
                proportional to 1/MINRGP, the user should be aware that
                he might be trading in precision when he decreases MINRGP.
          =-3:  Problem in DLARRB when refining a single eigenvalue
                after the Rayleigh correction was rejected.
          = 5:  The Rayleigh Quotient Iteration failed to converge to
                full accuracy in MAXITR steps.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 287 of file dlarrv.f.

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