1*> \brief \b SLARRF 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 SLARRF + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrf.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrf.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrf.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLARRF( 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* REAL CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 29* .. 30* .. Array Arguments .. 31* REAL 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 ), SLARRF 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 115*> estimate of the spectral diameter obtained from the 116*> Gerschgorin intervals 117*> \endverbatim 118*> 119*> \param[in] CLGAPL 120*> \verbatim 121*> CLGAPL is REAL 122*> \endverbatim 123*> 124*> \param[in] CLGAPR 125*> \verbatim 126*> CLGAPR is REAL 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 REAL 135*> The minimum pivot allowed in the Sturm sequence. 136*> \endverbatim 137*> 138*> \param[out] SIGMA 139*> \verbatim 140*> SIGMA is REAL 141*> The shift used to form L(+) D(+) L(+)^T. 142*> \endverbatim 143*> 144*> \param[out] DPLUS 145*> \verbatim 146*> DPLUS is REAL 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 REAL 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 REAL 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 November 2015 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 SLARRF( 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.6.0) -- 197* -- LAPACK is a software package provided by Univ. of Tennessee, -- 198* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 199* November 2015 200* 201* .. Scalar Arguments .. 202 INTEGER CLSTRT, CLEND, INFO, N 203 REAL CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 204* .. 205* .. Array Arguments .. 206 REAL D( * ), DPLUS( * ), L( * ), LD( * ), 207 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 208* .. 209* 210* ===================================================================== 211* 212* .. Parameters .. 213 REAL MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO 214 PARAMETER ( ONE = 1.0E0, TWO = 2.0E0, 215 $ QUART = 0.25E0, 216 $ MAXGROWTH1 = 8.E0, 217 $ MAXGROWTH2 = 8.E0 ) 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 REAL 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 SISNAN 230 REAL SLAMCH 231 EXTERNAL SISNAN, SLAMCH 232* .. 233* .. External Subroutines .. 234 EXTERNAL SCOPY 235* .. 236* .. Intrinsic Functions .. 237 INTRINSIC ABS 238* .. 239* .. Executable Statements .. 240* 241 INFO = 0 242 FACT = REAL(2**KTRYMAX) 243 EPS = SLAMCH( '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* Setting NOFAIL to .FALSE. for quick fix for bug 113 261 NOFAIL = .FALSE. 262* 263 264* Compute the average gap length of the cluster 265 CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT) 266 AVGAP = CLWDTH / REAL(CLEND-CLSTRT) 267 MINGAP = MIN(CLGAPL, CLGAPR) 268* Initial values for shifts to both ends of cluster 269 LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT ) 270 RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND ) 271 272* Use a small fudge to make sure that we really shift to the outside 273 LSIGMA = LSIGMA - ABS(LSIGMA)* TWO * EPS 274 RSIGMA = RSIGMA + ABS(RSIGMA)* TWO * EPS 275 276* Compute upper bounds for how much to back off the initial shifts 277 LDMAX = QUART * MINGAP + TWO * PIVMIN 278 RDMAX = QUART * MINGAP + TWO * PIVMIN 279 280 LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT 281 RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT 282* 283* Initialize the record of the best representation found 284* 285 S = SLAMCH( 'S' ) 286 SMLGROWTH = ONE / S 287 FAIL = REAL(N-1)*MINGAP/(SPDIAM*EPS) 288 FAIL2 = REAL(N-1)*MINGAP/(SPDIAM*SQRT(EPS)) 289 BESTSHIFT = LSIGMA 290* 291* while (KTRY <= KTRYMAX) 292 KTRY = 0 293 GROWTHBOUND = MAXGROWTH1*SPDIAM 294 295 5 CONTINUE 296 SAWNAN1 = .FALSE. 297 SAWNAN2 = .FALSE. 298* Ensure that we do not back off too much of the initial shifts 299 LDELTA = MIN(LDMAX,LDELTA) 300 RDELTA = MIN(RDMAX,RDELTA) 301 302* Compute the element growth when shifting to both ends of the cluster 303* accept the shift if there is no element growth at one of the two ends 304 305* Left end 306 S = -LSIGMA 307 DPLUS( 1 ) = D( 1 ) + S 308 IF(ABS(DPLUS(1)).LT.PIVMIN) THEN 309 DPLUS(1) = -PIVMIN 310* Need to set SAWNAN1 because refined RRR test should not be used 311* in this case 312 SAWNAN1 = .TRUE. 313 ENDIF 314 MAX1 = ABS( DPLUS( 1 ) ) 315 DO 6 I = 1, N - 1 316 LPLUS( I ) = LD( I ) / DPLUS( I ) 317 S = S*LPLUS( I )*L( I ) - LSIGMA 318 DPLUS( I+1 ) = D( I+1 ) + S 319 IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN 320 DPLUS(I+1) = -PIVMIN 321* Need to set SAWNAN1 because refined RRR test should not be used 322* in this case 323 SAWNAN1 = .TRUE. 324 ENDIF 325 MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) ) 326 6 CONTINUE 327 SAWNAN1 = SAWNAN1 .OR. SISNAN( MAX1 ) 328 329 IF( FORCER .OR. 330 $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN 331 SIGMA = LSIGMA 332 SHIFT = SLEFT 333 GOTO 100 334 ENDIF 335 336* Right end 337 S = -RSIGMA 338 WORK( 1 ) = D( 1 ) + S 339 IF(ABS(WORK(1)).LT.PIVMIN) THEN 340 WORK(1) = -PIVMIN 341* Need to set SAWNAN2 because refined RRR test should not be used 342* in this case 343 SAWNAN2 = .TRUE. 344 ENDIF 345 MAX2 = ABS( WORK( 1 ) ) 346 DO 7 I = 1, N - 1 347 WORK( N+I ) = LD( I ) / WORK( I ) 348 S = S*WORK( N+I )*L( I ) - RSIGMA 349 WORK( I+1 ) = D( I+1 ) + S 350 IF(ABS(WORK(I+1)).LT.PIVMIN) THEN 351 WORK(I+1) = -PIVMIN 352* Need to set SAWNAN2 because refined RRR test should not be used 353* in this case 354 SAWNAN2 = .TRUE. 355 ENDIF 356 MAX2 = MAX( MAX2,ABS(WORK(I+1)) ) 357 7 CONTINUE 358 SAWNAN2 = SAWNAN2 .OR. SISNAN( MAX2 ) 359 360 IF( FORCER .OR. 361 $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN 362 SIGMA = RSIGMA 363 SHIFT = SRIGHT 364 GOTO 100 365 ENDIF 366* If we are at this point, both shifts led to too much element growth 367 368* Record the better of the two shifts (provided it didn't lead to NaN) 369 IF(SAWNAN1.AND.SAWNAN2) THEN 370* both MAX1 and MAX2 are NaN 371 GOTO 50 372 ELSE 373 IF( .NOT.SAWNAN1 ) THEN 374 INDX = 1 375 IF(MAX1.LE.SMLGROWTH) THEN 376 SMLGROWTH = MAX1 377 BESTSHIFT = LSIGMA 378 ENDIF 379 ENDIF 380 IF( .NOT.SAWNAN2 ) THEN 381 IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2 382 IF(MAX2.LE.SMLGROWTH) THEN 383 SMLGROWTH = MAX2 384 BESTSHIFT = RSIGMA 385 ENDIF 386 ENDIF 387 ENDIF 388 389* If we are here, both the left and the right shift led to 390* element growth. If the element growth is moderate, then 391* we may still accept the representation, if it passes a 392* refined test for RRR. This test supposes that no NaN occurred. 393* Moreover, we use the refined RRR test only for isolated clusters. 394 IF((CLWDTH.LT.MINGAP/REAL(128)) .AND. 395 $ (MIN(MAX1,MAX2).LT.FAIL2) 396 $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN 397 DORRR1 = .TRUE. 398 ELSE 399 DORRR1 = .FALSE. 400 ENDIF 401 TRYRRR1 = .TRUE. 402 IF( TRYRRR1 .AND. DORRR1 ) THEN 403 IF(INDX.EQ.1) THEN 404 TMP = ABS( DPLUS( N ) ) 405 ZNM2 = ONE 406 PROD = ONE 407 OLDP = ONE 408 DO 15 I = N-1, 1, -1 409 IF( PROD .LE. EPS ) THEN 410 PROD = 411 $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP 412 ELSE 413 PROD = PROD*ABS(WORK(N+I)) 414 END IF 415 OLDP = PROD 416 ZNM2 = ZNM2 + PROD**2 417 TMP = MAX( TMP, ABS( DPLUS( I ) * PROD )) 418 15 CONTINUE 419 RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 420 IF (RRR1.LE.MAXGROWTH2) THEN 421 SIGMA = LSIGMA 422 SHIFT = SLEFT 423 GOTO 100 424 ENDIF 425 ELSE IF(INDX.EQ.2) THEN 426 TMP = ABS( WORK( N ) ) 427 ZNM2 = ONE 428 PROD = ONE 429 OLDP = ONE 430 DO 16 I = N-1, 1, -1 431 IF( PROD .LE. EPS ) THEN 432 PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP 433 ELSE 434 PROD = PROD*ABS(LPLUS(I)) 435 END IF 436 OLDP = PROD 437 ZNM2 = ZNM2 + PROD**2 438 TMP = MAX( TMP, ABS( WORK( I ) * PROD )) 439 16 CONTINUE 440 RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 441 IF (RRR2.LE.MAXGROWTH2) THEN 442 SIGMA = RSIGMA 443 SHIFT = SRIGHT 444 GOTO 100 445 ENDIF 446 END IF 447 ENDIF 448 449 50 CONTINUE 450 451 IF (KTRY.LT.KTRYMAX) THEN 452* If we are here, both shifts failed also the RRR test. 453* Back off to the outside 454 LSIGMA = MAX( LSIGMA - LDELTA, 455 $ LSIGMA - LDMAX) 456 RSIGMA = MIN( RSIGMA + RDELTA, 457 $ RSIGMA + RDMAX ) 458 LDELTA = TWO * LDELTA 459 RDELTA = TWO * RDELTA 460 KTRY = KTRY + 1 461 GOTO 5 462 ELSE 463* None of the representations investigated satisfied our 464* criteria. Take the best one we found. 465 IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN 466 LSIGMA = BESTSHIFT 467 RSIGMA = BESTSHIFT 468 FORCER = .TRUE. 469 GOTO 5 470 ELSE 471 INFO = 1 472 RETURN 473 ENDIF 474 END IF 475 476 100 CONTINUE 477 IF (SHIFT.EQ.SLEFT) THEN 478 ELSEIF (SHIFT.EQ.SRIGHT) THEN 479* store new L and D back into DPLUS, LPLUS 480 CALL SCOPY( N, WORK, 1, DPLUS, 1 ) 481 CALL SCOPY( N-1, WORK(N+1), 1, LPLUS, 1 ) 482 ENDIF 483 484 RETURN 485* 486* End of SLARRF 487* 488 END 489