SUBROUTINE DSTECT( N, A, B, SHIFT, NUM ) * * -- LAPACK test routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER N, NUM DOUBLE PRECISION SHIFT * .. * .. Array Arguments .. DOUBLE PRECISION A( * ), B( * ) * .. * * Purpose * ======= * * DSTECT counts the number NUM of eigenvalues of a tridiagonal * matrix T which are less than or equal to SHIFT. T has * diagonal entries A(1), ... , A(N), and offdiagonal entries * B(1), ..., B(N-1). * 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 dimension of the tridiagonal matrix T. * * A (input) DOUBLE PRECISION array, dimension (N) * The diagonal entries of the tridiagonal matrix T. * * B (input) DOUBLE PRECISION array, dimension (N-1) * The offdiagonal entries of the tridiagonal matrix T. * * SHIFT (input) DOUBLE PRECISION * The shift, used as described under Purpose. * * NUM (output) INTEGER * The number of eigenvalues of T less than or equal * to SHIFT. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, THREE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, THREE = 3.0D0 ) * .. * .. Local Scalars .. INTEGER I DOUBLE PRECISION M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP, $ TOM, U, UNFL * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. * .. Executable Statements .. * * Get machine constants * UNFL = DLAMCH( 'Safe minimum' ) OVFL = DLAMCH( 'Overflow' ) * * Find largest entry * MX = ABS( A( 1 ) ) DO 10 I = 1, N - 1 MX = MAX( MX, ABS( A( I+1 ) ), ABS( B( I ) ) ) 10 CONTINUE * * Handle easy cases, including zero matrix * IF( SHIFT.GE.THREE*MX ) THEN NUM = N RETURN END IF IF( SHIFT.LT.-THREE*MX ) THEN NUM = 0 RETURN END IF * * Compute scale factors as in Kahan's report * At this point, MX .NE. 0 so we can divide by it * SUN = SQRT( UNFL ) SSUN = SQRT( SUN ) SOV = SQRT( OVFL ) TOM = SSUN*SOV IF( MX.LE.ONE ) THEN M1 = ONE / MX M2 = TOM ELSE M1 = ONE M2 = TOM / MX END IF * * Begin counting * NUM = 0 SSHIFT = ( SHIFT*M1 )*M2 U = ( A( 1 )*M1 )*M2 - SSHIFT IF( U.LE.SUN ) THEN IF( U.LE.ZERO ) THEN NUM = NUM + 1 IF( U.GT.-SUN ) $ U = -SUN ELSE U = SUN END IF END IF DO 20 I = 2, N TMP = ( B( I-1 )*M1 )*M2 U = ( ( A( I )*M1 )*M2-TMP*( TMP / U ) ) - SSHIFT IF( U.LE.SUN ) THEN IF( U.LE.ZERO ) THEN NUM = NUM + 1 IF( U.GT.-SUN ) $ U = -SUN ELSE U = SUN END IF END IF 20 CONTINUE RETURN * * End of DSTECT * END