295 SUBROUTINE dlarre( 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 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
311 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
313 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
314 $ w( * ),werr( * ), wgap( * ), work( * )
320 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
321 $ maxgrowth, one, pert, two, zero
322 parameter( zero = 0.0d0, one = 1.0d0,
323 $ two = 2.0d0, four=4.0d0,
326 $ half = one/two, fourth = one/four, fac= half,
327 $ maxgrowth = 64.0d0, fudge = 2.0d0 )
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 DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
338 $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
339 $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
349 DOUBLE PRECISION DLAMCH
350 EXTERNAL dlamch, 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 = dlamch(
'S' )
389 IF( (irange.EQ.allrng).OR.
390 $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
391 $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) )
THEN
420 IF( eabs .GE. emax )
THEN
424 gers( 2*i-1) = d(i) - tmp1
425 gl = min( gl, gers( 2*i - 1))
426 gers( 2*i ) = d(i) + tmp1
427 gu = max( gu, gers(2*i) )
431 pivmin = safmin * max( one, emax**2 )
437 CALL dlarra( n, d, e, e2, spltol, spdiam,
438 $ nsplit, isplit, iinfo )
446 usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
448 IF( (irange.EQ.allrng) .AND. (.NOT. forceb) )
THEN
459 CALL dlarrd( range,
'B', n, vl, vu, il, iu, gers,
460 $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
461 $ mm, w, werr, vl, vu, iblock, indexw,
462 $ work, iwork, iinfo )
463 IF( iinfo.NE.0 )
THEN
481 DO 170 jblk = 1, nsplit
482 iend = isplit( jblk )
483 in = iend - ibegin + 1
487 IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
488 $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
489 $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
515 DO 15 i = ibegin , iend
516 gl = min( gers( 2*i-1 ), gl )
517 gu = max( gers( 2*i ), gu )
521 IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) )
THEN
525 IF( iblock(i).EQ.jblk )
THEN
542 usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
543 wend = wbegin + mb - 1
548 DO 30 i = wbegin, wend - 1
549 wgap( i ) = max( zero,
550 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
552 wgap( wend ) = max( zero,
553 $ vu - sigma - (w( wend )+werr( wend )))
555 indl = indexw(wbegin)
556 indu = indexw( wend )
559 IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd)
THEN
562 CALL dlarrk( in, 1, gl, gu, d(ibegin),
563 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
564 IF( iinfo.NE.0 )
THEN
568 isleft = max(gl, tmp - tmp1
569 $ - hndrd * eps* abs(tmp - tmp1))
571 CALL dlarrk( in, in, gl, gu, d(ibegin),
572 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
573 IF( iinfo.NE.0 )
THEN
577 isrght = min(gu, tmp + tmp1
578 $ + hndrd * eps * abs(tmp + tmp1))
580 spdiam = isrght - isleft
584 isleft = max(gl, w(wbegin) - werr(wbegin)
585 $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
586 isrght = min(gu,w(wend) + werr(wend)
587 $ + hndrd * eps * abs(w(wend)+ werr(wend)))
599 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
607 wend = wbegin + mb - 1
609 s1 = isleft + fourth * spdiam
610 s2 = isrght - fourth * spdiam
616 s1 = isleft + fourth * spdiam
617 s2 = isrght - fourth * spdiam
619 tmp = min(isrght,vu) - max(isleft,vl)
620 s1 = max(isleft,vl) + fourth * tmp
621 s2 = min(isrght,vu) - fourth * tmp
627 CALL dlarrc(
'T', in, s1, s2, d(ibegin),
628 $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
634 ELSEIF( cnt1 - indl .GE. indu - cnt2 )
THEN
635 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
636 sigma = max(isleft,gl)
637 ELSEIF( usedqd )
THEN
644 sigma = max(isleft,vl)
648 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
649 sigma = min(isrght,gu)
650 ELSEIF( usedqd )
THEN
657 sigma = min(isrght,vu)
671 tau = spdiam*eps*n + two*pivmin
672 tau = max( tau,two*eps*abs(sigma) )
675 clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
676 avgap = abs(clwdth / dble(wend-wbegin))
677 IF( sgndef.EQ.one )
THEN
678 tau = half*max(wgap(wbegin),avgap)
679 tau = max(tau,werr(wbegin))
681 tau = half*max(wgap(wend-1),avgap)
682 tau = max(tau,werr(wend))
689 DO 80 idum = 1, maxtry
693 dpivot = d( ibegin ) - sigma
695 dmax = abs( work(1) )
698 work( 2*in+i ) = one / work( i )
699 tmp = e( j )*work( 2*in+i )
701 dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
703 dmax = max( dmax, abs(dpivot) )
707 IF( dmax .GT. maxgrowth*spdiam )
THEN
712 IF( usedqd .AND. .NOT.norep )
THEN
716 tmp = sgndef*work( i )
717 IF( tmp.LT.zero ) norep = .true.
724 IF( idum.EQ.maxtry-1 )
THEN
725 IF( sgndef.EQ.one )
THEN
728 $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
731 $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
734 sigma = sigma - sgndef * tau
753 CALL dcopy( in, work, 1, d( ibegin ), 1 )
754 CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
767 CALL dlarnv(2, iseed, 2*in-1, work(1))
769 d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
770 e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
772 d(iend) = d(iend)*(one+eps*four*work(in))
782 IF ( .NOT.usedqd )
THEN
790 werr(j) = werr(j) + abs(w(j)) * eps
794 DO 135 i = ibegin, iend-1
795 work( i ) = d( i ) * e( i )**2
798 CALL dlarrb(in, d(ibegin), work(ibegin),
799 $ indl, indu, rtol1, rtol2, indl-1,
800 $ w(wbegin), wgap(wbegin), werr(wbegin),
801 $ work( 2*n+1 ), iwork, pivmin, spdiam,
803 IF( iinfo .NE. 0 )
THEN
809 wgap( wend ) = max( zero,
810 $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
811 DO 138 i = indl, indu
828 rtol = log(dble(in)) * four * eps
831 work( 2*i-1 ) = abs( d( j ) )
832 work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
835 work( 2*in-1 ) = abs( d( iend ) )
837 CALL dlasq2( in, work, iinfo )
838 IF( iinfo .NE. 0 )
THEN
847 IF( work( i ).LT.zero )
THEN
853 IF( sgndef.GT.zero )
THEN
854 DO 150 i = indl, indu
856 w( m ) = work( in-i+1 )
861 DO 160 i = indl, indu
869 DO 165 i = m - mb + 1, m
871 werr( i ) = rtol * abs( w(i) )
873 DO 166 i = m - mb + 1, m - 1
875 wgap( i ) = max( zero,
876 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
878 wgap( m ) = max( zero,
879 $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
subroutine dlarrk(N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
subroutine dlarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
DLARRA computes the splitting points with the specified threshold.
subroutine dlarrd(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)
DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
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.
subroutine dlarre(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)
DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
subroutine dlarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.