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