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*> \date June 2016 178* 179*> \ingroup OTHERauxiliary 180* 181*> \par Contributors: 182* ================== 183*> 184*> Beresford Parlett, University of California, Berkeley, USA \n 185*> Jim Demmel, University of California, Berkeley, USA \n 186*> Inderjit Dhillon, University of Texas, Austin, USA \n 187*> Osni Marques, LBNL/NERSC, USA \n 188*> Christof Voemel, University of California, Berkeley, USA 189* 190* ===================================================================== 191 SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND, 192 $ W, WGAP, WERR, 193 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, 194 $ DPLUS, LPLUS, WORK, INFO ) 195* 196* -- LAPACK auxiliary routine (version 3.7.1) -- 197* -- LAPACK is a software package provided by Univ. of Tennessee, -- 198* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 199* June 2016 200* 201* .. Scalar Arguments .. 202 INTEGER CLSTRT, CLEND, INFO, N 203 DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 204* .. 205* .. Array Arguments .. 206 DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ), 207 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 208* .. 209* 210* ===================================================================== 211* 212* .. Parameters .. 213 DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO 214 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, FOUR = 4.0D0, 215 $ QUART = 0.25D0, 216 $ MAXGROWTH1 = 8.D0, 217 $ MAXGROWTH2 = 8.D0 ) 218* .. 219* .. Local Scalars .. 220 LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1 221 INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT 222 PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 ) 223 DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL, 224 $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA, 225 $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX, 226 $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2 227* .. 228* .. External Functions .. 229 LOGICAL DISNAN 230 DOUBLE PRECISION DLAMCH 231 EXTERNAL DISNAN, DLAMCH 232* .. 233* .. External Subroutines .. 234 EXTERNAL DCOPY 235* .. 236* .. Intrinsic Functions .. 237 INTRINSIC ABS 238* .. 239* .. Executable Statements .. 240* 241 INFO = 0 242* 243* Quick return if possible 244* 245 IF( N.LE.0 ) THEN 246 RETURN 247 END IF 248* 249 FACT = DBLE(2**KTRYMAX) 250 EPS = DLAMCH( 'Precision' ) 251 SHIFT = 0 252 FORCER = .FALSE. 253 254 255* Note that we cannot guarantee that for any of the shifts tried, 256* the factorization has a small or even moderate element growth. 257* There could be Ritz values at both ends of the cluster and despite 258* backing off, there are examples where all factorizations tried 259* (in IEEE mode, allowing zero pivots & infinities) have INFINITE 260* element growth. 261* For this reason, we should use PIVMIN in this subroutine so that at 262* least the L D L^T factorization exists. It can be checked afterwards 263* whether the element growth caused bad residuals/orthogonality. 264 265* Decide whether the code should accept the best among all 266* representations despite large element growth or signal INFO=1 267* Setting NOFAIL to .FALSE. for quick fix for bug 113 268 NOFAIL = .FALSE. 269* 270 271* Compute the average gap length of the cluster 272 CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT) 273 AVGAP = CLWDTH / DBLE(CLEND-CLSTRT) 274 MINGAP = MIN(CLGAPL, CLGAPR) 275* Initial values for shifts to both ends of cluster 276 LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT ) 277 RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND ) 278 279* Use a small fudge to make sure that we really shift to the outside 280 LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS 281 RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS 282 283* Compute upper bounds for how much to back off the initial shifts 284 LDMAX = QUART * MINGAP + TWO * PIVMIN 285 RDMAX = QUART * MINGAP + TWO * PIVMIN 286 287 LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT 288 RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT 289* 290* Initialize the record of the best representation found 291* 292 S = DLAMCH( 'S' ) 293 SMLGROWTH = ONE / S 294 FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS) 295 FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS)) 296 BESTSHIFT = LSIGMA 297* 298* while (KTRY <= KTRYMAX) 299 KTRY = 0 300 GROWTHBOUND = MAXGROWTH1*SPDIAM 301 302 5 CONTINUE 303 SAWNAN1 = .FALSE. 304 SAWNAN2 = .FALSE. 305* Ensure that we do not back off too much of the initial shifts 306 LDELTA = MIN(LDMAX,LDELTA) 307 RDELTA = MIN(RDMAX,RDELTA) 308 309* Compute the element growth when shifting to both ends of the cluster 310* accept the shift if there is no element growth at one of the two ends 311 312* Left end 313 S = -LSIGMA 314 DPLUS( 1 ) = D( 1 ) + S 315 IF(ABS(DPLUS(1)).LT.PIVMIN) THEN 316 DPLUS(1) = -PIVMIN 317* Need to set SAWNAN1 because refined RRR test should not be used 318* in this case 319 SAWNAN1 = .TRUE. 320 ENDIF 321 MAX1 = ABS( DPLUS( 1 ) ) 322 DO 6 I = 1, N - 1 323 LPLUS( I ) = LD( I ) / DPLUS( I ) 324 S = S*LPLUS( I )*L( I ) - LSIGMA 325 DPLUS( I+1 ) = D( I+1 ) + S 326 IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN 327 DPLUS(I+1) = -PIVMIN 328* Need to set SAWNAN1 because refined RRR test should not be used 329* in this case 330 SAWNAN1 = .TRUE. 331 ENDIF 332 MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) ) 333 6 CONTINUE 334 SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 ) 335 336 IF( FORCER .OR. 337 $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN 338 SIGMA = LSIGMA 339 SHIFT = SLEFT 340 GOTO 100 341 ENDIF 342 343* Right end 344 S = -RSIGMA 345 WORK( 1 ) = D( 1 ) + S 346 IF(ABS(WORK(1)).LT.PIVMIN) THEN 347 WORK(1) = -PIVMIN 348* Need to set SAWNAN2 because refined RRR test should not be used 349* in this case 350 SAWNAN2 = .TRUE. 351 ENDIF 352 MAX2 = ABS( WORK( 1 ) ) 353 DO 7 I = 1, N - 1 354 WORK( N+I ) = LD( I ) / WORK( I ) 355 S = S*WORK( N+I )*L( I ) - RSIGMA 356 WORK( I+1 ) = D( I+1 ) + S 357 IF(ABS(WORK(I+1)).LT.PIVMIN) THEN 358 WORK(I+1) = -PIVMIN 359* Need to set SAWNAN2 because refined RRR test should not be used 360* in this case 361 SAWNAN2 = .TRUE. 362 ENDIF 363 MAX2 = MAX( MAX2,ABS(WORK(I+1)) ) 364 7 CONTINUE 365 SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 ) 366 367 IF( FORCER .OR. 368 $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN 369 SIGMA = RSIGMA 370 SHIFT = SRIGHT 371 GOTO 100 372 ENDIF 373* If we are at this point, both shifts led to too much element growth 374 375* Record the better of the two shifts (provided it didn't lead to NaN) 376 IF(SAWNAN1.AND.SAWNAN2) THEN 377* both MAX1 and MAX2 are NaN 378 GOTO 50 379 ELSE 380 IF( .NOT.SAWNAN1 ) THEN 381 INDX = 1 382 IF(MAX1.LE.SMLGROWTH) THEN 383 SMLGROWTH = MAX1 384 BESTSHIFT = LSIGMA 385 ENDIF 386 ENDIF 387 IF( .NOT.SAWNAN2 ) THEN 388 IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2 389 IF(MAX2.LE.SMLGROWTH) THEN 390 SMLGROWTH = MAX2 391 BESTSHIFT = RSIGMA 392 ENDIF 393 ENDIF 394 ENDIF 395 396* If we are here, both the left and the right shift led to 397* element growth. If the element growth is moderate, then 398* we may still accept the representation, if it passes a 399* refined test for RRR. This test supposes that no NaN occurred. 400* Moreover, we use the refined RRR test only for isolated clusters. 401 IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND. 402 $ (MIN(MAX1,MAX2).LT.FAIL2) 403 $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN 404 DORRR1 = .TRUE. 405 ELSE 406 DORRR1 = .FALSE. 407 ENDIF 408 TRYRRR1 = .TRUE. 409 IF( TRYRRR1 .AND. DORRR1 ) THEN 410 IF(INDX.EQ.1) THEN 411 TMP = ABS( DPLUS( N ) ) 412 ZNM2 = ONE 413 PROD = ONE 414 OLDP = ONE 415 DO 15 I = N-1, 1, -1 416 IF( PROD .LE. EPS ) THEN 417 PROD = 418 $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP 419 ELSE 420 PROD = PROD*ABS(WORK(N+I)) 421 END IF 422 OLDP = PROD 423 ZNM2 = ZNM2 + PROD**2 424 TMP = MAX( TMP, ABS( DPLUS( I ) * PROD )) 425 15 CONTINUE 426 RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 427 IF (RRR1.LE.MAXGROWTH2) THEN 428 SIGMA = LSIGMA 429 SHIFT = SLEFT 430 GOTO 100 431 ENDIF 432 ELSE IF(INDX.EQ.2) THEN 433 TMP = ABS( WORK( N ) ) 434 ZNM2 = ONE 435 PROD = ONE 436 OLDP = ONE 437 DO 16 I = N-1, 1, -1 438 IF( PROD .LE. EPS ) THEN 439 PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP 440 ELSE 441 PROD = PROD*ABS(LPLUS(I)) 442 END IF 443 OLDP = PROD 444 ZNM2 = ZNM2 + PROD**2 445 TMP = MAX( TMP, ABS( WORK( I ) * PROD )) 446 16 CONTINUE 447 RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 448 IF (RRR2.LE.MAXGROWTH2) THEN 449 SIGMA = RSIGMA 450 SHIFT = SRIGHT 451 GOTO 100 452 ENDIF 453 END IF 454 ENDIF 455 456 50 CONTINUE 457 458 IF (KTRY.LT.KTRYMAX) THEN 459* If we are here, both shifts failed also the RRR test. 460* Back off to the outside 461 LSIGMA = MAX( LSIGMA - LDELTA, 462 $ LSIGMA - LDMAX) 463 RSIGMA = MIN( RSIGMA + RDELTA, 464 $ RSIGMA + RDMAX ) 465 LDELTA = TWO * LDELTA 466 RDELTA = TWO * RDELTA 467 KTRY = KTRY + 1 468 GOTO 5 469 ELSE 470* None of the representations investigated satisfied our 471* criteria. Take the best one we found. 472 IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN 473 LSIGMA = BESTSHIFT 474 RSIGMA = BESTSHIFT 475 FORCER = .TRUE. 476 GOTO 5 477 ELSE 478 INFO = 1 479 RETURN 480 ENDIF 481 END IF 482 483 100 CONTINUE 484 IF (SHIFT.EQ.SLEFT) THEN 485 ELSEIF (SHIFT.EQ.SRIGHT) THEN 486* store new L and D back into DPLUS, LPLUS 487 CALL DCOPY( N, WORK, 1, DPLUS, 1 ) 488 CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 ) 489 ENDIF 490 491 RETURN 492* 493* End of DLARRF 494* 495 END 496