1*> \brief \b DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARRK + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrk.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrk.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrk.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARRK( N, IW, GL, GU, 22* D, E2, PIVMIN, RELTOL, W, WERR, INFO) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, IW, N 26* DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION D( * ), E2( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DLARRK computes one eigenvalue of a symmetric tridiagonal 39*> matrix T to suitable accuracy. This is an auxiliary code to be 40*> called from DSTEMR. 41*> 42*> To avoid overflow, the matrix must be scaled so that its 43*> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest 44*> accuracy, it should not be much smaller than that. 45*> 46*> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 47*> Matrix", Report CS41, Computer Science Dept., Stanford 48*> University, July 21, 1966. 49*> \endverbatim 50* 51* Arguments: 52* ========== 53* 54*> \param[in] N 55*> \verbatim 56*> N is INTEGER 57*> The order of the tridiagonal matrix T. N >= 0. 58*> \endverbatim 59*> 60*> \param[in] IW 61*> \verbatim 62*> IW is INTEGER 63*> The index of the eigenvalues to be returned. 64*> \endverbatim 65*> 66*> \param[in] GL 67*> \verbatim 68*> GL is DOUBLE PRECISION 69*> \endverbatim 70*> 71*> \param[in] GU 72*> \verbatim 73*> GU is DOUBLE PRECISION 74*> An upper and a lower bound on the eigenvalue. 75*> \endverbatim 76*> 77*> \param[in] D 78*> \verbatim 79*> D is DOUBLE PRECISION array, dimension (N) 80*> The n diagonal elements of the tridiagonal matrix T. 81*> \endverbatim 82*> 83*> \param[in] E2 84*> \verbatim 85*> E2 is DOUBLE PRECISION array, dimension (N-1) 86*> The (n-1) squared off-diagonal elements of the tridiagonal matrix T. 87*> \endverbatim 88*> 89*> \param[in] PIVMIN 90*> \verbatim 91*> PIVMIN is DOUBLE PRECISION 92*> The minimum pivot allowed in the Sturm sequence for T. 93*> \endverbatim 94*> 95*> \param[in] RELTOL 96*> \verbatim 97*> RELTOL is DOUBLE PRECISION 98*> The minimum relative width of an interval. When an interval 99*> is narrower than RELTOL times the larger (in 100*> magnitude) endpoint, then it is considered to be 101*> sufficiently small, i.e., converged. Note: this should 102*> always be at least radix*machine epsilon. 103*> \endverbatim 104*> 105*> \param[out] W 106*> \verbatim 107*> W is DOUBLE PRECISION 108*> \endverbatim 109*> 110*> \param[out] WERR 111*> \verbatim 112*> WERR is DOUBLE PRECISION 113*> The error bound on the corresponding eigenvalue approximation 114*> in W. 115*> \endverbatim 116*> 117*> \param[out] INFO 118*> \verbatim 119*> INFO is INTEGER 120*> = 0: Eigenvalue converged 121*> = -1: Eigenvalue did NOT converge 122*> \endverbatim 123* 124*> \par Internal Parameters: 125* ========================= 126*> 127*> \verbatim 128*> FUDGE DOUBLE PRECISION, default = 2 129*> A "fudge factor" to widen the Gershgorin intervals. 130*> \endverbatim 131* 132* Authors: 133* ======== 134* 135*> \author Univ. of Tennessee 136*> \author Univ. of California Berkeley 137*> \author Univ. of Colorado Denver 138*> \author NAG Ltd. 139* 140*> \ingroup OTHERauxiliary 141* 142* ===================================================================== 143 SUBROUTINE DLARRK( N, IW, GL, GU, 144 $ D, E2, PIVMIN, RELTOL, W, WERR, INFO) 145* 146* -- LAPACK auxiliary routine -- 147* -- LAPACK is a software package provided by Univ. of Tennessee, -- 148* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 149* 150* .. Scalar Arguments .. 151 INTEGER INFO, IW, N 152 DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR 153* .. 154* .. Array Arguments .. 155 DOUBLE PRECISION D( * ), E2( * ) 156* .. 157* 158* ===================================================================== 159* 160* .. Parameters .. 161 DOUBLE PRECISION FUDGE, HALF, TWO, ZERO 162 PARAMETER ( HALF = 0.5D0, TWO = 2.0D0, 163 $ FUDGE = TWO, ZERO = 0.0D0 ) 164* .. 165* .. Local Scalars .. 166 INTEGER I, IT, ITMAX, NEGCNT 167 DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1, 168 $ TMP2, TNORM 169* .. 170* .. External Functions .. 171 DOUBLE PRECISION DLAMCH 172 EXTERNAL DLAMCH 173* .. 174* .. Intrinsic Functions .. 175 INTRINSIC ABS, INT, LOG, MAX 176* .. 177* .. Executable Statements .. 178* 179* Quick return if possible 180* 181 IF( N.LE.0 ) THEN 182 INFO = 0 183 RETURN 184 END IF 185* 186* Get machine constants 187 EPS = DLAMCH( 'P' ) 188 189 TNORM = MAX( ABS( GL ), ABS( GU ) ) 190 RTOLI = RELTOL 191 ATOLI = FUDGE*TWO*PIVMIN 192 193 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 194 $ LOG( TWO ) ) + 2 195 196 INFO = -1 197 198 LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN 199 RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN 200 IT = 0 201 202 10 CONTINUE 203* 204* Check if interval converged or maximum number of iterations reached 205* 206 TMP1 = ABS( RIGHT - LEFT ) 207 TMP2 = MAX( ABS(RIGHT), ABS(LEFT) ) 208 IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN 209 INFO = 0 210 GOTO 30 211 ENDIF 212 IF(IT.GT.ITMAX) 213 $ GOTO 30 214 215* 216* Count number of negative pivots for mid-point 217* 218 IT = IT + 1 219 MID = HALF * (LEFT + RIGHT) 220 NEGCNT = 0 221 TMP1 = D( 1 ) - MID 222 IF( ABS( TMP1 ).LT.PIVMIN ) 223 $ TMP1 = -PIVMIN 224 IF( TMP1.LE.ZERO ) 225 $ NEGCNT = NEGCNT + 1 226* 227 DO 20 I = 2, N 228 TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID 229 IF( ABS( TMP1 ).LT.PIVMIN ) 230 $ TMP1 = -PIVMIN 231 IF( TMP1.LE.ZERO ) 232 $ NEGCNT = NEGCNT + 1 233 20 CONTINUE 234 235 IF(NEGCNT.GE.IW) THEN 236 RIGHT = MID 237 ELSE 238 LEFT = MID 239 ENDIF 240 GOTO 10 241 242 30 CONTINUE 243* 244* Converged or maximum number of iterations reached 245* 246 W = HALF * (LEFT + RIGHT) 247 WERR = HALF * ABS( RIGHT - LEFT ) 248 249 RETURN 250* 251* End of DLARRK 252* 253 END 254