283 SUBROUTINE clarrv( 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 )
312 parameter ( czero = ( 0.0e0, 0.0e0 ) )
313 REAL ZERO, ONE, TWO, THREE, FOUR, HALF
314 parameter ( zero = 0.0e0, one = 1.0e0,
315 $ two = 2.0e0, three = 3.0e0,
316 $ four = 4.0e0, half = 0.5e0)
319 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ
320 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1,
321 $ iindc2, iindr, iindwk, iinfo, im, in, indeig,
322 $ indld, indlld, indwrk, isupmn, isupmx, iter,
323 $ itmp1, j, jblk, k, miniwsize, minwsize, nclus,
324 $ ndepth, negcnt, newcls, newfst, newftt, newlst,
325 $ newsiz, offset, oldcls, oldfst, oldien, oldlst,
326 $ oldncl, p, parity, q, wbegin, wend, windex,
327 $ windmn, windpl, zfrom, zto, zusedl, zusedu,
329 INTEGER INDIN1, INDIN2
330 REAL 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
344 INTRINSIC abs,
REAL, MAX, MIN
388 zusedw = zusedu - zusedl + 1
391 CALL claset(
'Full', n, zusedw, czero, czero,
394 eps = slamch(
'Precision' )
400 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
419 DO 170 jblk = 1, iblock( m )
420 iend = isplit( jblk )
427 IF( iblock( wend+1 ).EQ.jblk )
THEN
432 IF( wend.LT.wbegin )
THEN
435 ELSEIF( (wend.LT.dol).OR.(wbegin.GT.dou) )
THEN
442 gl = gers( 2*ibegin-1 )
443 gu = gers( 2*ibegin )
444 DO 20 i = ibegin+1 , iend
445 gl = min( gers( 2*i-1 ), gl )
446 gu = max( gers( 2*i ), gu )
453 in = iend - ibegin + 1
455 im = wend - wbegin + 1
458 IF( ibegin.EQ.iend )
THEN
460 z( ibegin, wbegin ) = cmplx( one, zero )
461 isuppz( 2*wbegin-1 ) = ibegin
462 isuppz( 2*wbegin ) = ibegin
463 w( wbegin ) = w( wbegin ) + sigma
464 work( wbegin ) = w( wbegin )
476 CALL scopy( im, w( wbegin ), 1,
477 $ work( wbegin ), 1 )
482 w(wbegin+i-1) = w(wbegin+i-1)+sigma
493 iwork( iindc1+1 ) = 1
494 iwork( iindc1+2 ) = im
503 IF( idone.LT.im )
THEN
505 IF( ndepth.GT.m )
THEN
516 IF( parity.EQ.0 )
THEN
529 oldfst = iwork( j-1 )
531 IF( ndepth.GT.0 )
THEN
537 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
540 j = wbegin + oldfst - 1
542 IF(wbegin+oldfst-1.LT.dol)
THEN
545 ELSEIF(wbegin+oldfst-1.GT.dou)
THEN
549 j = wbegin + oldfst - 1
553 d( ibegin+k-1 ) =
REAL( Z( IBEGIN+K-1,
$ J )
554 REAL( Z( IBEGIN+K-1,
$ J+1 )
556 d( iend ) =
REAL( Z( IEND, J ) )
557 sigma =
REAL( Z( IEND, J+1 ) )
560 CALL claset(
'Full', in, 2, czero, czero,
561 $ z( ibegin, j), ldz )
565 DO 50 j = ibegin, iend-1
567 work( indld-1+j ) = tmp
568 work( indlld-1+j ) = tmp*l( j )
571 IF( ndepth.GT.0 )
THEN
574 p = indexw( wbegin-1+oldfst )
575 q = indexw( wbegin-1+oldlst )
579 offset = indexw( wbegin ) - 1
582 CALL slarrb( in, d( ibegin ),
583 $ work(indlld+ibegin-1),
584 $ p, q, rtol1, rtol2, offset,
585 $ work(wbegin),wgap(wbegin),werr(wbegin),
586 $ work( indwrk ), iwork( iindwk ),
587 $ pivmin, spdiam, in, iinfo )
588 IF( iinfo.NE.0 )
THEN
599 IF( oldfst.GT.1)
THEN
600 wgap( wbegin+oldfst-2 ) =
601 $ max(wgap(wbegin+oldfst-2),
602 $ w(wbegin+oldfst-1)-werr(wbegin+oldfst-1)
603 $ - w(wbegin+oldfst-2)-werr(wbegin+oldfst-2) )
605 IF( wbegin + oldlst -1 .LT. wend )
THEN
606 wgap( wbegin+oldlst-1 ) =
607 $ max(wgap(wbegin+oldlst-1),
608 $ w(wbegin+oldlst)-werr(wbegin+oldlst)
609 $ - w(wbegin+oldlst-1)-werr(wbegin+oldlst-1) )
613 DO 53 j=oldfst,oldlst
614 w(wbegin+j-1) = work(wbegin+j-1)+sigma
620 DO 140 j = oldfst, oldlst
621 IF( j.EQ.oldlst )
THEN
625 ELSE IF ( wgap( wbegin + j -1).GE.
626 $ minrgp* abs( work(wbegin + j -1) ) )
THEN
637 newsiz = newlst - newfst + 1
641 IF((dol.EQ.1).AND.(dou.EQ.m))
THEN
644 newftt = wbegin + newfst - 1
646 IF(wbegin+newfst-1.LT.dol)
THEN
649 ELSEIF(wbegin+newfst-1.GT.dou)
THEN
653 newftt = wbegin + newfst - 1
657 IF( newsiz.GT.1)
THEN
672 IF( newfst.EQ.1 )
THEN
674 $ w(wbegin)-werr(wbegin) - vl )
676 lgap = wgap( wbegin+newfst-2 )
678 rgap = wgap( wbegin+newlst-1 )
687 p = indexw( wbegin-1+newfst )
689 p = indexw( wbegin-1+newlst )
691 offset = indexw( wbegin ) - 1
692 CALL slarrb( in, d(ibegin),
693 $ work( indlld+ibegin-1 ),p,p,
694 $ rqtol, rqtol, offset,
695 $ work(wbegin),wgap(wbegin),
696 $ werr(wbegin),work( indwrk ),
697 $ iwork( iindwk ), pivmin, spdiam,
701 IF((wbegin+newlst-1.LT.dol).OR.
702 $ (wbegin+newfst-1.GT.dou))
THEN
710 idone = idone + newlst - newfst + 1
718 CALL slarrf( in, d( ibegin ), l( ibegin ),
719 $ work(indld+ibegin-1),
720 $ newfst, newlst, work(wbegin),
721 $ wgap(wbegin), werr(wbegin),
722 $ spdiam, lgap, rgap, pivmin, tau,
723 $ work( indin1 ), work( indin2 ),
724 $ work( indwrk ), iinfo )
729 z( ibegin+k-1, newftt ) =
730 $ cmplx( work( indin1+k-1 ), zero )
731 z( ibegin+k-1, newftt+1 ) =
732 $ cmplx( work( indin2+k-1 ), zero )
735 $ cmplx( work( indin1+in-1 ), zero )
736 IF( iinfo.EQ.0 )
THEN
740 z( iend, newftt+1 ) = cmplx( ssigma, zero )
743 DO 116 k = newfst, newlst
745 $ three*eps*abs(work(wbegin+k-1))
746 work( wbegin + k - 1 ) =
747 $ work( wbegin + k - 1) - tau
749 $ four*eps*abs(work(wbegin+k-1))
751 werr( wbegin + k - 1 ) =
752 $ werr( wbegin + k - 1 ) + fudge
764 iwork( k-1 ) = newfst
776 tol = four * log(
REAL(in)) * eps
779 windex = wbegin + k - 1
780 windmn = max(windex - 1,1)
781 windpl = min(windex + 1,m)
782 lambda = work( windex )
785 IF((windex.LT.dol).OR.
786 $ (windex.GT.dou))
THEN
792 left = work( windex ) - werr( windex )
793 right = work( windex ) + werr( windex )
794 indeig = indexw( windex )
809 lgap = eps*max(abs(left),abs(right))
819 rgap = eps*max(abs(left),abs(right))
823 gap = min( lgap, rgap )
824 IF(( k .EQ. 1).OR.(k .EQ. im))
THEN
839 savgap = wgap(windex)
856 itmp1 = iwork( iindr+windex )
857 offset = indexw( wbegin ) - 1
858 CALL slarrb( in, d(ibegin),
859 $ work(indlld+ibegin-1),indeig,indeig,
860 $ zero, two*eps, offset,
861 $ work(wbegin),wgap(wbegin),
862 $ werr(wbegin),work( indwrk ),
863 $ iwork( iindwk ), pivmin, spdiam,
865 IF( iinfo.NE.0 )
THEN
869 lambda = work( windex )
872 iwork( iindr+windex ) = 0
875 CALL clar1v( in, 1, in, lambda, d( ibegin ),
876 $ l( ibegin ), work(indld+ibegin-1),
877 $ work(indlld+ibegin-1),
878 $ pivmin, gaptol, z( ibegin, windex ),
879 $ .NOT.usedbs, negcnt, ztz, mingma,
880 $ iwork( iindr+windex ), isuppz( 2*windex-1 ),
881 $ nrminv, resid, rqcorr, work( indwrk ) )
885 ELSEIF(resid.LT.bstres)
THEN
889 isupmn = min(isupmn,isuppz( 2*windex-1 ))
890 isupmx = max(isupmx,isuppz( 2*windex ))
902 IF( resid.GT.tol*gap .AND. abs( rqcorr ).GT.
903 $ rqtol*abs( lambda ) .AND. .NOT. usedbs)
908 IF(indeig.LE.negcnt)
THEN
917 IF( ( rqcorr*sgndef.GE.zero )
918 $ .AND.( lambda + rqcorr.LE. right)
919 $ .AND.( lambda + rqcorr.GE. left)
923 IF(sgndef.EQ.one)
THEN
942 $ half * (right + left)
945 lambda = lambda + rqcorr
948 $ half * (right-left)
952 IF(right-left.LT.rqtol*abs(lambda))
THEN
957 ELSEIF( iter.LT.maxitr )
THEN
959 ELSEIF( iter.EQ.maxitr )
THEN
968 IF(usedrq .AND. usedbs .AND.
969 $ bstres.LE.resid)
THEN
975 CALL clar1v( in, 1, in, lambda,
976 $ d( ibegin ), l( ibegin ),
977 $ work(indld+ibegin-1),
978 $ work(indlld+ibegin-1),
979 $ pivmin, gaptol, z( ibegin, windex ),
980 $ .NOT.usedbs, negcnt, ztz, mingma,
981 $ iwork( iindr+windex ),
982 $ isuppz( 2*windex-1 ),
983 $ nrminv, resid, rqcorr, work( indwrk ) )
985 work( windex ) = lambda
990 isuppz( 2*windex-1 ) = isuppz( 2*windex-1 )+oldien
991 isuppz( 2*windex ) = isuppz( 2*windex )+oldien
992 zfrom = isuppz( 2*windex-1 )
993 zto = isuppz( 2*windex )
994 isupmn = isupmn + oldien
995 isupmx = isupmx + oldien
997 IF(isupmn.LT.zfrom)
THEN
998 DO 122 ii = isupmn,zfrom-1
999 z( ii, windex ) = zero
1002 IF(isupmx.GT.zto)
THEN
1003 DO 123 ii = zto+1,isupmx
1004 z( ii, windex ) = zero
1007 CALL csscal( zto-zfrom+1, nrminv,
1008 $ z( zfrom, windex ), 1 )
1011 w( windex ) = lambda+sigma
1020 wgap( windmn ) = max( wgap(windmn),
1021 $ w(windex)-werr(windex)
1022 $ - w(windmn)-werr(windmn) )
1024 IF( windex.LT.wend )
THEN
1025 wgap( windex ) = max( savgap,
1026 $ w( windpl )-werr( windpl )
1027 $ - w( windex )-werr( windex) )
1052 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 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...
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 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...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine csscal(N, SA, CX, INCX)
CSSCAL