dlarre()

subroutine dlarre	(	character	range,
		integer	n,
		double precision	vl,
		double precision	vu,
		integer	il,
		integer	iu,
		double precision, dimension( * )	d,
		double precision, dimension( * )	e,
		double precision, dimension( * )	e2,
		double precision	rtol1,
		double precision	rtol2,
		double precision	spltol,
		integer	nsplit,
		integer, dimension( * )	isplit,
		integer	m,
		double precision, dimension( * )	w,
		double precision, dimension( * )	werr,
		double precision, dimension( * )	wgap,
		integer, dimension( * )	iblock,
		integer, dimension( * )	indexw,
		double precision, dimension( * )	gers,
		double precision	pivmin,
		double precision, dimension( * )	work,
		integer, dimension( * )	iwork,
		integer	info )

DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.

Download DLARRE + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> To find the desired eigenvalues of a given real symmetric
!> tridiagonal matrix T, DLARRE sets any  off-diagonal
!> elements to zero, and for each unreduced block T_i, it finds
!> (a) a suitable shift at one end of the block's spectrum,
!> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
!> (c) eigenvalues of each L_i D_i L_i^T.
!> The representations and eigenvalues found are then used by
!> DSTEMR to compute the eigenvectors of T.
!> The accuracy varies depending on whether bisection is used to
!> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
!> compute all and then discard any unwanted one.
!> As an added benefit, DLARRE also outputs the n
!> Gerschgorin intervals for the matrices L_i D_i L_i^T.
!> 

Parameters

