295 SUBROUTINE slarre( RANGE, N, VL, VU, IL, IU, D, E, E2,
296 $ rtol1, rtol2, spltol, nsplit, isplit, m,
297 $ w, werr, wgap, iblock, indexw, gers, pivmin,
298 $ work, iwork, info )
307 INTEGER IL, INFO, IU, M, N, NSPLIT
308 REAL PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
311 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
313 REAL D( * ), E( * ), E2( * ), GERS( * ),
314 $ w( * ),werr( * ), wgap( * ), work( * )
320 REAL FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
321 $ maxgrowth, one, pert, two, zero
322 parameter( zero = 0.0e0, one = 1.0e0,
323 $ two = 2.0e0, four=4.0e0,
326 $ half = one/two, fourth = one/four, fac= half,
327 $ maxgrowth = 64.0e0, fudge = 2.0e0 )
328 INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
329 parameter( maxtry = 6, allrng = 1, indrng = 2,
333 LOGICAL FORCEB, NOREP, USEDQD
334 INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
335 $ in, indl, indu, irange, j, jblk, mb, mm,
337 REAL AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
338 $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
339 $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
350 EXTERNAL slamch, lsame
358 INTRINSIC abs, max, min
369 IF( lsame( range,
'A' ) )
THEN
371 ELSE IF( lsame( range,
'V' ) )
THEN
373 ELSE IF( lsame( range,
'I' ) )
THEN
380 safmin = slamch(
'S' )
389 bsrtol = sqrt(eps)*(0.5e-3)
393 IF( (irange.EQ.allrng).OR.
394 $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
395 $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) )
THEN
424 IF( eabs .GE. emax )
THEN
428 gers( 2*i-1) = d(i) - tmp1
429 gl = min( gl, gers( 2*i - 1))
430 gers( 2*i ) = d(i) + tmp1
431 gu = max( gu, gers(2*i) )
435 pivmin = safmin * max( one, emax**2 )
441 CALL slarra( n, d, e, e2, spltol, spdiam,
442 $ nsplit, isplit, iinfo )
450 usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
452 IF( (irange.EQ.allrng) .AND. (.NOT. forceb) )
THEN
463 CALL slarrd( range,
'B', n, vl, vu, il, iu, gers,
464 $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
465 $ mm, w, werr, vl, vu, iblock, indexw,
466 $ work, iwork, iinfo )
467 IF( iinfo.NE.0 )
THEN
485 DO 170 jblk = 1, nsplit
486 iend = isplit( jblk )
487 in = iend - ibegin + 1
491 IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
492 $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
493 $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
519 DO 15 i = ibegin , iend
520 gl = min( gers( 2*i-1 ), gl )
521 gu = max( gers( 2*i ), gu )
525 IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) )
THEN
529 IF( iblock(i).EQ.jblk )
THEN
546 usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
547 wend = wbegin + mb - 1
552 DO 30 i = wbegin, wend - 1
553 wgap( i ) = max( zero,
554 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
556 wgap( wend ) = max( zero,
557 $ vu - sigma - (w( wend )+werr( wend )))
559 indl = indexw(wbegin)
560 indu = indexw( wend )
563 IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd)
THEN
566 CALL slarrk( in, 1, gl, gu, d(ibegin),
567 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
568 IF( iinfo.NE.0 )
THEN
572 isleft = max(gl, tmp - tmp1
573 $ - hndrd * eps* abs(tmp - tmp1))
575 CALL slarrk( in, in, gl, gu, d(ibegin),
576 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
577 IF( iinfo.NE.0 )
THEN
581 isrght = min(gu, tmp + tmp1
582 $ + hndrd * eps * abs(tmp + tmp1))
584 spdiam = isrght - isleft
588 isleft = max(gl, w(wbegin) - werr(wbegin)
589 $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
590 isrght = min(gu,w(wend) + werr(wend)
591 $ + hndrd * eps * abs(w(wend)+ werr(wend)))
603 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
611 wend = wbegin + mb - 1
613 s1 = isleft + fourth * spdiam
614 s2 = isrght - fourth * spdiam
620 s1 = isleft + fourth * spdiam
621 s2 = isrght - fourth * spdiam
623 tmp = min(isrght,vu) - max(isleft,vl)
624 s1 = max(isleft,vl) + fourth * tmp
625 s2 = min(isrght,vu) - fourth * tmp
631 CALL slarrc(
'T', in, s1, s2, d(ibegin),
632 $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
638 ELSEIF( cnt1 - indl .GE. indu - cnt2 )
THEN
639 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
640 sigma = max(isleft,gl)
641 ELSEIF( usedqd )
THEN
648 sigma = max(isleft,vl)
652 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
653 sigma = min(isrght,gu)
654 ELSEIF( usedqd )
THEN
661 sigma = min(isrght,vu)
675 tau = spdiam*eps*n + two*pivmin
676 tau = max( tau,two*eps*abs(sigma) )
679 clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
680 avgap = abs(clwdth /
REAL(wend-wbegin))
681 IF( sgndef.EQ.one )
THEN
682 tau = half*max(wgap(wbegin),avgap)
683 tau = max(tau,werr(wbegin))
685 tau = half*max(wgap(wend-1),avgap)
686 tau = max(tau,werr(wend))
693 DO 80 idum = 1, maxtry
697 dpivot = d( ibegin ) - sigma
699 dmax = abs( work(1) )
702 work( 2*in+i ) = one / work( i )
703 tmp = e( j )*work( 2*in+i )
705 dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
707 dmax = max( dmax, abs(dpivot) )
711 IF( dmax .GT. maxgrowth*spdiam )
THEN
716 IF( usedqd .AND. .NOT.norep )
THEN
720 tmp = sgndef*work( i )
721 IF( tmp.LT.zero ) norep = .true.
728 IF( idum.EQ.maxtry-1 )
THEN
729 IF( sgndef.EQ.one )
THEN
732 $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
735 $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
738 sigma = sigma - sgndef * tau
757 CALL scopy( in, work, 1, d( ibegin ), 1 )
758 CALL scopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
771 CALL slarnv(2, iseed, 2*in-1, work(1))
773 d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
774 e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
776 d(iend) = d(iend)*(one+eps*four*work(in))
786 IF ( .NOT.usedqd )
THEN
794 werr(j) = werr(j) + abs(w(j)) * eps
798 DO 135 i = ibegin, iend-1
799 work( i ) = d( i ) * e( i )**2
802 CALL slarrb(in, d(ibegin), work(ibegin),
803 $ indl, indu, rtol1, rtol2, indl-1,
804 $ w(wbegin), wgap(wbegin), werr(wbegin),
805 $ work( 2*n+1 ), iwork, pivmin, spdiam,
807 IF( iinfo .NE. 0 )
THEN
813 wgap( wend ) = max( zero,
814 $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
815 DO 138 i = indl, indu
832 rtol = log(
REAL(in)) * four * eps
835 work( 2*i-1 ) = abs( d( j ) )
836 work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
839 work( 2*in-1 ) = abs( d( iend ) )
841 CALL slasq2( in, work, iinfo )
842 IF( iinfo .NE. 0 )
THEN
851 IF( work( i ).LT.zero )
THEN
857 IF( sgndef.GT.zero )
THEN
858 DO 150 i = indl, indu
860 w( m ) = work( in-i+1 )
865 DO 160 i = indl, indu
873 DO 165 i = m - mb + 1, m
875 werr( i ) = rtol * abs( w(i) )
877 DO 166 i = m - mb + 1, m - 1
879 wgap( i ) = max( zero,
880 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
882 wgap( m ) = max( zero,
883 $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
subroutine slasq2(N, Z, INFO)
SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
subroutine slarrd(RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
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 slarrk(N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
SLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
subroutine slarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
SLARRA computes the splitting points with the specified threshold.
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine slarre(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
subroutine slarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.