 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ 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: