LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 *>
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 DOUBLE PRECISION
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 DOUBLE PRECISION
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*DLAMCH('S'), not zero.
138 *> \endverbatim
139 *>
140 *> \param[in] D
141 *> \verbatim
142 *> D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
169 *> On exit, the first M elements of W will contain the
170 *> eigenvalues. (DSTEBZ 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. (DSTEBZ 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 DOUBLE PRECISION 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 DOUBLE PRECISION, 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 DOUBLE PRECISION, 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 dstebz( 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  DOUBLE PRECISION ABSTOL, VL, VU
285 * ..
286 * .. Array Arguments ..
287  INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
288  DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  DOUBLE PRECISION ZERO, ONE, TWO, HALF
295  PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0,
296  $ half = 1.0d0 / two )
297  DOUBLE PRECISION FUDGE, RELFAC
298  PARAMETER ( FUDGE = 2.1d0, relfac = 2.0d0 )
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  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
316  EXTERNAL lsame, ilaenv, dlamch
317 * ..
318 * .. External Subroutines ..
319  EXTERNAL dlaebz, 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 )
360  $ info = -5
361  ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1, n ) ) )
362  $ THEN
363  info = -6
364  ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n, il ) .OR. iu.GT.n ) )
365  $ THEN
366  info = -7
367  END IF
368 *
369  IF( info.NE.0 ) THEN
370  CALL xerbla( 'DSTEBZ', -info )
371  RETURN
372  END IF
373 *
374 * Initialize error flags
375 *
376  info = 0
377  ncnvrg = .false.
378  toofew = .false.
379 *
380 * Quick return if possible
381 *
382  m = 0
383  IF( n.EQ.0 )
384  $ RETURN
385 *
386 * Simplifications:
387 *
388  IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n )
389  $ irange = 1
390 *
391 * Get machine constants
392 * NB is the minimum vector length for vector bisection, or 0
393 * if only scalar is to be done.
394 *
395  safemn = dlamch( 'S' )
396  ulp = dlamch( 'P' )
397  rtoli = ulp*relfac
398  nb = ilaenv( 1, 'DSTEBZ', ' ', n, -1, -1, -1 )
399  IF( nb.LE.1 )
400  $ nb = 0
401 *
402 * Special Case when N=1
403 *
404  IF( n.EQ.1 ) THEN
405  nsplit = 1
406  isplit( 1 ) = 1
407  IF( irange.EQ.2 .AND. ( vl.GE.d( 1 ) .OR. vu.LT.d( 1 ) ) ) THEN
408  m = 0
409  ELSE
410  w( 1 ) = d( 1 )
411  iblock( 1 ) = 1
412  m = 1
413  END IF
414  RETURN
415  END IF
416 *
417 * Compute Splitting Points
418 *
419  nsplit = 1
420  work( n ) = zero
421  pivmin = one
422 *
423  DO 10 j = 2, n
424  tmp1 = e( j-1 )**2
425  IF( abs( d( j )*d( j-1 ) )*ulp**2+safemn.GT.tmp1 ) THEN
426  isplit( nsplit ) = j - 1
427  nsplit = nsplit + 1
428  work( j-1 ) = zero
429  ELSE
430  work( j-1 ) = tmp1
431  pivmin = max( pivmin, tmp1 )
432  END IF
433  10 CONTINUE
434  isplit( nsplit ) = n
435  pivmin = pivmin*safemn
436 *
437 * Compute Interval and ATOLI
438 *
439  IF( irange.EQ.3 ) THEN
440 *
441 * RANGE='I': Compute the interval containing eigenvalues
442 * IL through IU.
443 *
444 * Compute Gershgorin interval for entire (split) matrix
445 * and use it as the initial interval
446 *
447  gu = d( 1 )
448  gl = d( 1 )
449  tmp1 = zero
450 *
451  DO 20 j = 1, n - 1
452  tmp2 = sqrt( work( j ) )
453  gu = max( gu, d( j )+tmp1+tmp2 )
454  gl = min( gl, d( j )-tmp1-tmp2 )
455  tmp1 = tmp2
456  20 CONTINUE
457 *
458  gu = max( gu, d( n )+tmp1 )
459  gl = min( gl, d( n )-tmp1 )
460  tnorm = max( abs( gl ), abs( gu ) )
461  gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
462  gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
463 *
464 * Compute Iteration parameters
465 *
466  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
467  $ log( two ) ) + 2
468  IF( abstol.LE.zero ) THEN
469  atoli = ulp*tnorm
470  ELSE
471  atoli = abstol
472  END IF
473 *
474  work( n+1 ) = gl
475  work( n+2 ) = gl
476  work( n+3 ) = gu
477  work( n+4 ) = gu
478  work( n+5 ) = gl
479  work( n+6 ) = gu
480  iwork( 1 ) = -1
481  iwork( 2 ) = -1
482  iwork( 3 ) = n + 1
483  iwork( 4 ) = n + 1
484  iwork( 5 ) = il - 1
485  iwork( 6 ) = iu
486 *
487  CALL dlaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin, d, e,
488  $ work, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
489  $ iwork, w, iblock, iinfo )
490 *
491  IF( iwork( 6 ).EQ.iu ) THEN
492  wl = work( n+1 )
493  wlu = work( n+3 )
494  nwl = iwork( 1 )
495  wu = work( n+4 )
496  wul = work( n+2 )
497  nwu = iwork( 4 )
498  ELSE
499  wl = work( n+2 )
500  wlu = work( n+4 )
501  nwl = iwork( 2 )
502  wu = work( n+3 )
503  wul = work( n+1 )
504  nwu = iwork( 3 )
505  END IF
506 *
507  IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
508  info = 4
509  RETURN
510  END IF
511  ELSE
512 *
513 * RANGE='A' or 'V' -- Set ATOLI
514 *
515  tnorm = max( abs( d( 1 ) )+abs( e( 1 ) ),
516  $ abs( d( n ) )+abs( e( n-1 ) ) )
517 *
518  DO 30 j = 2, n - 1
519  tnorm = max( tnorm, abs( d( j ) )+abs( e( j-1 ) )+
520  $ abs( e( j ) ) )
521  30 CONTINUE
522 *
523  IF( abstol.LE.zero ) THEN
524  atoli = ulp*tnorm
525  ELSE
526  atoli = abstol
527  END IF
528 *
529  IF( irange.EQ.2 ) THEN
530  wl = vl
531  wu = vu
532  ELSE
533  wl = zero
534  wu = zero
535  END IF
536  END IF
537 *
538 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
539 * NWL accumulates the number of eigenvalues .le. WL,
540 * NWU accumulates the number of eigenvalues .le. WU
541 *
542  m = 0
543  iend = 0
544  info = 0
545  nwl = 0
546  nwu = 0
547 *
548  DO 70 jb = 1, nsplit
549  ioff = iend
550  ibegin = ioff + 1
551  iend = isplit( jb )
552  in = iend - ioff
553 *
554  IF( in.EQ.1 ) THEN
555 *
556 * Special Case -- IN=1
557 *
558  IF( irange.EQ.1 .OR. wl.GE.d( ibegin )-pivmin )
559  $ nwl = nwl + 1
560  IF( irange.EQ.1 .OR. wu.GE.d( ibegin )-pivmin )
561  $ nwu = nwu + 1
562  IF( irange.EQ.1 .OR. ( wl.LT.d( ibegin )-pivmin .AND. wu.GE.
563  $ d( ibegin )-pivmin ) ) THEN
564  m = m + 1
565  w( m ) = d( ibegin )
566  iblock( m ) = jb
567  END IF
568  ELSE
569 *
570 * General Case -- IN > 1
571 *
572 * Compute Gershgorin Interval
573 * and use it as the initial interval
574 *
575  gu = d( ibegin )
576  gl = d( ibegin )
577  tmp1 = zero
578 *
579  DO 40 j = ibegin, iend - 1
580  tmp2 = abs( e( j ) )
581  gu = max( gu, d( j )+tmp1+tmp2 )
582  gl = min( gl, d( j )-tmp1-tmp2 )
583  tmp1 = tmp2
584  40 CONTINUE
585 *
586  gu = max( gu, d( iend )+tmp1 )
587  gl = min( gl, d( iend )-tmp1 )
588  bnorm = max( abs( gl ), abs( gu ) )
589  gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
590  gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
591 *
592 * Compute ATOLI for the current submatrix
593 *
594  IF( abstol.LE.zero ) THEN
595  atoli = ulp*max( abs( gl ), abs( gu ) )
596  ELSE
597  atoli = abstol
598  END IF
599 *
600  IF( irange.GT.1 ) THEN
601  IF( gu.LT.wl ) THEN
602  nwl = nwl + in
603  nwu = nwu + in
604  GO TO 70
605  END IF
606  gl = max( gl, wl )
607  gu = min( gu, wu )
608  IF( gl.GE.gu )
609  $ GO TO 70
610  END IF
611 *
612 * Set Up Initial Interval
613 *
614  work( n+1 ) = gl
615  work( n+in+1 ) = gu
616  CALL dlaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
617  $ d( ibegin ), e( ibegin ), work( ibegin ),
618  $ idumma, work( n+1 ), work( n+2*in+1 ), im,
619  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
620 *
621  nwl = nwl + iwork( 1 )
622  nwu = nwu + iwork( in+1 )
623  iwoff = m - iwork( 1 )
624 *
625 * Compute Eigenvalues
626 *
627  itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
628  $ log( two ) ) + 2
629  CALL dlaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
630  $ d( ibegin ), e( ibegin ), work( ibegin ),
631  $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
632  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
633 *
634 * Copy Eigenvalues Into W and IBLOCK
635 * Use -JB for block number for unconverged eigenvalues.
636 *
637  DO 60 j = 1, iout
638  tmp1 = half*( work( j+n )+work( j+in+n ) )
639 *
640 * Flag non-convergence.
641 *
642  IF( j.GT.iout-iinfo ) THEN
643  ncnvrg = .true.
644  ib = -jb
645  ELSE
646  ib = jb
647  END IF
648  DO 50 je = iwork( j ) + 1 + iwoff,
649  $ iwork( j+in ) + iwoff
650  w( je ) = tmp1
651  iblock( je ) = ib
652  50 CONTINUE
653  60 CONTINUE
654 *
655  m = m + im
656  END IF
657  70 CONTINUE
658 *
659 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
660 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
661 *
662  IF( irange.EQ.3 ) THEN
663  im = 0
664  idiscl = il - 1 - nwl
665  idiscu = nwu - iu
666 *
667  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
668  DO 80 je = 1, m
669  IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
670  idiscl = idiscl - 1
671  ELSE IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
672  idiscu = idiscu - 1
673  ELSE
674  im = im + 1
675  w( im ) = w( je )
676  iblock( im ) = iblock( je )
677  END IF
678  80 CONTINUE
679  m = im
680  END IF
681  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
682 *
683 * Code to deal with effects of bad arithmetic:
684 * Some low eigenvalues to be discarded are not in (WL,WLU],
685 * or high eigenvalues to be discarded are not in (WUL,WU]
686 * so just kill off the smallest IDISCL/largest IDISCU
687 * eigenvalues, by simply finding the smallest/largest
688 * eigenvalue(s).
689 *
690 * (If N(w) is monotone non-decreasing, this should never
691 * happen.)
692 *
693  IF( idiscl.GT.0 ) THEN
694  wkill = wu
695  DO 100 jdisc = 1, idiscl
696  iw = 0
697  DO 90 je = 1, m
698  IF( iblock( je ).NE.0 .AND.
699  $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
700  iw = je
701  wkill = w( je )
702  END IF
703  90 CONTINUE
704  iblock( iw ) = 0
705  100 CONTINUE
706  END IF
707  IF( idiscu.GT.0 ) THEN
708 *
709  wkill = wl
710  DO 120 jdisc = 1, idiscu
711  iw = 0
712  DO 110 je = 1, m
713  IF( iblock( je ).NE.0 .AND.
714  $ ( w( je ).GT.wkill .OR. iw.EQ.0 ) ) THEN
715  iw = je
716  wkill = w( je )
717  END IF
718  110 CONTINUE
719  iblock( iw ) = 0
720  120 CONTINUE
721  END IF
722  im = 0
723  DO 130 je = 1, m
724  IF( iblock( je ).NE.0 ) THEN
725  im = im + 1
726  w( im ) = w( je )
727  iblock( im ) = iblock( je )
728  END IF
729  130 CONTINUE
730  m = im
731  END IF
732  IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
733  toofew = .true.
734  END IF
735  END IF
736 *
737 * If ORDER='B', do nothing -- the eigenvalues are already sorted
738 * by block.
739 * If ORDER='E', sort the eigenvalues from smallest to largest
740 *
741  IF( iorder.EQ.1 .AND. nsplit.GT.1 ) THEN
742  DO 150 je = 1, m - 1
743  ie = 0
744  tmp1 = w( je )
745  DO 140 j = je + 1, m
746  IF( w( j ).LT.tmp1 ) THEN
747  ie = j
748  tmp1 = w( j )
749  END IF
750  140 CONTINUE
751 *
752  IF( ie.NE.0 ) THEN
753  itmp1 = iblock( ie )
754  w( ie ) = w( je )
755  iblock( ie ) = iblock( je )
756  w( je ) = tmp1
757  iblock( je ) = itmp1
758  END IF
759  150 CONTINUE
760  END IF
761 *
762  info = 0
763  IF( ncnvrg )
764  $ info = info + 1
765  IF( toofew )
766  $ info = info + 2
767  RETURN
768 *
769 * End of DSTEBZ
770 *
771  END
dstebz
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:275
dlaebz
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:321
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62