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