1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
|
*> \brief \b SSTECT
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE SSTECT( N, A, B, SHIFT, NUM )
*
* .. Scalar Arguments ..
* INTEGER N, NUM
* REAL SHIFT
* ..
* .. Array Arguments ..
* REAL A( * ), B( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> SSTECT 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
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The dimension of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is REAL array, dimension (N)
*> The diagonal entries of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*> B is REAL array, dimension (N-1)
*> The offdiagonal entries of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] SHIFT
*> \verbatim
*> SHIFT is REAL
*> The shift, used as described under Purpose.
*> \endverbatim
*>
*> \param[out] NUM
*> \verbatim
*> NUM is INTEGER
*> The number of eigenvalues of T less than or equal
*> to SHIFT.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup single_eig
*
* =====================================================================
SUBROUTINE SSTECT( N, A, B, SHIFT, NUM )
*
* -- LAPACK test routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2011
*
* .. Scalar Arguments ..
INTEGER N, NUM
REAL SHIFT
* ..
* .. Array Arguments ..
REAL A( * ), B( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE, THREE
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, THREE = 3.0E0 )
* ..
* .. Local Scalars ..
INTEGER I
REAL M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
$ TOM, U, UNFL
* ..
* .. External Functions ..
REAL SLAMCH
EXTERNAL SLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT
* ..
* .. Executable Statements ..
*
* Get machine constants
*
UNFL = SLAMCH( 'Safe minimum' )
OVFL = SLAMCH( '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 SSTECT
*
END
|