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