283 SUBROUTINE slarrv( N, VL, VU, D, L, PIVMIN,
284 $ isplit, m, dol, dou, minrgp,
285 $ rtol1, rtol2, w, werr, wgap,
286 $ iblock, indexw, gers, z, ldz, isuppz,
287 $ work, iwork, info )
295 INTEGER DOL, DOU, INFO, LDZ, M, N
296 REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU
299 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ),
300 $ isuppz( * ), iwork( * )
301 REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ),
302 $ wgap( * ), work( * )
310 parameter ( maxitr = 10 )
311 REAL ZERO, ONE, TWO, THREE, FOUR, HALF
312 parameter ( zero = 0.0e0, one = 1.0e0,
313 $ two = 2.0e0, three = 3.0e0,
314 $ four = 4.0e0, half = 0.5e0)
317 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
318 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
319 $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
320 $ indld, indlld, indwrk, isupmn, isupmx, iter,
321 $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
322 $ ndepth, negcnt, newcls, newfst, newftt, newlst,
323 $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
324 $ oldncl, p, parity, q, wbegin, wend, windex,
325 $ windmn, windpl, zfrom, zto, zusedl, zusedu,
327 REAL 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
341 INTRINSIC abs,
REAL, MAX, MIN
382 zusedw = zusedu - zusedl + 1
385 CALL slaset(
'Full', n, zusedw, zero, zero,
388 eps = slamch(
'Precision' )
394 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
413 DO 170 jblk = 1, iblock( m )
414 iend = isplit( jblk )
421 IF( iblock( wend+1 ).EQ.jblk )
THEN
426 IF( wend.LT.wbegin )
THEN
429 ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) )
THEN
436 gl = gers( 2*ibegin-1 )
437 gu = gers( 2*ibegin )
438 DO 20 i = ibegin+1 , iend
439 gl = min( gers( 2*i-1 ), gl )
440 gu = max( gers( 2*i ), gu )
447 in = iend - ibegin + 1
449 im = wend - wbegin + 1
452 IF( ibegin.EQ.iend )
THEN
454 z( ibegin, wbegin ) = one
455 isuppz( 2*wbegin-1 ) = ibegin
456 isuppz( 2*wbegin ) = ibegin
457 w( wbegin ) = w( wbegin ) + sigma
458 work( wbegin ) = w( wbegin )
470 CALL scopy( im, w( wbegin ), 1,
471 $ work( wbegin ), 1 )
476 w(wbegin+i-1) = w(wbegin+i-1)+sigma
487 iwork( iindc1+1 ) = 1
488 iwork( iindc1+2 ) = im
497 IF( idone.LT.im )
THEN
499 IF( ndepth.GT.m )
THEN
510 IF( parity.EQ.0 )
THEN
523 oldfst = iwork( j-1 )
525 IF( ndepth.GT.0 )
THEN
531 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
534 j = wbegin + oldfst - 1
536 IF(wbegin+oldfst-1.LT.dol)
THEN
539 ELSEIF(wbegin+oldfst-1.GT.dou)
THEN
543 j = wbegin + oldfst - 1
546 CALL scopy( in, z( ibegin, j ), 1, d( ibegin ), 1 )
547 CALL scopy( in-1, z( ibegin, j+1 ), 1, l( ibegin ),
549 sigma = z( iend, j+1 )
552 CALL slaset(
'Full', in, 2, zero, zero,
553 $ z( ibegin, j), ldz )
557 DO 50 j = ibegin, iend-1
559 work( indld-1+j ) = tmp
560 work( indlld-1+j ) = tmp*l( j )
563 IF( ndepth.GT.0 )
THEN
566 p = indexw( wbegin-1+oldfst )
567 q = indexw( wbegin-1+oldlst )
571 offset = indexw( wbegin ) - 1
574 CALL slarrb( in, d( ibegin ),
575 $ work(indlld+ibegin-1),
576 $ p, q, rtol1, rtol2, offset,
577 $ work(wbegin),wgap(wbegin),werr(wbegin),
578 $ work( indwrk ), iwork( iindwk ),
579 $ pivmin, spdiam, in, iinfo )
580 IF( iinfo.NE.0 )
THEN
591 IF( oldfst.GT.1)
THEN
592 wgap( wbegin+oldfst-2 ) =
593 $ max(wgap(wbegin+oldfst-2),
594 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
595 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
597 IF( wbegin + oldlst -1 .LT. wend )
THEN
598 wgap( wbegin+oldlst-1 ) =
599 $ max(wgap(wbegin+oldlst-1),
600 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
601 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
605 DO 53 j=oldfst,oldlst
606 w(wbegin+j-1) = work(wbegin+j-1)+sigma
612 DO 140 j = oldfst, oldlst
613 IF( j.EQ.oldlst )
THEN
617 ELSE IF ( wgap( wbegin + j -1).GE.
618 $ minrgp* abs( work(wbegin + j -1) ) )
THEN
629 newsiz = newlst - newfst + 1
633 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
636 newftt = wbegin + newfst - 1
638 IF(wbegin+newfst-1.LT.dol)
THEN
641 ELSEIF(wbegin+newfst-1.GT.dou)
THEN
645 newftt = wbegin + newfst - 1
649 IF( newsiz.GT.1)
THEN
664 IF( newfst.EQ.1 )
THEN
666 $ w(wbegin)-werr(wbegin) - vl )
668 lgap = wgap( wbegin+newfst-2 )
670 rgap = wgap( wbegin+newlst-1 )
679 p = indexw( wbegin-1+newfst )
681 p = indexw( wbegin-1+newlst )
683 offset = indexw( wbegin ) - 1
684 CALL slarrb( in, d(ibegin),
685 $ work( indlld+ibegin-1 ),p,p,
686 $ rqtol, rqtol, offset,
687 $ work(wbegin),wgap(wbegin),
688 $ werr(wbegin),work( indwrk ),
689 $ iwork( iindwk ), pivmin, spdiam,
693 IF((wbegin+newlst-1.LT.dol).OR.
694 $ (wbegin+newfst-1.GT.dou))
THEN
702 idone = idone + newlst - newfst + 1
710 CALL slarrf( in, d( ibegin ), l( ibegin ),
711 $ work(indld+ibegin-1),
712 $ newfst, newlst, work(wbegin),
713 $ wgap(wbegin), werr(wbegin),
714 $ spdiam, lgap, rgap, pivmin, tau,
715 $ z(ibegin, newftt),z(ibegin, newftt+1),
716 $ work( indwrk ), iinfo )
717 IF( iinfo.EQ.0 )
THEN
721 z( iend, newftt+1 ) = ssigma
724 DO 116 k = newfst, newlst
726 $ three*eps*abs(work(wbegin+k-1))
727 work( wbegin + k - 1 ) =
728 $ work( wbegin + k - 1) - tau
730 $ four*eps*abs(work(wbegin+k-1))
732 werr( wbegin + k - 1 ) =
733 $ werr( wbegin + k - 1 ) + fudge
745 iwork( k-1 ) = newfst
757 tol = four * log(
REAL(in)) * eps
760 windex = wbegin + k - 1
761 windmn = max(windex - 1,1)
762 windpl = min(windex + 1,m)
763 lambda = work( windex )
766 IF((windex.LT.dol).OR.
767 $ (windex.GT.dou))
THEN
773 left = work( windex ) - werr( windex )
774 right = work( windex ) + werr( windex )
775 indeig = indexw( windex )
790 lgap = eps*max(abs(left),abs(right))
800 rgap = eps*max(abs(left),abs(right))
804 gap = min( lgap, rgap )
805 IF(( k .EQ. 1).OR.(k .EQ. im))
THEN
820 savgap = wgap(windex)
837 itmp1 = iwork( iindr+windex )
838 offset = indexw( wbegin ) - 1
839 CALL slarrb( in, d(ibegin),
840 $ work(indlld+ibegin-1),indeig,indeig,
841 $ zero, two*eps, offset,
842 $ work(wbegin),wgap(wbegin),
843 $ werr(wbegin),work( indwrk ),
844 $ iwork( iindwk ), pivmin, spdiam,
846 IF( iinfo.NE.0 )
THEN
850 lambda = work( windex )
853 iwork( iindr+windex ) = 0
856 CALL slar1v( in, 1, in, lambda, d( ibegin ),
857 $ l( ibegin ), work(indld+ibegin-1),
858 $ work(indlld+ibegin-1),
859 $ pivmin, gaptol, z( ibegin, windex ),
860 $ .NOT.usedbs, negcnt, ztz, mingma,
861 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
862 $ nrminv, resid, rqcorr, work( indwrk ) )
866 ELSEIF(resid.LT.bstres)
THEN
870 isupmn = min(isupmn,isuppz( 2*windex-1 ))
871 isupmx = max(isupmx,isuppz( 2*windex ))
883 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
884 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
889 IF(indeig.LE.negcnt)
THEN
898 IF( ( rqcorr*sgndef.GE.zero )
899 $ .AND.( lambda + rqcorr.LE. right)
900 $ .AND.( lambda + rqcorr.GE. left)
904 IF(sgndef.EQ.one)
THEN
923 $ half * (right + left)
926 lambda = lambda + rqcorr
929 $ half * (right-left)
933 IF(right-left.LT.rqtol*abs(lambda))
THEN
938 ELSEIF( iter.LT.maxitr )
THEN
940 ELSEIF( iter.EQ.maxitr )
THEN
949 IF(usedrq .AND. usedbs .AND.
950 $ bstres.LE.resid)
THEN
956 CALL slar1v( in, 1, in, lambda,
957 $ d( ibegin ), l( ibegin ),
958 $ work(indld+ibegin-1),
959 $ work(indlld+ibegin-1),
960 $ pivmin, gaptol, z( ibegin, windex ),
961 $ .NOT.usedbs, negcnt, ztz, mingma,
962 $ iwork( iindr+windex ),
963 $ isuppz( 2*windex-1 ),
964 $ nrminv, resid, rqcorr, work( indwrk ) )
966 work( windex ) = lambda
971 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
972 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
973 zfrom = isuppz( 2*windex-1 )
974 zto = isuppz( 2*windex )
975 isupmn = isupmn + oldien
976 isupmx = isupmx + oldien
978 IF(isupmn.LT.zfrom)
THEN
979 DO 122 ii = isupmn,zfrom-1
980 z( ii, windex ) = zero
983 IF(isupmx.GT.zto)
THEN
984 DO 123 ii = zto+1,isupmx
985 z( ii, windex ) = zero
988 CALL sscal( zto-zfrom+1, nrminv,
989 $ z( zfrom, windex ), 1 )
992 w( windex ) = lambda+sigma
1001 wgap( windmn ) = max( wgap(windmn),
1002 $ w(windex)-werr(windex)
1003 $ - w(windmn)-werr(windmn) )
1005 IF( windex.LT.wend )
THEN
1006 wgap( windex ) = max( savgap,
1007 $ w( windpl )-werr( windpl )
1008 $ - w( windex )-werr( windex) )
subroutine slar1v(N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)
SLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
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.
subroutine slarrv(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)
SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
subroutine sscal(N, SA, SX, INCX)
SSCAL
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...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY