pdlaebz()

subroutine pdlaebz	(	integer	ijob,
		integer	n,
		integer	mmax,
		integer	minp,
		double precision	abstol,
		double precision	reltol,
		double precision	pivmin,
		double precision, dimension( * )	d,
		integer, dimension( * )	nval,
		double precision, dimension( * )	intvl,
		integer, dimension( * )	intvlct,
		integer	mout,
		double precision	lsave,
		integer	ieflag,
		integer	info
	)

Definition at line 883 of file pdstebz.f.

*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*
*     .. Scalar Arguments ..
      INTEGER            IEFLAG, IJOB, INFO, MINP, MMAX, MOUT, N
      DOUBLE PRECISION   ABSTOL, LSAVE, PIVMIN, RELTOL
*     ..
*     .. Array Arguments ..
      INTEGER            INTVLCT( * ), NVAL( * )
      DOUBLE PRECISION   D( * ), INTVL( * )
*     ..
*
*  Purpose
*  =======
*
*  PDLAEBZ contains the iteration loop which computes the eigenvalues
*  contained in the input intervals [ INTVL(2*j-1), INTVL(2*j) ] where
*  j = 1,...,MINP. It uses and computes the function N(w), which is
*  the count of eigenvalues of a symmetric tridiagonal matrix less than
*  or equal to its argument w.
*
*  This is a ScaLAPACK internal subroutine and arguments are not
*  checked for unreasonable values.
*
*  Arguments
*  =========
*
*  IJOB    (input) INTEGER
*          Specifies the computation done by PDLAEBZ
*          = 0 : Find an interval with desired values of N(w) at the
*                endpoints of the interval.
*          = 1 : Find a floating point number contained in the initial
*                interval with a desired value of N(w).
*          = 2 : Perform bisection iteration to find eigenvalues of T.
*
*  N       (input) INTEGER
*          The order of the tridiagonal matrix T. N >= 1.
*
*  MMAX    (input) INTEGER
*          The maximum number of intervals that may be generated. If
*          more than MMAX intervals are generated, then PDLAEBZ will
*          quit with INFO = MMAX+1.
*
*  MINP    (input) INTEGER
*          The initial number of intervals. MINP <= MMAX.
*
*  ABSTOL  (input) DOUBLE PRECISION
*          The minimum (absolute) width of an interval. When an interval
*          is narrower than ABSTOL, or than RELTOL times the larger (in
*          magnitude) endpoint, then it is considered to be sufficiently
*          small, i.e., converged.
*          This must be at least zero.
*
*  RELTOL  (input) DOUBLE PRECISION
*          The minimum relative width of an interval. When an interval
*          is narrower than ABSTOL, or than RELTOL times the larger (in
*          magnitude) endpoint, then it is considered to be sufficiently
*          small, i.e., converged.
*          Note : This should be at least radix*machine epsilon.
*
*  PIVMIN  (input) DOUBLE PRECISION
*          The minimum absolute of a "pivot" in the "paranoid"
*          implementation of the Sturm sequence loop. This must be at
*          least max_j |e(j)^2| *safe_min, and at least safe_min, where
*          safe_min is at least the smallest number that can divide 1.0
*          without overflow.
*          See PDLAPDCT for the "paranoid" implementation of the Sturm
*          sequence loop.
*
*  D       (input) DOUBLE PRECISION array, dimension (2*N - 1)
*          Contains the diagonals and the squares of the off-diagonal
*          elements of the tridiagonal matrix T. These elements are
*          assumed to be interleaved in memory for better cache
*          performance. The diagonal entries of T are in the entries
*          D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
*          entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
*          matrix must be scaled so that its largest entry is no greater
*          than overflow**(1/2) * underflow**(1/4) in absolute value,
*          and for greatest accuracy, it should not be much smaller
*          than that.
*
*  NVAL    (input/output) INTEGER array, dimension (4)
*          If IJOB = 0, the desired values of N(w) are in NVAL(1) and
*                       NVAL(2).
*          If IJOB = 1, NVAL(2) is the desired value of N(w).
*          If IJOB = 2, not referenced.
*          This array will, in general, be reordered on output.
*
*  INTVL   (input/output) DOUBLE PRECISION array, dimension (2*MMAX)
*          The endpoints of the intervals. INTVL(2*j-1) is the left
*          endpoint of the j-th interval, and INTVL(2*j) is the right
*          endpoint of the j-th interval. The input intervals will,
*          in general, be modified, split and reordered by the
*          calculation.
*          On input, INTVL contains the MINP input intervals.
*          On output, INTVL contains the converged intervals.
*
*  INTVLCT (input/output) INTEGER array, dimension (2*MMAX)
*          The counts at the endpoints of the intervals. INTVLCT(2*j-1)
*          is the count at the left endpoint of the j-th interval, i.e.,
*          the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
*          count at the right endpoint of the j-th interval.
*          On input, INTVLCT contains the counts at the endpoints of
*          the MINP input intervals.
*          On output, INTVLCT contains the counts at the endpoints of
*          the converged intervals.
*
*  MOUT    (output) INTEGER
*          The number of intervals output.
*
*  LSAVE   (output) DOUBLE PRECISION
*          If IJOB = 0 or 2, not referenced.
*          If IJOB = 1, this is the largest floating point number
*          encountered which has count N(w) = NVAL(1).
*
*  IEFLAG  (input) INTEGER
*          A flag which indicates whether N(w) should be speeded up by
*          exploiting IEEE Arithmetic.
*
*  INFO    (output) INTEGER
*          = 0        : All intervals converged.
*          = 1 - MMAX : The last INFO intervals did not converge.
*          = MMAX + 1 : More than MMAX intervals were generated.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          abs, int, log, max, min
*     ..
*     .. External Subroutines ..
      EXTERNAL           pdlaecv, pdlaiectb, pdlaiectl, pdlapdct
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, TWO, HALF
      parameter( zero = 0.0d+0, two = 2.0d+0,
     $                   half = 1.0d+0 / two )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ITMAX, J, K, KF, KL, KLNEW, L, LCNT, LREQ,
     $                   NALPHA, NBETA, NMID, RCNT, RREQ
      DOUBLE PRECISION   ALPHA, BETA, MID