[in]	RANGE	!> RANGE is CHARACTER*1 !> = 'A': () all eigenvalues will be found. !> = 'V': () all eigenvalues in the half-open interval !> (VL, VU] will be found. !> = 'I': () the IL-th through IU-th eigenvalues (of the !> entire matrix) will be found. !>
[in]	N	!> N is INTEGER !> The order of the matrix. N > 0. !>
[in,out]	VL	!> VL is DOUBLE PRECISION !> If RANGE='V', the lower bound for the eigenvalues. !> Eigenvalues less than or equal to VL, or greater than VU, !> will not be returned. VL < VU. !> If RANGE='I' or ='A', DLARRE computes bounds on the desired !> part of the spectrum. !>
[in,out]	VU	!> VU is DOUBLE PRECISION !> If RANGE='V', the upper bound for the eigenvalues. !> Eigenvalues less than or equal to VL, or greater than VU, !> will not be returned. VL < VU. !> If RANGE='I' or ='A', DLARRE computes bounds on the desired !> part of the spectrum. !>
[in]	IL	!> IL is INTEGER !> If RANGE='I', the index of the !> smallest eigenvalue to be returned. !> 1 <= IL <= IU <= N. !>
[in]	IU	!> IU is INTEGER !> If RANGE='I', the index of the !> largest eigenvalue to be returned. !> 1 <= IL <= IU <= N. !>
[in,out]	D	!> D is DOUBLE PRECISION array, dimension (N) !> On entry, the N diagonal elements of the tridiagonal !> matrix T. !> On exit, the N diagonal elements of the diagonal !> matrices D_i. !>
[in,out]	E	!> E is DOUBLE PRECISION array, dimension (N) !> On entry, the first (N-1) entries contain the subdiagonal !> elements of the tridiagonal matrix T; E(N) need not be set. !> On exit, E contains the subdiagonal elements of the unit !> bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), !> 1 <= I <= NSPLIT, contain the base points sigma_i on output. !>
[in,out]	E2	!> E2 is DOUBLE PRECISION array, dimension (N) !> On entry, the first (N-1) entries contain the SQUARES of the !> subdiagonal elements of the tridiagonal matrix T; !> E2(N) need not be set. !> On exit, the entries E2( ISPLIT( I ) ), !> 1 <= I <= NSPLIT, have been set to zero !>
[in]	RTOL1	!> RTOL1 is DOUBLE PRECISION !>
[in]	RTOL2	!> RTOL2 is DOUBLE PRECISION !> Parameters for bisection. !> An interval [LEFT,RIGHT] has converged if !> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) !>
[in]	SPLTOL	!> SPLTOL is DOUBLE PRECISION !> The threshold for splitting. !>
[out]	NSPLIT	!> NSPLIT is INTEGER !> The number of blocks T splits into. 1 <= NSPLIT <= N. !>
[out]	ISPLIT	!> ISPLIT is INTEGER array, dimension (N) !> The splitting points, at which T breaks up into blocks. !> The first block consists of rows/columns 1 to ISPLIT(1), !> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), !> etc., and the NSPLIT-th consists of rows/columns !> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. !>
[out]	M	!> M is INTEGER !> The total number of eigenvalues (of all L_i D_i L_i^T) !> found. !>
[out]	W	!> W is DOUBLE PRECISION array, dimension (N) !> The first M elements contain the eigenvalues. The !> eigenvalues of each of the blocks, L_i D_i L_i^T, are !> sorted in ascending order ( DLARRE may use the !> remaining N-M elements as workspace). !>
[out]	WERR	!> WERR is DOUBLE PRECISION array, dimension (N) !> The error bound on the corresponding eigenvalue in W. !>
[out]	WGAP	!> WGAP is DOUBLE PRECISION array, dimension (N) !> The separation from the right neighbor eigenvalue in W. !> The gap is only with respect to the eigenvalues of the same block !> as each block has its own representation tree. !> Exception: at the right end of a block we store the left gap !>
[out]	IBLOCK	!> IBLOCK is INTEGER array, dimension (N) !> The indices of the blocks (submatrices) associated with the !> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue !> W(i) belongs to the first block from the top, =2 if W(i) !> belongs to the second block, etc. !>
[out]	INDEXW	!> INDEXW is INTEGER array, dimension (N) !> The indices of the eigenvalues within each block (submatrix); !> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the !> i-th eigenvalue W(i) is the 10-th eigenvalue in block 2 !>
[out]	GERS	!> GERS is DOUBLE PRECISION array, dimension (2*N) !> The N Gerschgorin intervals (the i-th Gerschgorin interval !> is (GERS(2*i-1), GERS(2*i)). !>
[out]	PIVMIN	!> PIVMIN is DOUBLE PRECISION !> The minimum pivot in the Sturm sequence for T. !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension (6*N) !> Workspace. !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (5*N) !> Workspace. !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> > 0: A problem occurred in DLARRE. !> < 0: One of the called subroutines signaled an internal problem. !> Needs inspection of the corresponding parameter IINFO !> for further information. !> !> =-1: Problem in DLARRD. !> = 2: No base representation could be found in MAXTRY iterations. !> Increasing MAXTRY and recompilation might be a remedy. !> =-3: Problem in DLARRB when computing the refined root !> representation for DLASQ2. !> =-4: Problem in DLARRB when preforming bisection on the !> desired part of the spectrum. !> =-5: Problem in DLASQ2. !> =-6: Problem in DLASQ2. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  The base representations are required to suffer very little
!>  element growth and consequently define all their eigenvalues to
!>  high relative accuracy.
!> 

Contributors:

Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 299 of file dlarre.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          RANGE
      INTEGER            IL, INFO, IU, M, N, NSPLIT
      DOUBLE PRECISION  PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
     $                   INDEXW( * )
      DOUBLE PRECISION   D( * ), E( * ), E2( * ), GERS( * ),
     $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
     $                   MAXGROWTH, ONE, PERT, TWO, ZERO
      parameter( zero = 0.0d0, one = 1.0d0,
     $                     two = 2.0d0, four=4.0d0,
     $                     hndrd = 100.0d0,
     $                     pert = 8.0d0,
     $                     half = one/two, fourth = one/four, fac= half,
     $                     maxgrowth = 64.0d0, fudge = 2.0d0 )
      INTEGER            MAXTRY, ALLRNG, INDRNG, VALRNG
      parameter( maxtry = 6, allrng = 1, indrng = 2,
     $                     valrng = 3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            FORCEB, NOREP, USEDQD
      INTEGER            CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
     $                   IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
     $                   WBEGIN, WEND
      DOUBLE PRECISION   AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
     $                   EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
     $                   RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
     $                   TAU, TMP, TMP1
 
 
*     ..
*     .. Local Arrays ..
      INTEGER            ISEED( 4 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION            DLAMCH
      EXTERNAL           dlamch, lsame
 
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlarnv, dlarra, dlarrb, dlarrc,
     $                   dlarrd,
     $                   dlasq2, dlarrk
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
 
*     ..
*     .. Executable Statements ..
*
 
      info = 0
      nsplit = 0
      m = 0
*
*     Quick return if possible
*
      IF( n.LE.0 ) THEN
         RETURN
      END IF
*
*     Decode RANGE
*
      IF( lsame( range, 'A' ) ) THEN
         irange = allrng
      ELSE IF( lsame( range, 'V' ) ) THEN
         irange = valrng
      ELSE IF( lsame( range, 'I' ) ) THEN
         irange = indrng
      END IF
 
*     Get machine constants
      safmin = dlamch( 'S' )
      eps = dlamch( 'P' )
 
*     Set parameters
      rtl = sqrt(eps)
      bsrtol = sqrt(eps)
 
*     Treat case of 1x1 matrix for quick return
      IF( n.EQ.1 ) THEN
         IF( (irange.EQ.allrng).OR.
     $       ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
     $       ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
            m = 1
            w(1) = d(1)
*           The computation error of the eigenvalue is zero
            werr(1) = zero
            wgap(1) = zero
            iblock( 1 ) = 1
            indexw( 1 ) = 1
            gers(1) = d( 1 )
            gers(2) = d( 1 )
         ENDIF
*        store the shift for the initial RRR, which is zero in this case
         e(1) = zero
         RETURN
      END IF
 
*     General case: tridiagonal matrix of order > 1
*
*     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
*     Compute maximum off-diagonal entry and pivmin.
      gl = d(1)
      gu = d(1)
      eold = zero
      emax = zero
      e(n) = zero
      DO 5 i = 1,n
         werr(i) = zero
         wgap(i) = zero
         eabs = abs( e(i) )
         IF( eabs .GE. emax ) THEN
            emax = eabs
         END IF
         tmp1 = eabs + eold
         gers( 2*i-1) = d(i) - tmp1
         gl =  min( gl, gers( 2*i - 1))
         gers( 2*i ) = d(i) + tmp1
         gu = max( gu, gers(2*i) )
         eold  = eabs
 5    CONTINUE
*     The minimum pivot allowed in the Sturm sequence for T
      pivmin = safmin * max( one, emax**2 )
*     Compute spectral diameter. The Gerschgorin bounds give an
*     estimate that is wrong by at most a factor of SQRT(2)
      spdiam = gu - gl
 
*     Compute splitting points
      CALL dlarra( n, d, e, e2, spltol, spdiam,
     $                    nsplit, isplit, iinfo )
 
*     Can force use of bisection instead of faster DQDS.
*     Option left in the code for future multisection work.
      forceb = .false.
 
*     Initialize USEDQD, DQDS should be used for ALLRNG unless someone
*     explicitly wants bisection.
      usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
 
      IF( (irange.EQ.allrng) .AND. (.NOT. forceb) ) THEN
*        Set interval [VL,VU] that contains all eigenvalues
         vl = gl
         vu = gu
      ELSE
*        We call DLARRD to find crude approximations to the eigenvalues
*        in the desired range. In case IRANGE = INDRNG, we also obtain the
*        interval (VL,VU] that contains all the wanted eigenvalues.
*        An interval [LEFT,RIGHT] has converged if
*        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
*        DLARRD needs a WORK of size 4*N, IWORK of size 3*N
         CALL dlarrd( range, 'B', n, vl, vu, il, iu, gers,
     $                    bsrtol, d, e, e2, pivmin, nsplit, isplit,
     $                    mm, w, werr, vl, vu, iblock, indexw,
     $                    work, iwork, iinfo )
         IF( iinfo.NE.0 ) THEN
            info = -1
            RETURN
         ENDIF
*        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
         DO 14 i = mm+1,n
            w( i ) = zero
            werr( i ) = zero
            iblock( i ) = 0
            indexw( i ) = 0
 14      CONTINUE
      END IF
 
 
***
*     Loop over unreduced blocks
      ibegin = 1
      wbegin = 1
      DO 170 jblk = 1, nsplit
         iend = isplit( jblk )
         in = iend - ibegin + 1
 
*        1 X 1 block
         IF( in.EQ.1 ) THEN
            IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
     $         ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
     $        .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
     $        ) THEN
               m = m + 1
               w( m ) = d( ibegin )
               werr(m) = zero
*              The gap for a single block doesn't matter for the later
*              algorithm and is assigned an arbitrary large value
               wgap(m) = zero
               iblock( m ) = jblk
               indexw( m ) = 1
               wbegin = wbegin + 1
            ENDIF
*           E( IEND ) holds the shift for the initial RRR
            e( iend ) = zero
            ibegin = iend + 1
            GO TO 170
         END IF
*
*        Blocks of size larger than 1x1
*
*        E( IEND ) will hold the shift for the initial RRR, for now set it =0
         e( iend ) = zero
*
*        Find local outer bounds GL,GU for the block
         gl = d(ibegin)
         gu = d(ibegin)
         DO 15 i = ibegin , iend
            gl = min( gers( 2*i-1 ), gl )
            gu = max( gers( 2*i ), gu )
 15      CONTINUE
         spdiam = gu - gl
 
         IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) ) THEN
*           Count the number of eigenvalues in the current block.
            mb = 0
            DO 20 i = wbegin,mm
               IF( iblock(i).EQ.jblk ) THEN
                  mb = mb+1
               ELSE
                  GOTO 21
               ENDIF
 20         CONTINUE
 21         CONTINUE
 
            IF( mb.EQ.0) THEN
*              No eigenvalue in the current block lies in the desired range
*              E( IEND ) holds the shift for the initial RRR
               e( iend ) = zero
               ibegin = iend + 1
               GO TO 170
            ELSE
 
*              Decide whether dqds or bisection is more efficient
               usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
               wend = wbegin + mb - 1
*              Calculate gaps for the current block
*              In later stages, when representations for individual
*              eigenvalues are different, we use SIGMA = E( IEND ).
               sigma = zero
               DO 30 i = wbegin, wend - 1
                  wgap( i ) = max( zero,
     $                        w(i+1)-werr(i+1) - (w(i)+werr(i)) )
 30            CONTINUE
               wgap( wend ) = max( zero,
     $                     vu - sigma - (w( wend )+werr( wend )))
*              Find local index of the first and last desired evalue.
               indl = indexw(wbegin)
               indu = indexw( wend )
            ENDIF
         ENDIF
         IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd) THEN
