*> \brief \b SSTEBZ * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download SSTEBZ + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, * M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, * INFO ) * * .. Scalar Arguments .. * CHARACTER ORDER, RANGE * INTEGER IL, INFO, IU, M, N, NSPLIT * REAL ABSTOL, VL, VU * .. * .. Array Arguments .. * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) * REAL D( * ), E( * ), W( * ), WORK( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] RANGE *> \verbatim *> 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. *> \endverbatim *> *> \param[in] ORDER *> \verbatim *> 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. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the tridiagonal matrix T. N >= 0. *> \endverbatim *> *> \param[in] VL *> \verbatim *> 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'. *> \endverbatim *> *> \param[in] VU *> \verbatim *> 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'. *> \endverbatim *> *> \param[in] IL *> \verbatim *> 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'. *> \endverbatim *> *> \param[in] IU *> \verbatim *> 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'. *> \endverbatim *> *> \param[in] ABSTOL *> \verbatim *> 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. *> \endverbatim *> *> \param[in] D *> \verbatim *> D is REAL array, dimension (N) *> The n diagonal elements of the tridiagonal matrix T. *> \endverbatim *> *> \param[in] E *> \verbatim *> E is REAL array, dimension (N-1) *> The (n-1) off-diagonal elements of the tridiagonal matrix T. *> \endverbatim *> *> \param[out] M *> \verbatim *> M is INTEGER *> The actual number of eigenvalues found. 0 <= M <= N. *> (See also the description of INFO=2,3.) *> \endverbatim *> *> \param[out] NSPLIT *> \verbatim *> NSPLIT is INTEGER *> The number of diagonal blocks in the matrix T. *> 1 <= NSPLIT <= N. *> \endverbatim *> *> \param[out] W *> \verbatim *> 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.) *> \endverbatim *> *> \param[out] IBLOCK *> \verbatim *> 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.) *> \endverbatim *> *> \param[out] ISPLIT *> \verbatim *> 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.) *> \endverbatim *> *> \param[out] WORK *> \verbatim *> WORK is REAL array, dimension (4*N) *> \endverbatim *> *> \param[out] IWORK *> \verbatim *> IWORK is INTEGER array, dimension (3*N) *> \endverbatim *> *> \param[out] INFO *> \verbatim *> 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. *> \endverbatim * *> \par Internal Parameters: * ========================= *> *> \verbatim *> 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. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup auxOTHERcomputational * * ===================================================================== SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, $ INFO ) * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * * .. Scalar Arguments .. CHARACTER ORDER, RANGE INTEGER IL, INFO, IU, M, N, NSPLIT REAL ABSTOL, VL, VU * .. * .. Array Arguments .. INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) REAL D( * ), E( * ), W( * ), WORK( * ) * .. * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE, TWO, HALF PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, $ HALF = 1.0E0 / TWO ) REAL FUDGE, RELFAC PARAMETER ( FUDGE = 2.1E0, RELFAC = 2.0E0 ) * .. * .. Local Scalars .. LOGICAL NCNVRG, TOOFEW INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX, $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL, $ NWU REAL ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN, $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL * .. * .. Local Arrays .. INTEGER IDUMMA( 1 ) * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV REAL SLAMCH EXTERNAL LSAME, ILAENV, SLAMCH * .. * .. External Subroutines .. EXTERNAL SLAEBZ, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT * .. * .. Executable Statements .. * INFO = 0 * * Decode RANGE * IF( LSAME( RANGE, 'A' ) ) THEN IRANGE = 1 ELSE IF( LSAME( RANGE, 'V' ) ) THEN IRANGE = 2 ELSE IF( LSAME( RANGE, 'I' ) ) THEN IRANGE = 3 ELSE IRANGE = 0 END IF * * Decode ORDER * IF( LSAME( ORDER, 'B' ) ) THEN IORDER = 2 ELSE IF( LSAME( ORDER, 'E' ) ) THEN IORDER = 1 ELSE IORDER = 0 END IF * * Check for Errors * IF( IRANGE.LE.0 ) THEN INFO = -1 ELSE IF( IORDER.LE.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( IRANGE.EQ.2 ) THEN IF( VL.GE.VU ) INFO = -5 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) $ THEN INFO = -6 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) $ THEN INFO = -7 END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'SSTEBZ', -INFO ) RETURN END IF * * Initialize error flags * INFO = 0 NCNVRG = .FALSE. TOOFEW = .FALSE. * * Quick return if possible * M = 0 IF( N.EQ.0 ) $ RETURN * * Simplifications: * IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) $ IRANGE = 1 * * Get machine constants * NB is the minimum vector length for vector bisection, or 0 * if only scalar is to be done. * SAFEMN = SLAMCH( 'S' ) ULP = SLAMCH( 'P' ) RTOLI = ULP*RELFAC NB = ILAENV( 1, 'SSTEBZ', ' ', N, -1, -1, -1 ) IF( NB.LE.1 ) $ NB = 0 * * Special Case when N=1 * IF( N.EQ.1 ) THEN NSPLIT = 1 ISPLIT( 1 ) = 1 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN M = 0 ELSE W( 1 ) = D( 1 ) IBLOCK( 1 ) = 1 M = 1 END IF RETURN END IF * * Compute Splitting Points * NSPLIT = 1 WORK( N ) = ZERO PIVMIN = ONE * DO 10 J = 2, N TMP1 = E( J-1 )**2 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN ISPLIT( NSPLIT ) = J - 1 NSPLIT = NSPLIT + 1 WORK( J-1 ) = ZERO ELSE WORK( J-1 ) = TMP1 PIVMIN = MAX( PIVMIN, TMP1 ) END IF 10 CONTINUE ISPLIT( NSPLIT ) = N PIVMIN = PIVMIN*SAFEMN * * Compute Interval and ATOLI * IF( IRANGE.EQ.3 ) THEN * * RANGE='I': Compute the interval containing eigenvalues * IL through IU. * * Compute Gershgorin interval for entire (split) matrix * and use it as the initial interval * GU = D( 1 ) GL = D( 1 ) TMP1 = ZERO * DO 20 J = 1, N - 1 TMP2 = SQRT( WORK( J ) ) GU = MAX( GU, D( J )+TMP1+TMP2 ) GL = MIN( GL, D( J )-TMP1-TMP2 ) TMP1 = TMP2 20 CONTINUE * GU = MAX( GU, D( N )+TMP1 ) GL = MIN( GL, D( N )-TMP1 ) TNORM = MAX( ABS( GL ), ABS( GU ) ) GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN * * Compute Iteration parameters * ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 IF( ABSTOL.LE.ZERO ) THEN ATOLI = ULP*TNORM ELSE ATOLI = ABSTOL END IF * WORK( N+1 ) = GL WORK( N+2 ) = GL WORK( N+3 ) = GU WORK( N+4 ) = GU WORK( N+5 ) = GL WORK( N+6 ) = GU IWORK( 1 ) = -1 IWORK( 2 ) = -1 IWORK( 3 ) = N + 1 IWORK( 4 ) = N + 1 IWORK( 5 ) = IL - 1 IWORK( 6 ) = IU * CALL SLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, $ IWORK, W, IBLOCK, IINFO ) * IF( IWORK( 6 ).EQ.IU ) THEN WL = WORK( N+1 ) WLU = WORK( N+3 ) NWL = IWORK( 1 ) WU = WORK( N+4 ) WUL = WORK( N+2 ) NWU = IWORK( 4 ) ELSE WL = WORK( N+2 ) WLU = WORK( N+4 ) NWL = IWORK( 2 ) WU = WORK( N+3 ) WUL = WORK( N+1 ) NWU = IWORK( 3 ) END IF * IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN INFO = 4 RETURN END IF ELSE * * RANGE='A' or 'V' -- Set ATOLI * TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), $ ABS( D( N ) )+ABS( E( N-1 ) ) ) * DO 30 J = 2, N - 1 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ $ ABS( E( J ) ) ) 30 CONTINUE * IF( ABSTOL.LE.ZERO ) THEN ATOLI = ULP*TNORM ELSE ATOLI = ABSTOL END IF * IF( IRANGE.EQ.2 ) THEN WL = VL WU = VU ELSE WL = ZERO WU = ZERO END IF END IF * * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. * NWL accumulates the number of eigenvalues .le. WL, * NWU accumulates the number of eigenvalues .le. WU * M = 0 IEND = 0 INFO = 0 NWL = 0 NWU = 0 * DO 70 JB = 1, NSPLIT IOFF = IEND IBEGIN = IOFF + 1 IEND = ISPLIT( JB ) IN = IEND - IOFF * IF( IN.EQ.1 ) THEN * * Special Case -- IN=1 * IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) $ NWL = NWL + 1 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) $ NWU = NWU + 1 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. $ D( IBEGIN )-PIVMIN ) ) THEN M = M + 1 W( M ) = D( IBEGIN ) IBLOCK( M ) = JB END IF ELSE * * General Case -- IN > 1 * * Compute Gershgorin Interval * and use it as the initial interval * GU = D( IBEGIN ) GL = D( IBEGIN ) TMP1 = ZERO * DO 40 J = IBEGIN, IEND - 1 TMP2 = ABS( E( J ) ) GU = MAX( GU, D( J )+TMP1+TMP2 ) GL = MIN( GL, D( J )-TMP1-TMP2 ) TMP1 = TMP2 40 CONTINUE * GU = MAX( GU, D( IEND )+TMP1 ) GL = MIN( GL, D( IEND )-TMP1 ) BNORM = MAX( ABS( GL ), ABS( GU ) ) GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN * * Compute ATOLI for the current submatrix * IF( ABSTOL.LE.ZERO ) THEN ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) ELSE ATOLI = ABSTOL END IF * IF( IRANGE.GT.1 ) THEN IF( GU.LT.WL ) THEN NWL = NWL + IN NWU = NWU + IN GO TO 70 END IF GL = MAX( GL, WL ) GU = MIN( GU, WU ) IF( GL.GE.GU ) $ GO TO 70 END IF * * Set Up Initial Interval * WORK( N+1 ) = GL WORK( N+IN+1 ) = GU CALL SLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) * NWL = NWL + IWORK( 1 ) NWU = NWU + IWORK( IN+1 ) IWOFF = M - IWORK( 1 ) * * Compute Eigenvalues * ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 CALL SLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) * * Copy Eigenvalues Into W and IBLOCK * Use -JB for block number for unconverged eigenvalues. * DO 60 J = 1, IOUT TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) * * Flag non-convergence. * IF( J.GT.IOUT-IINFO ) THEN NCNVRG = .TRUE. IB = -JB ELSE IB = JB END IF DO 50 JE = IWORK( J ) + 1 + IWOFF, $ IWORK( J+IN ) + IWOFF W( JE ) = TMP1 IBLOCK( JE ) = IB 50 CONTINUE 60 CONTINUE * M = M + IM END IF 70 CONTINUE * * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. * IF( IRANGE.EQ.3 ) THEN IM = 0 IDISCL = IL - 1 - NWL IDISCU = NWU - IU * IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN DO 80 JE = 1, M IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN IDISCL = IDISCL - 1 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN IDISCU = IDISCU - 1 ELSE IM = IM + 1 W( IM ) = W( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 80 CONTINUE M = IM END IF IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN * * Code to deal with effects of bad arithmetic: * Some low eigenvalues to be discarded are not in (WL,WLU], * or high eigenvalues to be discarded are not in (WUL,WU] * so just kill off the smallest IDISCL/largest IDISCU * eigenvalues, by simply finding the smallest/largest * eigenvalue(s). * * (If N(w) is monotone non-decreasing, this should never * happen.) * IF( IDISCL.GT.0 ) THEN WKILL = WU DO 100 JDISC = 1, IDISCL IW = 0 DO 90 JE = 1, M IF( IBLOCK( JE ).NE.0 .AND. $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN IW = JE WKILL = W( JE ) END IF 90 CONTINUE IBLOCK( IW ) = 0 100 CONTINUE END IF IF( IDISCU.GT.0 ) THEN * WKILL = WL DO 120 JDISC = 1, IDISCU IW = 0 DO 110 JE = 1, M IF( IBLOCK( JE ).NE.0 .AND. $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN IW = JE WKILL = W( JE ) END IF 110 CONTINUE IBLOCK( IW ) = 0 120 CONTINUE END IF IM = 0 DO 130 JE = 1, M IF( IBLOCK( JE ).NE.0 ) THEN IM = IM + 1 W( IM ) = W( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 130 CONTINUE M = IM END IF IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN TOOFEW = .TRUE. END IF END IF * * If ORDER='B', do nothing -- the eigenvalues are already sorted * by block. * If ORDER='E', sort the eigenvalues from smallest to largest * IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN DO 150 JE = 1, M - 1 IE = 0 TMP1 = W( JE ) DO 140 J = JE + 1, M IF( W( J ).LT.TMP1 ) THEN IE = J TMP1 = W( J ) END IF 140 CONTINUE * IF( IE.NE.0 ) THEN ITMP1 = IBLOCK( IE ) W( IE ) = W( JE ) IBLOCK( IE ) = IBLOCK( JE ) W( JE ) = TMP1 IBLOCK( JE ) = ITMP1 END IF 150 CONTINUE END IF * INFO = 0 IF( NCNVRG ) $ INFO = INFO + 1 IF( TOOFEW ) $ INFO = INFO + 2 RETURN * * End of SSTEBZ * END