*     ..
*     .. Executable Statements ..
*
      kf = 1
      kl = minp + 1
      info = 0
      IF( intvl( 2 )-intvl( 1 ).LE.zero ) THEN
         info = minp
         mout = kf
         RETURN
      END IF
      IF( ijob.EQ.0 ) THEN
*
*        Check if some input intervals have "converged"
*
         CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
     $                 max( abstol, pivmin ), reltol )
         IF( kf.GE.kl )
     $      GO TO 60
*
*        Compute upper bound on number of iterations needed
*
         itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
     $           log( pivmin ) ) / log( two ) ) + 2
*
*        Iteration Loop
*
         DO 20 i = 1, itmax
            klnew = kl
            DO 10 j = kf, kl - 1
               k = 2*j
*
*           Bisect the interval and find the count at that point
*
               mid = half*( intvl( k-1 )+intvl( k ) )
               IF( ieflag.EQ.0 ) THEN
                  CALL pdlapdct( mid, n, d, pivmin, nmid )
               ELSE IF( ieflag.EQ.1 ) THEN
                  CALL pdlaiectb( mid, n, d, nmid )
               ELSE
                  CALL pdlaiectl( mid, n, d, nmid )
               END IF
               lreq = nval( k-1 )
               rreq = nval( k )
               IF( kl.EQ.1 )
     $            nmid = min( intvlct( k ),
     $                   max( intvlct( k-1 ), nmid ) )
               IF( nmid.LE.nval( k-1 ) ) THEN
                  intvl( k-1 ) = mid
                  intvlct( k-1 ) = nmid
               END IF
               IF( nmid.GE.nval( k ) ) THEN
                  intvl( k ) = mid
                  intvlct( k ) = nmid
               END IF
               IF( nmid.GT.lreq .AND. nmid.LT.rreq ) THEN
                  l = 2*klnew
                  intvl( l-1 ) = mid
                  intvl( l ) = intvl( k )
                  intvlct( l-1 ) = nval( k )
                  intvlct( l ) = intvlct( k )
                  intvl( k ) = mid
                  intvlct( k ) = nval( k-1 )
                  nval( l-1 ) = nval( k )
                  nval( l ) = nval( l-1 )
                  nval( k ) = nval( k-1 )
                  klnew = klnew + 1
               END IF
   10       CONTINUE
            kl = klnew
            CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
     $                    max( abstol, pivmin ), reltol )
            IF( kf.GE.kl )
     $         GO TO 60
   20    CONTINUE
      ELSE IF( ijob.EQ.1 ) THEN
         alpha = intvl( 1 )
         beta = intvl( 2 )
         nalpha = intvlct( 1 )
         nbeta = intvlct( 2 )
         lsave = alpha
         lreq = nval( 1 )
         rreq = nval( 2 )
   30    CONTINUE
         IF( nbeta.NE.rreq .AND. beta-alpha.GT.
     $       max( abstol, reltol*max( abs( alpha ), abs( beta ) ) ) )
     $        THEN
*
*           Bisect the interval and find the count at that point
*
            mid = half*( alpha+beta )
            IF( ieflag.EQ.0 ) THEN
               CALL pdlapdct( mid, n, d, pivmin, nmid )
            ELSE IF( ieflag.EQ.1 ) THEN
               CALL pdlaiectb( mid, n, d, nmid )
            ELSE
               CALL pdlaiectl( mid, n, d, nmid )
            END IF
            nmid = min( nbeta, max( nalpha, nmid ) )
            IF( nmid.GE.rreq ) THEN
               beta = mid
               nbeta = nmid
            ELSE
               alpha = mid
               nalpha = nmid
               IF( nmid.EQ.lreq )
     $            lsave = alpha
            END IF
            GO TO 30
         END IF
         kl = kf
         intvl( 1 ) = alpha
         intvl( 2 ) = beta
         intvlct( 1 ) = nalpha
         intvlct( 2 ) = nbeta
      ELSE IF( ijob.EQ.2 ) THEN
*
*        Check if some input intervals have "converged"
*
         CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
     $                 max( abstol, pivmin ), reltol )
         IF( kf.GE.kl )
     $      GO TO 60
*
*        Compute upper bound on number of iterations needed
*
         itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
     $           log( pivmin ) ) / log( two ) ) + 2
*
*        Iteration Loop
*
         DO 50 i = 1, itmax
            klnew = kl
            DO 40 j = kf, kl - 1
               k = 2*j
               mid = half*( intvl( k-1 )+intvl( k ) )
               IF( ieflag.EQ.0 ) THEN
                  CALL pdlapdct( mid, n, d, pivmin, nmid )
               ELSE IF( ieflag.EQ.1 ) THEN
                  CALL pdlaiectb( mid, n, d, nmid )
               ELSE
                  CALL pdlaiectl( mid, n, d, nmid )
               END IF
               lcnt = intvlct( k-1 )
               rcnt = intvlct( k )
               nmid = min( rcnt, max( lcnt, nmid ) )
*
*              Form New Interval(s)
*
               IF( nmid.EQ.lcnt ) THEN
                  intvl( k-1 ) = mid
               ELSE IF( nmid.EQ.rcnt ) THEN
                  intvl( k ) = mid
               ELSE IF( klnew.LT.mmax+1 ) THEN
                  l = 2*klnew
                  intvl( l-1 ) = mid
                  intvl( l ) = intvl( k )
                  intvlct( l-1 ) = nmid
                  intvlct( l ) = intvlct( k )
                  intvl( k ) = mid
                  intvlct( k ) = nmid
                  klnew = klnew + 1
               ELSE
                  info = mmax + 1
                  RETURN
               END IF
   40       CONTINUE
            kl = klnew
            CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
     $                    max( abstol, pivmin ), reltol )
            IF( kf.GE.kl )
     $         GO TO 60
   50    CONTINUE
      END IF
   60 CONTINUE
      info = max( kl-kf, 0 )
      mout = kl - 1
      RETURN
*
*     End of PDLAEBZ
*

Here is the call graph for this function:

Here is the caller graph for this function: