LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ sstebz()

 subroutine sstebz ( character RANGE, character ORDER, integer N, real VL, real VU, integer IL, integer IU, real ABSTOL, real, dimension( * ) D, real, dimension( * ) E, integer M, integer NSPLIT, real, dimension( * ) W, integer, dimension( * ) IBLOCK, integer, dimension( * ) ISPLIT, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO )

SSTEBZ

Purpose:
``` SSTEBZ computes the eigenvalues of a symmetric tridiagonal
matrix T.  The user may ask for all eigenvalues, all eigenvalues
in the half-open interval (VL, VU], or the IL-th through IU-th
eigenvalues.

To avoid overflow, the matrix must be scaled so that its
largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
accuracy, it should not be much smaller than that.

See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
Matrix", Report CS41, Computer Science Dept., Stanford
University, July 21, 1966.```
Parameters
 [in] RANGE ``` RANGE is CHARACTER*1 = 'A': ("All") all eigenvalues will be found. = 'V': ("Value") all eigenvalues in the half-open interval (VL, VU] will be found. = 'I': ("Index") the IL-th through IU-th eigenvalues (of the entire matrix) will be found.``` [in] ORDER ``` ORDER is CHARACTER*1 = 'B': ("By Block") the eigenvalues will be grouped by split-off block (see IBLOCK, ISPLIT) and ordered from smallest to largest within the block. = 'E': ("Entire matrix") the eigenvalues for the entire matrix will be ordered from smallest to largest.``` [in] N ``` N is INTEGER The order of the tridiagonal matrix T. N >= 0.``` [in] VL ``` VL is REAL If RANGE='V', the lower bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.``` [in] VU ``` VU is REAL If RANGE='V', the upper bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.``` [in] IL ``` IL is INTEGER If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.``` [in] IU ``` IU is INTEGER If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.``` [in] ABSTOL ``` ABSTOL is REAL The absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less. If ABSTOL is less than or equal to zero, then ULP*|T| will be used, where |T| means the 1-norm of T. Eigenvalues will be computed most accurately when ABSTOL is set to twice the underflow threshold 2*SLAMCH('S'), not zero.``` [in] D ``` D is REAL array, dimension (N) The n diagonal elements of the tridiagonal matrix T.``` [in] E ``` E is REAL array, dimension (N-1) The (n-1) off-diagonal elements of the tridiagonal matrix T.``` [out] M ``` M is INTEGER The actual number of eigenvalues found. 0 <= M <= N. (See also the description of INFO=2,3.)``` [out] NSPLIT ``` NSPLIT is INTEGER The number of diagonal blocks in the matrix T. 1 <= NSPLIT <= N.``` [out] W ``` W is REAL array, dimension (N) On exit, the first M elements of W will contain the eigenvalues. (SSTEBZ may use the remaining N-M elements as workspace.)``` [out] IBLOCK ``` IBLOCK is INTEGER array, dimension (N) At each row/column j where E(j) is zero or small, the matrix T is considered to split into a block diagonal matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which block (from 1 to the number of blocks) the eigenvalue W(i) belongs. (SSTEBZ may use the remaining N-M elements as workspace.)``` [out] ISPLIT ``` ISPLIT is INTEGER array, dimension (N) The splitting points, at which T breaks up into submatrices. The first submatrix consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. (Only the first NSPLIT elements will actually be used, but since the user cannot know a priori what value NSPLIT will have, N words must be reserved for ISPLIT.)``` [out] WORK ` WORK is REAL array, dimension (4*N)` [out] IWORK ` IWORK is INTEGER array, dimension (3*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: some or all of the eigenvalues failed to converge or were not computed: =1 or 3: Bisection failed to converge for some eigenvalues; these eigenvalues are flagged by a negative block number. The effect is that the eigenvalues may not be as accurate as the absolute and relative tolerances. This is generally caused by unexpectedly inaccurate arithmetic. =2 or 3: RANGE='I' only: Not all of the eigenvalues IL:IU were found. Effect: M < IU+1-IL Cause: non-monotonic arithmetic, causing the Sturm sequence to be non-monotonic. Cure: recalculate, using RANGE='A', and pick out eigenvalues IL:IU. In some cases, increasing the PARAMETER "FUDGE" may make things work. = 4: RANGE='I', and the Gershgorin interval initially used was too small. No eigenvalues were computed. Probable cause: your machine has sloppy floating-point arithmetic. Cure: Increase the PARAMETER "FUDGE", recompile, and try again.```
Internal Parameters:
```  RELFAC  REAL, default = 2.0e0
The relative tolerance.  An interval (a,b] lies within
"relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),
where "ulp" is the machine precision (distance from 1 to
the next larger floating point number.)

FUDGE   REAL, default = 2
A "fudge factor" to widen the Gershgorin intervals.  Ideally,
a value of 1 should work, but on machines with sloppy
arithmetic, this needs to be larger.  The default for
publicly released versions should be large enough to handle
the worst machine around.  Note that this has no effect
on accuracy of the solution.```

Definition at line 270 of file sstebz.f.

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*
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: