LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dstebz.f
Go to the documentation of this file.
1 *> \brief \b DSTEBZ
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSTEBZ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
22 * M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER ORDER, RANGE
27 * INTEGER IL, INFO, IU, M, N, NSPLIT
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
32 * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DSTEBZ computes the eigenvalues of a symmetric tridiagonal
42 *> matrix T. The user may ask for all eigenvalues, all eigenvalues
43 *> in the half-open interval (VL, VU], or the IL-th through IU-th
44 *> eigenvalues.
45 *>
46 *> To avoid overflow, the matrix must be scaled so that its
47 *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
48 *> accuracy, it should not be much smaller than that.
49 *>
50 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
51 *> Matrix", Report CS41, Computer Science Dept., Stanford
52 *> University, July 21, 1966.
53 *> \endverbatim
54 *
55 * Arguments:
56 * ==========
57 *
58 *> \param[in] RANGE
59 *> \verbatim
60 *> RANGE is CHARACTER*1
61 *> = 'A': ("All") all eigenvalues will be found.
62 *> = 'V': ("Value") all eigenvalues in the half-open interval
63 *> (VL, VU] will be found.
64 *> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
65 *> entire matrix) will be found.
66 *> \endverbatim
67 *>
68 *> \param[in] ORDER
69 *> \verbatim
70 *> ORDER is CHARACTER*1
71 *> = 'B': ("By Block") the eigenvalues will be grouped by
72 *> split-off block (see IBLOCK, ISPLIT) and
73 *> ordered from smallest to largest within
74 *> the block.
75 *> = 'E': ("Entire matrix")
76 *> the eigenvalues for the entire matrix
77 *> will be ordered from smallest to
78 *> largest.
79 *> \endverbatim
80 *>
81 *> \param[in] N
82 *> \verbatim
83 *> N is INTEGER
84 *> The order of the tridiagonal matrix T. N >= 0.
85 *> \endverbatim
86 *>
87 *> \param[in] VL
88 *> \verbatim
89 *> VL is DOUBLE PRECISION
90 *> \endverbatim
91 *>
92 *> \param[in] VU
93 *> \verbatim
94 *> VU is DOUBLE PRECISION
95 *>
96 *> If RANGE='V', the lower and upper bounds of the interval to
97 *> be searched for eigenvalues. Eigenvalues less than or equal
98 *> to VL, or greater than VU, will not be returned. VL < VU.
99 *> Not referenced if RANGE = 'A' or 'I'.
100 *> \endverbatim
101 *>
102 *> \param[in] IL
103 *> \verbatim
104 *> IL is INTEGER
105 *> \endverbatim
106 *>
107 *> \param[in] IU
108 *> \verbatim
109 *> IU is INTEGER
110 *>
111 *> If RANGE='I', the indices (in ascending order) of the
112 *> smallest and largest eigenvalues to be returned.
113 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
114 *> Not referenced if RANGE = 'A' or 'V'.
115 *> \endverbatim
116 *>
117 *> \param[in] ABSTOL
118 *> \verbatim
119 *> ABSTOL is DOUBLE PRECISION
120 *> The absolute tolerance for the eigenvalues. An eigenvalue
121 *> (or cluster) is considered to be located if it has been
122 *> determined to lie in an interval whose width is ABSTOL or
123 *> less. If ABSTOL is less than or equal to zero, then ULP*|T|
124 *> will be used, where |T| means the 1-norm of T.
125 *>
126 *> Eigenvalues will be computed most accurately when ABSTOL is
127 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
128 *> \endverbatim
129 *>
130 *> \param[in] D
131 *> \verbatim
132 *> D is DOUBLE PRECISION array, dimension (N)
133 *> The n diagonal elements of the tridiagonal matrix T.
134 *> \endverbatim
135 *>
136 *> \param[in] E
137 *> \verbatim
138 *> E is DOUBLE PRECISION array, dimension (N-1)
139 *> The (n-1) off-diagonal elements of the tridiagonal matrix T.
140 *> \endverbatim
141 *>
142 *> \param[out] M
143 *> \verbatim
144 *> M is INTEGER
145 *> The actual number of eigenvalues found. 0 <= M <= N.
146 *> (See also the description of INFO=2,3.)
147 *> \endverbatim
148 *>
149 *> \param[out] NSPLIT
150 *> \verbatim
151 *> NSPLIT is INTEGER
152 *> The number of diagonal blocks in the matrix T.
153 *> 1 <= NSPLIT <= N.
154 *> \endverbatim
155 *>
156 *> \param[out] W
157 *> \verbatim
158 *> W is DOUBLE PRECISION array, dimension (N)
159 *> On exit, the first M elements of W will contain the
160 *> eigenvalues. (DSTEBZ may use the remaining N-M elements as
161 *> workspace.)
162 *> \endverbatim
163 *>
164 *> \param[out] IBLOCK
165 *> \verbatim
166 *> IBLOCK is INTEGER array, dimension (N)
167 *> At each row/column j where E(j) is zero or small, the
168 *> matrix T is considered to split into a block diagonal
169 *> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
170 *> block (from 1 to the number of blocks) the eigenvalue W(i)
171 *> belongs. (DSTEBZ may use the remaining N-M elements as
172 *> workspace.)
173 *> \endverbatim
174 *>
175 *> \param[out] ISPLIT
176 *> \verbatim
177 *> ISPLIT is INTEGER array, dimension (N)
178 *> The splitting points, at which T breaks up into submatrices.
179 *> The first submatrix consists of rows/columns 1 to ISPLIT(1),
180 *> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
181 *> etc., and the NSPLIT-th consists of rows/columns
182 *> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
183 *> (Only the first NSPLIT elements will actually be used, but
184 *> since the user cannot know a priori what value NSPLIT will
185 *> have, N words must be reserved for ISPLIT.)
186 *> \endverbatim
187 *>
188 *> \param[out] WORK
189 *> \verbatim
190 *> WORK is DOUBLE PRECISION array, dimension (4*N)
191 *> \endverbatim
192 *>
193 *> \param[out] IWORK
194 *> \verbatim
195 *> IWORK is INTEGER array, dimension (3*N)
196 *> \endverbatim
197 *>
198 *> \param[out] INFO
199 *> \verbatim
200 *> INFO is INTEGER
201 *> = 0: successful exit
202 *> < 0: if INFO = -i, the i-th argument had an illegal value
203 *> > 0: some or all of the eigenvalues failed to converge or
204 *> were not computed:
205 *> =1 or 3: Bisection failed to converge for some
206 *> eigenvalues; these eigenvalues are flagged by a
207 *> negative block number. The effect is that the
208 *> eigenvalues may not be as accurate as the
209 *> absolute and relative tolerances. This is
210 *> generally caused by unexpectedly inaccurate
211 *> arithmetic.
212 *> =2 or 3: RANGE='I' only: Not all of the eigenvalues
213 *> IL:IU were found.
214 *> Effect: M < IU+1-IL
215 *> Cause: non-monotonic arithmetic, causing the
216 *> Sturm sequence to be non-monotonic.
217 *> Cure: recalculate, using RANGE='A', and pick
218 *> out eigenvalues IL:IU. In some cases,
219 *> increasing the PARAMETER "FUDGE" may
220 *> make things work.
221 *> = 4: RANGE='I', and the Gershgorin interval
222 *> initially used was too small. No eigenvalues
223 *> were computed.
224 *> Probable cause: your machine has sloppy
225 *> floating-point arithmetic.
226 *> Cure: Increase the PARAMETER "FUDGE",
227 *> recompile, and try again.
228 *> \endverbatim
229 *
230 *> \par Internal Parameters:
231 * =========================
232 *>
233 *> \verbatim
234 *> RELFAC DOUBLE PRECISION, default = 2.0e0
235 *> The relative tolerance. An interval (a,b] lies within
236 *> "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
237 *> where "ulp" is the machine precision (distance from 1 to
238 *> the next larger floating point number.)
239 *>
240 *> FUDGE DOUBLE PRECISION, default = 2
241 *> A "fudge factor" to widen the Gershgorin intervals. Ideally,
242 *> a value of 1 should work, but on machines with sloppy
243 *> arithmetic, this needs to be larger. The default for
244 *> publicly released versions should be large enough to handle
245 *> the worst machine around. Note that this has no effect
246 *> on accuracy of the solution.
247 *> \endverbatim
248 *
249 * Authors:
250 * ========
251 *
252 *> \author Univ. of Tennessee
253 *> \author Univ. of California Berkeley
254 *> \author Univ. of Colorado Denver
255 *> \author NAG Ltd.
256 *
257 *> \date November 2011
258 *
259 *> \ingroup auxOTHERcomputational
260 *
261 * =====================================================================
262  SUBROUTINE dstebz( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
263  $ m, nsplit, w, iblock, isplit, work, iwork,
264  $ info )
265 *
266 * -- LAPACK computational routine (version 3.4.0) --
267 * -- LAPACK is a software package provided by Univ. of Tennessee, --
268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269 * November 2011
270 *
271 * .. Scalar Arguments ..
272  CHARACTER order, range
273  INTEGER il, info, iu, m, n, nsplit
274  DOUBLE PRECISION abstol, vl, vu
275 * ..
276 * .. Array Arguments ..
277  INTEGER iblock( * ), isplit( * ), iwork( * )
278  DOUBLE PRECISION d( * ), e( * ), w( * ), work( * )
279 * ..
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284  DOUBLE PRECISION zero, one, two, half
285  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
286  $ half = 1.0d0 / two )
287  DOUBLE PRECISION fudge, relfac
288  parameter( fudge = 2.1d0, relfac = 2.0d0 )
289 * ..
290 * .. Local Scalars ..
291  LOGICAL ncnvrg, toofew
292  INTEGER ib, ibegin, idiscl, idiscu, ie, iend, iinfo,
293  $ im, in, ioff, iorder, iout, irange, itmax,
294  $ itmp1, iw, iwoff, j, jb, jdisc, je, nb, nwl,
295  $ nwu
296  DOUBLE PRECISION atoli, bnorm, gl, gu, pivmin, rtoli, safemn,
297  $ tmp1, tmp2, tnorm, ulp, wkill, wl, wlu, wu, wul
298 * ..
299 * .. Local Arrays ..
300  INTEGER idumma( 1 )
301 * ..
302 * .. External Functions ..
303  LOGICAL lsame
304  INTEGER ilaenv
305  DOUBLE PRECISION dlamch
306  EXTERNAL lsame, ilaenv, dlamch
307 * ..
308 * .. External Subroutines ..
309  EXTERNAL dlaebz, xerbla
310 * ..
311 * .. Intrinsic Functions ..
312  INTRINSIC abs, int, log, max, min, sqrt
313 * ..
314 * .. Executable Statements ..
315 *
316  info = 0
317 *
318 * Decode RANGE
319 *
320  IF( lsame( range, 'A' ) ) THEN
321  irange = 1
322  ELSE IF( lsame( range, 'V' ) ) THEN
323  irange = 2
324  ELSE IF( lsame( range, 'I' ) ) THEN
325  irange = 3
326  ELSE
327  irange = 0
328  END IF
329 *
330 * Decode ORDER
331 *
332  IF( lsame( order, 'B' ) ) THEN
333  iorder = 2
334  ELSE IF( lsame( order, 'E' ) ) THEN
335  iorder = 1
336  ELSE
337  iorder = 0
338  END IF
339 *
340 * Check for Errors
341 *
342  IF( irange.LE.0 ) THEN
343  info = -1
344  ELSE IF( iorder.LE.0 ) THEN
345  info = -2
346  ELSE IF( n.LT.0 ) THEN
347  info = -3
348  ELSE IF( irange.EQ.2 ) THEN
349  IF( vl.GE.vu )
350  $ info = -5
351  ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1, n ) ) )
352  $ THEN
353  info = -6
354  ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n, il ) .OR. iu.GT.n ) )
355  $ THEN
356  info = -7
357  END IF
358 *
359  IF( info.NE.0 ) THEN
360  CALL xerbla( 'DSTEBZ', -info )
361  return
362  END IF
363 *
364 * Initialize error flags
365 *
366  info = 0
367  ncnvrg = .false.
368  toofew = .false.
369 *
370 * Quick return if possible
371 *
372  m = 0
373  IF( n.EQ.0 )
374  $ return
375 *
376 * Simplifications:
377 *
378  IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n )
379  $ irange = 1
380 *
381 * Get machine constants
382 * NB is the minimum vector length for vector bisection, or 0
383 * if only scalar is to be done.
384 *
385  safemn = dlamch( 'S' )
386  ulp = dlamch( 'P' )
387  rtoli = ulp*relfac
388  nb = ilaenv( 1, 'DSTEBZ', ' ', n, -1, -1, -1 )
389  IF( nb.LE.1 )
390  $ nb = 0
391 *
392 * Special Case when N=1
393 *
394  IF( n.EQ.1 ) THEN
395  nsplit = 1
396  isplit( 1 ) = 1
397  IF( irange.EQ.2 .AND. ( vl.GE.d( 1 ) .OR. vu.LT.d( 1 ) ) ) THEN
398  m = 0
399  ELSE
400  w( 1 ) = d( 1 )
401  iblock( 1 ) = 1
402  m = 1
403  END IF
404  return
405  END IF
406 *
407 * Compute Splitting Points
408 *
409  nsplit = 1
410  work( n ) = zero
411  pivmin = one
412 *
413  DO 10 j = 2, n
414  tmp1 = e( j-1 )**2
415  IF( abs( d( j )*d( j-1 ) )*ulp**2+safemn.GT.tmp1 ) THEN
416  isplit( nsplit ) = j - 1
417  nsplit = nsplit + 1
418  work( j-1 ) = zero
419  ELSE
420  work( j-1 ) = tmp1
421  pivmin = max( pivmin, tmp1 )
422  END IF
423  10 continue
424  isplit( nsplit ) = n
425  pivmin = pivmin*safemn
426 *
427 * Compute Interval and ATOLI
428 *
429  IF( irange.EQ.3 ) THEN
430 *
431 * RANGE='I': Compute the interval containing eigenvalues
432 * IL through IU.
433 *
434 * Compute Gershgorin interval for entire (split) matrix
435 * and use it as the initial interval
436 *
437  gu = d( 1 )
438  gl = d( 1 )
439  tmp1 = zero
440 *
441  DO 20 j = 1, n - 1
442  tmp2 = sqrt( work( j ) )
443  gu = max( gu, d( j )+tmp1+tmp2 )
444  gl = min( gl, d( j )-tmp1-tmp2 )
445  tmp1 = tmp2
446  20 continue
447 *
448  gu = max( gu, d( n )+tmp1 )
449  gl = min( gl, d( n )-tmp1 )
450  tnorm = max( abs( gl ), abs( gu ) )
451  gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
452  gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
453 *
454 * Compute Iteration parameters
455 *
456  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
457  $ log( two ) ) + 2
458  IF( abstol.LE.zero ) THEN
459  atoli = ulp*tnorm
460  ELSE
461  atoli = abstol
462  END IF
463 *
464  work( n+1 ) = gl
465  work( n+2 ) = gl
466  work( n+3 ) = gu
467  work( n+4 ) = gu
468  work( n+5 ) = gl
469  work( n+6 ) = gu
470  iwork( 1 ) = -1
471  iwork( 2 ) = -1
472  iwork( 3 ) = n + 1
473  iwork( 4 ) = n + 1
474  iwork( 5 ) = il - 1
475  iwork( 6 ) = iu
476 *
477  CALL dlaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin, d, e,
478  $ work, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
479  $ iwork, w, iblock, iinfo )
480 *
481  IF( iwork( 6 ).EQ.iu ) THEN
482  wl = work( n+1 )
483  wlu = work( n+3 )
484  nwl = iwork( 1 )
485  wu = work( n+4 )
486  wul = work( n+2 )
487  nwu = iwork( 4 )
488  ELSE
489  wl = work( n+2 )
490  wlu = work( n+4 )
491  nwl = iwork( 2 )
492  wu = work( n+3 )
493  wul = work( n+1 )
494  nwu = iwork( 3 )
495  END IF
496 *
497  IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
498  info = 4
499  return
500  END IF
501  ELSE
502 *
503 * RANGE='A' or 'V' -- Set ATOLI
504 *
505  tnorm = max( abs( d( 1 ) )+abs( e( 1 ) ),
506  $ abs( d( n ) )+abs( e( n-1 ) ) )
507 *
508  DO 30 j = 2, n - 1
509  tnorm = max( tnorm, abs( d( j ) )+abs( e( j-1 ) )+
510  $ abs( e( j ) ) )
511  30 continue
512 *
513  IF( abstol.LE.zero ) THEN
514  atoli = ulp*tnorm
515  ELSE
516  atoli = abstol
517  END IF
518 *
519  IF( irange.EQ.2 ) THEN
520  wl = vl
521  wu = vu
522  ELSE
523  wl = zero
524  wu = zero
525  END IF
526  END IF
527 *
528 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
529 * NWL accumulates the number of eigenvalues .le. WL,
530 * NWU accumulates the number of eigenvalues .le. WU
531 *
532  m = 0
533  iend = 0
534  info = 0
535  nwl = 0
536  nwu = 0
537 *
538  DO 70 jb = 1, nsplit
539  ioff = iend
540  ibegin = ioff + 1
541  iend = isplit( jb )
542  in = iend - ioff
543 *
544  IF( in.EQ.1 ) THEN
545 *
546 * Special Case -- IN=1
547 *
548  IF( irange.EQ.1 .OR. wl.GE.d( ibegin )-pivmin )
549  $ nwl = nwl + 1
550  IF( irange.EQ.1 .OR. wu.GE.d( ibegin )-pivmin )
551  $ nwu = nwu + 1
552  IF( irange.EQ.1 .OR. ( wl.LT.d( ibegin )-pivmin .AND. wu.GE.
553  $ d( ibegin )-pivmin ) ) THEN
554  m = m + 1
555  w( m ) = d( ibegin )
556  iblock( m ) = jb
557  END IF
558  ELSE
559 *
560 * General Case -- IN > 1
561 *
562 * Compute Gershgorin Interval
563 * and use it as the initial interval
564 *
565  gu = d( ibegin )
566  gl = d( ibegin )
567  tmp1 = zero
568 *
569  DO 40 j = ibegin, iend - 1
570  tmp2 = abs( e( j ) )
571  gu = max( gu, d( j )+tmp1+tmp2 )
572  gl = min( gl, d( j )-tmp1-tmp2 )
573  tmp1 = tmp2
574  40 continue
575 *
576  gu = max( gu, d( iend )+tmp1 )
577  gl = min( gl, d( iend )-tmp1 )
578  bnorm = max( abs( gl ), abs( gu ) )
579  gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
580  gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
581 *
582 * Compute ATOLI for the current submatrix
583 *
584  IF( abstol.LE.zero ) THEN
585  atoli = ulp*max( abs( gl ), abs( gu ) )
586  ELSE
587  atoli = abstol
588  END IF
589 *
590  IF( irange.GT.1 ) THEN
591  IF( gu.LT.wl ) THEN
592  nwl = nwl + in
593  nwu = nwu + in
594  go to 70
595  END IF
596  gl = max( gl, wl )
597  gu = min( gu, wu )
598  IF( gl.GE.gu )
599  $ go to 70
600  END IF
601 *
602 * Set Up Initial Interval
603 *
604  work( n+1 ) = gl
605  work( n+in+1 ) = gu
606  CALL dlaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
607  $ d( ibegin ), e( ibegin ), work( ibegin ),
608  $ idumma, work( n+1 ), work( n+2*in+1 ), im,
609  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
610 *
611  nwl = nwl + iwork( 1 )
612  nwu = nwu + iwork( in+1 )
613  iwoff = m - iwork( 1 )
614 *
615 * Compute Eigenvalues
616 *
617  itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
618  $ log( two ) ) + 2
619  CALL dlaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
620  $ d( ibegin ), e( ibegin ), work( ibegin ),
621  $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
622  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
623 *
624 * Copy Eigenvalues Into W and IBLOCK
625 * Use -JB for block number for unconverged eigenvalues.
626 *
627  DO 60 j = 1, iout
628  tmp1 = half*( work( j+n )+work( j+in+n ) )
629 *
630 * Flag non-convergence.
631 *
632  IF( j.GT.iout-iinfo ) THEN
633  ncnvrg = .true.
634  ib = -jb
635  ELSE
636  ib = jb
637  END IF
638  DO 50 je = iwork( j ) + 1 + iwoff,
639  $ iwork( j+in ) + iwoff
640  w( je ) = tmp1
641  iblock( je ) = ib
642  50 continue
643  60 continue
644 *
645  m = m + im
646  END IF
647  70 continue
648 *
649 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
650 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
651 *
652  IF( irange.EQ.3 ) THEN
653  im = 0
654  idiscl = il - 1 - nwl
655  idiscu = nwu - iu
656 *
657  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
658  DO 80 je = 1, m
659  IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
660  idiscl = idiscl - 1
661  ELSE IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
662  idiscu = idiscu - 1
663  ELSE
664  im = im + 1
665  w( im ) = w( je )
666  iblock( im ) = iblock( je )
667  END IF
668  80 continue
669  m = im
670  END IF
671  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
672 *
673 * Code to deal with effects of bad arithmetic:
674 * Some low eigenvalues to be discarded are not in (WL,WLU],
675 * or high eigenvalues to be discarded are not in (WUL,WU]
676 * so just kill off the smallest IDISCL/largest IDISCU
677 * eigenvalues, by simply finding the smallest/largest
678 * eigenvalue(s).
679 *
680 * (If N(w) is monotone non-decreasing, this should never
681 * happen.)
682 *
683  IF( idiscl.GT.0 ) THEN
684  wkill = wu
685  DO 100 jdisc = 1, idiscl
686  iw = 0
687  DO 90 je = 1, m
688  IF( iblock( je ).NE.0 .AND.
689  $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
690  iw = je
691  wkill = w( je )
692  END IF
693  90 continue
694  iblock( iw ) = 0
695  100 continue
696  END IF
697  IF( idiscu.GT.0 ) THEN
698 *
699  wkill = wl
700  DO 120 jdisc = 1, idiscu
701  iw = 0
702  DO 110 je = 1, m
703  IF( iblock( je ).NE.0 .AND.
704  $ ( w( je ).GT.wkill .OR. iw.EQ.0 ) ) THEN
705  iw = je
706  wkill = w( je )
707  END IF
708  110 continue
709  iblock( iw ) = 0
710  120 continue
711  END IF
712  im = 0
713  DO 130 je = 1, m
714  IF( iblock( je ).NE.0 ) THEN
715  im = im + 1
716  w( im ) = w( je )
717  iblock( im ) = iblock( je )
718  END IF
719  130 continue
720  m = im
721  END IF
722  IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
723  toofew = .true.
724  END IF
725  END IF
726 *
727 * If ORDER='B', do nothing -- the eigenvalues are already sorted
728 * by block.
729 * If ORDER='E', sort the eigenvalues from smallest to largest
730 *
731  IF( iorder.EQ.1 .AND. nsplit.GT.1 ) THEN
732  DO 150 je = 1, m - 1
733  ie = 0
734  tmp1 = w( je )
735  DO 140 j = je + 1, m
736  IF( w( j ).LT.tmp1 ) THEN
737  ie = j
738  tmp1 = w( j )
739  END IF
740  140 continue
741 *
742  IF( ie.NE.0 ) THEN
743  itmp1 = iblock( ie )
744  w( ie ) = w( je )
745  iblock( ie ) = iblock( je )
746  w( je ) = tmp1
747  iblock( je ) = itmp1
748  END IF
749  150 continue
750  END IF
751 *
752  info = 0
753  IF( ncnvrg )
754  $ info = info + 1
755  IF( toofew )
756  $ info = info + 2
757  return
758 *
759 * End of DSTEBZ
760 *
761  END