DLAEBZ - contain the iteration loops which compute and use
       the function N(w), which is the count of eigenvalues of  a
       symmetric  tridiagonal  matrix T less than or equal to its
       argument w


SYNOPSIS

       SUBROUTINE DLAEBZ( IJOB, NITMAX,  N,  MMAX,  MINP,  NBMIN,
                          ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL,
                          AB, C, MOUT, NAB, WORK, IWORK, INFO )

           INTEGER        IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN,
                          NITMAX

           DOUBLE         PRECISION ABSTOL, PIVMIN, RELTOL

           INTEGER        IWORK( * ), NAB( MMAX, * ), NVAL( * )

           DOUBLE         PRECISION  AB(  MMAX, * ), C( * ), D( *
                          ), E( * ), E2( * ), WORK( * )


PURPOSE

       DLAEBZ contains the iteration loops which compute and  use
       the  function N(w), which is the count of eigenvalues of a
       symmetric tridiagonal matrix T less than or equal  to  its
       argument w. It performs a choice of two types of loops:

       IJOB=1, followed by
       IJOB=2:  It takes as input a list of intervals and returns
       a list of
               sufficiently small intervals whose union  contains
       the same
               eigenvalues  as  the  union of the original inter­
       vals.
               The   input   intervals   are   (AB(j,1),AB(j,2)],
       j=1,...,MINP.
               The output interval (AB(j,1),AB(j,2)] will contain
               eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <=  j
       <= MOUT.

       IJOB=3: It performs a binary search in each input interval
               (AB(j,1),AB(j,2)] for a point  w(j)  such that
               N(w(j))=NVAL(j), and uses  C(j)  as  the  starting
       point of
               the search.  If such a w(j) is found, then on out­
       put
               AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then
       on output
               (AB(j,1),AB(j,2)]  will  be  a small interval con­
       taining the
               point where N(w)  jumps  through  NVAL(j),  unless
       that point
               lies outside the initial interval.

       vals, i.e., of the form  (a,b] , which  includes   b   but
       not  a .

       To  avoid  underflow,  the matrix should be scaled so that
       its largest element is no greater than  overflow**(1/2)  *
       underflow**(1/4)  in  absolute  value.  To assure the most
       accurate computation  of  small  eigenvalues,  the  matrix
       should be scaled to be
       not much smaller than that, either.

       See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiag­
       onal Matrix", Report CS41, Computer Science  Dept.,  Stan­
       ford
       University, July 21, 1966

       Note:  the  arguments  are,  in general, *not* checked for
       unreasonable values.


