SUBROUTINE SLARRK( N, IW, GL, GU, \$ D, E2, PIVMIN, RELTOL, W, WERR, INFO) IMPLICIT NONE * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, IW, N REAL PIVMIN, RELTOL, GL, GU, W, WERR * .. * .. Array Arguments .. REAL D( * ), E2( * ) * .. * * Purpose * ======= * * SLARRK computes one eigenvalue of a symmetric tridiagonal * matrix T to suitable accuracy. This is an auxiliary code to be * called from SSTEMR. * * 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. * * Arguments * ========= * * N (input) INTEGER * The order of the tridiagonal matrix T. N >= 0. * * IW (input) INTEGER * The index of the eigenvalues to be returned. * * GL (input) REAL * GU (input) REAL * An upper and a lower bound on the eigenvalue. * * D (input) REAL array, dimension (N) * The n diagonal elements of the tridiagonal matrix T. * * E2 (input) REAL array, dimension (N-1) * The (n-1) squared off-diagonal elements of the tridiagonal matrix T. * * PIVMIN (input) REAL * The minimum pivot allowed in the Sturm sequence for T. * * RELTOL (input) REAL * The minimum relative width of an interval. When an interval * is narrower than RELTOL times the larger (in * magnitude) endpoint, then it is considered to be * sufficiently small, i.e., converged. Note: this should * always be at least radix*machine epsilon. * * W (output) REAL * * WERR (output) REAL * The error bound on the corresponding eigenvalue approximation * in W. * * INFO (output) INTEGER * = 0: Eigenvalue converged * = -1: Eigenvalue did NOT converge * * Internal Parameters * =================== * * FUDGE REAL , default = 2 * A "fudge factor" to widen the Gershgorin intervals. * * ===================================================================== * * .. Parameters .. REAL FUDGE, HALF, TWO, ZERO PARAMETER ( HALF = 0.5E0, TWO = 2.0E0, \$ FUDGE = TWO, ZERO = 0.0E0 ) * .. * .. Local Scalars .. INTEGER I, IT, ITMAX, NEGCNT REAL ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1, \$ TMP2, TNORM * .. * .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX * .. * .. Executable Statements .. * * Get machine constants EPS = SLAMCH( 'P' ) TNORM = MAX( ABS( GL ), ABS( GU ) ) RTOLI = RELTOL ATOLI = FUDGE*TWO*PIVMIN ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / \$ LOG( TWO ) ) + 2 INFO = -1 LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN IT = 0 10 CONTINUE * * Check if interval converged or maximum number of iterations reached * TMP1 = ABS( RIGHT - LEFT ) TMP2 = MAX( ABS(RIGHT), ABS(LEFT) ) IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN INFO = 0 GOTO 30 ENDIF IF(IT.GT.ITMAX) \$ GOTO 30 * * Count number of negative pivots for mid-point * IT = IT + 1 MID = HALF * (LEFT + RIGHT) NEGCNT = 0 TMP1 = D( 1 ) - MID IF( ABS( TMP1 ).LT.PIVMIN ) \$ TMP1 = -PIVMIN IF( TMP1.LE.ZERO ) \$ NEGCNT = NEGCNT + 1 * DO 20 I = 2, N TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID IF( ABS( TMP1 ).LT.PIVMIN ) \$ TMP1 = -PIVMIN IF( TMP1.LE.ZERO ) \$ NEGCNT = NEGCNT + 1 20 CONTINUE IF(NEGCNT.GE.IW) THEN RIGHT = MID ELSE LEFT = MID ENDIF GOTO 10 30 CONTINUE * * Converged or maximum number of iterations reached * W = HALF * (LEFT + RIGHT) WERR = HALF * ABS( RIGHT - LEFT ) RETURN * * End of SLARRK * END