280 SUBROUTINE clarrv( 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 REAL minrgp, pivmin, rtol1, rtol2, vl, vu
296 INTEGER iblock( * ), indexw( * ), isplit( * ),
297 $ isuppz( * ), iwork( * )
298 REAL d( * ), gers( * ), l( * ), w( * ), werr( * ),
299 $ wgap( * ), work( * )
307 parameter( maxitr = 10 )
309 parameter( czero = ( 0.0e0, 0.0e0 ) )
310 REAL zero, one, two, three, four, half
311 parameter( zero = 0.0e0, one = 1.0e0,
312 $ two = 2.0e0, three = 3.0e0,
313 $ four = 4.0e0, half = 0.5e0)
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,
326 INTEGER indin1, indin2
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
384 zusedw = zusedu - zusedl + 1
387 CALL
claset(
'Full', n, zusedw, czero, czero,
390 eps =
slamch(
'Precision' )
396 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
415 DO 170 jblk = 1, iblock( m )
416 iend = isplit( jblk )
423 IF( iblock( wend+1 ).EQ.jblk )
THEN
428 IF( wend.LT.wbegin )
THEN
431 ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) )
THEN
438 gl = gers( 2*ibegin-1 )
439 gu = gers( 2*ibegin )
440 DO 20 i = ibegin+1 , iend
441 gl = min( gers( 2*i-1 ), gl )
442 gu = max( gers( 2*i ), gu )
449 in = iend - ibegin + 1
451 im = wend - wbegin + 1
454 IF( ibegin.EQ.iend )
THEN
456 z( ibegin, wbegin ) = cmplx( one, zero )
457 isuppz( 2*wbegin-1 ) = ibegin
458 isuppz( 2*wbegin ) = ibegin
459 w( wbegin ) = w( wbegin ) + sigma
460 work( wbegin ) = w( wbegin )
472 CALL
scopy( im, w( wbegin ), 1,
473 $ work( wbegin ), 1 )
478 w(wbegin+i-1) = w(wbegin+i-1)+sigma
489 iwork( iindc1+1 ) = 1
490 iwork( iindc1+2 ) = im
499 IF( idone.LT.im )
THEN
501 IF( ndepth.GT.m )
THEN
512 IF( parity.EQ.0 )
THEN
525 oldfst = iwork(
j-1 )
527 IF( ndepth.GT.0 )
THEN
533 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
536 j = wbegin + oldfst - 1
538 IF(wbegin+oldfst-1.LT.dol)
THEN
541 ELSEIF(wbegin+oldfst-1.GT.dou)
THEN
545 j = wbegin + oldfst - 1
549 d( ibegin+k-1 ) =
REAL( Z( IBEGIN+K-1,
$ J ) )
550 l( ibegin+k-1 ) =
REAL( Z( IBEGIN+K-1,
$ J+1 ) )
552 d( iend ) =
REAL( Z( IEND, J ) )
553 sigma =
REAL( Z( IEND, J+1 ) )
556 CALL
claset(
'Full', in, 2, czero, czero,
557 $ z( ibegin,
j), ldz )
561 DO 50
j = ibegin, iend-1
563 work( indld-1+
j ) = tmp
564 work( indlld-1+
j ) = tmp*l(
j )
567 IF( ndepth.GT.0 )
THEN
570 p = indexw( wbegin-1+oldfst )
571 q = indexw( wbegin-1+oldlst )
575 offset = indexw( wbegin ) - 1
578 CALL
slarrb( in, d( ibegin ),
579 $ work(indlld+ibegin-1),
580 $ p, q, rtol1, rtol2, offset,
581 $ work(wbegin),wgap(wbegin),werr(wbegin),
582 $ work( indwrk ), iwork( iindwk ),
583 $ pivmin, spdiam, in, iinfo )
584 IF( iinfo.NE.0 )
THEN
595 IF( oldfst.GT.1)
THEN
596 wgap( wbegin+oldfst-2 ) =
597 $ max(wgap(wbegin+oldfst-2),
598 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
599 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
601 IF( wbegin + oldlst -1 .LT. wend )
THEN
602 wgap( wbegin+oldlst-1 ) =
603 $ max(wgap(wbegin+oldlst-1),
604 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
605 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
609 DO 53
j=oldfst,oldlst
610 w(wbegin+
j-1) = work(wbegin+
j-1)+sigma
616 DO 140
j = oldfst, oldlst
617 IF(
j.EQ.oldlst )
THEN
621 ELSE IF ( wgap( wbegin +
j -1).GE.
622 $ minrgp* abs( work(wbegin +
j -1) ) )
THEN
633 newsiz = newlst - newfst + 1
637 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
640 newftt = wbegin + newfst - 1
642 IF(wbegin+newfst-1.LT.dol)
THEN
645 ELSEIF(wbegin+newfst-1.GT.dou)
THEN
649 newftt = wbegin + newfst - 1
653 IF( newsiz.GT.1)
THEN
668 IF( newfst.EQ.1 )
THEN
670 $ w(wbegin)-werr(wbegin) - vl )
672 lgap = wgap( wbegin+newfst-2 )
674 rgap = wgap( wbegin+newlst-1 )
683 p = indexw( wbegin-1+newfst )
685 p = indexw( wbegin-1+newlst )
687 offset = indexw( wbegin ) - 1
688 CALL
slarrb( in, d(ibegin),
689 $ work( indlld+ibegin-1 ),p,p,
690 $ rqtol, rqtol, offset,
691 $ work(wbegin),wgap(wbegin),
692 $ werr(wbegin),work( indwrk ),
693 $ iwork( iindwk ), pivmin, spdiam,
697 IF((wbegin+newlst-1.LT.dol).OR.
698 $ (wbegin+newfst-1.GT.dou))
THEN
706 idone = idone + newlst - newfst + 1
714 CALL
slarrf( in, d( ibegin ), l( ibegin ),
715 $ work(indld+ibegin-1),
716 $ newfst, newlst, work(wbegin),
717 $ wgap(wbegin), werr(wbegin),
718 $ spdiam, lgap, rgap, pivmin, tau,
719 $ work( indin1 ), work( indin2 ),
720 $ work( indwrk ), iinfo )
725 z( ibegin+k-1, newftt ) =
726 $ cmplx( work( indin1+k-1 ), zero )
727 z( ibegin+k-1, newftt+1 ) =
728 $ cmplx( work( indin2+k-1 ), zero )
731 $ cmplx( work( indin1+in-1 ), zero )
732 IF( iinfo.EQ.0 )
THEN
736 z( iend, newftt+1 ) = cmplx( ssigma, zero )
739 DO 116 k = newfst, newlst
741 $ three*eps*abs(work(wbegin+k-1))
742 work( wbegin + k - 1 ) =
743 $ work( wbegin + k - 1) - tau
745 $ four*eps*abs(work(wbegin+k-1))
747 werr( wbegin + k - 1 ) =
748 $ werr( wbegin + k - 1 ) + fudge
760 iwork( k-1 ) = newfst
772 tol = four * log(
REAL(in)) * eps
775 windex = wbegin + k - 1
776 windmn = max(windex - 1,1)
777 windpl = min(windex + 1,m)
778 lambda = work( windex )
781 IF((windex.LT.dol).OR.
782 $ (windex.GT.dou))
THEN
788 left = work( windex ) - werr( windex )
789 right = work( windex ) + werr( windex )
790 indeig = indexw( windex )
805 lgap = eps*max(abs(left),abs(right))
815 rgap = eps*max(abs(left),abs(right))
819 gap = min( lgap, rgap )
820 IF(( k .EQ. 1).OR.(k .EQ. im))
THEN
835 savgap = wgap(windex)
852 itmp1 = iwork( iindr+windex )
853 offset = indexw( wbegin ) - 1
854 CALL
slarrb( in, d(ibegin),
855 $ work(indlld+ibegin-1),indeig,indeig,
856 $ zero, two*eps, offset,
857 $ work(wbegin),wgap(wbegin),
858 $ werr(wbegin),work( indwrk ),
859 $ iwork( iindwk ), pivmin, spdiam,
861 IF( iinfo.NE.0 )
THEN
865 lambda = work( windex )
868 iwork( iindr+windex ) = 0
871 CALL
clar1v( in, 1, in, lambda, d( ibegin ),
872 $ l( ibegin ), work(indld+ibegin-1),
873 $ work(indlld+ibegin-1),
874 $ pivmin, gaptol, z( ibegin, windex ),
875 $ .NOT.usedbs, negcnt, ztz, mingma,
876 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
877 $ nrminv, resid, rqcorr, work( indwrk ) )
881 ELSEIF(resid.LT.bstres)
THEN
885 isupmn = min(isupmn,isuppz( 2*windex-1 ))
886 isupmx = max(isupmx,isuppz( 2*windex ))
898 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
899 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
904 IF(indeig.LE.negcnt)
THEN
913 IF( ( rqcorr*sgndef.GE.zero )
914 $ .AND.( lambda + rqcorr.LE. right)
915 $ .AND.( lambda + rqcorr.GE. left)
919 IF(sgndef.EQ.one)
THEN
938 $ half * (right + left)
941 lambda = lambda + rqcorr
944 $ half * (right-left)
948 IF(right-left.LT.rqtol*abs(lambda))
THEN
953 ELSEIF( iter.LT.maxitr )
THEN
955 ELSEIF( iter.EQ.maxitr )
THEN
964 IF(usedrq .AND. usedbs .AND.
965 $ bstres.LE.resid)
THEN
971 CALL
clar1v( in, 1, in, lambda,
972 $ d( ibegin ), l( ibegin ),
973 $ work(indld+ibegin-1),
974 $ work(indlld+ibegin-1),
975 $ pivmin, gaptol, z( ibegin, windex ),
976 $ .NOT.usedbs, negcnt, ztz, mingma,
977 $ iwork( iindr+windex ),
978 $ isuppz( 2*windex-1 ),
979 $ nrminv, resid, rqcorr, work( indwrk ) )
981 work( windex ) = lambda
986 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
987 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
988 zfrom = isuppz( 2*windex-1 )
989 zto = isuppz( 2*windex )
990 isupmn = isupmn + oldien
991 isupmx = isupmx + oldien
993 IF(isupmn.LT.zfrom)
THEN
994 DO 122 ii = isupmn,zfrom-1
995 z( ii, windex ) = zero
998 IF(isupmx.GT.zto)
THEN
999 DO 123 ii = zto+1,isupmx
1000 z( ii, windex ) = zero
1003 CALL
csscal( zto-zfrom+1, nrminv,
1004 $ z( zfrom, windex ), 1 )
1007 w( windex ) = lambda+sigma
1016 wgap( windmn ) = max( wgap(windmn),
1017 $ w(windex)-werr(windex)
1018 $ - w(windmn)-werr(windmn) )
1020 IF( windex.LT.wend )
THEN
1021 wgap( windex ) = max( savgap,
1022 $ w( windpl )-werr( windpl )
1023 $ - w( windex )-werr( windex) )
1048 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...
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...
REAL function slamch(CMACH)
SLAMCH
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 scopy(N, SX, INCX, SY, INCY)
SCOPY
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 ...
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...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine csscal(N, SA, CX, INCX)
CSSCAL