LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slaebz.f
Go to the documentation of this file.
1 *> \brief \b 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.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLAEBZ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaebz.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaebz.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaebz.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
22 * RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
23 * NAB, WORK, IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
27 * REAL ABSTOL, PIVMIN, RELTOL
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
31 * REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
32 * $ WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> SLAEBZ contains the iteration loops which compute and use the
42 *> function N(w), which is the count of eigenvalues of a symmetric
43 *> tridiagonal matrix T less than or equal to its argument w. It
44 *> performs a choice of two types of loops:
45 *>
46 *> IJOB=1, followed by
47 *> IJOB=2: It takes as input a list of intervals and returns a list of
48 *> sufficiently small intervals whose union contains the same
49 *> eigenvalues as the union of the original intervals.
50 *> The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
51 *> The output interval (AB(j,1),AB(j,2)] will contain
52 *> eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
53 *>
54 *> IJOB=3: It performs a binary search in each input interval
55 *> (AB(j,1),AB(j,2)] for a point w(j) such that
56 *> N(w(j))=NVAL(j), and uses C(j) as the starting point of
57 *> the search. If such a w(j) is found, then on output
58 *> AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output
59 *> (AB(j,1),AB(j,2)] will be a small interval containing the
60 *> point where N(w) jumps through NVAL(j), unless that point
61 *> lies outside the initial interval.
62 *>
63 *> Note that the intervals are in all cases half-open intervals,
64 *> i.e., of the form (a,b] , which includes b but not a .
65 *>
66 *> To avoid underflow, the matrix should be scaled so that its largest
67 *> element is no greater than overflow**(1/2) * underflow**(1/4)
68 *> in absolute value. To assure the most accurate computation
69 *> of small eigenvalues, the matrix should be scaled to be
70 *> not much smaller than that, either.
71 *>
72 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
73 *> Matrix", Report CS41, Computer Science Dept., Stanford
74 *> University, July 21, 1966
75 *>
76 *> Note: the arguments are, in general, *not* checked for unreasonable
77 *> values.
78 *> \endverbatim
79 *
80 * Arguments:
81 * ==========
82 *
83 *> \param[in] IJOB
84 *> \verbatim
85 *> IJOB is INTEGER
86 *> Specifies what is to be done:
87 *> = 1: Compute NAB for the initial intervals.
88 *> = 2: Perform bisection iteration to find eigenvalues of T.
89 *> = 3: Perform bisection iteration to invert N(w), i.e.,
90 *> to find a point which has a specified number of
91 *> eigenvalues of T to its left.
92 *> Other values will cause SLAEBZ to return with INFO=-1.
93 *> \endverbatim
94 *>
95 *> \param[in] NITMAX
96 *> \verbatim
97 *> NITMAX is INTEGER
98 *> The maximum number of "levels" of bisection to be
99 *> performed, i.e., an interval of width W will not be made
100 *> smaller than 2^(-NITMAX) * W. If not all intervals
101 *> have converged after NITMAX iterations, then INFO is set
102 *> to the number of non-converged intervals.
103 *> \endverbatim
104 *>
105 *> \param[in] N
106 *> \verbatim
107 *> N is INTEGER
108 *> The dimension n of the tridiagonal matrix T. It must be at
109 *> least 1.
110 *> \endverbatim
111 *>
112 *> \param[in] MMAX
113 *> \verbatim
114 *> MMAX is INTEGER
115 *> The maximum number of intervals. If more than MMAX intervals
116 *> are generated, then SLAEBZ will quit with INFO=MMAX+1.
117 *> \endverbatim
118 *>
119 *> \param[in] MINP
120 *> \verbatim
121 *> MINP is INTEGER
122 *> The initial number of intervals. It may not be greater than
123 *> MMAX.
124 *> \endverbatim
125 *>
126 *> \param[in] NBMIN
127 *> \verbatim
128 *> NBMIN is INTEGER
129 *> The smallest number of intervals that should be processed
130 *> using a vector loop. If zero, then only the scalar loop
131 *> will be used.
132 *> \endverbatim
133 *>
134 *> \param[in] ABSTOL
135 *> \verbatim
136 *> ABSTOL is REAL
137 *> The minimum (absolute) width of an interval. When an
138 *> interval is narrower than ABSTOL, or than RELTOL times the
139 *> larger (in magnitude) endpoint, then it is considered to be
140 *> sufficiently small, i.e., converged. This must be at least
141 *> zero.
142 *> \endverbatim
143 *>
144 *> \param[in] RELTOL
145 *> \verbatim
146 *> RELTOL is REAL
147 *> The minimum relative width of an interval. When an interval
148 *> is narrower than ABSTOL, or than RELTOL times the larger (in
149 *> magnitude) endpoint, then it is considered to be
150 *> sufficiently small, i.e., converged. Note: this should
151 *> always be at least radix*machine epsilon.
152 *> \endverbatim
153 *>
154 *> \param[in] PIVMIN
155 *> \verbatim
156 *> PIVMIN is REAL
157 *> The minimum absolute value of a "pivot" in the Sturm
158 *> sequence loop.
159 *> This must be at least max |e(j)**2|*safe_min and at
160 *> least safe_min, where safe_min is at least
161 *> the smallest number that can divide one without overflow.
162 *> \endverbatim
163 *>
164 *> \param[in] D
165 *> \verbatim
166 *> D is REAL array, dimension (N)
167 *> The diagonal elements of the tridiagonal matrix T.
168 *> \endverbatim
169 *>
170 *> \param[in] E
171 *> \verbatim
172 *> E is REAL array, dimension (N)
173 *> The offdiagonal elements of the tridiagonal matrix T in
174 *> positions 1 through N-1. E(N) is arbitrary.
175 *> \endverbatim
176 *>
177 *> \param[in] E2
178 *> \verbatim
179 *> E2 is REAL array, dimension (N)
180 *> The squares of the offdiagonal elements of the tridiagonal
181 *> matrix T. E2(N) is ignored.
182 *> \endverbatim
183 *>
184 *> \param[in,out] NVAL
185 *> \verbatim
186 *> NVAL is INTEGER array, dimension (MINP)
187 *> If IJOB=1 or 2, not referenced.
188 *> If IJOB=3, the desired values of N(w). The elements of NVAL
189 *> will be reordered to correspond with the intervals in AB.
190 *> Thus, NVAL(j) on output will not, in general be the same as
191 *> NVAL(j) on input, but it will correspond with the interval
192 *> (AB(j,1),AB(j,2)] on output.
193 *> \endverbatim
194 *>
195 *> \param[in,out] AB
196 *> \verbatim
197 *> AB is REAL array, dimension (MMAX,2)
198 *> The endpoints of the intervals. AB(j,1) is a(j), the left
199 *> endpoint of the j-th interval, and AB(j,2) is b(j), the
200 *> right endpoint of the j-th interval. The input intervals
201 *> will, in general, be modified, split, and reordered by the
202 *> calculation.
203 *> \endverbatim
204 *>
205 *> \param[in,out] C
206 *> \verbatim
207 *> C is REAL array, dimension (MMAX)
208 *> If IJOB=1, ignored.
209 *> If IJOB=2, workspace.
210 *> If IJOB=3, then on input C(j) should be initialized to the
211 *> first search point in the binary search.
212 *> \endverbatim
213 *>
214 *> \param[out] MOUT
215 *> \verbatim
216 *> MOUT is INTEGER
217 *> If IJOB=1, the number of eigenvalues in the intervals.
218 *> If IJOB=2 or 3, the number of intervals output.
219 *> If IJOB=3, MOUT will equal MINP.
220 *> \endverbatim
221 *>
222 *> \param[in,out] NAB
223 *> \verbatim
224 *> NAB is INTEGER array, dimension (MMAX,2)
225 *> If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
226 *> If IJOB=2, then on input, NAB(i,j) should be set. It must
227 *> satisfy the condition:
228 *> N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
229 *> which means that in interval i only eigenvalues
230 *> NAB(i,1)+1,...,NAB(i,2) will be considered. Usually,
231 *> NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with
232 *> IJOB=1.
233 *> On output, NAB(i,j) will contain
234 *> max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
235 *> the input interval that the output interval
236 *> (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
237 *> the input values of NAB(k,1) and NAB(k,2).
238 *> If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
239 *> unless N(w) > NVAL(i) for all search points w , in which
240 *> case NAB(i,1) will not be modified, i.e., the output
241 *> value will be the same as the input value (modulo
242 *> reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
243 *> for all search points w , in which case NAB(i,2) will
244 *> not be modified. Normally, NAB should be set to some
245 *> distinctive value(s) before SLAEBZ is called.
246 *> \endverbatim
247 *>
248 *> \param[out] WORK
249 *> \verbatim
250 *> WORK is REAL array, dimension (MMAX)
251 *> Workspace.
252 *> \endverbatim
253 *>
254 *> \param[out] IWORK
255 *> \verbatim
256 *> IWORK is INTEGER array, dimension (MMAX)
257 *> Workspace.
258 *> \endverbatim
259 *>
260 *> \param[out] INFO
261 *> \verbatim
262 *> INFO is INTEGER
263 *> = 0: All intervals converged.
264 *> = 1--MMAX: The last INFO intervals did not converge.
265 *> = MMAX+1: More than MMAX intervals were generated.
266 *> \endverbatim
267 *
268 * Authors:
269 * ========
270 *
271 *> \author Univ. of Tennessee
272 *> \author Univ. of California Berkeley
273 *> \author Univ. of Colorado Denver
274 *> \author NAG Ltd.
275 *
276 *> \date September 2012
277 *
278 *> \ingroup auxOTHERauxiliary
279 *
280 *> \par Further Details:
281 * =====================
282 *>
283 *> \verbatim
284 *>
285 *> This routine is intended to be called only by other LAPACK
286 *> routines, thus the interface is less user-friendly. It is intended
287 *> for two purposes:
288 *>
289 *> (a) finding eigenvalues. In this case, SLAEBZ should have one or
290 *> more initial intervals set up in AB, and SLAEBZ should be called
291 *> with IJOB=1. This sets up NAB, and also counts the eigenvalues.
292 *> Intervals with no eigenvalues would usually be thrown out at
293 *> this point. Also, if not all the eigenvalues in an interval i
294 *> are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
295 *> For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
296 *> eigenvalue. SLAEBZ is then called with IJOB=2 and MMAX
297 *> no smaller than the value of MOUT returned by the call with
298 *> IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
299 *> through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
300 *> tolerance specified by ABSTOL and RELTOL.
301 *>
302 *> (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
303 *> In this case, start with a Gershgorin interval (a,b). Set up
304 *> AB to contain 2 search intervals, both initially (a,b). One
305 *> NVAL element should contain f-1 and the other should contain l
306 *> , while C should contain a and b, resp. NAB(i,1) should be -1
307 *> and NAB(i,2) should be N+1, to flag an error if the desired
308 *> interval does not lie in (a,b). SLAEBZ is then called with
309 *> IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
310 *> j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
311 *> if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
312 *> >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
313 *> N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
314 *> w(l-r)=...=w(l+k) are handled similarly.
315 *> \endverbatim
316 *>
317 * =====================================================================
318  SUBROUTINE slaebz( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
319  $ reltol, pivmin, d, e, e2, nval, ab, c, mout,
320  $ nab, work, iwork, info )
321 *
322 * -- LAPACK auxiliary routine (version 3.4.2) --
323 * -- LAPACK is a software package provided by Univ. of Tennessee, --
324 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
325 * September 2012
326 *
327 * .. Scalar Arguments ..
328  INTEGER ijob, info, minp, mmax, mout, n, nbmin, nitmax
329  REAL abstol, pivmin, reltol
330 * ..
331 * .. Array Arguments ..
332  INTEGER iwork( * ), nab( mmax, * ), nval( * )
333  REAL ab( mmax, * ), c( * ), d( * ), e( * ), e2( * ),
334  $ work( * )
335 * ..
336 *
337 * =====================================================================
338 *
339 * .. Parameters ..
340  REAL zero, two, half
341  parameter( zero = 0.0e0, two = 2.0e0,
342  $ half = 1.0e0 / two )
343 * ..
344 * .. Local Scalars ..
345  INTEGER itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl,
346  $ klnew
347  REAL tmp1, tmp2
348 * ..
349 * .. Intrinsic Functions ..
350  INTRINSIC abs, max, min
351 * ..
352 * .. Executable Statements ..
353 *
354 * Check for Errors
355 *
356  info = 0
357  IF( ijob.LT.1 .OR. ijob.GT.3 ) THEN
358  info = -1
359  return
360  END IF
361 *
362 * Initialize NAB
363 *
364  IF( ijob.EQ.1 ) THEN
365 *
366 * Compute the number of eigenvalues in the initial intervals.
367 *
368  mout = 0
369  DO 30 ji = 1, minp
370  DO 20 jp = 1, 2
371  tmp1 = d( 1 ) - ab( ji, jp )
372  IF( abs( tmp1 ).LT.pivmin )
373  $ tmp1 = -pivmin
374  nab( ji, jp ) = 0
375  IF( tmp1.LE.zero )
376  $ nab( ji, jp ) = 1
377 *
378  DO 10 j = 2, n
379  tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
380  IF( abs( tmp1 ).LT.pivmin )
381  $ tmp1 = -pivmin
382  IF( tmp1.LE.zero )
383  $ nab( ji, jp ) = nab( ji, jp ) + 1
384  10 continue
385  20 continue
386  mout = mout + nab( ji, 2 ) - nab( ji, 1 )
387  30 continue
388  return
389  END IF
390 *
391 * Initialize for loop
392 *
393 * KF and KL have the following meaning:
394 * Intervals 1,...,KF-1 have converged.
395 * Intervals KF,...,KL still need to be refined.
396 *
397  kf = 1
398  kl = minp
399 *
400 * If IJOB=2, initialize C.
401 * If IJOB=3, use the user-supplied starting point.
402 *
403  IF( ijob.EQ.2 ) THEN
404  DO 40 ji = 1, minp
405  c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
406  40 continue
407  END IF
408 *
409 * Iteration loop
410 *
411  DO 130 jit = 1, nitmax
412 *
413 * Loop over intervals
414 *
415  IF( kl-kf+1.GE.nbmin .AND. nbmin.GT.0 ) THEN
416 *
417 * Begin of Parallel Version of the loop
418 *
419  DO 60 ji = kf, kl
420 *
421 * Compute N(c), the number of eigenvalues less than c
422 *
423  work( ji ) = d( 1 ) - c( ji )
424  iwork( ji ) = 0
425  IF( work( ji ).LE.pivmin ) THEN
426  iwork( ji ) = 1
427  work( ji ) = min( work( ji ), -pivmin )
428  END IF
429 *
430  DO 50 j = 2, n
431  work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
432  IF( work( ji ).LE.pivmin ) THEN
433  iwork( ji ) = iwork( ji ) + 1
434  work( ji ) = min( work( ji ), -pivmin )
435  END IF
436  50 continue
437  60 continue
438 *
439  IF( ijob.LE.2 ) THEN
440 *
441 * IJOB=2: Choose all intervals containing eigenvalues.
442 *
443  klnew = kl
444  DO 70 ji = kf, kl
445 *
446 * Insure that N(w) is monotone
447 *
448  iwork( ji ) = min( nab( ji, 2 ),
449  $ max( nab( ji, 1 ), iwork( ji ) ) )
450 *
451 * Update the Queue -- add intervals if both halves
452 * contain eigenvalues.
453 *
454  IF( iwork( ji ).EQ.nab( ji, 2 ) ) THEN
455 *
456 * No eigenvalue in the upper interval:
457 * just use the lower interval.
458 *
459  ab( ji, 2 ) = c( ji )
460 *
461  ELSE IF( iwork( ji ).EQ.nab( ji, 1 ) ) THEN
462 *
463 * No eigenvalue in the lower interval:
464 * just use the upper interval.
465 *
466  ab( ji, 1 ) = c( ji )
467  ELSE
468  klnew = klnew + 1
469  IF( klnew.LE.mmax ) THEN
470 *
471 * Eigenvalue in both intervals -- add upper to
472 * queue.
473 *
474  ab( klnew, 2 ) = ab( ji, 2 )
475  nab( klnew, 2 ) = nab( ji, 2 )
476  ab( klnew, 1 ) = c( ji )
477  nab( klnew, 1 ) = iwork( ji )
478  ab( ji, 2 ) = c( ji )
479  nab( ji, 2 ) = iwork( ji )
480  ELSE
481  info = mmax + 1
482  END IF
483  END IF
484  70 continue
485  IF( info.NE.0 )
486  $ return
487  kl = klnew
488  ELSE
489 *
490 * IJOB=3: Binary search. Keep only the interval containing
491 * w s.t. N(w) = NVAL
492 *
493  DO 80 ji = kf, kl
494  IF( iwork( ji ).LE.nval( ji ) ) THEN
495  ab( ji, 1 ) = c( ji )
496  nab( ji, 1 ) = iwork( ji )
497  END IF
498  IF( iwork( ji ).GE.nval( ji ) ) THEN
499  ab( ji, 2 ) = c( ji )
500  nab( ji, 2 ) = iwork( ji )
501  END IF
502  80 continue
503  END IF
504 *
505  ELSE
506 *
507 * End of Parallel Version of the loop
508 *
509 * Begin of Serial Version of the loop
510 *
511  klnew = kl
512  DO 100 ji = kf, kl
513 *
514 * Compute N(w), the number of eigenvalues less than w
515 *
516  tmp1 = c( ji )
517  tmp2 = d( 1 ) - tmp1
518  itmp1 = 0
519  IF( tmp2.LE.pivmin ) THEN
520  itmp1 = 1
521  tmp2 = min( tmp2, -pivmin )
522  END IF
523 *
524  DO 90 j = 2, n
525  tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
526  IF( tmp2.LE.pivmin ) THEN
527  itmp1 = itmp1 + 1
528  tmp2 = min( tmp2, -pivmin )
529  END IF
530  90 continue
531 *
532  IF( ijob.LE.2 ) THEN
533 *
534 * IJOB=2: Choose all intervals containing eigenvalues.
535 *
536 * Insure that N(w) is monotone
537 *
538  itmp1 = min( nab( ji, 2 ),
539  $ max( nab( ji, 1 ), itmp1 ) )
540 *
541 * Update the Queue -- add intervals if both halves
542 * contain eigenvalues.
543 *
544  IF( itmp1.EQ.nab( ji, 2 ) ) THEN
545 *
546 * No eigenvalue in the upper interval:
547 * just use the lower interval.
548 *
549  ab( ji, 2 ) = tmp1
550 *
551  ELSE IF( itmp1.EQ.nab( ji, 1 ) ) THEN
552 *
553 * No eigenvalue in the lower interval:
554 * just use the upper interval.
555 *
556  ab( ji, 1 ) = tmp1
557  ELSE IF( klnew.LT.mmax ) THEN
558 *
559 * Eigenvalue in both intervals -- add upper to queue.
560 *
561  klnew = klnew + 1
562  ab( klnew, 2 ) = ab( ji, 2 )
563  nab( klnew, 2 ) = nab( ji, 2 )
564  ab( klnew, 1 ) = tmp1
565  nab( klnew, 1 ) = itmp1
566  ab( ji, 2 ) = tmp1
567  nab( ji, 2 ) = itmp1
568  ELSE
569  info = mmax + 1
570  return
571  END IF
572  ELSE
573 *
574 * IJOB=3: Binary search. Keep only the interval
575 * containing w s.t. N(w) = NVAL
576 *
577  IF( itmp1.LE.nval( ji ) ) THEN
578  ab( ji, 1 ) = tmp1
579  nab( ji, 1 ) = itmp1
580  END IF
581  IF( itmp1.GE.nval( ji ) ) THEN
582  ab( ji, 2 ) = tmp1
583  nab( ji, 2 ) = itmp1
584  END IF
585  END IF
586  100 continue
587  kl = klnew
588 *
589  END IF
590 *
591 * Check for convergence
592 *
593  kfnew = kf
594  DO 110 ji = kf, kl
595  tmp1 = abs( ab( ji, 2 )-ab( ji, 1 ) )
596  tmp2 = max( abs( ab( ji, 2 ) ), abs( ab( ji, 1 ) ) )
597  IF( tmp1.LT.max( abstol, pivmin, reltol*tmp2 ) .OR.
598  $ nab( ji, 1 ).GE.nab( ji, 2 ) ) THEN
599 *
600 * Converged -- Swap with position KFNEW,
601 * then increment KFNEW
602 *
603  IF( ji.GT.kfnew ) THEN
604  tmp1 = ab( ji, 1 )
605  tmp2 = ab( ji, 2 )
606  itmp1 = nab( ji, 1 )
607  itmp2 = nab( ji, 2 )
608  ab( ji, 1 ) = ab( kfnew, 1 )
609  ab( ji, 2 ) = ab( kfnew, 2 )
610  nab( ji, 1 ) = nab( kfnew, 1 )
611  nab( ji, 2 ) = nab( kfnew, 2 )
612  ab( kfnew, 1 ) = tmp1
613  ab( kfnew, 2 ) = tmp2
614  nab( kfnew, 1 ) = itmp1
615  nab( kfnew, 2 ) = itmp2
616  IF( ijob.EQ.3 ) THEN
617  itmp1 = nval( ji )
618  nval( ji ) = nval( kfnew )
619  nval( kfnew ) = itmp1
620  END IF
621  END IF
622  kfnew = kfnew + 1
623  END IF
624  110 continue
625  kf = kfnew
626 *
627 * Choose Midpoints
628 *
629  DO 120 ji = kf, kl
630  c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
631  120 continue
632 *
633 * If no more intervals to refine, quit.
634 *
635  IF( kf.GT.kl )
636  $ go to 140
637  130 continue
638 *
639 * Converged
640 *
641  140 continue
642  info = max( kl+1-kf, 0 )
643  mout = kl
644 *
645  return
646 *
647 * End of SLAEBZ
648 *
649  END