diff options
Diffstat (limited to 'SRC/slarrf.f')
| -rw-r--r-- | SRC/slarrf.f | 373 |
1 files changed, 373 insertions, 0 deletions
diff --git a/SRC/slarrf.f b/SRC/slarrf.f new file mode 100644 index 00000000..529e4e70 --- /dev/null +++ b/SRC/slarrf.f @@ -0,0 +1,373 @@ + SUBROUTINE SLARRF( N, D, L, LD, CLSTRT, CLEND, + $ W, WGAP, WERR, + $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, + $ DPLUS, LPLUS, WORK, INFO ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +** +* .. Scalar Arguments .. + INTEGER CLSTRT, CLEND, INFO, N + REAL CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM +* .. +* .. Array Arguments .. + REAL D( * ), DPLUS( * ), L( * ), LD( * ), + $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* Given the initial representation L D L^T and its cluster of close +* eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... +* W( CLEND ), SLARRF finds a new relatively robust representation +* L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the +* eigenvalues of L(+) D(+) L(+)^T is relatively isolated. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix (subblock, if the matrix splitted). +* +* D (input) REAL array, dimension (N) +* The N diagonal elements of the diagonal matrix D. +* +* L (input) REAL array, dimension (N-1) +* The (N-1) subdiagonal elements of the unit bidiagonal +* matrix L. +* +* LD (input) REAL array, dimension (N-1) +* The (N-1) elements L(i)*D(i). +* +* CLSTRT (input) INTEGER +* The index of the first eigenvalue in the cluster. +* +* CLEND (input) INTEGER +* The index of the last eigenvalue in the cluster. +* +* W (input) REAL array, dimension >= (CLEND-CLSTRT+1) +* The eigenvalue APPROXIMATIONS of L D L^T in ascending order. +* W( CLSTRT ) through W( CLEND ) form the cluster of relatively +* close eigenalues. +* +* WGAP (input/output) REAL array, dimension >= (CLEND-CLSTRT+1) +* The separation from the right neighbor eigenvalue in W. +* +* WERR (input) REAL array, dimension >= (CLEND-CLSTRT+1) +* WERR contain the semiwidth of the uncertainty +* interval of the corresponding eigenvalue APPROXIMATION in W +* +* SPDIAM (input) estimate of the spectral diameter obtained from the +* Gerschgorin intervals +* +* CLGAPL, CLGAPR (input) absolute gap on each end of the cluster. +* Set by the calling routine to protect against shifts too close +* to eigenvalues outside the cluster. +* +* PIVMIN (input) DOUBLE PRECISION +* The minimum pivot allowed in the Sturm sequence. +* +* SIGMA (output) REAL +* The shift used to form L(+) D(+) L(+)^T. +* +* DPLUS (output) REAL array, dimension (N) +* The N diagonal elements of the diagonal matrix D(+). +* +* LPLUS (output) REAL array, dimension (N-1) +* The first (N-1) elements of LPLUS contain the subdiagonal +* elements of the unit bidiagonal matrix L(+). +* +* WORK (workspace) REAL array, dimension (2*N) +* Workspace. +* +* Further Details +* =============== +* +* Based on contributions by +* Beresford Parlett, University of California, Berkeley, USA +* Jim Demmel, University of California, Berkeley, USA +* Inderjit Dhillon, University of Texas, Austin, USA +* Osni Marques, LBNL/NERSC, USA +* Christof Voemel, University of California, Berkeley, USA +* +* ===================================================================== +* +* .. Parameters .. + REAL FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO, + $ ZERO + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, + $ FOUR = 4.0E0, QUART = 0.25E0, + $ MAXGROWTH1 = 8.E0, + $ MAXGROWTH2 = 8.E0 ) +* .. +* .. Local Scalars .. + LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1 + INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT + PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 ) + REAL AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL, + $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA, + $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX, + $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2 +* .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL SCOPY +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* + INFO = 0 + FACT = REAL(2**KTRYMAX) + EPS = SLAMCH( 'Precision' ) + SHIFT = 0 + FORCER = .FALSE. + + +* Note that we cannot guarantee that for any of the shifts tried, +* the factorization has a small or even moderate element growth. +* There could be Ritz values at both ends of the cluster and despite +* backing off, there are examples where all factorizations tried +* (in IEEE mode, allowing zero pivots & infinities) have INFINITE +* element growth. +* For this reason, we should use PIVMIN in this subroutine so that at +* least the L D L^T factorization exists. It can be checked afterwards +* whether the element growth caused bad residuals/orthogonality. + +* Decide whether the code should accept the best among all +* representations despite large element growth or signal INFO=1 + NOFAIL = .TRUE. +* + +* Compute the average gap length of the cluster + CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT) + AVGAP = CLWDTH / REAL(CLEND-CLSTRT) + MINGAP = MIN(CLGAPL, CLGAPR) +* Initial values for shifts to both ends of cluster + LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT ) + RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND ) + +* Use a small fudge to make sure that we really shift to the outside + LSIGMA = LSIGMA - ABS(LSIGMA)* TWO * EPS + RSIGMA = RSIGMA + ABS(RSIGMA)* TWO * EPS + +* Compute upper bounds for how much to back off the initial shifts + LDMAX = QUART * MINGAP + TWO * PIVMIN + RDMAX = QUART * MINGAP + TWO * PIVMIN + + LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT + RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT +* +* Initialize the record of the best representation found +* + S = SLAMCH( 'S' ) + SMLGROWTH = ONE / S + FAIL = REAL(N-1)*MINGAP/(SPDIAM*EPS) + FAIL2 = REAL(N-1)*MINGAP/(SPDIAM*SQRT(EPS)) + BESTSHIFT = LSIGMA +* +* while (KTRY <= KTRYMAX) + KTRY = 0 + GROWTHBOUND = MAXGROWTH1*SPDIAM + + 5 CONTINUE + SAWNAN1 = .FALSE. + SAWNAN2 = .FALSE. +* Ensure that we do not back off too much of the initial shifts + LDELTA = MIN(LDMAX,LDELTA) + RDELTA = MIN(RDMAX,RDELTA) + +* Compute the element growth when shifting to both ends of the cluster +* accept the shift if there is no element growth at one of the two ends + +* Left end + S = -LSIGMA + DPLUS( 1 ) = D( 1 ) + S + IF(ABS(DPLUS(1)).LT.PIVMIN) THEN + DPLUS(1) = -PIVMIN +* Need to set SAWNAN1 because refined RRR test should not be used +* in this case + SAWNAN1 = .TRUE. + ENDIF + MAX1 = ABS( DPLUS( 1 ) ) + DO 6 I = 1, N - 1 + LPLUS( I ) = LD( I ) / DPLUS( I ) + S = S*LPLUS( I )*L( I ) - LSIGMA + DPLUS( I+1 ) = D( I+1 ) + S + IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN + DPLUS(I+1) = -PIVMIN +* Need to set SAWNAN1 because refined RRR test should not be used +* in this case + SAWNAN1 = .TRUE. + ENDIF + MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) ) + 6 CONTINUE + SAWNAN1 = SAWNAN1 .OR. SISNAN( MAX1 ) + + IF( FORCER .OR. + $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN + SIGMA = LSIGMA + SHIFT = SLEFT + GOTO 100 + ENDIF + +* Right end + S = -RSIGMA + WORK( 1 ) = D( 1 ) + S + IF(ABS(WORK(1)).LT.PIVMIN) THEN + WORK(1) = -PIVMIN +* Need to set SAWNAN2 because refined RRR test should not be used +* in this case + SAWNAN2 = .TRUE. + ENDIF + MAX2 = ABS( WORK( 1 ) ) + DO 7 I = 1, N - 1 + WORK( N+I ) = LD( I ) / WORK( I ) + S = S*WORK( N+I )*L( I ) - RSIGMA + WORK( I+1 ) = D( I+1 ) + S + IF(ABS(WORK(I+1)).LT.PIVMIN) THEN + WORK(I+1) = -PIVMIN +* Need to set SAWNAN2 because refined RRR test should not be used +* in this case + SAWNAN2 = .TRUE. + ENDIF + MAX2 = MAX( MAX2,ABS(WORK(I+1)) ) + 7 CONTINUE + SAWNAN2 = SAWNAN2 .OR. SISNAN( MAX2 ) + + IF( FORCER .OR. + $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN + SIGMA = RSIGMA + SHIFT = SRIGHT + GOTO 100 + ENDIF +* If we are at this point, both shifts led to too much element growth + +* Record the better of the two shifts (provided it didn't lead to NaN) + IF(SAWNAN1.AND.SAWNAN2) THEN +* both MAX1 and MAX2 are NaN + GOTO 50 + ELSE + IF( .NOT.SAWNAN1 ) THEN + INDX = 1 + IF(MAX1.LE.SMLGROWTH) THEN + SMLGROWTH = MAX1 + BESTSHIFT = LSIGMA + ENDIF + ENDIF + IF( .NOT.SAWNAN2 ) THEN + IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2 + IF(MAX2.LE.SMLGROWTH) THEN + SMLGROWTH = MAX2 + BESTSHIFT = RSIGMA + ENDIF + ENDIF + ENDIF + +* If we are here, both the left and the right shift led to +* element growth. If the element growth is moderate, then +* we may still accept the representation, if it passes a +* refined test for RRR. This test supposes that no NaN occurred. +* Moreover, we use the refined RRR test only for isolated clusters. + IF((CLWDTH.LT.MINGAP/REAL(128)) .AND. + $ (MIN(MAX1,MAX2).LT.FAIL2) + $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN + DORRR1 = .TRUE. + ELSE + DORRR1 = .FALSE. + ENDIF + TRYRRR1 = .TRUE. + IF( TRYRRR1 .AND. DORRR1 ) THEN + IF(INDX.EQ.1) THEN + TMP = ABS( DPLUS( N ) ) + ZNM2 = ONE + PROD = ONE + OLDP = ONE + DO 15 I = N-1, 1, -1 + IF( PROD .LE. EPS ) THEN + PROD = + $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP + ELSE + PROD = PROD*ABS(WORK(N+I)) + END IF + OLDP = PROD + ZNM2 = ZNM2 + PROD**2 + TMP = MAX( TMP, ABS( DPLUS( I ) * PROD )) + 15 CONTINUE + RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) ) + IF (RRR1.LE.MAXGROWTH2) THEN + SIGMA = LSIGMA + SHIFT = SLEFT + GOTO 100 + ENDIF + ELSE IF(INDX.EQ.2) THEN + TMP = ABS( WORK( N ) ) + ZNM2 = ONE + PROD = ONE + OLDP = ONE + DO 16 I = N-1, 1, -1 + IF( PROD .LE. EPS ) THEN + PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP + ELSE + PROD = PROD*ABS(LPLUS(I)) + END IF + OLDP = PROD + ZNM2 = ZNM2 + PROD**2 + TMP = MAX( TMP, ABS( WORK( I ) * PROD )) + 16 CONTINUE + RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) ) + IF (RRR2.LE.MAXGROWTH2) THEN + SIGMA = RSIGMA + SHIFT = SRIGHT + GOTO 100 + ENDIF + END IF + ENDIF + + 50 CONTINUE + + IF (KTRY.LT.KTRYMAX) THEN +* If we are here, both shifts failed also the RRR test. +* Back off to the outside + LSIGMA = MAX( LSIGMA - LDELTA, + $ LSIGMA - LDMAX) + RSIGMA = MIN( RSIGMA + RDELTA, + $ RSIGMA + RDMAX ) + LDELTA = TWO * LDELTA + RDELTA = TWO * RDELTA + KTRY = KTRY + 1 + GOTO 5 + ELSE +* None of the representations investigated satisfied our +* criteria. Take the best one we found. + IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN + LSIGMA = BESTSHIFT + RSIGMA = BESTSHIFT + FORCER = .TRUE. + GOTO 5 + ELSE + INFO = 1 + RETURN + ENDIF + END IF + + 100 CONTINUE + IF (SHIFT.EQ.SLEFT) THEN + ELSEIF (SHIFT.EQ.SRIGHT) THEN +* store new L and D back into DPLUS, LPLUS + CALL SCOPY( N, WORK, 1, DPLUS, 1 ) + CALL SCOPY( N-1, WORK(N+1), 1, LPLUS, 1 ) + ENDIF + + RETURN +* +* End of SLARRF +* + END |