280 SUBROUTINE zlarrv( 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 COMPLEX*16 z( ldz, * )
307 parameter( maxitr = 10 )
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)
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 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
341 INTRINSIC abs, dble, max, min
384 zusedw = zusedu - zusedl + 1
387 CALL
zlaset(
'Full', n, zusedw, czero, czero,
390 eps =
dlamch(
'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 ) = dcmplx( 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
dcopy( 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 ) = dble( z( ibegin+k-1,
551 l( ibegin+k-1 ) = dble( z( ibegin+k-1,
554 d( iend ) = dble( z( iend,
j ) )
555 sigma = dble( z( iend,
j+1 ) )
558 CALL
zlaset(
'Full', in, 2, czero, czero,
559 $ z( ibegin,
j), ldz )
563 DO 50
j = ibegin, iend-1
565 work( indld-1+
j ) = tmp
566 work( indlld-1+
j ) = tmp*l(
j )
569 IF( ndepth.GT.0 )
THEN
572 p = indexw( wbegin-1+oldfst )
573 q = indexw( wbegin-1+oldlst )
577 offset = indexw( wbegin ) - 1
580 CALL
dlarrb( in, d( ibegin ),
581 $ work(indlld+ibegin-1),
582 $ p, q, rtol1, rtol2, offset,
583 $ work(wbegin),wgap(wbegin),werr(wbegin),
584 $ work( indwrk ), iwork( iindwk ),
585 $ pivmin, spdiam, in, iinfo )
586 IF( iinfo.NE.0 )
THEN
597 IF( oldfst.GT.1)
THEN
598 wgap( wbegin+oldfst-2 ) =
599 $ max(wgap(wbegin+oldfst-2),
600 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
601 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
603 IF( wbegin + oldlst -1 .LT. wend )
THEN
604 wgap( wbegin+oldlst-1 ) =
605 $ max(wgap(wbegin+oldlst-1),
606 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
607 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
611 DO 53
j=oldfst,oldlst
612 w(wbegin+
j-1) = work(wbegin+
j-1)+sigma
618 DO 140
j = oldfst, oldlst
619 IF(
j.EQ.oldlst )
THEN
623 ELSE IF ( wgap( wbegin +
j -1).GE.
624 $ minrgp* abs( work(wbegin +
j -1) ) )
THEN
635 newsiz = newlst - newfst + 1
639 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
642 newftt = wbegin + newfst - 1
644 IF(wbegin+newfst-1.LT.dol)
THEN
647 ELSEIF(wbegin+newfst-1.GT.dou)
THEN
651 newftt = wbegin + newfst - 1
655 IF( newsiz.GT.1)
THEN
670 IF( newfst.EQ.1 )
THEN
672 $ w(wbegin)-werr(wbegin) - vl )
674 lgap = wgap( wbegin+newfst-2 )
676 rgap = wgap( wbegin+newlst-1 )
685 p = indexw( wbegin-1+newfst )
687 p = indexw( wbegin-1+newlst )
689 offset = indexw( wbegin ) - 1
690 CALL
dlarrb( in, d(ibegin),
691 $ work( indlld+ibegin-1 ),p,p,
692 $ rqtol, rqtol, offset,
693 $ work(wbegin),wgap(wbegin),
694 $ werr(wbegin),work( indwrk ),
695 $ iwork( iindwk ), pivmin, spdiam,
699 IF((wbegin+newlst-1.LT.dol).OR.
700 $ (wbegin+newfst-1.GT.dou))
THEN
708 idone = idone + newlst - newfst + 1
716 CALL
dlarrf( in, d( ibegin ), l( ibegin ),
717 $ work(indld+ibegin-1),
718 $ newfst, newlst, work(wbegin),
719 $ wgap(wbegin), werr(wbegin),
720 $ spdiam, lgap, rgap, pivmin, tau,
721 $ work( indin1 ), work( indin2 ),
722 $ work( indwrk ), iinfo )
727 z( ibegin+k-1, newftt ) =
728 $ dcmplx( work( indin1+k-1 ), zero )
729 z( ibegin+k-1, newftt+1 ) =
730 $ dcmplx( work( indin2+k-1 ), zero )
733 $ dcmplx( work( indin1+in-1 ), zero )
734 IF( iinfo.EQ.0 )
THEN
738 z( iend, newftt+1 ) = dcmplx( ssigma, zero )
741 DO 116 k = newfst, newlst
743 $ three*eps*abs(work(wbegin+k-1))
744 work( wbegin + k - 1 ) =
745 $ work( wbegin + k - 1) - tau
747 $ four*eps*abs(work(wbegin+k-1))
749 werr( wbegin + k - 1 ) =
750 $ werr( wbegin + k - 1 ) + fudge
762 iwork( k-1 ) = newfst
774 tol = four * log(dble(in)) * eps
777 windex = wbegin + k - 1
778 windmn = max(windex - 1,1)
779 windpl = min(windex + 1,m)
780 lambda = work( windex )
783 IF((windex.LT.dol).OR.
784 $ (windex.GT.dou))
THEN
790 left = work( windex ) - werr( windex )
791 right = work( windex ) + werr( windex )
792 indeig = indexw( windex )
807 lgap = eps*max(abs(left),abs(right))
817 rgap = eps*max(abs(left),abs(right))
821 gap = min( lgap, rgap )
822 IF(( k .EQ. 1).OR.(k .EQ. im))
THEN
837 savgap = wgap(windex)
854 itmp1 = iwork( iindr+windex )
855 offset = indexw( wbegin ) - 1
856 CALL
dlarrb( in, d(ibegin),
857 $ work(indlld+ibegin-1),indeig,indeig,
858 $ zero, two*eps, offset,
859 $ work(wbegin),wgap(wbegin),
860 $ werr(wbegin),work( indwrk ),
861 $ iwork( iindwk ), pivmin, spdiam,
863 IF( iinfo.NE.0 )
THEN
867 lambda = work( windex )
870 iwork( iindr+windex ) = 0
873 CALL
zlar1v( in, 1, in, lambda, d( ibegin ),
874 $ l( ibegin ), work(indld+ibegin-1),
875 $ work(indlld+ibegin-1),
876 $ pivmin, gaptol, z( ibegin, windex ),
877 $ .NOT.usedbs, negcnt, ztz, mingma,
878 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
879 $ nrminv, resid, rqcorr, work( indwrk ) )
883 ELSEIF(resid.LT.bstres)
THEN
887 isupmn = min(isupmn,isuppz( 2*windex-1 ))
888 isupmx = max(isupmx,isuppz( 2*windex ))
900 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
901 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
906 IF(indeig.LE.negcnt)
THEN
915 IF( ( rqcorr*sgndef.GE.zero )
916 $ .AND.( lambda + rqcorr.LE. right)
917 $ .AND.( lambda + rqcorr.GE. left)
921 IF(sgndef.EQ.one)
THEN
940 $ half * (right + left)
943 lambda = lambda + rqcorr
946 $ half * (right-left)
950 IF(right-left.LT.rqtol*abs(lambda))
THEN
955 ELSEIF( iter.LT.maxitr )
THEN
957 ELSEIF( iter.EQ.maxitr )
THEN
966 IF(usedrq .AND. usedbs .AND.
967 $ bstres.LE.resid)
THEN
973 CALL
zlar1v( in, 1, in, lambda,
974 $ d( ibegin ), l( ibegin ),
975 $ work(indld+ibegin-1),
976 $ work(indlld+ibegin-1),
977 $ pivmin, gaptol, z( ibegin, windex ),
978 $ .NOT.usedbs, negcnt, ztz, mingma,
979 $ iwork( iindr+windex ),
980 $ isuppz( 2*windex-1 ),
981 $ nrminv, resid, rqcorr, work( indwrk ) )
983 work( windex ) = lambda
988 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
989 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
990 zfrom = isuppz( 2*windex-1 )
991 zto = isuppz( 2*windex )
992 isupmn = isupmn + oldien
993 isupmx = isupmx + oldien
995 IF(isupmn.LT.zfrom)
THEN
996 DO 122 ii = isupmn,zfrom-1
997 z( ii, windex ) = zero
1000 IF(isupmx.GT.zto)
THEN
1001 DO 123 ii = zto+1,isupmx
1002 z( ii, windex ) = zero
1005 CALL
zdscal( zto-zfrom+1, nrminv,
1006 $ z( zfrom, windex ), 1 )
1009 w( windex ) = lambda+sigma
1018 wgap( windmn ) = max( wgap(windmn),
1019 $ w(windex)-werr(windex)
1020 $ - w(windmn)-werr(windmn) )
1022 IF( windex.LT.wend )
THEN
1023 wgap( windex ) = max( savgap,
1024 $ w( windpl )-werr( windpl )
1025 $ - w( windex )-werr( windex) )
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...
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
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 zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 ...
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH