303 SUBROUTINE dlarre( RANGE, N, VL, VU, IL, IU, D, E, E2,
304 $ rtol1, rtol2, spltol, nsplit, isplit, m,
305 $ w, werr, wgap, iblock, indexw, gers, pivmin,
306 $ work, iwork, info )
315 INTEGER IL, INFO, IU, M, N, NSPLIT
316 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
319 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
321 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
322 $ w( * ),werr( * ), wgap( * ), work( * )
328 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
329 $ maxgrowth, one, pert, two, zero
330 parameter ( zero = 0.0d0, one = 1.0d0,
331 $ two = 2.0d0, four=4.0d0,
334 $ half = one/two, fourth = one/four, fac= half,
335 $ maxgrowth = 64.0d0, fudge = 2.0d0 )
336 INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
337 parameter ( maxtry = 6, allrng = 1, indrng = 2,
341 LOGICAL FORCEB, NOREP, USEDQD
342 INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
343 $ in, indl, indu, irange, j, jblk, mb, mm,
345 DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
346 $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
347 $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
357 DOUBLE PRECISION DLAMCH
358 EXTERNAL dlamch, lsame
366 INTRINSIC abs, max, min
377 IF( lsame( range,
'A' ) )
THEN
379 ELSE IF( lsame( range,
'V' ) )
THEN
381 ELSE IF( lsame( range,
'I' ) )
THEN
388 safmin = dlamch(
'S' )
397 IF( (irange.EQ.allrng).OR.
398 $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
399 $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) )
THEN
428 IF( eabs .GE. emax )
THEN
432 gers( 2*i-1) = d(i) - tmp1
433 gl = min( gl, gers( 2*i - 1))
434 gers( 2*i ) = d(i) + tmp1
435 gu = max( gu, gers(2*i) )
439 pivmin = safmin * max( one, emax**2 )
445 CALL dlarra( n, d, e, e2, spltol, spdiam,
446 $ nsplit, isplit, iinfo )
454 usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
456 IF( (irange.EQ.allrng) .AND. (.NOT. forceb) )
THEN
467 CALL dlarrd( range,
'B', n, vl, vu, il, iu, gers,
468 $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
469 $ mm, w, werr, vl, vu, iblock, indexw,
470 $ work, iwork, iinfo )
471 IF( iinfo.NE.0 )
THEN
489 DO 170 jblk = 1, nsplit
490 iend = isplit( jblk )
491 in = iend - ibegin + 1
495 IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
496 $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
497 $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
523 DO 15 i = ibegin , iend
524 gl = min( gers( 2*i-1 ), gl )
525 gu = max( gers( 2*i ), gu )
529 IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) )
THEN
533 IF( iblock(i).EQ.jblk )
THEN
550 usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
551 wend = wbegin + mb - 1
556 DO 30 i = wbegin, wend - 1
557 wgap( i ) = max( zero,
558 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
560 wgap( wend ) = max( zero,
561 $ vu - sigma - (w( wend )+werr( wend )))
563 indl = indexw(wbegin)
564 indu = indexw( wend )
567 IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd)
THEN
570 CALL dlarrk( in, 1, gl, gu, d(ibegin),
571 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
572 IF( iinfo.NE.0 )
THEN
576 isleft = max(gl, tmp - tmp1
577 $ - hndrd * eps* abs(tmp - tmp1))
579 CALL dlarrk( in, in, gl, gu, d(ibegin),
580 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
581 IF( iinfo.NE.0 )
THEN
585 isrght = min(gu, tmp + tmp1
586 $ + hndrd * eps * abs(tmp + tmp1))
588 spdiam = isrght - isleft
592 isleft = max(gl, w(wbegin) - werr(wbegin)
593 $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
594 isrght = min(gu,w(wend) + werr(wend)
595 $ + hndrd * eps * abs(w(wend)+ werr(wend)))
607 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
615 wend = wbegin + mb - 1
617 s1 = isleft + fourth * spdiam
618 s2 = isrght - fourth * spdiam
624 s1 = isleft + fourth * spdiam
625 s2 = isrght - fourth * spdiam
627 tmp = min(isrght,vu) - max(isleft,vl)
628 s1 = max(isleft,vl) + fourth * tmp
629 s2 = min(isrght,vu) - fourth * tmp
635 CALL dlarrc(
'T', in, s1, s2, d(ibegin),
636 $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
642 ELSEIF( cnt1 - indl .GE. indu - cnt2 )
THEN
643 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
644 sigma = max(isleft,gl)
645 ELSEIF( usedqd )
THEN
652 sigma = max(isleft,vl)
656 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
657 sigma = min(isrght,gu)
658 ELSEIF( usedqd )
THEN
665 sigma = min(isrght,vu)
679 tau = spdiam*eps*n + two*pivmin
680 tau = max( tau,two*eps*abs(sigma) )
683 clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
684 avgap = abs(clwdth / dble(wend-wbegin))
685 IF( sgndef.EQ.one )
THEN
686 tau = half*max(wgap(wbegin),avgap)
687 tau = max(tau,werr(wbegin))
689 tau = half*max(wgap(wend-1),avgap)
690 tau = max(tau,werr(wend))
697 DO 80 idum = 1, maxtry
701 dpivot = d( ibegin ) - sigma
703 dmax = abs( work(1) )
706 work( 2*in+i ) = one / work( i )
707 tmp = e( j )*work( 2*in+i )
709 dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
711 dmax = max( dmax, abs(dpivot) )
715 IF( dmax .GT. maxgrowth*spdiam )
THEN
720 IF( usedqd .AND. .NOT.norep )
THEN
724 tmp = sgndef*work( i )
725 IF( tmp.LT.zero ) norep = .true.
732 IF( idum.EQ.maxtry-1 )
THEN
733 IF( sgndef.EQ.one )
THEN
736 $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
739 $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
742 sigma = sigma - sgndef * tau
761 CALL dcopy( in, work, 1, d( ibegin ), 1 )
762 CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
775 CALL dlarnv(2, iseed, 2*in-1, work(1))
777 d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
778 e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
780 d(iend) = d(iend)*(one+eps*four*work(in))
790 IF ( .NOT.usedqd )
THEN
798 werr(j) = werr(j) + abs(w(j)) * eps
802 DO 135 i = ibegin, iend-1
803 work( i ) = d( i ) * e( i )**2
806 CALL dlarrb(in, d(ibegin), work(ibegin),
807 $ indl, indu, rtol1, rtol2, indl-1,
808 $ w(wbegin), wgap(wbegin), werr(wbegin),
809 $ work( 2*n+1 ), iwork, pivmin, spdiam,
811 IF( iinfo .NE. 0 )
THEN
817 wgap( wend ) = max( zero,
818 $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
819 DO 138 i = indl, indu
836 rtol = log(dble(in)) * four * eps
839 work( 2*i-1 ) = abs( d( j ) )
840 work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
843 work( 2*in-1 ) = abs( d( iend ) )
845 CALL dlasq2( in, work, iinfo )
846 IF( iinfo .NE. 0 )
THEN
855 IF( work( i ).LT.zero )
THEN
861 IF( sgndef.GT.zero )
THEN
862 DO 150 i = indl, indu
864 w( m ) = work( in-i+1 )
869 DO 160 i = indl, indu
877 DO 165 i = m - mb + 1, m
879 werr( i ) = rtol * abs( w(i) )
881 DO 166 i = m - mb + 1, m - 1
883 wgap( i ) = max( zero,
884 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
886 wgap( m ) = max( zero,
887 $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
subroutine dlarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
DLARRA computes the splitting points with the specified threshold.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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 dlarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
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 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.