 LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 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.```
Date
June 2016

Definition at line 275 of file sstebz.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: