Loading...
Searching...
No Matches

## slaebz()

subroutine slaebz | ( | integer | ijob,
|  | integer | nitmax,
|  | integer | n,
|  | integer | mmax,
|  | integer | minp,
|  | integer | nbmin,
|  | real | abstol,
|  | real | reltol,
|  | real | pivmin,
|  | real, dimension( * ) | d,
|  | real, dimension( * ) | e,
|  | real, dimension( * ) | e2,
|  | integer, dimension( * ) | nval,
|  | real, dimension( mmax, * ) | ab,
|  | real, dimension( * ) | c,
|  | integer | mout,
|  | integer, dimension( mmax, * ) | nab,
|  | real, dimension( * ) | work,
|  | integer, dimension( * ) | iwork,
|  | integer | info )

SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz.

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

Purpose:

      !>
      !> SLAEBZ 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 intervals.
      !>         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 output
      !>         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 containing the
      !>         point where N(w) jumps through NVAL(j), unless that point
      !>         lies outside the initial interval.
      !>
      !> Note that the intervals are in all cases half-open intervals,
      !> 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 , Report CS41, Computer Science Dept., Stanford
      !> University, July 21, 1966
      !>
      !> Note: the arguments are, in general, *not* checked for unreasonable
      !> values.
      !>

Parameters

  [in]

  IJOB

      !>          IJOB is INTEGER
      !>          Specifies what is to be done:
      !>          = 1:  Compute NAB for the initial intervals.
      !>          = 2:  Perform bisection iteration to find eigenvalues 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 SLAEBZ to return with INFO=-1.
      !>

  [in]

  NITMAX

      !>          NITMAX is INTEGER
      !>          The maximum number of  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.
      !>

  [in]

  N

      !>          N is INTEGER
      !>          The dimension n of the tridiagonal matrix T.  It must be at
      !>          least 1.
      !>

  [in]

  MMAX

      !>          MMAX is INTEGER
      !>          The maximum number of intervals.  If more than MMAX intervals
      !>          are generated, then SLAEBZ will quit with INFO=MMAX+1.
      !>

  [in]

  MINP

      !>          MINP is INTEGER
      !>          The initial number of intervals.  It may not be greater than
      !>          MMAX.
      !>

  [in]

  NBMIN

      !>          NBMIN is INTEGER
      !>          The smallest number of intervals that should be processed
      !>          using a vector loop.  If zero, then only the scalar loop
      !>          will be used.
      !>

  [in]

  ABSTOL

      !>          ABSTOL is REAL
      !>          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.
      !>

  [in]

  RELTOL

      !>          RELTOL is REAL
      !>          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
      !>          always be at least radix*machine epsilon.
      !>

  [in]

  PIVMIN

      !>          PIVMIN is REAL
      !>          The minimum absolute value of a  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.
      !>

  [in]

  D

      !>          D is REAL array, dimension (N)
      !>          The diagonal elements of the tridiagonal matrix T.
      !>

  [in]

  E

      !>          E is REAL array, dimension (N)
      !>          The offdiagonal elements of the tridiagonal matrix T in
      !>          positions 1 through N-1.  E(N) is arbitrary.
      !>

  [in]

  E2

      !>          E2 is REAL array, dimension (N)
      !>          The squares of the offdiagonal elements of the tridiagonal
      !>          matrix T.  E2(N) is ignored.
      !>

  [in,out]

  NVAL

      !>          NVAL is 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 correspond with the interval
      !>          (AB(j,1),AB(j,2)] on output.
      !>

  [in,out]

  AB

      !>          AB is REAL 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
      !>          calculation.
      !>

  [in,out]

  C

      !>          C is REAL array, dimension (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.
      !>

  [out]

  MOUT

      !>          MOUT is INTEGER
      !>          If IJOB=1, the number of eigenvalues in the intervals.
      !>          If IJOB=2 or 3, the number of intervals output.
      !>          If IJOB=3, MOUT will equal MINP.
      !>

  [in,out]

  NAB

      !>          NAB is 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.  Usually,
      !>             NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with
      !>             IJOB=1.
      !>             On output, NAB(i,j) will contain
      !>             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 SLAEBZ is called.
      !>

  [out]

  WORK

      !>          WORK is REAL array, dimension (MMAX)
      !>          Workspace.
      !>

  [out]

  IWORK

      !>          IWORK is INTEGER array, dimension (MMAX)
      !>          Workspace.
      !>

  [out]

  INFO

      !>          INFO is INTEGER
      !>          = 0:       All intervals converged.
      !>          = 1--MMAX: The last INFO intervals did not converge.
      !>          = MMAX+1:  More than MMAX intervals were generated.
      !>

Author
  Univ. of Tennessee
  Univ. of California Berkeley
  Univ. of Colorado Denver
  NAG Ltd.

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, SLAEBZ should have one or
      !>      more initial intervals set up in AB, and SLAEBZ 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.  SLAEBZ 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).  SLAEBZ 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.
      !>