LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slaebz.f
Go to the documentation of this file.
00001 *> \brief \b SLAEBZ
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLAEBZ + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaebz.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaebz.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaebz.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
00022 *                          RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
00023 *                          NAB, WORK, IWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
00027 *       REAL               ABSTOL, PIVMIN, RELTOL
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
00031 *       REAL               AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
00032 *      $                   WORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> SLAEBZ contains the iteration loops which compute and use the
00042 *> function N(w), which is the count of eigenvalues of a symmetric
00043 *> tridiagonal matrix T less than or equal to its argument  w.  It
00044 *> performs a choice of two types of loops:
00045 *>
00046 *> IJOB=1, followed by
00047 *> IJOB=2: It takes as input a list of intervals and returns a list of
00048 *>         sufficiently small intervals whose union contains the same
00049 *>         eigenvalues as the union of the original intervals.
00050 *>         The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
00051 *>         The output interval (AB(j,1),AB(j,2)] will contain
00052 *>         eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
00053 *>
00054 *> IJOB=3: It performs a binary search in each input interval
00055 *>         (AB(j,1),AB(j,2)] for a point  w(j)  such that
00056 *>         N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
00057 *>         the search.  If such a w(j) is found, then on output
00058 *>         AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
00059 *>         (AB(j,1),AB(j,2)] will be a small interval containing the
00060 *>         point where N(w) jumps through NVAL(j), unless that point
00061 *>         lies outside the initial interval.
00062 *>
00063 *> Note that the intervals are in all cases half-open intervals,
00064 *> i.e., of the form  (a,b] , which includes  b  but not  a .
00065 *>
00066 *> To avoid underflow, the matrix should be scaled so that its largest
00067 *> element is no greater than  overflow**(1/2) * underflow**(1/4)
00068 *> in absolute value.  To assure the most accurate computation
00069 *> of small eigenvalues, the matrix should be scaled to be
00070 *> not much smaller than that, either.
00071 *>
00072 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00073 *> Matrix", Report CS41, Computer Science Dept., Stanford
00074 *> University, July 21, 1966
00075 *>
00076 *> Note: the arguments are, in general, *not* checked for unreasonable
00077 *> values.
00078 *> \endverbatim
00079 *
00080 *  Arguments:
00081 *  ==========
00082 *
00083 *> \param[in] IJOB
00084 *> \verbatim
00085 *>          IJOB is INTEGER
00086 *>          Specifies what is to be done:
00087 *>          = 1:  Compute NAB for the initial intervals.
00088 *>          = 2:  Perform bisection iteration to find eigenvalues of T.
00089 *>          = 3:  Perform bisection iteration to invert N(w), i.e.,
00090 *>                to find a point which has a specified number of
00091 *>                eigenvalues of T to its left.
00092 *>          Other values will cause SLAEBZ to return with INFO=-1.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] NITMAX
00096 *> \verbatim
00097 *>          NITMAX is INTEGER
00098 *>          The maximum number of "levels" of bisection to be
00099 *>          performed, i.e., an interval of width W will not be made
00100 *>          smaller than 2^(-NITMAX) * W.  If not all intervals
00101 *>          have converged after NITMAX iterations, then INFO is set
00102 *>          to the number of non-converged intervals.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] N
00106 *> \verbatim
00107 *>          N is INTEGER
00108 *>          The dimension n of the tridiagonal matrix T.  It must be at
00109 *>          least 1.
00110 *> \endverbatim
00111 *>
00112 *> \param[in] MMAX
00113 *> \verbatim
00114 *>          MMAX is INTEGER
00115 *>          The maximum number of intervals.  If more than MMAX intervals
00116 *>          are generated, then SLAEBZ will quit with INFO=MMAX+1.
00117 *> \endverbatim
00118 *>
00119 *> \param[in] MINP
00120 *> \verbatim
00121 *>          MINP is INTEGER
00122 *>          The initial number of intervals.  It may not be greater than
00123 *>          MMAX.
00124 *> \endverbatim
00125 *>
00126 *> \param[in] NBMIN
00127 *> \verbatim
00128 *>          NBMIN is INTEGER
00129 *>          The smallest number of intervals that should be processed
00130 *>          using a vector loop.  If zero, then only the scalar loop
00131 *>          will be used.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] ABSTOL
00135 *> \verbatim
00136 *>          ABSTOL is REAL
00137 *>          The minimum (absolute) width of an interval.  When an
00138 *>          interval is narrower than ABSTOL, or than RELTOL times the
00139 *>          larger (in magnitude) endpoint, then it is considered to be
00140 *>          sufficiently small, i.e., converged.  This must be at least
00141 *>          zero.
00142 *> \endverbatim
00143 *>
00144 *> \param[in] RELTOL
00145 *> \verbatim
00146 *>          RELTOL is REAL
00147 *>          The minimum relative width of an interval.  When an interval
00148 *>          is narrower than ABSTOL, or than RELTOL times the larger (in
00149 *>          magnitude) endpoint, then it is considered to be
00150 *>          sufficiently small, i.e., converged.  Note: this should
00151 *>          always be at least radix*machine epsilon.
00152 *> \endverbatim
00153 *>
00154 *> \param[in] PIVMIN
00155 *> \verbatim
00156 *>          PIVMIN is REAL
00157 *>          The minimum absolute value of a "pivot" in the Sturm
00158 *>          sequence loop.
00159 *>          This must be at least  max |e(j)**2|*safe_min  and at
00160 *>          least safe_min, where safe_min is at least
00161 *>          the smallest number that can divide one without overflow.
00162 *> \endverbatim
00163 *>
00164 *> \param[in] D
00165 *> \verbatim
00166 *>          D is REAL array, dimension (N)
00167 *>          The diagonal elements of the tridiagonal matrix T.
00168 *> \endverbatim
00169 *>
00170 *> \param[in] E
00171 *> \verbatim
00172 *>          E is REAL array, dimension (N)
00173 *>          The offdiagonal elements of the tridiagonal matrix T in
00174 *>          positions 1 through N-1.  E(N) is arbitrary.
00175 *> \endverbatim
00176 *>
00177 *> \param[in] E2
00178 *> \verbatim
00179 *>          E2 is REAL array, dimension (N)
00180 *>          The squares of the offdiagonal elements of the tridiagonal
00181 *>          matrix T.  E2(N) is ignored.
00182 *> \endverbatim
00183 *>
00184 *> \param[in,out] NVAL
00185 *> \verbatim
00186 *>          NVAL is INTEGER array, dimension (MINP)
00187 *>          If IJOB=1 or 2, not referenced.
00188 *>          If IJOB=3, the desired values of N(w).  The elements of NVAL
00189 *>          will be reordered to correspond with the intervals in AB.
00190 *>          Thus, NVAL(j) on output will not, in general be the same as
00191 *>          NVAL(j) on input, but it will correspond with the interval
00192 *>          (AB(j,1),AB(j,2)] on output.
00193 *> \endverbatim
00194 *>
00195 *> \param[in,out] AB
00196 *> \verbatim
00197 *>          AB is REAL array, dimension (MMAX,2)
00198 *>          The endpoints of the intervals.  AB(j,1) is  a(j), the left
00199 *>          endpoint of the j-th interval, and AB(j,2) is b(j), the
00200 *>          right endpoint of the j-th interval.  The input intervals
00201 *>          will, in general, be modified, split, and reordered by the
00202 *>          calculation.
00203 *> \endverbatim
00204 *>
00205 *> \param[in,out] C
00206 *> \verbatim
00207 *>          C is REAL array, dimension (MMAX)
00208 *>          If IJOB=1, ignored.
00209 *>          If IJOB=2, workspace.
00210 *>          If IJOB=3, then on input C(j) should be initialized to the
00211 *>          first search point in the binary search.
00212 *> \endverbatim
00213 *>
00214 *> \param[out] MOUT
00215 *> \verbatim
00216 *>          MOUT is INTEGER
00217 *>          If IJOB=1, the number of eigenvalues in the intervals.
00218 *>          If IJOB=2 or 3, the number of intervals output.
00219 *>          If IJOB=3, MOUT will equal MINP.
00220 *> \endverbatim
00221 *>
00222 *> \param[in,out] NAB
00223 *> \verbatim
00224 *>          NAB is INTEGER array, dimension (MMAX,2)
00225 *>          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
00226 *>          If IJOB=2, then on input, NAB(i,j) should be set.  It must
00227 *>             satisfy the condition:
00228 *>             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
00229 *>             which means that in interval i only eigenvalues
00230 *>             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
00231 *>             NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with
00232 *>             IJOB=1.
00233 *>             On output, NAB(i,j) will contain
00234 *>             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
00235 *>             the input interval that the output interval
00236 *>             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
00237 *>             the input values of NAB(k,1) and NAB(k,2).
00238 *>          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
00239 *>             unless N(w) > NVAL(i) for all search points  w , in which
00240 *>             case NAB(i,1) will not be modified, i.e., the output
00241 *>             value will be the same as the input value (modulo
00242 *>             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
00243 *>             for all search points  w , in which case NAB(i,2) will
00244 *>             not be modified.  Normally, NAB should be set to some
00245 *>             distinctive value(s) before SLAEBZ is called.
00246 *> \endverbatim
00247 *>
00248 *> \param[out] WORK
00249 *> \verbatim
00250 *>          WORK is REAL array, dimension (MMAX)
00251 *>          Workspace.
00252 *> \endverbatim
00253 *>
00254 *> \param[out] IWORK
00255 *> \verbatim
00256 *>          IWORK is INTEGER array, dimension (MMAX)
00257 *>          Workspace.
00258 *> \endverbatim
00259 *>
00260 *> \param[out] INFO
00261 *> \verbatim
00262 *>          INFO is INTEGER
00263 *>          = 0:       All intervals converged.
00264 *>          = 1--MMAX: The last INFO intervals did not converge.
00265 *>          = MMAX+1:  More than MMAX intervals were generated.
00266 *> \endverbatim
00267 *
00268 *  Authors:
00269 *  ========
00270 *
00271 *> \author Univ. of Tennessee 
00272 *> \author Univ. of California Berkeley 
00273 *> \author Univ. of Colorado Denver 
00274 *> \author NAG Ltd. 
00275 *
00276 *> \date November 2011
00277 *
00278 *> \ingroup auxOTHERauxiliary
00279 *
00280 *> \par Further Details:
00281 *  =====================
00282 *>
00283 *> \verbatim
00284 *>
00285 *>      This routine is intended to be called only by other LAPACK
00286 *>  routines, thus the interface is less user-friendly.  It is intended
00287 *>  for two purposes:
00288 *>
00289 *>  (a) finding eigenvalues.  In this case, SLAEBZ should have one or
00290 *>      more initial intervals set up in AB, and SLAEBZ should be called
00291 *>      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
00292 *>      Intervals with no eigenvalues would usually be thrown out at
00293 *>      this point.  Also, if not all the eigenvalues in an interval i
00294 *>      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
00295 *>      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
00296 *>      eigenvalue.  SLAEBZ is then called with IJOB=2 and MMAX
00297 *>      no smaller than the value of MOUT returned by the call with
00298 *>      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
00299 *>      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
00300 *>      tolerance specified by ABSTOL and RELTOL.
00301 *>
00302 *>  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
00303 *>      In this case, start with a Gershgorin interval  (a,b).  Set up
00304 *>      AB to contain 2 search intervals, both initially (a,b).  One
00305 *>      NVAL element should contain  f-1  and the other should contain  l
00306 *>      , while C should contain a and b, resp.  NAB(i,1) should be -1
00307 *>      and NAB(i,2) should be N+1, to flag an error if the desired
00308 *>      interval does not lie in (a,b).  SLAEBZ is then called with
00309 *>      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
00310 *>      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
00311 *>      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
00312 *>      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
00313 *>      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
00314 *>      w(l-r)=...=w(l+k) are handled similarly.
00315 *> \endverbatim
00316 *>
00317 *  =====================================================================
00318       SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
00319      $                   RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
00320      $                   NAB, WORK, IWORK, INFO )
00321 *
00322 *  -- LAPACK auxiliary routine (version 3.4.0) --
00323 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00324 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00325 *     November 2011
00326 *
00327 *     .. Scalar Arguments ..
00328       INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
00329       REAL               ABSTOL, PIVMIN, RELTOL
00330 *     ..
00331 *     .. Array Arguments ..
00332       INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
00333       REAL               AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
00334      $                   WORK( * )
00335 *     ..
00336 *
00337 *  =====================================================================
00338 *
00339 *     .. Parameters ..
00340       REAL               ZERO, TWO, HALF
00341       PARAMETER          ( ZERO = 0.0E0, TWO = 2.0E0,
00342      $                   HALF = 1.0E0 / TWO )
00343 *     ..
00344 *     .. Local Scalars ..
00345       INTEGER            ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
00346      $                   KLNEW
00347       REAL               TMP1, TMP2
00348 *     ..
00349 *     .. Intrinsic Functions ..
00350       INTRINSIC          ABS, MAX, MIN
00351 *     ..
00352 *     .. Executable Statements ..
00353 *
00354 *     Check for Errors
00355 *
00356       INFO = 0
00357       IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
00358          INFO = -1
00359          RETURN
00360       END IF
00361 *
00362 *     Initialize NAB
00363 *
00364       IF( IJOB.EQ.1 ) THEN
00365 *
00366 *        Compute the number of eigenvalues in the initial intervals.
00367 *
00368          MOUT = 0
00369          DO 30 JI = 1, MINP
00370             DO 20 JP = 1, 2
00371                TMP1 = D( 1 ) - AB( JI, JP )
00372                IF( ABS( TMP1 ).LT.PIVMIN )
00373      $            TMP1 = -PIVMIN
00374                NAB( JI, JP ) = 0
00375                IF( TMP1.LE.ZERO )
00376      $            NAB( JI, JP ) = 1
00377 *
00378                DO 10 J = 2, N
00379                   TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
00380                   IF( ABS( TMP1 ).LT.PIVMIN )
00381      $               TMP1 = -PIVMIN
00382                   IF( TMP1.LE.ZERO )
00383      $               NAB( JI, JP ) = NAB( JI, JP ) + 1
00384    10          CONTINUE
00385    20       CONTINUE
00386             MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
00387    30    CONTINUE
00388          RETURN
00389       END IF
00390 *
00391 *     Initialize for loop
00392 *
00393 *     KF and KL have the following meaning:
00394 *        Intervals 1,...,KF-1 have converged.
00395 *        Intervals KF,...,KL  still need to be refined.
00396 *
00397       KF = 1
00398       KL = MINP
00399 *
00400 *     If IJOB=2, initialize C.
00401 *     If IJOB=3, use the user-supplied starting point.
00402 *
00403       IF( IJOB.EQ.2 ) THEN
00404          DO 40 JI = 1, MINP
00405             C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
00406    40    CONTINUE
00407       END IF
00408 *
00409 *     Iteration loop
00410 *
00411       DO 130 JIT = 1, NITMAX
00412 *
00413 *        Loop over intervals
00414 *
00415          IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
00416 *
00417 *           Begin of Parallel Version of the loop
00418 *
00419             DO 60 JI = KF, KL
00420 *
00421 *              Compute N(c), the number of eigenvalues less than c
00422 *
00423                WORK( JI ) = D( 1 ) - C( JI )
00424                IWORK( JI ) = 0
00425                IF( WORK( JI ).LE.PIVMIN ) THEN
00426                   IWORK( JI ) = 1
00427                   WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
00428                END IF
00429 *
00430                DO 50 J = 2, N
00431                   WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
00432                   IF( WORK( JI ).LE.PIVMIN ) THEN
00433                      IWORK( JI ) = IWORK( JI ) + 1
00434                      WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
00435                   END IF
00436    50          CONTINUE
00437    60       CONTINUE
00438 *
00439             IF( IJOB.LE.2 ) THEN
00440 *
00441 *              IJOB=2: Choose all intervals containing eigenvalues.
00442 *
00443                KLNEW = KL
00444                DO 70 JI = KF, KL
00445 *
00446 *                 Insure that N(w) is monotone
00447 *
00448                   IWORK( JI ) = MIN( NAB( JI, 2 ),
00449      $                          MAX( NAB( JI, 1 ), IWORK( JI ) ) )
00450 *
00451 *                 Update the Queue -- add intervals if both halves
00452 *                 contain eigenvalues.
00453 *
00454                   IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
00455 *
00456 *                    No eigenvalue in the upper interval:
00457 *                    just use the lower interval.
00458 *
00459                      AB( JI, 2 ) = C( JI )
00460 *
00461                   ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
00462 *
00463 *                    No eigenvalue in the lower interval:
00464 *                    just use the upper interval.
00465 *
00466                      AB( JI, 1 ) = C( JI )
00467                   ELSE
00468                      KLNEW = KLNEW + 1
00469                      IF( KLNEW.LE.MMAX ) THEN
00470 *
00471 *                       Eigenvalue in both intervals -- add upper to
00472 *                       queue.
00473 *
00474                         AB( KLNEW, 2 ) = AB( JI, 2 )
00475                         NAB( KLNEW, 2 ) = NAB( JI, 2 )
00476                         AB( KLNEW, 1 ) = C( JI )
00477                         NAB( KLNEW, 1 ) = IWORK( JI )
00478                         AB( JI, 2 ) = C( JI )
00479                         NAB( JI, 2 ) = IWORK( JI )
00480                      ELSE
00481                         INFO = MMAX + 1
00482                      END IF
00483                   END IF
00484    70          CONTINUE
00485                IF( INFO.NE.0 )
00486      $            RETURN
00487                KL = KLNEW
00488             ELSE
00489 *
00490 *              IJOB=3: Binary search.  Keep only the interval containing
00491 *                      w   s.t. N(w) = NVAL
00492 *
00493                DO 80 JI = KF, KL
00494                   IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
00495                      AB( JI, 1 ) = C( JI )
00496                      NAB( JI, 1 ) = IWORK( JI )
00497                   END IF
00498                   IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
00499                      AB( JI, 2 ) = C( JI )
00500                      NAB( JI, 2 ) = IWORK( JI )
00501                   END IF
00502    80          CONTINUE
00503             END IF
00504 *
00505          ELSE
00506 *
00507 *           End of Parallel Version of the loop
00508 *
00509 *           Begin of Serial Version of the loop
00510 *
00511             KLNEW = KL
00512             DO 100 JI = KF, KL
00513 *
00514 *              Compute N(w), the number of eigenvalues less than w
00515 *
00516                TMP1 = C( JI )
00517                TMP2 = D( 1 ) - TMP1
00518                ITMP1 = 0
00519                IF( TMP2.LE.PIVMIN ) THEN
00520                   ITMP1 = 1
00521                   TMP2 = MIN( TMP2, -PIVMIN )
00522                END IF
00523 *
00524                DO 90 J = 2, N
00525                   TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
00526                   IF( TMP2.LE.PIVMIN ) THEN
00527                      ITMP1 = ITMP1 + 1
00528                      TMP2 = MIN( TMP2, -PIVMIN )
00529                   END IF
00530    90          CONTINUE
00531 *
00532                IF( IJOB.LE.2 ) THEN
00533 *
00534 *                 IJOB=2: Choose all intervals containing eigenvalues.
00535 *
00536 *                 Insure that N(w) is monotone
00537 *
00538                   ITMP1 = MIN( NAB( JI, 2 ),
00539      $                    MAX( NAB( JI, 1 ), ITMP1 ) )
00540 *
00541 *                 Update the Queue -- add intervals if both halves
00542 *                 contain eigenvalues.
00543 *
00544                   IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
00545 *
00546 *                    No eigenvalue in the upper interval:
00547 *                    just use the lower interval.
00548 *
00549                      AB( JI, 2 ) = TMP1
00550 *
00551                   ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
00552 *
00553 *                    No eigenvalue in the lower interval:
00554 *                    just use the upper interval.
00555 *
00556                      AB( JI, 1 ) = TMP1
00557                   ELSE IF( KLNEW.LT.MMAX ) THEN
00558 *
00559 *                    Eigenvalue in both intervals -- add upper to queue.
00560 *
00561                      KLNEW = KLNEW + 1
00562                      AB( KLNEW, 2 ) = AB( JI, 2 )
00563                      NAB( KLNEW, 2 ) = NAB( JI, 2 )
00564                      AB( KLNEW, 1 ) = TMP1
00565                      NAB( KLNEW, 1 ) = ITMP1
00566                      AB( JI, 2 ) = TMP1
00567                      NAB( JI, 2 ) = ITMP1
00568                   ELSE
00569                      INFO = MMAX + 1
00570                      RETURN
00571                   END IF
00572                ELSE
00573 *
00574 *                 IJOB=3: Binary search.  Keep only the interval
00575 *                         containing  w  s.t. N(w) = NVAL
00576 *
00577                   IF( ITMP1.LE.NVAL( JI ) ) THEN
00578                      AB( JI, 1 ) = TMP1
00579                      NAB( JI, 1 ) = ITMP1
00580                   END IF
00581                   IF( ITMP1.GE.NVAL( JI ) ) THEN
00582                      AB( JI, 2 ) = TMP1
00583                      NAB( JI, 2 ) = ITMP1
00584                   END IF
00585                END IF
00586   100       CONTINUE
00587             KL = KLNEW
00588 *
00589          END IF
00590 *
00591 *        Check for convergence
00592 *
00593          KFNEW = KF
00594          DO 110 JI = KF, KL
00595             TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
00596             TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
00597             IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
00598      $          NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
00599 *
00600 *              Converged -- Swap with position KFNEW,
00601 *                           then increment KFNEW
00602 *
00603                IF( JI.GT.KFNEW ) THEN
00604                   TMP1 = AB( JI, 1 )
00605                   TMP2 = AB( JI, 2 )
00606                   ITMP1 = NAB( JI, 1 )
00607                   ITMP2 = NAB( JI, 2 )
00608                   AB( JI, 1 ) = AB( KFNEW, 1 )
00609                   AB( JI, 2 ) = AB( KFNEW, 2 )
00610                   NAB( JI, 1 ) = NAB( KFNEW, 1 )
00611                   NAB( JI, 2 ) = NAB( KFNEW, 2 )
00612                   AB( KFNEW, 1 ) = TMP1
00613                   AB( KFNEW, 2 ) = TMP2
00614                   NAB( KFNEW, 1 ) = ITMP1
00615                   NAB( KFNEW, 2 ) = ITMP2
00616                   IF( IJOB.EQ.3 ) THEN
00617                      ITMP1 = NVAL( JI )
00618                      NVAL( JI ) = NVAL( KFNEW )
00619                      NVAL( KFNEW ) = ITMP1
00620                   END IF
00621                END IF
00622                KFNEW = KFNEW + 1
00623             END IF
00624   110    CONTINUE
00625          KF = KFNEW
00626 *
00627 *        Choose Midpoints
00628 *
00629          DO 120 JI = KF, KL
00630             C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
00631   120    CONTINUE
00632 *
00633 *        If no more intervals to refine, quit.
00634 *
00635          IF( KF.GT.KL )
00636      $      GO TO 140
00637   130 CONTINUE
00638 *
00639 *     Converged
00640 *
00641   140 CONTINUE
00642       INFO = MAX( KL+1-KF, 0 )
00643       MOUT = KL
00644 *
00645       RETURN
00646 *
00647 *     End of SLAEBZ
00648 *
00649       END
 All Files Functions