LAPACK 3.3.0

dlaebz.f

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