ARGUMENTS

       IJOB    (input) INTEGER
               Specifies what is to be done:
               = 1:  Compute NAB for the initial intervals.
               = 2:  Perform bisection iteration to  find  eigen­
               values of T.
               =  3:  Perform bisection iteration to invert N(w),
               i.e., to find a point which has a specified number
               of  eigenvalues  of  T  to its left.  Other values
               will cause DLAEBZ to return with INFO=-1.

       NITMAX  (input) INTEGER
               The maximum number of "levels" of bisection to  be
               performed,  i.e.,  an interval of width W will not
               be made smaller than 2^(-NITMAX) * W.  If not  all
               intervals  have converged after NITMAX iterations,
               then INFO is set to the  number  of  non-converged
               intervals.

       N       (input) INTEGER
               The  dimension  n of the tridiagonal matrix T.  It
               must be at least 1.

       MMAX    (input) INTEGER
               The maximum number of  intervals.   If  more  than
               MMAX  intervals  are  generated,  then DLAEBZ will
               quit with INFO=MMAX+1.

       MINP    (input) INTEGER
               The initial number of intervals.  It  may  not  be
               greater than MMAX.

       NBMIN   (input) INTEGER
               The  smallest  number  of intervals that should be

               the scalar loop will be used.

       ABSTOL  (input) DOUBLE PRECISION
               The minimum (absolute) width of an interval.  When
               an interval is narrower than ABSTOL, or than  REL­
               TOL 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  REL­
               TOL times the larger (in magnitude) endpoint, then
               it is considered to be sufficiently  small,  i.e.,
               converged.   Note:  this should always be at least
               radix*machine epsilon.

       PIVMIN  (input) DOUBLE PRECISION
               The minimum absolute value of  a  "pivot"  in  the
               Sturm sequence loop.  This *must* be at least  max
               |e(j)**2| * safe_min  and at least safe_min, where
               safe_min  is at least the smallest number that can
               divide one without overflow.

       D       (input) DOUBLE PRECISION array, dimension (N)
               The diagonal elements of the tridiagonal matrix T.

       E       (input) DOUBLE PRECISION array, dimension (N)
               The offdiagonal elements of the tridiagonal matrix
               T in positions 1 through N-1.  E(N) is  arbitrary.

       E2      (input) DOUBLE PRECISION array, dimension (N)
               The  squares  of  the  offdiagonal elements of the
               tridiagonal matrix T.  E2(N) is ignored.

       NVAL    (input/output) INTEGER array, dimension (MINP)
               If IJOB=1 or 2, not referenced.   If  IJOB=3,  the
               desired values of N(w).  The elements of NVAL will
               be reordered to correspond with the  intervals  in
               AB.   Thus, NVAL(j) on output will not, in general
               be the same as NVAL(j) on input, but it will  cor­
               respond  with  the  interval  (AB(j,1),AB(j,2)] on
               output.

       AB      (input/output) DOUBLE PRECISION array, dimension
               (MMAX,2)
               The endpoints of the intervals.  AB(j,1) is  a(j),
               the  left  endpoint  of  the  j-th  interval,  and
               AB(j,2) is b(j), the right endpoint  of  the  j-th
               interval.   The  input intervals will, in general,
               be modified, split, and reordered by the  calcula­
               tion.

               (MMAX)
               If IJOB=1, ignored.   If  IJOB=2,  workspace.   If
               IJOB=3,  then  on input C(j) should be initialized
               to the first search point in the binary search.

       MOUT    (output) INTEGER
               If IJOB=1, the number of eigenvalues in the inter­
               vals.   If  IJOB=2  or  3, the number of intervals
               output.  If IJOB=3, MOUT will equal MINP.

       NAB     (input/output) INTEGER array, dimension (MMAX,2)
               If IJOB=1, then on output NAB(i,j) will be set  to
               N(AB(i,j)).   If  IJOB=2,  then on input, NAB(i,j)
               should be set.  It  must  satisfy  the  condition:
               N(AB(i,1))  <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
               which means that in interval  i  only  eigenvalues
               NAB(i,1)+1,...,NAB(i,2)  will be considered.  Usu­
               ally, NAB(i,j)=N(AB(i,j)), from a previous call to
               DLAEBZ with IJOB=1.  On output, NAB(i,j) will con­
               tain max(na(k),min(nb(k),N(AB(i,j)))), where k  is
               the  index  of  the input interval that the output
               interval (AB(j,1),AB(j,2)] came  from,  and  na(k)
               and nb(k) are the the input values of NAB(k,1) and
               NAB(k,2).  If IJOB=3,  then  on  output,  NAB(i,j)
               contains N(AB(i,j)), unless N(w) > NVAL(i) for all
               search points  w , in which case NAB(i,1) will not
               be  modified,  i.e.,  the output value will be the
               same as the input value (modulo reorderings -- see
               NVAL  and  AB),  or  unless N(w) < NVAL(i) for all
               search points  w , in which case NAB(i,2) will not
               be  modified.  Normally, NAB should be set to some
               distinctive value(s) before DLAEBZ is called.

       WORK    (workspace) DOUBLE PRECISION array, dimension
               (MMAX)
               Workspace.

       IWORK   (workspace) INTEGER array, dimension (MMAX)
               Workspace.

       INFO    (output) INTEGER
               = 0:       All intervals converged.
               =  1--MMAX:  The  last INFO intervals did not con­
               verge.
               = MMAX+1:  More than MMAX  intervals  were  gener­
               ated.


FURTHER DETAILS

           This  routine  is  intended to be called only by other
       LAPACK routines, thus the interface is less user-friendly.
       It is intended for two purposes:

       (a) finding eigenvalues.  In this case, DLAEBZ should have

           more initial intervals set up in AB, and DLAEBZ should
       be called
           with  IJOB=1.   This  sets up NAB, and also counts the
       eigenvalues.
           Intervals with no eigenvalues would usually be  thrown
       out at
           this  point.   Also,  if not all the eigenvalues in an
       interval i
           are desired, NAB(i,1) can  be  increased  or  NAB(i,2)
       decreased.
           For   example,  set  NAB(i,1)=NAB(i,2)-1  to  get  the
       largest
           eigenvalue.  DLAEBZ is then  called  with  IJOB=2  and
       MMAX
           no smaller than the value of MOUT returned by the call
       with
           IJOB=1.   After  this   (IJOB=2)   call,   eigenvalues
       NAB(i,1)+1
           through   NAB(i,2)   are   approximately  AB(i,1)  (or
       AB(i,2)) to the
           tolerance specified by ABSTOL and RELTOL.

       (b) finding an  interval  (a',b']  containing  eigenvalues
       w(f),...,w(l).
           In this case, start with a Gershgorin interval  (a,b).
       Set up
           AB to  contain  2  search  intervals,  both  initially
       (a,b).  One
           NVAL element should contain  f-1  and the other should
       contain  l
           , while C should contain  a  and  b,  resp.   NAB(i,1)
       should be -1
           and  NAB(i,2)  should  be N+1, to flag an error if the
       desired
           interval does not lie in (a,b).  DLAEBZ is then called
       with
           IJOB=3.   On  exit,  if w(f-1) < w(f), then one of the
       intervals --
           j     --     will     have     AB(j,1)=AB(j,2)     and
       NAB(j,1)=NAB(j,2)=f-1, while
           if, to the specified tolerance, w(f-k)=...=w(f+r), k >
       0 and r
           >=    0,    then    the     interval     will     have
       N(AB(j,1))=NAB(j,1)=f-k and
           N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
           w(l-r)=...=w(l+k) are handled similarly.


Man(1) output converted with man2html