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