1*> \brief \b DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARRF + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrf.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrf.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrf.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND, 22* W, WGAP, WERR, 23* SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, 24* DPLUS, LPLUS, WORK, INFO ) 25* 26* .. Scalar Arguments .. 27* INTEGER CLSTRT, CLEND, INFO, N 28* DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 29* .. 30* .. Array Arguments .. 31* DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ), 32* $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> Given the initial representation L D L^T and its cluster of close 42*> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... 43*> W( CLEND ), DLARRF finds a new relatively robust representation 44*> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the 45*> eigenvalues of L(+) D(+) L(+)^T is relatively isolated. 46*> \endverbatim 47* 48* Arguments: 49* ========== 50* 51*> \param[in] N 52*> \verbatim 53*> N is INTEGER 54*> The order of the matrix (subblock, if the matrix split). 55*> \endverbatim 56*> 57*> \param[in] D 58*> \verbatim 59*> D is DOUBLE PRECISION array, dimension (N) 60*> The N diagonal elements of the diagonal matrix D. 61*> \endverbatim 62*> 63*> \param[in] L 64*> \verbatim 65*> L is DOUBLE PRECISION array, dimension (N-1) 66*> The (N-1) subdiagonal elements of the unit bidiagonal 67*> matrix L. 68*> \endverbatim 69*> 70*> \param[in] LD 71*> \verbatim 72*> LD is DOUBLE PRECISION array, dimension (N-1) 73*> The (N-1) elements L(i)*D(i). 74*> \endverbatim 75*> 76*> \param[in] CLSTRT 77*> \verbatim 78*> CLSTRT is INTEGER 79*> The index of the first eigenvalue in the cluster. 80*> \endverbatim 81*> 82*> \param[in] CLEND 83*> \verbatim 84*> CLEND is INTEGER 85*> The index of the last eigenvalue in the cluster. 86*> \endverbatim 87*> 88*> \param[in] W 89*> \verbatim 90*> W is DOUBLE PRECISION array, dimension 91*> dimension is >= (CLEND-CLSTRT+1) 92*> The eigenvalue APPROXIMATIONS of L D L^T in ascending order. 93*> W( CLSTRT ) through W( CLEND ) form the cluster of relatively 94*> close eigenalues. 95*> \endverbatim 96*> 97*> \param[in,out] WGAP 98*> \verbatim 99*> WGAP is DOUBLE PRECISION array, dimension 100*> dimension is >= (CLEND-CLSTRT+1) 101*> The separation from the right neighbor eigenvalue in W. 102*> \endverbatim 103*> 104*> \param[in] WERR 105*> \verbatim 106*> WERR is DOUBLE PRECISION array, dimension 107*> dimension is >= (CLEND-CLSTRT+1) 108*> WERR contain the semiwidth of the uncertainty 109*> interval of the corresponding eigenvalue APPROXIMATION in W 110*> \endverbatim 111*> 112*> \param[in] SPDIAM 113*> \verbatim 114*> SPDIAM is DOUBLE PRECISION 115*> estimate of the spectral diameter obtained from the 116*> Gerschgorin intervals 117*> \endverbatim 118*> 119*> \param[in] CLGAPL 120*> \verbatim 121*> CLGAPL is DOUBLE PRECISION 122*> \endverbatim 123*> 124*> \param[in] CLGAPR 125*> \verbatim 126*> CLGAPR is DOUBLE PRECISION 127*> absolute gap on each end of the cluster. 128*> Set by the calling routine to protect against shifts too close 129*> to eigenvalues outside the cluster. 130*> \endverbatim 131*> 132*> \param[in] PIVMIN 133*> \verbatim 134*> PIVMIN is DOUBLE PRECISION 135*> The minimum pivot allowed in the Sturm sequence. 136*> \endverbatim 137*> 138*> \param[out] SIGMA 139*> \verbatim 140*> SIGMA is DOUBLE PRECISION 141*> The shift used to form L(+) D(+) L(+)^T. 142*> \endverbatim 143*> 144*> \param[out] DPLUS 145*> \verbatim 146*> DPLUS is DOUBLE PRECISION array, dimension (N) 147*> The N diagonal elements of the diagonal matrix D(+). 148*> \endverbatim 149*> 150*> \param[out] LPLUS 151*> \verbatim 152*> LPLUS is DOUBLE PRECISION array, dimension (N-1) 153*> The first (N-1) elements of LPLUS contain the subdiagonal 154*> elements of the unit bidiagonal matrix L(+). 155*> \endverbatim 156*> 157*> \param[out] WORK 158*> \verbatim 159*> WORK is DOUBLE PRECISION array, dimension (2*N) 160*> Workspace. 161*> \endverbatim 162*> 163*> \param[out] INFO 164*> \verbatim 165*> INFO is INTEGER 166*> Signals processing OK (=0) or failure (=1) 167*> \endverbatim 168* 169* Authors: 170* ======== 171* 172*> \author Univ. of Tennessee 173*> \author Univ. of California Berkeley 174*> \author Univ. of Colorado Denver 175*> \author NAG Ltd. 176* 177*> \ingroup OTHERauxiliary 178* 179*> \par Contributors: 180* ================== 181*> 182*> Beresford Parlett, University of California, Berkeley, USA \n 183*> Jim Demmel, University of California, Berkeley, USA \n 184*> Inderjit Dhillon, University of Texas, Austin, USA \n 185*> Osni Marques, LBNL/NERSC, USA \n 186*> Christof Voemel, University of California, Berkeley, USA 187* 188* ===================================================================== 189 SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND, 190 $ W, WGAP, WERR, 191 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, 192 $ DPLUS, LPLUS, WORK, INFO ) 193* 194* -- LAPACK auxiliary routine -- 195* -- LAPACK is a software package provided by Univ. of Tennessee, -- 196* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 197* 198* .. Scalar Arguments .. 199 INTEGER CLSTRT, CLEND, INFO, N 200 DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 201* .. 202* .. Array Arguments .. 203 DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ), 204 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 205* .. 206* 207* ===================================================================== 208* 209* .. Parameters .. 210 DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO 211 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, FOUR = 4.0D0, 212 $ QUART = 0.25D0, 213 $ MAXGROWTH1 = 8.D0, 214 $ MAXGROWTH2 = 8.D0 ) 215* .. 216* .. Local Scalars .. 217 LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1 218 INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT 219 PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 ) 220 DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL, 221 $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA, 222 $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX, 223 $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2 224* .. 225* .. External Functions .. 226 LOGICAL DISNAN 227 DOUBLE PRECISION DLAMCH 228 EXTERNAL DISNAN, DLAMCH 229* .. 230* .. External Subroutines .. 231 EXTERNAL DCOPY 232* .. 233* .. Intrinsic Functions .. 234 INTRINSIC ABS 235* .. 236* .. Executable Statements .. 237* 238 INFO = 0 239* 240* Quick return if possible 241* 242 IF( N.LE.0 ) THEN 243 RETURN 244 END IF 245* 246 FACT = DBLE(2**KTRYMAX) 247 EPS = DLAMCH( 'Precision' ) 248 SHIFT = 0 249 FORCER = .FALSE. 250 251 252* Note that we cannot guarantee that for any of the shifts tried, 253* the factorization has a small or even moderate element growth. 254* There could be Ritz values at both ends of the cluster and despite 255* backing off, there are examples where all factorizations tried 256* (in IEEE mode, allowing zero pivots & infinities) have INFINITE 257* element growth. 258* For this reason, we should use PIVMIN in this subroutine so that at 259* least the L D L^T factorization exists. It can be checked afterwards 260* whether the element growth caused bad residuals/orthogonality. 261 262* Decide whether the code should accept the best among all 263* representations despite large element growth or signal INFO=1 264* Setting NOFAIL to .FALSE. for quick fix for bug 113 265 NOFAIL = .FALSE. 266* 267 268* Compute the average gap length of the cluster 269 CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT) 270 AVGAP = CLWDTH / DBLE(CLEND-CLSTRT) 271 MINGAP = MIN(CLGAPL, CLGAPR) 272* Initial values for shifts to both ends of cluster 273 LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT ) 274 RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND ) 275 276* Use a small fudge to make sure that we really shift to the outside 277 LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS 278 RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS 279 280* Compute upper bounds for how much to back off the initial shifts 281 LDMAX = QUART * MINGAP + TWO * PIVMIN 282 RDMAX = QUART * MINGAP + TWO * PIVMIN 283 284 LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT 285 RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT 286* 287* Initialize the record of the best representation found 288* 289 S = DLAMCH( 'S' ) 290 SMLGROWTH = ONE / S 291 FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS) 292 FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS)) 293 BESTSHIFT = LSIGMA 294* 295* while (KTRY <= KTRYMAX) 296 KTRY = 0 297 GROWTHBOUND = MAXGROWTH1*SPDIAM 298 299 5 CONTINUE 300 SAWNAN1 = .FALSE. 301 SAWNAN2 = .FALSE. 302* Ensure that we do not back off too much of the initial shifts 303 LDELTA = MIN(LDMAX,LDELTA) 304 RDELTA = MIN(RDMAX,RDELTA) 305 306* Compute the element growth when shifting to both ends of the cluster 307* accept the shift if there is no element growth at one of the two ends 308 309* Left end 310 S = -LSIGMA 311 DPLUS( 1 ) = D( 1 ) + S 312 IF(ABS(DPLUS(1)).LT.PIVMIN) THEN 313 DPLUS(1) = -PIVMIN 314* Need to set SAWNAN1 because refined RRR test should not be used 315* in this case 316 SAWNAN1 = .TRUE. 317 ENDIF 318 MAX1 = ABS( DPLUS( 1 ) ) 319 DO 6 I = 1, N - 1 320 LPLUS( I ) = LD( I ) / DPLUS( I ) 321 S = S*LPLUS( I )*L( I ) - LSIGMA 322 DPLUS( I+1 ) = D( I+1 ) + S 323 IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN 324 DPLUS(I+1) = -PIVMIN 325* Need to set SAWNAN1 because refined RRR test should not be used 326* in this case 327 SAWNAN1 = .TRUE. 328 ENDIF 329 MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) ) 330 6 CONTINUE 331 SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 ) 332 333 IF( FORCER .OR. 334 $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN 335 SIGMA = LSIGMA 336 SHIFT = SLEFT 337 GOTO 100 338 ENDIF 339 340* Right end 341 S = -RSIGMA 342 WORK( 1 ) = D( 1 ) + S 343 IF(ABS(WORK(1)).LT.PIVMIN) THEN 344 WORK(1) = -PIVMIN 345* Need to set SAWNAN2 because refined RRR test should not be used 346* in this case 347 SAWNAN2 = .TRUE. 348 ENDIF 349 MAX2 = ABS( WORK( 1 ) ) 350 DO 7 I = 1, N - 1 351 WORK( N+I ) = LD( I ) / WORK( I ) 352 S = S*WORK( N+I )*L( I ) - RSIGMA 353 WORK( I+1 ) = D( I+1 ) + S 354 IF(ABS(WORK(I+1)).LT.PIVMIN) THEN 355 WORK(I+1) = -PIVMIN 356* Need to set SAWNAN2 because refined RRR test should not be used 357* in this case 358 SAWNAN2 = .TRUE. 359 ENDIF 360 MAX2 = MAX( MAX2,ABS(WORK(I+1)) ) 361 7 CONTINUE 362 SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 ) 363 364 IF( FORCER .OR. 365 $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN 366 SIGMA = RSIGMA 367 SHIFT = SRIGHT 368 GOTO 100 369 ENDIF 370* If we are at this point, both shifts led to too much element growth 371 372* Record the better of the two shifts (provided it didn't lead to NaN) 373 IF(SAWNAN1.AND.SAWNAN2) THEN 374* both MAX1 and MAX2 are NaN 375 GOTO 50 376 ELSE 377 IF( .NOT.SAWNAN1 ) THEN 378 INDX = 1 379 IF(MAX1.LE.SMLGROWTH) THEN 380 SMLGROWTH = MAX1 381 BESTSHIFT = LSIGMA 382 ENDIF 383 ENDIF 384 IF( .NOT.SAWNAN2 ) THEN 385 IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2 386 IF(MAX2.LE.SMLGROWTH) THEN 387 SMLGROWTH = MAX2 388 BESTSHIFT = RSIGMA 389 ENDIF 390 ENDIF 391 ENDIF 392 393* If we are here, both the left and the right shift led to 394* element growth. If the element growth is moderate, then 395* we may still accept the representation, if it passes a 396* refined test for RRR. This test supposes that no NaN occurred. 397* Moreover, we use the refined RRR test only for isolated clusters. 398 IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND. 399 $ (MIN(MAX1,MAX2).LT.FAIL2) 400 $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN 401 DORRR1 = .TRUE. 402 ELSE 403 DORRR1 = .FALSE. 404 ENDIF 405 TRYRRR1 = .TRUE. 406 IF( TRYRRR1 .AND. DORRR1 ) THEN 407 IF(INDX.EQ.1) THEN 408 TMP = ABS( DPLUS( N ) ) 409 ZNM2 = ONE 410 PROD = ONE 411 OLDP = ONE 412 DO 15 I = N-1, 1, -1 413 IF( PROD .LE. EPS ) THEN 414 PROD = 415 $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP 416 ELSE 417 PROD = PROD*ABS(WORK(N+I)) 418 END IF 419 OLDP = PROD 420 ZNM2 = ZNM2 + PROD**2 421 TMP = MAX( TMP, ABS( DPLUS( I ) * PROD )) 422 15 CONTINUE 423 RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 424 IF (RRR1.LE.MAXGROWTH2) THEN 425 SIGMA = LSIGMA 426 SHIFT = SLEFT 427 GOTO 100 428 ENDIF 429 ELSE IF(INDX.EQ.2) THEN 430 TMP = ABS( WORK( N ) ) 431 ZNM2 = ONE 432 PROD = ONE 433 OLDP = ONE 434 DO 16 I = N-1, 1, -1 435 IF( PROD .LE. EPS ) THEN 436 PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP 437 ELSE 438 PROD = PROD*ABS(LPLUS(I)) 439 END IF 440 OLDP = PROD 441 ZNM2 = ZNM2 + PROD**2 442 TMP = MAX( TMP, ABS( WORK( I ) * PROD )) 443 16 CONTINUE 444 RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 445 IF (RRR2.LE.MAXGROWTH2) THEN 446 SIGMA = RSIGMA 447 SHIFT = SRIGHT 448 GOTO 100 449 ENDIF 450 END IF 451 ENDIF 452 453 50 CONTINUE 454 455 IF (KTRY.LT.KTRYMAX) THEN 456* If we are here, both shifts failed also the RRR test. 457* Back off to the outside 458 LSIGMA = MAX( LSIGMA - LDELTA, 459 $ LSIGMA - LDMAX) 460 RSIGMA = MIN( RSIGMA + RDELTA, 461 $ RSIGMA + RDMAX ) 462 LDELTA = TWO * LDELTA 463 RDELTA = TWO * RDELTA 464 KTRY = KTRY + 1 465 GOTO 5 466 ELSE 467* None of the representations investigated satisfied our 468* criteria. Take the best one we found. 469 IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN 470 LSIGMA = BESTSHIFT 471 RSIGMA = BESTSHIFT 472 FORCER = .TRUE. 473 GOTO 5 474 ELSE 475 INFO = 1 476 RETURN 477 ENDIF 478 END IF 479 480 100 CONTINUE 481 IF (SHIFT.EQ.SLEFT) THEN 482 ELSEIF (SHIFT.EQ.SRIGHT) THEN 483* store new L and D back into DPLUS, LPLUS 484 CALL DCOPY( N, WORK, 1, DPLUS, 1 ) 485 CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 ) 486 ENDIF 487 488 RETURN 489* 490* End of DLARRF 491* 492 END 493