1 SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
2 $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
3 $ INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 * 8-18-00: Increase FUDGE factor for T3E (eca)
10 *
11 * .. Scalar Arguments ..
12 CHARACTER ORDER, RANGE
13 INTEGER IL, INFO, IU, M, N, NSPLIT
14 DOUBLE PRECISION ABSTOL, VL, VU
15 * ..
16 * .. Array Arguments ..
17 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
18 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSTEBZ computes the eigenvalues of a symmetric tridiagonal
25 * matrix T. The user may ask for all eigenvalues, all eigenvalues
26 * in the half-open interval (VL, VU], or the IL-th through IU-th
27 * eigenvalues.
28 *
29 * To avoid overflow, the matrix must be scaled so that its
30 * largest element is no greater than overflow**(1/2) *
31 * underflow**(1/4) in absolute value, and for greatest
32 * accuracy, it should not be much smaller than that.
33 *
34 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
35 * Matrix", Report CS41, Computer Science Dept., Stanford
36 * University, July 21, 1966.
37 *
38 * Arguments
39 * =========
40 *
41 * RANGE (input) CHARACTER*1
42 * = 'A': ("All") all eigenvalues will be found.
43 * = 'V': ("Value") all eigenvalues in the half-open interval
44 * (VL, VU] will be found.
45 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
46 * entire matrix) will be found.
47 *
48 * ORDER (input) CHARACTER*1
49 * = 'B': ("By Block") the eigenvalues will be grouped by
50 * split-off block (see IBLOCK, ISPLIT) and
51 * ordered from smallest to largest within
52 * the block.
53 * = 'E': ("Entire matrix")
54 * the eigenvalues for the entire matrix
55 * will be ordered from smallest to
56 * largest.
57 *
58 * N (input) INTEGER
59 * The order of the tridiagonal matrix T. N >= 0.
60 *
61 * VL (input) DOUBLE PRECISION
62 * VU (input) DOUBLE PRECISION
63 * If RANGE='V', the lower and upper bounds of the interval to
64 * be searched for eigenvalues. Eigenvalues less than or equal
65 * to VL, or greater than VU, will not be returned. VL < VU.
66 * Not referenced if RANGE = 'A' or 'I'.
67 *
68 * IL (input) INTEGER
69 * IU (input) INTEGER
70 * If RANGE='I', the indices (in ascending order) of the
71 * smallest and largest eigenvalues to be returned.
72 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
73 * Not referenced if RANGE = 'A' or 'V'.
74 *
75 * ABSTOL (input) DOUBLE PRECISION
76 * The absolute tolerance for the eigenvalues. An eigenvalue
77 * (or cluster) is considered to be located if it has been
78 * determined to lie in an interval whose width is ABSTOL or
79 * less. If ABSTOL is less than or equal to zero, then ULP*|T|
80 * will be used, where |T| means the 1-norm of T.
81 *
82 * Eigenvalues will be computed most accurately when ABSTOL is
83 * set to twice the underflow threshold 2*DLAMCH('S'), not zero.
84 *
85 * D (input) DOUBLE PRECISION array, dimension (N)
86 * The n diagonal elements of the tridiagonal matrix T.
87 *
88 * E (input) DOUBLE PRECISION array, dimension (N-1)
89 * The (n-1) off-diagonal elements of the tridiagonal matrix T.
90 *
91 * M (output) INTEGER
92 * The actual number of eigenvalues found. 0 <= M <= N.
93 * (See also the description of INFO=2,3.)
94 *
95 * NSPLIT (output) INTEGER
96 * The number of diagonal blocks in the matrix T.
97 * 1 <= NSPLIT <= N.
98 *
99 * W (output) DOUBLE PRECISION array, dimension (N)
100 * On exit, the first M elements of W will contain the
101 * eigenvalues. (DSTEBZ may use the remaining N-M elements as
102 * workspace.)
103 *
104 * IBLOCK (output) INTEGER array, dimension (N)
105 * At each row/column j where E(j) is zero or small, the
106 * matrix T is considered to split into a block diagonal
107 * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
108 * block (from 1 to the number of blocks) the eigenvalue W(i)
109 * belongs. (DSTEBZ may use the remaining N-M elements as
110 * workspace.)
111 *
112 * ISPLIT (output) INTEGER array, dimension (N)
113 * The splitting points, at which T breaks up into submatrices.
114 * The first submatrix consists of rows/columns 1 to ISPLIT(1),
115 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
116 * etc., and the NSPLIT-th consists of rows/columns
117 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
118 * (Only the first NSPLIT elements will actually be used, but
119 * since the user cannot know a priori what value NSPLIT will
120 * have, N words must be reserved for ISPLIT.)
121 *
122 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
123 *
124 * IWORK (workspace) INTEGER array, dimension (3*N)
125 *
126 * INFO (output) INTEGER
127 * = 0: successful exit
128 * < 0: if INFO = -i, the i-th argument had an illegal value
129 * > 0: some or all of the eigenvalues failed to converge or
130 * were not computed:
131 * =1 or 3: Bisection failed to converge for some
132 * eigenvalues; these eigenvalues are flagged by a
133 * negative block number. The effect is that the
134 * eigenvalues may not be as accurate as the
135 * absolute and relative tolerances. This is
136 * generally caused by unexpectedly inaccurate
137 * arithmetic.
138 * =2 or 3: RANGE='I' only: Not all of the eigenvalues
139 * IL:IU were found.
140 * Effect: M < IU+1-IL
141 * Cause: non-monotonic arithmetic, causing the
142 * Sturm sequence to be non-monotonic.
143 * Cure: recalculate, using RANGE='A', and pick
144 * out eigenvalues IL:IU. In some cases,
145 * increasing the PARAMETER "FUDGE" may
146 * make things work.
147 * = 4: RANGE='I', and the Gershgorin interval
148 * initially used was too small. No eigenvalues
149 * were computed.
150 * Probable cause: your machine has sloppy
151 * floating-point arithmetic.
152 * Cure: Increase the PARAMETER "FUDGE",
153 * recompile, and try again.
154 *
155 * Internal Parameters
156 * ===================
157 *
158 * RELFAC DOUBLE PRECISION, default = 2.0e0
159 * The relative tolerance. An interval (a,b] lies within
160 * "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
161 * where "ulp" is the machine precision (distance from 1 to
162 * the next larger floating point number.)
163 *
164 * FUDGE DOUBLE PRECISION, default = 2
165 * A "fudge factor" to widen the Gershgorin intervals. Ideally,
166 * a value of 1 should work, but on machines with sloppy
167 * arithmetic, this needs to be larger. The default for
168 * publicly released versions should be large enough to handle
169 * the worst machine around. Note that this has no effect
170 * on accuracy of the solution.
171 *
172 * =====================================================================
173 *
174 * .. Parameters ..
175 DOUBLE PRECISION ZERO, ONE, TWO, HALF
176 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
177 $ HALF = 1.0D0 / TWO )
178 DOUBLE PRECISION FUDGE, RELFAC
179 PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 )
180 * ..
181 * .. Local Scalars ..
182 LOGICAL NCNVRG, TOOFEW
183 INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
184 $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
185 $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
186 $ NWU
187 DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
188 $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
189 * ..
190 * .. Local Arrays ..
191 INTEGER IDUMMA( 1 )
192 * ..
193 * .. External Functions ..
194 LOGICAL LSAME
195 INTEGER ILAENV
196 DOUBLE PRECISION DLAMCH
197 EXTERNAL LSAME, ILAENV, DLAMCH
198 * ..
199 * .. External Subroutines ..
200 EXTERNAL DLAEBZ, XERBLA
201 * ..
202 * .. Intrinsic Functions ..
203 INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
204 * ..
205 * .. Executable Statements ..
206 *
207 INFO = 0
208 *
209 * Decode RANGE
210 *
211 IF( LSAME( RANGE, 'A' ) ) THEN
212 IRANGE = 1
213 ELSE IF( LSAME( RANGE, 'V' ) ) THEN
214 IRANGE = 2
215 ELSE IF( LSAME( RANGE, 'I' ) ) THEN
216 IRANGE = 3
217 ELSE
218 IRANGE = 0
219 END IF
220 *
221 * Decode ORDER
222 *
223 IF( LSAME( ORDER, 'B' ) ) THEN
224 IORDER = 2
225 ELSE IF( LSAME( ORDER, 'E' ) ) THEN
226 IORDER = 1
227 ELSE
228 IORDER = 0
229 END IF
230 *
231 * Check for Errors
232 *
233 IF( IRANGE.LE.0 ) THEN
234 INFO = -1
235 ELSE IF( IORDER.LE.0 ) THEN
236 INFO = -2
237 ELSE IF( N.LT.0 ) THEN
238 INFO = -3
239 ELSE IF( IRANGE.EQ.2 ) THEN
240 IF( VL.GE.VU )
241 $ INFO = -5
242 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
243 $ THEN
244 INFO = -6
245 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
246 $ THEN
247 INFO = -7
248 END IF
249 *
250 IF( INFO.NE.0 ) THEN
251 CALL XERBLA( 'DSTEBZ', -INFO )
252 RETURN
253 END IF
254 *
255 * Initialize error flags
256 *
257 INFO = 0
258 NCNVRG = .FALSE.
259 TOOFEW = .FALSE.
260 *
261 * Quick return if possible
262 *
263 M = 0
264 IF( N.EQ.0 )
265 $ RETURN
266 *
267 * Simplifications:
268 *
269 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
270 $ IRANGE = 1
271 *
272 * Get machine constants
273 * NB is the minimum vector length for vector bisection, or 0
274 * if only scalar is to be done.
275 *
276 SAFEMN = DLAMCH( 'S' )
277 ULP = DLAMCH( 'P' )
278 RTOLI = ULP*RELFAC
279 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
280 IF( NB.LE.1 )
281 $ NB = 0
282 *
283 * Special Case when N=1
284 *
285 IF( N.EQ.1 ) THEN
286 NSPLIT = 1
287 ISPLIT( 1 ) = 1
288 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
289 M = 0
290 ELSE
291 W( 1 ) = D( 1 )
292 IBLOCK( 1 ) = 1
293 M = 1
294 END IF
295 RETURN
296 END IF
297 *
298 * Compute Splitting Points
299 *
300 NSPLIT = 1
301 WORK( N ) = ZERO
302 PIVMIN = ONE
303 *
304 DO 10 J = 2, N
305 TMP1 = E( J-1 )**2
306 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
307 ISPLIT( NSPLIT ) = J - 1
308 NSPLIT = NSPLIT + 1
309 WORK( J-1 ) = ZERO
310 ELSE
311 WORK( J-1 ) = TMP1
312 PIVMIN = MAX( PIVMIN, TMP1 )
313 END IF
314 10 CONTINUE
315 ISPLIT( NSPLIT ) = N
316 PIVMIN = PIVMIN*SAFEMN
317 *
318 * Compute Interval and ATOLI
319 *
320 IF( IRANGE.EQ.3 ) THEN
321 *
322 * RANGE='I': Compute the interval containing eigenvalues
323 * IL through IU.
324 *
325 * Compute Gershgorin interval for entire (split) matrix
326 * and use it as the initial interval
327 *
328 GU = D( 1 )
329 GL = D( 1 )
330 TMP1 = ZERO
331 *
332 DO 20 J = 1, N - 1
333 TMP2 = SQRT( WORK( J ) )
334 GU = MAX( GU, D( J )+TMP1+TMP2 )
335 GL = MIN( GL, D( J )-TMP1-TMP2 )
336 TMP1 = TMP2
337 20 CONTINUE
338 *
339 GU = MAX( GU, D( N )+TMP1 )
340 GL = MIN( GL, D( N )-TMP1 )
341 TNORM = MAX( ABS( GL ), ABS( GU ) )
342 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
343 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
344 *
345 * Compute Iteration parameters
346 *
347 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
348 $ LOG( TWO ) ) + 2
349 IF( ABSTOL.LE.ZERO ) THEN
350 ATOLI = ULP*TNORM
351 ELSE
352 ATOLI = ABSTOL
353 END IF
354 *
355 WORK( N+1 ) = GL
356 WORK( N+2 ) = GL
357 WORK( N+3 ) = GU
358 WORK( N+4 ) = GU
359 WORK( N+5 ) = GL
360 WORK( N+6 ) = GU
361 IWORK( 1 ) = -1
362 IWORK( 2 ) = -1
363 IWORK( 3 ) = N + 1
364 IWORK( 4 ) = N + 1
365 IWORK( 5 ) = IL - 1
366 IWORK( 6 ) = IU
367 *
368 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
369 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
370 $ IWORK, W, IBLOCK, IINFO )
371 *
372 IF( IWORK( 6 ).EQ.IU ) THEN
373 WL = WORK( N+1 )
374 WLU = WORK( N+3 )
375 NWL = IWORK( 1 )
376 WU = WORK( N+4 )
377 WUL = WORK( N+2 )
378 NWU = IWORK( 4 )
379 ELSE
380 WL = WORK( N+2 )
381 WLU = WORK( N+4 )
382 NWL = IWORK( 2 )
383 WU = WORK( N+3 )
384 WUL = WORK( N+1 )
385 NWU = IWORK( 3 )
386 END IF
387 *
388 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
389 INFO = 4
390 RETURN
391 END IF
392 ELSE
393 *
394 * RANGE='A' or 'V' -- Set ATOLI
395 *
396 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
397 $ ABS( D( N ) )+ABS( E( N-1 ) ) )
398 *
399 DO 30 J = 2, N - 1
400 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
401 $ ABS( E( J ) ) )
402 30 CONTINUE
403 *
404 IF( ABSTOL.LE.ZERO ) THEN
405 ATOLI = ULP*TNORM
406 ELSE
407 ATOLI = ABSTOL
408 END IF
409 *
410 IF( IRANGE.EQ.2 ) THEN
411 WL = VL
412 WU = VU
413 ELSE
414 WL = ZERO
415 WU = ZERO
416 END IF
417 END IF
418 *
419 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
420 * NWL accumulates the number of eigenvalues .le. WL,
421 * NWU accumulates the number of eigenvalues .le. WU
422 *
423 M = 0
424 IEND = 0
425 INFO = 0
426 NWL = 0
427 NWU = 0
428 *
429 DO 70 JB = 1, NSPLIT
430 IOFF = IEND
431 IBEGIN = IOFF + 1
432 IEND = ISPLIT( JB )
433 IN = IEND - IOFF
434 *
435 IF( IN.EQ.1 ) THEN
436 *
437 * Special Case -- IN=1
438 *
439 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
440 $ NWL = NWL + 1
441 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
442 $ NWU = NWU + 1
443 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
444 $ D( IBEGIN )-PIVMIN ) ) THEN
445 M = M + 1
446 W( M ) = D( IBEGIN )
447 IBLOCK( M ) = JB
448 END IF
449 ELSE
450 *
451 * General Case -- IN > 1
452 *
453 * Compute Gershgorin Interval
454 * and use it as the initial interval
455 *
456 GU = D( IBEGIN )
457 GL = D( IBEGIN )
458 TMP1 = ZERO
459 *
460 DO 40 J = IBEGIN, IEND - 1
461 TMP2 = ABS( E( J ) )
462 GU = MAX( GU, D( J )+TMP1+TMP2 )
463 GL = MIN( GL, D( J )-TMP1-TMP2 )
464 TMP1 = TMP2
465 40 CONTINUE
466 *
467 GU = MAX( GU, D( IEND )+TMP1 )
468 GL = MIN( GL, D( IEND )-TMP1 )
469 BNORM = MAX( ABS( GL ), ABS( GU ) )
470 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
471 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
472 *
473 * Compute ATOLI for the current submatrix
474 *
475 IF( ABSTOL.LE.ZERO ) THEN
476 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
477 ELSE
478 ATOLI = ABSTOL
479 END IF
480 *
481 IF( IRANGE.GT.1 ) THEN
482 IF( GU.LT.WL ) THEN
483 NWL = NWL + IN
484 NWU = NWU + IN
485 GO TO 70
486 END IF
487 GL = MAX( GL, WL )
488 GU = MIN( GU, WU )
489 IF( GL.GE.GU )
490 $ GO TO 70
491 END IF
492 *
493 * Set Up Initial Interval
494 *
495 WORK( N+1 ) = GL
496 WORK( N+IN+1 ) = GU
497 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
498 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
499 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
500 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
501 *
502 NWL = NWL + IWORK( 1 )
503 NWU = NWU + IWORK( IN+1 )
504 IWOFF = M - IWORK( 1 )
505 *
506 * Compute Eigenvalues
507 *
508 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
509 $ LOG( TWO ) ) + 2
510 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
511 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
512 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
513 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
514 *
515 * Copy Eigenvalues Into W and IBLOCK
516 * Use -JB for block number for unconverged eigenvalues.
517 *
518 DO 60 J = 1, IOUT
519 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
520 *
521 * Flag non-convergence.
522 *
523 IF( J.GT.IOUT-IINFO ) THEN
524 NCNVRG = .TRUE.
525 IB = -JB
526 ELSE
527 IB = JB
528 END IF
529 DO 50 JE = IWORK( J ) + 1 + IWOFF,
530 $ IWORK( J+IN ) + IWOFF
531 W( JE ) = TMP1
532 IBLOCK( JE ) = IB
533 50 CONTINUE
534 60 CONTINUE
535 *
536 M = M + IM
537 END IF
538 70 CONTINUE
539 *
540 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
541 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
542 *
543 IF( IRANGE.EQ.3 ) THEN
544 IM = 0
545 IDISCL = IL - 1 - NWL
546 IDISCU = NWU - IU
547 *
548 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
549 DO 80 JE = 1, M
550 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
551 IDISCL = IDISCL - 1
552 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
553 IDISCU = IDISCU - 1
554 ELSE
555 IM = IM + 1
556 W( IM ) = W( JE )
557 IBLOCK( IM ) = IBLOCK( JE )
558 END IF
559 80 CONTINUE
560 M = IM
561 END IF
562 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
563 *
564 * Code to deal with effects of bad arithmetic:
565 * Some low eigenvalues to be discarded are not in (WL,WLU],
566 * or high eigenvalues to be discarded are not in (WUL,WU]
567 * so just kill off the smallest IDISCL/largest IDISCU
568 * eigenvalues, by simply finding the smallest/largest
569 * eigenvalue(s).
570 *
571 * (If N(w) is monotone non-decreasing, this should never
572 * happen.)
573 *
574 IF( IDISCL.GT.0 ) THEN
575 WKILL = WU
576 DO 100 JDISC = 1, IDISCL
577 IW = 0
578 DO 90 JE = 1, M
579 IF( IBLOCK( JE ).NE.0 .AND.
580 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
581 IW = JE
582 WKILL = W( JE )
583 END IF
584 90 CONTINUE
585 IBLOCK( IW ) = 0
586 100 CONTINUE
587 END IF
588 IF( IDISCU.GT.0 ) THEN
589 *
590 WKILL = WL
591 DO 120 JDISC = 1, IDISCU
592 IW = 0
593 DO 110 JE = 1, M
594 IF( IBLOCK( JE ).NE.0 .AND.
595 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
596 IW = JE
597 WKILL = W( JE )
598 END IF
599 110 CONTINUE
600 IBLOCK( IW ) = 0
601 120 CONTINUE
602 END IF
603 IM = 0
604 DO 130 JE = 1, M
605 IF( IBLOCK( JE ).NE.0 ) THEN
606 IM = IM + 1
607 W( IM ) = W( JE )
608 IBLOCK( IM ) = IBLOCK( JE )
609 END IF
610 130 CONTINUE
611 M = IM
612 END IF
613 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
614 TOOFEW = .TRUE.
615 END IF
616 END IF
617 *
618 * If ORDER='B', do nothing -- the eigenvalues are already sorted
619 * by block.
620 * If ORDER='E', sort the eigenvalues from smallest to largest
621 *
622 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
623 DO 150 JE = 1, M - 1
624 IE = 0
625 TMP1 = W( JE )
626 DO 140 J = JE + 1, M
627 IF( W( J ).LT.TMP1 ) THEN
628 IE = J
629 TMP1 = W( J )
630 END IF
631 140 CONTINUE
632 *
633 IF( IE.NE.0 ) THEN
634 ITMP1 = IBLOCK( IE )
635 W( IE ) = W( JE )
636 IBLOCK( IE ) = IBLOCK( JE )
637 W( JE ) = TMP1
638 IBLOCK( JE ) = ITMP1
639 END IF
640 150 CONTINUE
641 END IF
642 *
643 INFO = 0
644 IF( NCNVRG )
645 $ INFO = INFO + 1
646 IF( TOOFEW )
647 $ INFO = INFO + 2
648 RETURN
649 *
650 * End of DSTEBZ
651 *
652 END
2 $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
3 $ INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 * 8-18-00: Increase FUDGE factor for T3E (eca)
10 *
11 * .. Scalar Arguments ..
12 CHARACTER ORDER, RANGE
13 INTEGER IL, INFO, IU, M, N, NSPLIT
14 DOUBLE PRECISION ABSTOL, VL, VU
15 * ..
16 * .. Array Arguments ..
17 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
18 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSTEBZ computes the eigenvalues of a symmetric tridiagonal
25 * matrix T. The user may ask for all eigenvalues, all eigenvalues
26 * in the half-open interval (VL, VU], or the IL-th through IU-th
27 * eigenvalues.
28 *
29 * To avoid overflow, the matrix must be scaled so that its
30 * largest element is no greater than overflow**(1/2) *
31 * underflow**(1/4) in absolute value, and for greatest
32 * accuracy, it should not be much smaller than that.
33 *
34 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
35 * Matrix", Report CS41, Computer Science Dept., Stanford
36 * University, July 21, 1966.
37 *
38 * Arguments
39 * =========
40 *
41 * RANGE (input) CHARACTER*1
42 * = 'A': ("All") all eigenvalues will be found.
43 * = 'V': ("Value") all eigenvalues in the half-open interval
44 * (VL, VU] will be found.
45 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
46 * entire matrix) will be found.
47 *
48 * ORDER (input) CHARACTER*1
49 * = 'B': ("By Block") the eigenvalues will be grouped by
50 * split-off block (see IBLOCK, ISPLIT) and
51 * ordered from smallest to largest within
52 * the block.
53 * = 'E': ("Entire matrix")
54 * the eigenvalues for the entire matrix
55 * will be ordered from smallest to
56 * largest.
57 *
58 * N (input) INTEGER
59 * The order of the tridiagonal matrix T. N >= 0.
60 *
61 * VL (input) DOUBLE PRECISION
62 * VU (input) DOUBLE PRECISION
63 * If RANGE='V', the lower and upper bounds of the interval to
64 * be searched for eigenvalues. Eigenvalues less than or equal
65 * to VL, or greater than VU, will not be returned. VL < VU.
66 * Not referenced if RANGE = 'A' or 'I'.
67 *
68 * IL (input) INTEGER
69 * IU (input) INTEGER
70 * If RANGE='I', the indices (in ascending order) of the
71 * smallest and largest eigenvalues to be returned.
72 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
73 * Not referenced if RANGE = 'A' or 'V'.
74 *
75 * ABSTOL (input) DOUBLE PRECISION
76 * The absolute tolerance for the eigenvalues. An eigenvalue
77 * (or cluster) is considered to be located if it has been
78 * determined to lie in an interval whose width is ABSTOL or
79 * less. If ABSTOL is less than or equal to zero, then ULP*|T|
80 * will be used, where |T| means the 1-norm of T.
81 *
82 * Eigenvalues will be computed most accurately when ABSTOL is
83 * set to twice the underflow threshold 2*DLAMCH('S'), not zero.
84 *
85 * D (input) DOUBLE PRECISION array, dimension (N)
86 * The n diagonal elements of the tridiagonal matrix T.
87 *
88 * E (input) DOUBLE PRECISION array, dimension (N-1)
89 * The (n-1) off-diagonal elements of the tridiagonal matrix T.
90 *
91 * M (output) INTEGER
92 * The actual number of eigenvalues found. 0 <= M <= N.
93 * (See also the description of INFO=2,3.)
94 *
95 * NSPLIT (output) INTEGER
96 * The number of diagonal blocks in the matrix T.
97 * 1 <= NSPLIT <= N.
98 *
99 * W (output) DOUBLE PRECISION array, dimension (N)
100 * On exit, the first M elements of W will contain the
101 * eigenvalues. (DSTEBZ may use the remaining N-M elements as
102 * workspace.)
103 *
104 * IBLOCK (output) INTEGER array, dimension (N)
105 * At each row/column j where E(j) is zero or small, the
106 * matrix T is considered to split into a block diagonal
107 * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
108 * block (from 1 to the number of blocks) the eigenvalue W(i)
109 * belongs. (DSTEBZ may use the remaining N-M elements as
110 * workspace.)
111 *
112 * ISPLIT (output) INTEGER array, dimension (N)
113 * The splitting points, at which T breaks up into submatrices.
114 * The first submatrix consists of rows/columns 1 to ISPLIT(1),
115 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
116 * etc., and the NSPLIT-th consists of rows/columns
117 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
118 * (Only the first NSPLIT elements will actually be used, but
119 * since the user cannot know a priori what value NSPLIT will
120 * have, N words must be reserved for ISPLIT.)
121 *
122 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
123 *
124 * IWORK (workspace) INTEGER array, dimension (3*N)
125 *
126 * INFO (output) INTEGER
127 * = 0: successful exit
128 * < 0: if INFO = -i, the i-th argument had an illegal value
129 * > 0: some or all of the eigenvalues failed to converge or
130 * were not computed:
131 * =1 or 3: Bisection failed to converge for some
132 * eigenvalues; these eigenvalues are flagged by a
133 * negative block number. The effect is that the
134 * eigenvalues may not be as accurate as the
135 * absolute and relative tolerances. This is
136 * generally caused by unexpectedly inaccurate
137 * arithmetic.
138 * =2 or 3: RANGE='I' only: Not all of the eigenvalues
139 * IL:IU were found.
140 * Effect: M < IU+1-IL
141 * Cause: non-monotonic arithmetic, causing the
142 * Sturm sequence to be non-monotonic.
143 * Cure: recalculate, using RANGE='A', and pick
144 * out eigenvalues IL:IU. In some cases,
145 * increasing the PARAMETER "FUDGE" may
146 * make things work.
147 * = 4: RANGE='I', and the Gershgorin interval
148 * initially used was too small. No eigenvalues
149 * were computed.
150 * Probable cause: your machine has sloppy
151 * floating-point arithmetic.
152 * Cure: Increase the PARAMETER "FUDGE",
153 * recompile, and try again.
154 *
155 * Internal Parameters
156 * ===================
157 *
158 * RELFAC DOUBLE PRECISION, default = 2.0e0
159 * The relative tolerance. An interval (a,b] lies within
160 * "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
161 * where "ulp" is the machine precision (distance from 1 to
162 * the next larger floating point number.)
163 *
164 * FUDGE DOUBLE PRECISION, default = 2
165 * A "fudge factor" to widen the Gershgorin intervals. Ideally,
166 * a value of 1 should work, but on machines with sloppy
167 * arithmetic, this needs to be larger. The default for
168 * publicly released versions should be large enough to handle
169 * the worst machine around. Note that this has no effect
170 * on accuracy of the solution.
171 *
172 * =====================================================================
173 *
174 * .. Parameters ..
175 DOUBLE PRECISION ZERO, ONE, TWO, HALF
176 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
177 $ HALF = 1.0D0 / TWO )
178 DOUBLE PRECISION FUDGE, RELFAC
179 PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 )
180 * ..
181 * .. Local Scalars ..
182 LOGICAL NCNVRG, TOOFEW
183 INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
184 $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
185 $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
186 $ NWU
187 DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
188 $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
189 * ..
190 * .. Local Arrays ..
191 INTEGER IDUMMA( 1 )
192 * ..
193 * .. External Functions ..
194 LOGICAL LSAME
195 INTEGER ILAENV
196 DOUBLE PRECISION DLAMCH
197 EXTERNAL LSAME, ILAENV, DLAMCH
198 * ..
199 * .. External Subroutines ..
200 EXTERNAL DLAEBZ, XERBLA
201 * ..
202 * .. Intrinsic Functions ..
203 INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT
204 * ..
205 * .. Executable Statements ..
206 *
207 INFO = 0
208 *
209 * Decode RANGE
210 *
211 IF( LSAME( RANGE, 'A' ) ) THEN
212 IRANGE = 1
213 ELSE IF( LSAME( RANGE, 'V' ) ) THEN
214 IRANGE = 2
215 ELSE IF( LSAME( RANGE, 'I' ) ) THEN
216 IRANGE = 3
217 ELSE
218 IRANGE = 0
219 END IF
220 *
221 * Decode ORDER
222 *
223 IF( LSAME( ORDER, 'B' ) ) THEN
224 IORDER = 2
225 ELSE IF( LSAME( ORDER, 'E' ) ) THEN
226 IORDER = 1
227 ELSE
228 IORDER = 0
229 END IF
230 *
231 * Check for Errors
232 *
233 IF( IRANGE.LE.0 ) THEN
234 INFO = -1
235 ELSE IF( IORDER.LE.0 ) THEN
236 INFO = -2
237 ELSE IF( N.LT.0 ) THEN
238 INFO = -3
239 ELSE IF( IRANGE.EQ.2 ) THEN
240 IF( VL.GE.VU )
241 $ INFO = -5
242 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
243 $ THEN
244 INFO = -6
245 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
246 $ THEN
247 INFO = -7
248 END IF
249 *
250 IF( INFO.NE.0 ) THEN
251 CALL XERBLA( 'DSTEBZ', -INFO )
252 RETURN
253 END IF
254 *
255 * Initialize error flags
256 *
257 INFO = 0
258 NCNVRG = .FALSE.
259 TOOFEW = .FALSE.
260 *
261 * Quick return if possible
262 *
263 M = 0
264 IF( N.EQ.0 )
265 $ RETURN
266 *
267 * Simplifications:
268 *
269 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
270 $ IRANGE = 1
271 *
272 * Get machine constants
273 * NB is the minimum vector length for vector bisection, or 0
274 * if only scalar is to be done.
275 *
276 SAFEMN = DLAMCH( 'S' )
277 ULP = DLAMCH( 'P' )
278 RTOLI = ULP*RELFAC
279 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 )
280 IF( NB.LE.1 )
281 $ NB = 0
282 *
283 * Special Case when N=1
284 *
285 IF( N.EQ.1 ) THEN
286 NSPLIT = 1
287 ISPLIT( 1 ) = 1
288 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
289 M = 0
290 ELSE
291 W( 1 ) = D( 1 )
292 IBLOCK( 1 ) = 1
293 M = 1
294 END IF
295 RETURN
296 END IF
297 *
298 * Compute Splitting Points
299 *
300 NSPLIT = 1
301 WORK( N ) = ZERO
302 PIVMIN = ONE
303 *
304 DO 10 J = 2, N
305 TMP1 = E( J-1 )**2
306 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
307 ISPLIT( NSPLIT ) = J - 1
308 NSPLIT = NSPLIT + 1
309 WORK( J-1 ) = ZERO
310 ELSE
311 WORK( J-1 ) = TMP1
312 PIVMIN = MAX( PIVMIN, TMP1 )
313 END IF
314 10 CONTINUE
315 ISPLIT( NSPLIT ) = N
316 PIVMIN = PIVMIN*SAFEMN
317 *
318 * Compute Interval and ATOLI
319 *
320 IF( IRANGE.EQ.3 ) THEN
321 *
322 * RANGE='I': Compute the interval containing eigenvalues
323 * IL through IU.
324 *
325 * Compute Gershgorin interval for entire (split) matrix
326 * and use it as the initial interval
327 *
328 GU = D( 1 )
329 GL = D( 1 )
330 TMP1 = ZERO
331 *
332 DO 20 J = 1, N - 1
333 TMP2 = SQRT( WORK( J ) )
334 GU = MAX( GU, D( J )+TMP1+TMP2 )
335 GL = MIN( GL, D( J )-TMP1-TMP2 )
336 TMP1 = TMP2
337 20 CONTINUE
338 *
339 GU = MAX( GU, D( N )+TMP1 )
340 GL = MIN( GL, D( N )-TMP1 )
341 TNORM = MAX( ABS( GL ), ABS( GU ) )
342 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
343 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
344 *
345 * Compute Iteration parameters
346 *
347 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
348 $ LOG( TWO ) ) + 2
349 IF( ABSTOL.LE.ZERO ) THEN
350 ATOLI = ULP*TNORM
351 ELSE
352 ATOLI = ABSTOL
353 END IF
354 *
355 WORK( N+1 ) = GL
356 WORK( N+2 ) = GL
357 WORK( N+3 ) = GU
358 WORK( N+4 ) = GU
359 WORK( N+5 ) = GL
360 WORK( N+6 ) = GU
361 IWORK( 1 ) = -1
362 IWORK( 2 ) = -1
363 IWORK( 3 ) = N + 1
364 IWORK( 4 ) = N + 1
365 IWORK( 5 ) = IL - 1
366 IWORK( 6 ) = IU
367 *
368 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
369 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
370 $ IWORK, W, IBLOCK, IINFO )
371 *
372 IF( IWORK( 6 ).EQ.IU ) THEN
373 WL = WORK( N+1 )
374 WLU = WORK( N+3 )
375 NWL = IWORK( 1 )
376 WU = WORK( N+4 )
377 WUL = WORK( N+2 )
378 NWU = IWORK( 4 )
379 ELSE
380 WL = WORK( N+2 )
381 WLU = WORK( N+4 )
382 NWL = IWORK( 2 )
383 WU = WORK( N+3 )
384 WUL = WORK( N+1 )
385 NWU = IWORK( 3 )
386 END IF
387 *
388 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
389 INFO = 4
390 RETURN
391 END IF
392 ELSE
393 *
394 * RANGE='A' or 'V' -- Set ATOLI
395 *
396 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
397 $ ABS( D( N ) )+ABS( E( N-1 ) ) )
398 *
399 DO 30 J = 2, N - 1
400 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
401 $ ABS( E( J ) ) )
402 30 CONTINUE
403 *
404 IF( ABSTOL.LE.ZERO ) THEN
405 ATOLI = ULP*TNORM
406 ELSE
407 ATOLI = ABSTOL
408 END IF
409 *
410 IF( IRANGE.EQ.2 ) THEN
411 WL = VL
412 WU = VU
413 ELSE
414 WL = ZERO
415 WU = ZERO
416 END IF
417 END IF
418 *
419 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
420 * NWL accumulates the number of eigenvalues .le. WL,
421 * NWU accumulates the number of eigenvalues .le. WU
422 *
423 M = 0
424 IEND = 0
425 INFO = 0
426 NWL = 0
427 NWU = 0
428 *
429 DO 70 JB = 1, NSPLIT
430 IOFF = IEND
431 IBEGIN = IOFF + 1
432 IEND = ISPLIT( JB )
433 IN = IEND - IOFF
434 *
435 IF( IN.EQ.1 ) THEN
436 *
437 * Special Case -- IN=1
438 *
439 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
440 $ NWL = NWL + 1
441 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
442 $ NWU = NWU + 1
443 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
444 $ D( IBEGIN )-PIVMIN ) ) THEN
445 M = M + 1
446 W( M ) = D( IBEGIN )
447 IBLOCK( M ) = JB
448 END IF
449 ELSE
450 *
451 * General Case -- IN > 1
452 *
453 * Compute Gershgorin Interval
454 * and use it as the initial interval
455 *
456 GU = D( IBEGIN )
457 GL = D( IBEGIN )
458 TMP1 = ZERO
459 *
460 DO 40 J = IBEGIN, IEND - 1
461 TMP2 = ABS( E( J ) )
462 GU = MAX( GU, D( J )+TMP1+TMP2 )
463 GL = MIN( GL, D( J )-TMP1-TMP2 )
464 TMP1 = TMP2
465 40 CONTINUE
466 *
467 GU = MAX( GU, D( IEND )+TMP1 )
468 GL = MIN( GL, D( IEND )-TMP1 )
469 BNORM = MAX( ABS( GL ), ABS( GU ) )
470 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
471 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
472 *
473 * Compute ATOLI for the current submatrix
474 *
475 IF( ABSTOL.LE.ZERO ) THEN
476 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
477 ELSE
478 ATOLI = ABSTOL
479 END IF
480 *
481 IF( IRANGE.GT.1 ) THEN
482 IF( GU.LT.WL ) THEN
483 NWL = NWL + IN
484 NWU = NWU + IN
485 GO TO 70
486 END IF
487 GL = MAX( GL, WL )
488 GU = MIN( GU, WU )
489 IF( GL.GE.GU )
490 $ GO TO 70
491 END IF
492 *
493 * Set Up Initial Interval
494 *
495 WORK( N+1 ) = GL
496 WORK( N+IN+1 ) = GU
497 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
498 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
499 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
500 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
501 *
502 NWL = NWL + IWORK( 1 )
503 NWU = NWU + IWORK( IN+1 )
504 IWOFF = M - IWORK( 1 )
505 *
506 * Compute Eigenvalues
507 *
508 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
509 $ LOG( TWO ) ) + 2
510 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
511 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
512 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
513 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
514 *
515 * Copy Eigenvalues Into W and IBLOCK
516 * Use -JB for block number for unconverged eigenvalues.
517 *
518 DO 60 J = 1, IOUT
519 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
520 *
521 * Flag non-convergence.
522 *
523 IF( J.GT.IOUT-IINFO ) THEN
524 NCNVRG = .TRUE.
525 IB = -JB
526 ELSE
527 IB = JB
528 END IF
529 DO 50 JE = IWORK( J ) + 1 + IWOFF,
530 $ IWORK( J+IN ) + IWOFF
531 W( JE ) = TMP1
532 IBLOCK( JE ) = IB
533 50 CONTINUE
534 60 CONTINUE
535 *
536 M = M + IM
537 END IF
538 70 CONTINUE
539 *
540 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
541 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
542 *
543 IF( IRANGE.EQ.3 ) THEN
544 IM = 0
545 IDISCL = IL - 1 - NWL
546 IDISCU = NWU - IU
547 *
548 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
549 DO 80 JE = 1, M
550 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
551 IDISCL = IDISCL - 1
552 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
553 IDISCU = IDISCU - 1
554 ELSE
555 IM = IM + 1
556 W( IM ) = W( JE )
557 IBLOCK( IM ) = IBLOCK( JE )
558 END IF
559 80 CONTINUE
560 M = IM
561 END IF
562 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
563 *
564 * Code to deal with effects of bad arithmetic:
565 * Some low eigenvalues to be discarded are not in (WL,WLU],
566 * or high eigenvalues to be discarded are not in (WUL,WU]
567 * so just kill off the smallest IDISCL/largest IDISCU
568 * eigenvalues, by simply finding the smallest/largest
569 * eigenvalue(s).
570 *
571 * (If N(w) is monotone non-decreasing, this should never
572 * happen.)
573 *
574 IF( IDISCL.GT.0 ) THEN
575 WKILL = WU
576 DO 100 JDISC = 1, IDISCL
577 IW = 0
578 DO 90 JE = 1, M
579 IF( IBLOCK( JE ).NE.0 .AND.
580 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
581 IW = JE
582 WKILL = W( JE )
583 END IF
584 90 CONTINUE
585 IBLOCK( IW ) = 0
586 100 CONTINUE
587 END IF
588 IF( IDISCU.GT.0 ) THEN
589 *
590 WKILL = WL
591 DO 120 JDISC = 1, IDISCU
592 IW = 0
593 DO 110 JE = 1, M
594 IF( IBLOCK( JE ).NE.0 .AND.
595 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
596 IW = JE
597 WKILL = W( JE )
598 END IF
599 110 CONTINUE
600 IBLOCK( IW ) = 0
601 120 CONTINUE
602 END IF
603 IM = 0
604 DO 130 JE = 1, M
605 IF( IBLOCK( JE ).NE.0 ) THEN
606 IM = IM + 1
607 W( IM ) = W( JE )
608 IBLOCK( IM ) = IBLOCK( JE )
609 END IF
610 130 CONTINUE
611 M = IM
612 END IF
613 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
614 TOOFEW = .TRUE.
615 END IF
616 END IF
617 *
618 * If ORDER='B', do nothing -- the eigenvalues are already sorted
619 * by block.
620 * If ORDER='E', sort the eigenvalues from smallest to largest
621 *
622 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
623 DO 150 JE = 1, M - 1
624 IE = 0
625 TMP1 = W( JE )
626 DO 140 J = JE + 1, M
627 IF( W( J ).LT.TMP1 ) THEN
628 IE = J
629 TMP1 = W( J )
630 END IF
631 140 CONTINUE
632 *
633 IF( IE.NE.0 ) THEN
634 ITMP1 = IBLOCK( IE )
635 W( IE ) = W( JE )
636 IBLOCK( IE ) = IBLOCK( JE )
637 W( JE ) = TMP1
638 IBLOCK( JE ) = ITMP1
639 END IF
640 150 CONTINUE
641 END IF
642 *
643 INFO = 0
644 IF( NCNVRG )
645 $ INFO = INFO + 1
646 IF( TOOFEW )
647 $ INFO = INFO + 2
648 RETURN
649 *
650 * End of DSTEBZ
651 *
652 END