*           Case of DQDS
*           Find approximations to the extremal eigenvalues of the block
            CALL dlarrk( in, 1, gl, gu, d(ibegin),
     $               e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
            IF( iinfo.NE.0 ) THEN
               info = -1
               RETURN
            ENDIF
            isleft = max(gl, tmp - tmp1
     $               - hndrd * eps* abs(tmp - tmp1))
 
            CALL dlarrk( in, in, gl, gu, d(ibegin),
     $               e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
            IF( iinfo.NE.0 ) THEN
               info = -1
               RETURN
            ENDIF
            isrght = min(gu, tmp + tmp1
     $                 + hndrd * eps * abs(tmp + tmp1))
*           Improve the estimate of the spectral diameter
            spdiam = isrght - isleft
         ELSE
*           Case of bisection
*           Find approximations to the wanted extremal eigenvalues
            isleft = max(gl, w(wbegin) - werr(wbegin)
     $                  - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
            isrght = min(gu,w(wend) + werr(wend)
     $                  + hndrd * eps * abs(w(wend)+ werr(wend)))
         ENDIF
 
 
*        Decide whether the base representation for the current block
*        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
*        should be on the left or the right end of the current block.
*        The strategy is to shift to the end which is "more populated"
*        Furthermore, decide whether to use DQDS for the computation of
*        the eigenvalue approximations at the end of DLARRE or bisection.
*        dqds is chosen if all eigenvalues are desired or the number of
*        eigenvalues to be computed is large compared to the blocksize.
         IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
*           If all the eigenvalues have to be computed, we use dqd
            usedqd = .true.
*           INDL is the local index of the first eigenvalue to compute
            indl = 1
            indu = in
*           MB =  number of eigenvalues to compute
            mb = in
            wend = wbegin + mb - 1
*           Define 1/4 and 3/4 points of the spectrum
            s1 = isleft + fourth * spdiam
            s2 = isrght - fourth * spdiam
         ELSE
*           DLARRD has computed IBLOCK and INDEXW for each eigenvalue
*           approximation.
*           choose sigma
            IF( usedqd ) THEN
               s1 = isleft + fourth * spdiam
               s2 = isrght - fourth * spdiam
            ELSE
               tmp = min(isrght,vu) -  max(isleft,vl)
               s1 =  max(isleft,vl) + fourth * tmp
               s2 =  min(isrght,vu) - fourth * tmp
            ENDIF
         ENDIF
 
*        Compute the negcount at the 1/4 and 3/4 points
         IF(mb.GT.1) THEN
            CALL dlarrc( 'T', in, s1, s2, d(ibegin),
     $                    e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
         ENDIF
 
         IF(mb.EQ.1) THEN
            sigma = gl
            sgndef = one
         ELSEIF( cnt1 - indl .GE. indu - cnt2 ) THEN
            IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
               sigma = max(isleft,gl)
            ELSEIF( usedqd ) THEN
*              use Gerschgorin bound as shift to get pos def matrix
*              for dqds
               sigma = isleft
            ELSE
*              use approximation of the first desired eigenvalue of the
*              block as shift
               sigma = max(isleft,vl)
            ENDIF
            sgndef = one
         ELSE
            IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
               sigma = min(isrght,gu)
            ELSEIF( usedqd ) THEN
*              use Gerschgorin bound as shift to get neg def matrix
*              for dqds
               sigma = isrght
            ELSE
*              use approximation of the first desired eigenvalue of the
*              block as shift
               sigma = min(isrght,vu)
            ENDIF
            sgndef = -one
         ENDIF
 
 
*        An initial SIGMA has been chosen that will be used for computing
*        T - SIGMA I = L D L^T
*        Define the increment TAU of the shift in case the initial shift
*        needs to be refined to obtain a factorization with not too much
*        element growth.
         IF( usedqd ) THEN
*           The initial SIGMA was to the outer end of the spectrum
*           the matrix is definite and we need not retreat.
            tau = spdiam*eps*n + two*pivmin
            tau = max( tau,two*eps*abs(sigma) )
         ELSE
            IF(mb.GT.1) THEN
               clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
               avgap = abs(clwdth / dble(wend-wbegin))
               IF( sgndef.EQ.one ) THEN
                  tau = half*max(wgap(wbegin),avgap)
                  tau = max(tau,werr(wbegin))
               ELSE
                  tau = half*max(wgap(wend-1),avgap)
                  tau = max(tau,werr(wend))
               ENDIF
            ELSE
               tau = werr(wbegin)
            ENDIF
         ENDIF
*
         DO 80 idum = 1, maxtry
*           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
*           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
*           pivots in WORK(2*IN+1:3*IN)
            dpivot = d( ibegin ) - sigma
            work( 1 ) = dpivot
            dmax = abs( work(1) )
            j = ibegin
            DO 70 i = 1, in - 1
               work( 2*in+i ) = one / work( i )
               tmp = e( j )*work( 2*in+i )
               work( in+i ) = tmp
               dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
               work( i+1 ) = dpivot
               dmax = max( dmax, abs(dpivot) )
               j = j + 1
 70         CONTINUE
*           check for element growth
            IF( dmax .GT. maxgrowth*spdiam ) THEN
               norep = .true.
            ELSE
               norep = .false.
            ENDIF
            IF( usedqd .AND. .NOT.norep ) THEN
*              Ensure the definiteness of the representation
*              All entries of D (of L D L^T) must have the same sign
               DO 71 i = 1, in
                  tmp = sgndef*work( i )
                  IF( tmp.LT.zero ) norep = .true.
 71            CONTINUE
            ENDIF
            IF(norep) THEN
*              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
*              shift which makes the matrix definite. So we should end up
*              here really only in the case of IRANGE = VALRNG or INDRNG.
               IF( idum.EQ.maxtry-1 ) THEN
                  IF( sgndef.EQ.one ) THEN
*                    The fudged Gerschgorin shift should succeed
                     sigma =
     $                    gl - fudge*spdiam*eps*n - fudge*two*pivmin
                  ELSE
                     sigma =
     $                    gu + fudge*spdiam*eps*n + fudge*two*pivmin
                  END IF
               ELSE
                  sigma = sigma - sgndef * tau
                  tau = two * tau
               END IF
            ELSE
*              an initial RRR is found
               GO TO 83
            END IF
 80      CONTINUE
*        if the program reaches this point, no base representation could be
*        found in MAXTRY iterations.
         info = 2
         RETURN
 
 83      CONTINUE
*        At this point, we have found an initial base representation
*        T - SIGMA I = L D L^T with not too much element growth.
*        Store the shift.
         e( iend ) = sigma
*        Store D and L.
         CALL dcopy( in, work, 1, d( ibegin ), 1 )
         CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
 
 
         IF(mb.GT.1 ) THEN
*
*           Perturb each entry of the base representation by a small
*           (but random) relative amount to overcome difficulties with
*           glued matrices.
*
            DO 122 i = 1, 4
               iseed( i ) = 1
 122        CONTINUE
 
            CALL dlarnv(2, iseed, 2*in-1, work(1))
            DO 125 i = 1,in-1
               d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
               e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
 125        CONTINUE
            d(iend) = d(iend)*(one+eps*four*work(in))
*
         ENDIF
*
*        Don't update the Gerschgorin intervals because keeping track
*        of the updates would be too much work in DLARRV.
*        We update W instead and use it to locate the proper Gerschgorin
*        intervals.
 
*        Compute the required eigenvalues of L D L' by bisection or dqds
         IF ( .NOT.usedqd ) THEN
*           If DLARRD has been used, shift the eigenvalue approximations
*           according to their representation. This is necessary for
*           a uniform DLARRV since dqds computes eigenvalues of the
*           shifted representation. In DLARRV, W will always hold the
*           UNshifted eigenvalue approximation.
            DO 134 j=wbegin,wend
               w(j) = w(j) - sigma
               werr(j) = werr(j) + abs(w(j)) * eps
 134        CONTINUE
*           call DLARRB to reduce eigenvalue error of the approximations
*           from DLARRD
            DO 135 i = ibegin, iend-1
               work( i ) = d( i ) * e( i )**2
 135        CONTINUE
*           use bisection to find EV from INDL to INDU
            CALL dlarrb(in, d(ibegin), work(ibegin),
     $                  indl, indu, rtol1, rtol2, indl-1,
     $                  w(wbegin), wgap(wbegin), werr(wbegin),
     $                  work( 2*n+1 ), iwork, pivmin, spdiam,
     $                  in, iinfo )
            IF( iinfo .NE. 0 ) THEN
               info = -4
               RETURN
            END IF
*           DLARRB computes all gaps correctly except for the last one
*           Record distance to VU/GU
            wgap( wend ) = max( zero,
     $           ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
            DO 138 i = indl, indu
               m = m + 1
               iblock(m) = jblk
               indexw(m) = i
 138        CONTINUE
         ELSE
*           Call dqds to get all eigs (and then possibly delete unwanted
*           eigenvalues).
*           Note that dqds finds the eigenvalues of the L D L^T representation
*           of T to high relative accuracy. High relative accuracy
*           might be lost when the shift of the RRR is subtracted to obtain
*           the eigenvalues of T. However, T is not guaranteed to define its
*           eigenvalues to high relative accuracy anyway.
*           Set RTOL to the order of the tolerance used in DLASQ2
*           This is an ESTIMATED error, the worst case bound is 4*N*EPS
*           which is usually too large and requires unnecessary work to be
*           done by bisection when computing the eigenvectors
            rtol = log(dble(in)) * four * eps
            j = ibegin
            DO 140 i = 1, in - 1
               work( 2*i-1 ) = abs( d( j ) )
               work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
               j = j + 1
  140       CONTINUE
            work( 2*in-1 ) = abs( d( iend ) )
            work( 2*in ) = zero
            CALL dlasq2( in, work, iinfo )
            IF( iinfo .NE. 0 ) THEN
*              If IINFO = -5 then an index is part of a tight cluster
*              and should be changed. The index is in IWORK(1) and the
*              gap is in WORK(N+1)
               info = -5
               RETURN
            ELSE
*              Test that all eigenvalues are positive as expected
               DO 149 i = 1, in
                  IF( work( i ).LT.zero ) THEN
                     info = -6
                     RETURN
                  ENDIF
 149           CONTINUE
            END IF
            IF( sgndef.GT.zero ) THEN
               DO 150 i = indl, indu
                  m = m + 1
                  w( m ) = work( in-i+1 )
                  iblock( m ) = jblk
                  indexw( m ) = i
 150           CONTINUE
            ELSE
               DO 160 i = indl, indu
                  m = m + 1
                  w( m ) = -work( i )
                  iblock( m ) = jblk
                  indexw( m ) = i
 160           CONTINUE
            END IF
 
            DO 165 i = m - mb + 1, m
*              the value of RTOL below should be the tolerance in DLASQ2
               werr( i ) = rtol * abs( w(i) )
 165        CONTINUE
            DO 166 i = m - mb + 1, m - 1
*              compute the right gap between the intervals
               wgap( i ) = max( zero,
     $                          w(i+1)-werr(i+1) - (w(i)+werr(i)) )
 166        CONTINUE
            wgap( m ) = max( zero,
     $           ( vu-sigma ) - ( w( m ) + werr( m ) ) )
         END IF
*        proceed with next block
         ibegin = iend + 1
         wbegin = wend + 1
 170  CONTINUE
*
 
      RETURN
*
*     End of DLARRE
*

Here is the call graph for this function:

Here is the caller graph for this function: