1*> \brief \b DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARRJ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrj.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrj.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrj.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, 22* RTOL, OFFSET, W, WERR, WORK, IWORK, 23* PIVMIN, SPDIAM, INFO ) 24* 25* .. Scalar Arguments .. 26* INTEGER IFIRST, ILAST, INFO, N, OFFSET 27* DOUBLE PRECISION PIVMIN, RTOL, SPDIAM 28* .. 29* .. Array Arguments .. 30* INTEGER IWORK( * ) 31* DOUBLE PRECISION D( * ), E2( * ), W( * ), 32* $ WERR( * ), WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> Given the initial eigenvalue approximations of T, DLARRJ 42*> does bisection to refine the eigenvalues of T, 43*> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial 44*> guesses for these eigenvalues are input in W, the corresponding estimate 45*> of the error in these guesses in WERR. During bisection, intervals 46*> [left, right] are maintained by storing their mid-points and 47*> semi-widths in the arrays W and WERR respectively. 48*> \endverbatim 49* 50* Arguments: 51* ========== 52* 53*> \param[in] N 54*> \verbatim 55*> N is INTEGER 56*> The order of the matrix. 57*> \endverbatim 58*> 59*> \param[in] D 60*> \verbatim 61*> D is DOUBLE PRECISION array, dimension (N) 62*> The N diagonal elements of T. 63*> \endverbatim 64*> 65*> \param[in] E2 66*> \verbatim 67*> E2 is DOUBLE PRECISION array, dimension (N-1) 68*> The Squares of the (N-1) subdiagonal elements of T. 69*> \endverbatim 70*> 71*> \param[in] IFIRST 72*> \verbatim 73*> IFIRST is INTEGER 74*> The index of the first eigenvalue to be computed. 75*> \endverbatim 76*> 77*> \param[in] ILAST 78*> \verbatim 79*> ILAST is INTEGER 80*> The index of the last eigenvalue to be computed. 81*> \endverbatim 82*> 83*> \param[in] RTOL 84*> \verbatim 85*> RTOL is DOUBLE PRECISION 86*> Tolerance for the convergence of the bisection intervals. 87*> An interval [LEFT,RIGHT] has converged if 88*> RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|). 89*> \endverbatim 90*> 91*> \param[in] OFFSET 92*> \verbatim 93*> OFFSET is INTEGER 94*> Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET 95*> through ILAST-OFFSET elements of these arrays are to be used. 96*> \endverbatim 97*> 98*> \param[in,out] W 99*> \verbatim 100*> W is DOUBLE PRECISION array, dimension (N) 101*> On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are 102*> estimates of the eigenvalues of L D L^T indexed IFIRST through 103*> ILAST. 104*> On output, these estimates are refined. 105*> \endverbatim 106*> 107*> \param[in,out] WERR 108*> \verbatim 109*> WERR is DOUBLE PRECISION array, dimension (N) 110*> On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are 111*> the errors in the estimates of the corresponding elements in W. 112*> On output, these errors are refined. 113*> \endverbatim 114*> 115*> \param[out] WORK 116*> \verbatim 117*> WORK is DOUBLE PRECISION array, dimension (2*N) 118*> Workspace. 119*> \endverbatim 120*> 121*> \param[out] IWORK 122*> \verbatim 123*> IWORK is INTEGER array, dimension (2*N) 124*> Workspace. 125*> \endverbatim 126*> 127*> \param[in] PIVMIN 128*> \verbatim 129*> PIVMIN is DOUBLE PRECISION 130*> The minimum pivot in the Sturm sequence for T. 131*> \endverbatim 132*> 133*> \param[in] SPDIAM 134*> \verbatim 135*> SPDIAM is DOUBLE PRECISION 136*> The spectral diameter of T. 137*> \endverbatim 138*> 139*> \param[out] INFO 140*> \verbatim 141*> INFO is INTEGER 142*> Error flag. 143*> \endverbatim 144* 145* Authors: 146* ======== 147* 148*> \author Univ. of Tennessee 149*> \author Univ. of California Berkeley 150*> \author Univ. of Colorado Denver 151*> \author NAG Ltd. 152* 153*> \date June 2017 154* 155*> \ingroup OTHERauxiliary 156* 157*> \par Contributors: 158* ================== 159*> 160*> Beresford Parlett, University of California, Berkeley, USA \n 161*> Jim Demmel, University of California, Berkeley, USA \n 162*> Inderjit Dhillon, University of Texas, Austin, USA \n 163*> Osni Marques, LBNL/NERSC, USA \n 164*> Christof Voemel, University of California, Berkeley, USA 165* 166* ===================================================================== 167 SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, 168 $ RTOL, OFFSET, W, WERR, WORK, IWORK, 169 $ PIVMIN, SPDIAM, INFO ) 170* 171* -- LAPACK auxiliary routine (version 3.7.1) -- 172* -- LAPACK is a software package provided by Univ. of Tennessee, -- 173* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 174* June 2017 175* 176* .. Scalar Arguments .. 177 INTEGER IFIRST, ILAST, INFO, N, OFFSET 178 DOUBLE PRECISION PIVMIN, RTOL, SPDIAM 179* .. 180* .. Array Arguments .. 181 INTEGER IWORK( * ) 182 DOUBLE PRECISION D( * ), E2( * ), W( * ), 183 $ WERR( * ), WORK( * ) 184* .. 185* 186* ===================================================================== 187* 188* .. Parameters .. 189 DOUBLE PRECISION ZERO, ONE, TWO, HALF 190 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 191 $ HALF = 0.5D0 ) 192 INTEGER MAXITR 193* .. 194* .. Local Scalars .. 195 INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT, 196 $ OLNINT, P, PREV, SAVI1 197 DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH 198* 199* .. 200* .. Intrinsic Functions .. 201 INTRINSIC ABS, MAX 202* .. 203* .. Executable Statements .. 204* 205 INFO = 0 206* 207* Quick return if possible 208* 209 IF( N.LE.0 ) THEN 210 RETURN 211 END IF 212* 213 MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) / 214 $ LOG( TWO ) ) + 2 215* 216* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. 217* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while 218* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) 219* for an unconverged interval is set to the index of the next unconverged 220* interval, and is -1 or 0 for a converged interval. Thus a linked 221* list of unconverged intervals is set up. 222* 223 224 I1 = IFIRST 225 I2 = ILAST 226* The number of unconverged intervals 227 NINT = 0 228* The last unconverged interval found 229 PREV = 0 230 DO 75 I = I1, I2 231 K = 2*I 232 II = I - OFFSET 233 LEFT = W( II ) - WERR( II ) 234 MID = W(II) 235 RIGHT = W( II ) + WERR( II ) 236 WIDTH = RIGHT - MID 237 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 238 239* The following test prevents the test of converged intervals 240 IF( WIDTH.LT.RTOL*TMP ) THEN 241* This interval has already converged and does not need refinement. 242* (Note that the gaps might change through refining the 243* eigenvalues, however, they can only get bigger.) 244* Remove it from the list. 245 IWORK( K-1 ) = -1 246* Make sure that I1 always points to the first unconverged interval 247 IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1 248 IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1 249 ELSE 250* unconverged interval found 251 PREV = I 252* Make sure that [LEFT,RIGHT] contains the desired eigenvalue 253* 254* Do while( CNT(LEFT).GT.I-1 ) 255* 256 FAC = ONE 257 20 CONTINUE 258 CNT = 0 259 S = LEFT 260 DPLUS = D( 1 ) - S 261 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 262 DO 30 J = 2, N 263 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 264 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 265 30 CONTINUE 266 IF( CNT.GT.I-1 ) THEN 267 LEFT = LEFT - WERR( II )*FAC 268 FAC = TWO*FAC 269 GO TO 20 270 END IF 271* 272* Do while( CNT(RIGHT).LT.I ) 273* 274 FAC = ONE 275 50 CONTINUE 276 CNT = 0 277 S = RIGHT 278 DPLUS = D( 1 ) - S 279 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 280 DO 60 J = 2, N 281 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 282 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 283 60 CONTINUE 284 IF( CNT.LT.I ) THEN 285 RIGHT = RIGHT + WERR( II )*FAC 286 FAC = TWO*FAC 287 GO TO 50 288 END IF 289 NINT = NINT + 1 290 IWORK( K-1 ) = I + 1 291 IWORK( K ) = CNT 292 END IF 293 WORK( K-1 ) = LEFT 294 WORK( K ) = RIGHT 295 75 CONTINUE 296 297 298 SAVI1 = I1 299* 300* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals 301* and while (ITER.LT.MAXITR) 302* 303 ITER = 0 304 80 CONTINUE 305 PREV = I1 - 1 306 I = I1 307 OLNINT = NINT 308 309 DO 100 P = 1, OLNINT 310 K = 2*I 311 II = I - OFFSET 312 NEXT = IWORK( K-1 ) 313 LEFT = WORK( K-1 ) 314 RIGHT = WORK( K ) 315 MID = HALF*( LEFT + RIGHT ) 316 317* semiwidth of interval 318 WIDTH = RIGHT - MID 319 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 320 321 IF( ( WIDTH.LT.RTOL*TMP ) .OR. 322 $ (ITER.EQ.MAXITR) )THEN 323* reduce number of unconverged intervals 324 NINT = NINT - 1 325* Mark interval as converged. 326 IWORK( K-1 ) = 0 327 IF( I1.EQ.I ) THEN 328 I1 = NEXT 329 ELSE 330* Prev holds the last unconverged interval previously examined 331 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT 332 END IF 333 I = NEXT 334 GO TO 100 335 END IF 336 PREV = I 337* 338* Perform one bisection step 339* 340 CNT = 0 341 S = MID 342 DPLUS = D( 1 ) - S 343 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 344 DO 90 J = 2, N 345 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 346 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 347 90 CONTINUE 348 IF( CNT.LE.I-1 ) THEN 349 WORK( K-1 ) = MID 350 ELSE 351 WORK( K ) = MID 352 END IF 353 I = NEXT 354 355 100 CONTINUE 356 ITER = ITER + 1 357* do another loop if there are still unconverged intervals 358* However, in the last iteration, all intervals are accepted 359* since this is the best we can do. 360 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80 361* 362* 363* At this point, all the intervals have converged 364 DO 110 I = SAVI1, ILAST 365 K = 2*I 366 II = I - OFFSET 367* All intervals marked by '0' have been refined. 368 IF( IWORK( K-1 ).EQ.0 ) THEN 369 W( II ) = HALF*( WORK( K-1 )+WORK( K ) ) 370 WERR( II ) = WORK( K ) - W( II ) 371 END IF 372 110 CONTINUE 373* 374 375 RETURN 376* 377* End of DLARRJ 378* 379 END 380