280 SUBROUTINE dlarrv( N, VL, VU, D, L, PIVMIN,
281 $ isplit, m, dol, dou, minrgp,
282 $ rtol1, rtol2, w, werr, wgap,
283 $ iblock, indexw, gers, z, ldz, isuppz,
284 $ work, iwork, info )
292 INTEGER DOL, DOU, INFO, LDZ, M, N
293 DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
296 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
297 $ isuppz( * ), iwork( * )
298 DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
299 $ wgap( * ), work( * )
300 DOUBLE PRECISION Z( ldz, * )
307 parameter( maxitr = 10 )
308 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF
309 parameter( zero = 0.0d0, one = 1.0d0,
310 $ two = 2.0d0, three = 3.0d0,
311 $ four = 4.0d0, half = 0.5d0)
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,
324 DOUBLE PRECISION BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU,
325 $ lambda, left, lgap, mingma, nrminv, resid,
326 $ rgap, right, rqcorr, rqtol, savgap, sgndef,
327 $ sigma, spdiam, ssigma, tau, tmp, tol, ztz
330 DOUBLE PRECISION DLAMCH
338 INTRINSIC abs, dble, max, min
378 zusedw = zusedu - zusedl + 1
381 CALL dlaset(
'Full', n, zusedw, zero, zero,
384 eps = dlamch(
'Precision' )
390 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
409 DO 170 jblk = 1, iblock( m )
410 iend = isplit( jblk )
417 IF( iblock( wend+1 ).EQ.jblk )
THEN
422 IF( wend.LT.wbegin )
THEN
425 ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) )
THEN
432 gl = gers( 2*ibegin-1 )
433 gu = gers( 2*ibegin )
434 DO 20 i = ibegin+1 , iend
435 gl = min( gers( 2*i-1 ), gl )
436 gu = max( gers( 2*i ), gu )
443 in = iend - ibegin + 1
445 im = wend - wbegin + 1
448 IF( ibegin.EQ.iend )
THEN
450 z( ibegin, wbegin ) = one
451 isuppz( 2*wbegin-1 ) = ibegin
452 isuppz( 2*wbegin ) = ibegin
453 w( wbegin ) = w( wbegin ) + sigma
454 work( wbegin ) = w( wbegin )
466 CALL dcopy( im, w( wbegin ), 1,
467 $ work( wbegin ), 1 )
472 w(wbegin+i-1) = w(wbegin+i-1)+sigma
483 iwork( iindc1+1 ) = 1
484 iwork( iindc1+2 ) = im
493 IF( idone.LT.im )
THEN
495 IF( ndepth.GT.m )
THEN
506 IF( parity.EQ.0 )
THEN
519 oldfst = iwork( j-1 )
521 IF( ndepth.GT.0 )
THEN
527 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
530 j = wbegin + oldfst - 1
532 IF(wbegin+oldfst-1.LT.dol)
THEN
535 ELSEIF(wbegin+oldfst-1.GT.dou)
THEN
539 j = wbegin + oldfst - 1
542 CALL dcopy( in, z( ibegin, j ), 1, d( ibegin ), 1 )
543 CALL dcopy( in-1, z( ibegin, j+1 ), 1, l( ibegin ),
545 sigma = z( iend, j+1 )
548 CALL dlaset(
'Full', in, 2, zero, zero,
549 $ z( ibegin, j), ldz )
553 DO 50 j = ibegin, iend-1
555 work( indld-1+j ) = tmp
556 work( indlld-1+j ) = tmp*l( j )
559 IF( ndepth.GT.0 )
THEN
562 p = indexw( wbegin-1+oldfst )
563 q = indexw( wbegin-1+oldlst )
567 offset = indexw( wbegin ) - 1
570 CALL dlarrb( in, d( ibegin ),
571 $ work(indlld+ibegin-1),
572 $ p, q, rtol1, rtol2, offset,
573 $ work(wbegin),wgap(wbegin),werr(wbegin),
574 $ work( indwrk ), iwork( iindwk ),
575 $ pivmin, spdiam, in, iinfo )
576 IF( iinfo.NE.0 )
THEN
587 IF( oldfst.GT.1)
THEN
588 wgap( wbegin+oldfst-2 ) =
589 $ max(wgap(wbegin+oldfst-2),
590 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
591 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
593 IF( wbegin + oldlst -1 .LT. wend )
THEN
594 wgap( wbegin+oldlst-1 ) =
595 $ max(wgap(wbegin+oldlst-1),
596 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
597 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
601 DO 53 j=oldfst,oldlst
602 w(wbegin+j-1) = work(wbegin+j-1)+sigma
608 DO 140 j = oldfst, oldlst
609 IF( j.EQ.oldlst )
THEN
613 ELSE IF ( wgap( wbegin + j -1).GE.
614 $ minrgp* abs( work(wbegin + j -1) ) )
THEN
625 newsiz = newlst - newfst + 1
629 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
632 newftt = wbegin + newfst - 1
634 IF(wbegin+newfst-1.LT.dol)
THEN
637 ELSEIF(wbegin+newfst-1.GT.dou)
THEN
641 newftt = wbegin + newfst - 1
645 IF( newsiz.GT.1)
THEN
660 IF( newfst.EQ.1 )
THEN
662 $ w(wbegin)-werr(wbegin) - vl )
664 lgap = wgap( wbegin+newfst-2 )
666 rgap = wgap( wbegin+newlst-1 )
675 p = indexw( wbegin-1+newfst )
677 p = indexw( wbegin-1+newlst )
679 offset = indexw( wbegin ) - 1
680 CALL dlarrb( in, d(ibegin),
681 $ work( indlld+ibegin-1 ),p,p,
682 $ rqtol, rqtol, offset,
683 $ work(wbegin),wgap(wbegin),
684 $ werr(wbegin),work( indwrk ),
685 $ iwork( iindwk ), pivmin, spdiam,
689 IF((wbegin+newlst-1.LT.dol).OR.
690 $ (wbegin+newfst-1.GT.dou))
THEN
698 idone = idone + newlst - newfst + 1
706 CALL dlarrf( in, d( ibegin ), l( ibegin ),
707 $ work(indld+ibegin-1),
708 $ newfst, newlst, work(wbegin),
709 $ wgap(wbegin), werr(wbegin),
710 $ spdiam, lgap, rgap, pivmin, tau,
711 $ z(ibegin, newftt),z(ibegin, newftt+1),
712 $ work( indwrk ), iinfo )
713 IF( iinfo.EQ.0 )
THEN
717 z( iend, newftt+1 ) = ssigma
720 DO 116 k = newfst, newlst
722 $ three*eps*abs(work(wbegin+k-1))
723 work( wbegin + k - 1 ) =
724 $ work( wbegin + k - 1) - tau
726 $ four*eps*abs(work(wbegin+k-1))
728 werr( wbegin + k - 1 ) =
729 $ werr( wbegin + k - 1 ) + fudge
741 iwork( k-1 ) = newfst
753 tol = four * log(dble(in)) * eps
756 windex = wbegin + k - 1
757 windmn = max(windex - 1,1)
758 windpl = min(windex + 1,m)
759 lambda = work( windex )
762 IF((windex.LT.dol).OR.
763 $ (windex.GT.dou))
THEN
769 left = work( windex ) - werr( windex )
770 right = work( windex ) + werr( windex )
771 indeig = indexw( windex )
786 lgap = eps*max(abs(left),abs(right))
796 rgap = eps*max(abs(left),abs(right))
800 gap = min( lgap, rgap )
801 IF(( k .EQ. 1).OR.(k .EQ. im))
THEN
816 savgap = wgap(windex)
833 itmp1 = iwork( iindr+windex )
834 offset = indexw( wbegin ) - 1
835 CALL dlarrb( in, d(ibegin),
836 $ work(indlld+ibegin-1),indeig,indeig,
837 $ zero, two*eps, offset,
838 $ work(wbegin),wgap(wbegin),
839 $ werr(wbegin),work( indwrk ),
840 $ iwork( iindwk ), pivmin, spdiam,
842 IF( iinfo.NE.0 )
THEN
846 lambda = work( windex )
849 iwork( iindr+windex ) = 0
852 CALL dlar1v( in, 1, in, lambda, d( ibegin ),
853 $ l( ibegin ), work(indld+ibegin-1),
854 $ work(indlld+ibegin-1),
855 $ pivmin, gaptol, z( ibegin, windex ),
856 $ .NOT.usedbs, negcnt, ztz, mingma,
857 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
858 $ nrminv, resid, rqcorr, work( indwrk ) )
862 ELSEIF(resid.LT.bstres)
THEN
866 isupmn = min(isupmn,isuppz( 2*windex-1 ))
867 isupmx = max(isupmx,isuppz( 2*windex ))
879 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
880 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
885 IF(indeig.LE.negcnt)
THEN
894 IF( ( rqcorr*sgndef.GE.zero )
895 $ .AND.( lambda + rqcorr.LE. right)
896 $ .AND.( lambda + rqcorr.GE. left)
900 IF(sgndef.EQ.one)
THEN
919 $ half * (right + left)
922 lambda = lambda + rqcorr
925 $ half * (right-left)
929 IF(right-left.LT.rqtol*abs(lambda))
THEN
934 ELSEIF( iter.LT.maxitr )
THEN
936 ELSEIF( iter.EQ.maxitr )
THEN
945 IF(usedrq .AND. usedbs .AND.
946 $ bstres.LE.resid)
THEN
952 CALL dlar1v( in, 1, in, lambda,
953 $ d( ibegin ), l( ibegin ),
954 $ work(indld+ibegin-1),
955 $ work(indlld+ibegin-1),
956 $ pivmin, gaptol, z( ibegin, windex ),
957 $ .NOT.usedbs, negcnt, ztz, mingma,
958 $ iwork( iindr+windex ),
959 $ isuppz( 2*windex-1 ),
960 $ nrminv, resid, rqcorr, work( indwrk ) )
962 work( windex ) = lambda
967 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
968 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
969 zfrom = isuppz( 2*windex-1 )
970 zto = isuppz( 2*windex )
971 isupmn = isupmn + oldien
972 isupmx = isupmx + oldien
974 IF(isupmn.LT.zfrom)
THEN
975 DO 122 ii = isupmn,zfrom-1
976 z( ii, windex ) = zero
979 IF(isupmx.GT.zto)
THEN
980 DO 123 ii = zto+1,isupmx
981 z( ii, windex ) = zero
984 CALL dscal( zto-zfrom+1, nrminv,
985 $ z( zfrom, windex ), 1 )
988 w( windex ) = lambda+sigma
997 wgap( windmn ) = max( wgap(windmn),
998 $ w(windex)-werr(windex)
999 $ - w(windmn)-werr(windmn) )
1001 IF( windex.LT.wend )
THEN
1002 wgap( windex ) = max( savgap,
1003 $ w( windpl )-werr( windpl )
1004 $ - w( windex )-werr( windex) )
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...
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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...
subroutine dlarrv(N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, WORK, IWORK, INFO)
DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
subroutine dscal(N, DA, DX, INCX)
DSCAL
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.
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...