LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlaebz.f
Go to the documentation of this file.
1 *> \brief \b DLAEBZ 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 DLAEBZ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaebz.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaebz.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaebz.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAEBZ( 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 * DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
31 * DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
32 * $ WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DLAEBZ 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 DLAEBZ 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 DLAEBZ 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION array, dimension (N)
167 *> The diagonal elements of the tridiagonal matrix T.
168 *> \endverbatim
169 *>
170 *> \param[in] E
171 *> \verbatim
172 *> E is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLAEBZ 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 DLAEBZ is called.
246 *> \endverbatim
247 *>
248 *> \param[out] WORK
249 *> \verbatim
250 *> WORK is DOUBLE PRECISION 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 *> \ingroup OTHERauxiliary
277 *
278 *> \par Further Details:
279 * =====================
280 *>
281 *> \verbatim
282 *>
283 *> This routine is intended to be called only by other LAPACK
284 *> routines, thus the interface is less user-friendly. It is intended
285 *> for two purposes:
286 *>
287 *> (a) finding eigenvalues. In this case, DLAEBZ should have one or
288 *> more initial intervals set up in AB, and DLAEBZ should be called
289 *> with IJOB=1. This sets up NAB, and also counts the eigenvalues.
290 *> Intervals with no eigenvalues would usually be thrown out at
291 *> this point. Also, if not all the eigenvalues in an interval i
292 *> are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
293 *> For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
294 *> eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX
295 *> no smaller than the value of MOUT returned by the call with
296 *> IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
297 *> through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
298 *> tolerance specified by ABSTOL and RELTOL.
299 *>
300 *> (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
301 *> In this case, start with a Gershgorin interval (a,b). Set up
302 *> AB to contain 2 search intervals, both initially (a,b). One
303 *> NVAL element should contain f-1 and the other should contain l
304 *> , while C should contain a and b, resp. NAB(i,1) should be -1
305 *> and NAB(i,2) should be N+1, to flag an error if the desired
306 *> interval does not lie in (a,b). DLAEBZ is then called with
307 *> IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
308 *> j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
309 *> if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
310 *> >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
311 *> N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
312 *> w(l-r)=...=w(l+k) are handled similarly.
313 *> \endverbatim
314 *>
315 * =====================================================================
316  SUBROUTINE dlaebz( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
317  $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
318  $ NAB, WORK, IWORK, INFO )
319 *
320 * -- LAPACK auxiliary routine --
321 * -- LAPACK is a software package provided by Univ. of Tennessee, --
322 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
323 *
324 * .. Scalar Arguments ..
325  INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
326  DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
327 * ..
328 * .. Array Arguments ..
329  INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
330  DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
331  $ work( * )
332 * ..
333 *
334 * =====================================================================
335 *
336 * .. Parameters ..
337  DOUBLE PRECISION ZERO, TWO, HALF
338  PARAMETER ( ZERO = 0.0d0, two = 2.0d0,
339  $ half = 1.0d0 / two )
340 * ..
341 * .. Local Scalars ..
342  INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
343  $ KLNEW
344  DOUBLE PRECISION TMP1, TMP2
345 * ..
346 * .. Intrinsic Functions ..
347  INTRINSIC abs, max, min
348 * ..
349 * .. Executable Statements ..
350 *
351 * Check for Errors
352 *
353  info = 0
354  IF( ijob.LT.1 .OR. ijob.GT.3 ) THEN
355  info = -1
356  RETURN
357  END IF
358 *
359 * Initialize NAB
360 *
361  IF( ijob.EQ.1 ) THEN
362 *
363 * Compute the number of eigenvalues in the initial intervals.
364 *
365  mout = 0
366  DO 30 ji = 1, minp
367  DO 20 jp = 1, 2
368  tmp1 = d( 1 ) - ab( ji, jp )
369  IF( abs( tmp1 ).LT.pivmin )
370  $ tmp1 = -pivmin
371  nab( ji, jp ) = 0
372  IF( tmp1.LE.zero )
373  $ nab( ji, jp ) = 1
374 *
375  DO 10 j = 2, n
376  tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
377  IF( abs( tmp1 ).LT.pivmin )
378  $ tmp1 = -pivmin
379  IF( tmp1.LE.zero )
380  $ nab( ji, jp ) = nab( ji, jp ) + 1
381  10 CONTINUE
382  20 CONTINUE
383  mout = mout + nab( ji, 2 ) - nab( ji, 1 )
384  30 CONTINUE
385  RETURN
386  END IF
387 *
388 * Initialize for loop
389 *
390 * KF and KL have the following meaning:
391 * Intervals 1,...,KF-1 have converged.
392 * Intervals KF,...,KL still need to be refined.
393 *
394  kf = 1
395  kl = minp
396 *
397 * If IJOB=2, initialize C.
398 * If IJOB=3, use the user-supplied starting point.
399 *
400  IF( ijob.EQ.2 ) THEN
401  DO 40 ji = 1, minp
402  c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
403  40 CONTINUE
404  END IF
405 *
406 * Iteration loop
407 *
408  DO 130 jit = 1, nitmax
409 *
410 * Loop over intervals
411 *
412  IF( kl-kf+1.GE.nbmin .AND. nbmin.GT.0 ) THEN
413 *
414 * Begin of Parallel Version of the loop
415 *
416  DO 60 ji = kf, kl
417 *
418 * Compute N(c), the number of eigenvalues less than c
419 *
420  work( ji ) = d( 1 ) - c( ji )
421  iwork( ji ) = 0
422  IF( work( ji ).LE.pivmin ) THEN
423  iwork( ji ) = 1
424  work( ji ) = min( work( ji ), -pivmin )
425  END IF
426 *
427  DO 50 j = 2, n
428  work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
429  IF( work( ji ).LE.pivmin ) THEN
430  iwork( ji ) = iwork( ji ) + 1
431  work( ji ) = min( work( ji ), -pivmin )
432  END IF
433  50 CONTINUE
434  60 CONTINUE
435 *
436  IF( ijob.LE.2 ) THEN
437 *
438 * IJOB=2: Choose all intervals containing eigenvalues.
439 *
440  klnew = kl
441  DO 70 ji = kf, kl
442 *
443 * Insure that N(w) is monotone
444 *
445  iwork( ji ) = min( nab( ji, 2 ),
446  $ max( nab( ji, 1 ), iwork( ji ) ) )
447 *
448 * Update the Queue -- add intervals if both halves
449 * contain eigenvalues.
450 *
451  IF( iwork( ji ).EQ.nab( ji, 2 ) ) THEN
452 *
453 * No eigenvalue in the upper interval:
454 * just use the lower interval.
455 *
456  ab( ji, 2 ) = c( ji )
457 *
458  ELSE IF( iwork( ji ).EQ.nab( ji, 1 ) ) THEN
459 *
460 * No eigenvalue in the lower interval:
461 * just use the upper interval.
462 *
463  ab( ji, 1 ) = c( ji )
464  ELSE
465  klnew = klnew + 1
466  IF( klnew.LE.mmax ) THEN
467 *
468 * Eigenvalue in both intervals -- add upper to
469 * queue.
470 *
471  ab( klnew, 2 ) = ab( ji, 2 )
472  nab( klnew, 2 ) = nab( ji, 2 )
473  ab( klnew, 1 ) = c( ji )
474  nab( klnew, 1 ) = iwork( ji )
475  ab( ji, 2 ) = c( ji )
476  nab( ji, 2 ) = iwork( ji )
477  ELSE
478  info = mmax + 1
479  END IF
480  END IF
481  70 CONTINUE
482  IF( info.NE.0 )
483  $ RETURN
484  kl = klnew
485  ELSE
486 *
487 * IJOB=3: Binary search. Keep only the interval containing
488 * w s.t. N(w) = NVAL
489 *
490  DO 80 ji = kf, kl
491  IF( iwork( ji ).LE.nval( ji ) ) THEN
492  ab( ji, 1 ) = c( ji )
493  nab( ji, 1 ) = iwork( ji )
494  END IF
495  IF( iwork( ji ).GE.nval( ji ) ) THEN
496  ab( ji, 2 ) = c( ji )
497  nab( ji, 2 ) = iwork( ji )
498  END IF
499  80 CONTINUE
500  END IF
501 *
502  ELSE
503 *
504 * End of Parallel Version of the loop
505 *
506 * Begin of Serial Version of the loop
507 *
508  klnew = kl
509  DO 100 ji = kf, kl
510 *
511 * Compute N(w), the number of eigenvalues less than w
512 *
513  tmp1 = c( ji )
514  tmp2 = d( 1 ) - tmp1
515  itmp1 = 0
516  IF( tmp2.LE.pivmin ) THEN
517  itmp1 = 1
518  tmp2 = min( tmp2, -pivmin )
519  END IF
520 *
521  DO 90 j = 2, n
522  tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
523  IF( tmp2.LE.pivmin ) THEN
524  itmp1 = itmp1 + 1
525  tmp2 = min( tmp2, -pivmin )
526  END IF
527  90 CONTINUE
528 *
529  IF( ijob.LE.2 ) THEN
530 *
531 * IJOB=2: Choose all intervals containing eigenvalues.
532 *
533 * Insure that N(w) is monotone
534 *
535  itmp1 = min( nab( ji, 2 ),
536  $ max( nab( ji, 1 ), itmp1 ) )
537 *
538 * Update the Queue -- add intervals if both halves
539 * contain eigenvalues.
540 *
541  IF( itmp1.EQ.nab( ji, 2 ) ) THEN
542 *
543 * No eigenvalue in the upper interval:
544 * just use the lower interval.
545 *
546  ab( ji, 2 ) = tmp1
547 *
548  ELSE IF( itmp1.EQ.nab( ji, 1 ) ) THEN
549 *
550 * No eigenvalue in the lower interval:
551 * just use the upper interval.
552 *
553  ab( ji, 1 ) = tmp1
554  ELSE IF( klnew.LT.mmax ) THEN
555 *
556 * Eigenvalue in both intervals -- add upper to queue.
557 *
558  klnew = klnew + 1
559  ab( klnew, 2 ) = ab( ji, 2 )
560  nab( klnew, 2 ) = nab( ji, 2 )
561  ab( klnew, 1 ) = tmp1
562  nab( klnew, 1 ) = itmp1
563  ab( ji, 2 ) = tmp1
564  nab( ji, 2 ) = itmp1
565  ELSE
566  info = mmax + 1
567  RETURN
568  END IF
569  ELSE
570 *
571 * IJOB=3: Binary search. Keep only the interval
572 * containing w s.t. N(w) = NVAL
573 *
574  IF( itmp1.LE.nval( ji ) ) THEN
575  ab( ji, 1 ) = tmp1
576  nab( ji, 1 ) = itmp1
577  END IF
578  IF( itmp1.GE.nval( ji ) ) THEN
579  ab( ji, 2 ) = tmp1
580  nab( ji, 2 ) = itmp1
581  END IF
582  END IF
583  100 CONTINUE
584  kl = klnew
585 *
586  END IF
587 *
588 * Check for convergence
589 *
590  kfnew = kf
591  DO 110 ji = kf, kl
592  tmp1 = abs( ab( ji, 2 )-ab( ji, 1 ) )
593  tmp2 = max( abs( ab( ji, 2 ) ), abs( ab( ji, 1 ) ) )
594  IF( tmp1.LT.max( abstol, pivmin, reltol*tmp2 ) .OR.
595  $ nab( ji, 1 ).GE.nab( ji, 2 ) ) THEN
596 *
597 * Converged -- Swap with position KFNEW,
598 * then increment KFNEW
599 *
600  IF( ji.GT.kfnew ) THEN
601  tmp1 = ab( ji, 1 )
602  tmp2 = ab( ji, 2 )
603  itmp1 = nab( ji, 1 )
604  itmp2 = nab( ji, 2 )
605  ab( ji, 1 ) = ab( kfnew, 1 )
606  ab( ji, 2 ) = ab( kfnew, 2 )
607  nab( ji, 1 ) = nab( kfnew, 1 )
608  nab( ji, 2 ) = nab( kfnew, 2 )
609  ab( kfnew, 1 ) = tmp1
610  ab( kfnew, 2 ) = tmp2
611  nab( kfnew, 1 ) = itmp1
612  nab( kfnew, 2 ) = itmp2
613  IF( ijob.EQ.3 ) THEN
614  itmp1 = nval( ji )
615  nval( ji ) = nval( kfnew )
616  nval( kfnew ) = itmp1
617  END IF
618  END IF
619  kfnew = kfnew + 1
620  END IF
621  110 CONTINUE
622  kf = kfnew
623 *
624 * Choose Midpoints
625 *
626  DO 120 ji = kf, kl
627  c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
628  120 CONTINUE
629 *
630 * If no more intervals to refine, quit.
631 *
632  IF( kf.GT.kl )
633  $ GO TO 140
634  130 CONTINUE
635 *
636 * Converged
637 *
638  140 CONTINUE
639  info = max( kl+1-kf, 0 )
640  mout = kl
641 *
642  RETURN
643 *
644 * End of DLAEBZ
645 *
646  END
subroutine dlaebz(IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than ...
Definition: dlaebz.f:319