303 SUBROUTINE slarre( 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 REAL PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
319 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
321 REAL D( * ), E( * ), E2( * ), GERS( * ),
322 $ w( * ),werr( * ), wgap( * ), work( * )
328 REAL FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
329 $ maxgrowth, one, pert, two, zero
330 parameter ( zero = 0.0e0, one = 1.0e0,
331 $ two = 2.0e0, four=4.0e0,
334 $ half = one/two, fourth = one/four, fac= half,
335 $ maxgrowth = 64.0e0, fudge = 2.0e0 )
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 REAL AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
346 $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
347 $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
358 EXTERNAL slamch, 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 = slamch(
'S' )
397 bsrtol = sqrt(eps)*(0.5e-3)
401 IF( (irange.EQ.allrng).OR.
402 $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
403 $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) )
THEN
432 IF( eabs .GE. emax )
THEN
436 gers( 2*i-1) = d(i) - tmp1
437 gl = min( gl, gers( 2*i - 1))
438 gers( 2*i ) = d(i) + tmp1
439 gu = max( gu, gers(2*i) )
443 pivmin = safmin * max( one, emax**2 )
449 CALL slarra( n, d, e, e2, spltol, spdiam,
450 $ nsplit, isplit, iinfo )
458 usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
460 IF( (irange.EQ.allrng) .AND. (.NOT. forceb) )
THEN
471 CALL slarrd( range,
'B', n, vl, vu, il, iu, gers,
472 $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
473 $ mm, w, werr, vl, vu, iblock, indexw,
474 $ work, iwork, iinfo )
475 IF( iinfo.NE.0 )
THEN
493 DO 170 jblk = 1, nsplit
494 iend = isplit( jblk )
495 in = iend - ibegin + 1
499 IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
500 $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
501 $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
527 DO 15 i = ibegin , iend
528 gl = min( gers( 2*i-1 ), gl )
529 gu = max( gers( 2*i ), gu )
533 IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) )
THEN
537 IF( iblock(i).EQ.jblk )
THEN
554 usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
555 wend = wbegin + mb - 1
560 DO 30 i = wbegin, wend - 1
561 wgap( i ) = max( zero,
562 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
564 wgap( wend ) = max( zero,
565 $ vu - sigma - (w( wend )+werr( wend )))
567 indl = indexw(wbegin)
568 indu = indexw( wend )
571 IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd)
THEN
574 CALL slarrk( in, 1, gl, gu, d(ibegin),
575 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
576 IF( iinfo.NE.0 )
THEN
580 isleft = max(gl, tmp - tmp1
581 $ - hndrd * eps* abs(tmp - tmp1))
583 CALL slarrk( in, in, gl, gu, d(ibegin),
584 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
585 IF( iinfo.NE.0 )
THEN
589 isrght = min(gu, tmp + tmp1
590 $ + hndrd * eps * abs(tmp + tmp1))
592 spdiam = isrght - isleft
596 isleft = max(gl, w(wbegin) - werr(wbegin)
597 $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
598 isrght = min(gu,w(wend) + werr(wend)
599 $ + hndrd * eps * abs(w(wend)+ werr(wend)))
611 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
619 wend = wbegin + mb - 1
621 s1 = isleft + fourth * spdiam
622 s2 = isrght - fourth * spdiam
628 s1 = isleft + fourth * spdiam
629 s2 = isrght - fourth * spdiam
631 tmp = min(isrght,vu) - max(isleft,vl)
632 s1 = max(isleft,vl) + fourth * tmp
633 s2 = min(isrght,vu) - fourth * tmp
639 CALL slarrc(
'T', in, s1, s2, d(ibegin),
640 $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
646 ELSEIF( cnt1 - indl .GE. indu - cnt2 )
THEN
647 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
648 sigma = max(isleft,gl)
649 ELSEIF( usedqd )
THEN
656 sigma = max(isleft,vl)
660 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) )
THEN
661 sigma = min(isrght,gu)
662 ELSEIF( usedqd )
THEN
669 sigma = min(isrght,vu)
683 tau = spdiam*eps*n + two*pivmin
684 tau = max( tau,two*eps*abs(sigma) )
687 clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
688 avgap = abs(clwdth /
REAL(wend-wbegin))
689 IF( sgndef.EQ.one )
THEN
690 tau = half*max(wgap(wbegin),avgap)
691 tau = max(tau,werr(wbegin))
693 tau = half*max(wgap(wend-1),avgap)
694 tau = max(tau,werr(wend))
701 DO 80 idum = 1, maxtry
705 dpivot = d( ibegin ) - sigma
707 dmax = abs( work(1) )
710 work( 2*in+i ) = one / work( i )
711 tmp = e( j )*work( 2*in+i )
713 dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
715 dmax = max( dmax, abs(dpivot) )
719 IF( dmax .GT. maxgrowth*spdiam )
THEN
724 IF( usedqd .AND. .NOT.norep )
THEN
728 tmp = sgndef*work( i )
729 IF( tmp.LT.zero ) norep = .true.
736 IF( idum.EQ.maxtry-1 )
THEN
737 IF( sgndef.EQ.one )
THEN
740 $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
743 $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
746 sigma = sigma - sgndef * tau
765 CALL scopy( in, work, 1, d( ibegin ), 1 )
766 CALL scopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
779 CALL slarnv(2, iseed, 2*in-1, work(1))
781 d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
782 e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
784 d(iend) = d(iend)*(one+eps*four*work(in))
794 IF ( .NOT.usedqd )
THEN
802 werr(j) = werr(j) + abs(w(j)) * eps
806 DO 135 i = ibegin, iend-1
807 work( i ) = d( i ) * e( i )**2
810 CALL slarrb(in, d(ibegin), work(ibegin),
811 $ indl, indu, rtol1, rtol2, indl-1,
812 $ w(wbegin), wgap(wbegin), werr(wbegin),
813 $ work( 2*n+1 ), iwork, pivmin, spdiam,
815 IF( iinfo .NE. 0 )
THEN
821 wgap( wend ) = max( zero,
822 $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
823 DO 138 i = indl, indu
840 rtol = log(
REAL(in)) * four * eps
843 work( 2*i-1 ) = abs( d( j ) )
844 work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
847 work( 2*in-1 ) = abs( d( iend ) )
849 CALL slasq2( in, work, iinfo )
850 IF( iinfo .NE. 0 )
THEN
859 IF( work( i ).LT.zero )
THEN
865 IF( sgndef.GT.zero )
THEN
866 DO 150 i = indl, indu
868 w( m ) = work( in-i+1 )
873 DO 160 i = indl, indu
881 DO 165 i = m - mb + 1, m
883 werr( i ) = rtol * abs( w(i) )
885 DO 166 i = m - mb + 1, m - 1
887 wgap( i ) = max( zero,
888 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
890 wgap( m ) = max( zero,
891 $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
subroutine slarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
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 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 slarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
SLARRA computes the splitting points with the specified threshold.
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 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 scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine slasq2(N, Z, INFO)
SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...