|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE SLARRF( N, D, L, LD, CLSTRT, CLEND, 00002 $ W, WGAP, WERR, 00003 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, 00004 $ DPLUS, LPLUS, WORK, INFO ) 00005 * 00006 * -- LAPACK auxiliary routine (version 3.2.2) -- 00007 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00008 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00009 * June 2010 00010 ** 00011 * .. Scalar Arguments .. 00012 INTEGER CLSTRT, CLEND, INFO, N 00013 REAL CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 00014 * .. 00015 * .. Array Arguments .. 00016 REAL D( * ), DPLUS( * ), L( * ), LD( * ), 00017 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * Given the initial representation L D L^T and its cluster of close 00024 * eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... 00025 * W( CLEND ), SLARRF finds a new relatively robust representation 00026 * L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the 00027 * eigenvalues of L(+) D(+) L(+)^T is relatively isolated. 00028 * 00029 * Arguments 00030 * ========= 00031 * 00032 * N (input) INTEGER 00033 * The order of the matrix (subblock, if the matrix splitted). 00034 * 00035 * D (input) REAL array, dimension (N) 00036 * The N diagonal elements of the diagonal matrix D. 00037 * 00038 * L (input) REAL array, dimension (N-1) 00039 * The (N-1) subdiagonal elements of the unit bidiagonal 00040 * matrix L. 00041 * 00042 * LD (input) REAL array, dimension (N-1) 00043 * The (N-1) elements L(i)*D(i). 00044 * 00045 * CLSTRT (input) INTEGER 00046 * The index of the first eigenvalue in the cluster. 00047 * 00048 * CLEND (input) INTEGER 00049 * The index of the last eigenvalue in the cluster. 00050 * 00051 * W (input) REAL array, dimension 00052 * dimension is >= (CLEND-CLSTRT+1) 00053 * The eigenvalue APPROXIMATIONS of L D L^T in ascending order. 00054 * W( CLSTRT ) through W( CLEND ) form the cluster of relatively 00055 * close eigenalues. 00056 * 00057 * WGAP (input/output) REAL array, dimension 00058 * dimension is >= (CLEND-CLSTRT+1) 00059 * The separation from the right neighbor eigenvalue in W. 00060 * 00061 * WERR (input) REAL array, dimension 00062 * dimension is >= (CLEND-CLSTRT+1) 00063 * WERR contain the semiwidth of the uncertainty 00064 * interval of the corresponding eigenvalue APPROXIMATION in W 00065 * 00066 * SPDIAM (input) REAL 00067 * estimate of the spectral diameter obtained from the 00068 * Gerschgorin intervals 00069 * 00070 * CLGAPL (input) REAL 00071 * 00072 * CLGAPR (input) REAL 00073 * absolute gap on each end of the cluster. 00074 * Set by the calling routine to protect against shifts too close 00075 * to eigenvalues outside the cluster. 00076 * 00077 * PIVMIN (input) REAL 00078 * The minimum pivot allowed in the Sturm sequence. 00079 * 00080 * SIGMA (output) REAL 00081 * The shift used to form L(+) D(+) L(+)^T. 00082 * 00083 * DPLUS (output) REAL array, dimension (N) 00084 * The N diagonal elements of the diagonal matrix D(+). 00085 * 00086 * LPLUS (output) REAL array, dimension (N-1) 00087 * The first (N-1) elements of LPLUS contain the subdiagonal 00088 * elements of the unit bidiagonal matrix L(+). 00089 * 00090 * WORK (workspace) REAL array, dimension (2*N) 00091 * Workspace. 00092 * 00093 * Further Details 00094 * =============== 00095 * 00096 * Based on contributions by 00097 * Beresford Parlett, University of California, Berkeley, USA 00098 * Jim Demmel, University of California, Berkeley, USA 00099 * Inderjit Dhillon, University of Texas, Austin, USA 00100 * Osni Marques, LBNL/NERSC, USA 00101 * Christof Voemel, University of California, Berkeley, USA 00102 * 00103 * ===================================================================== 00104 * 00105 * .. Parameters .. 00106 REAL FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO, 00107 $ ZERO 00108 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 00109 $ FOUR = 4.0E0, QUART = 0.25E0, 00110 $ MAXGROWTH1 = 8.E0, 00111 $ MAXGROWTH2 = 8.E0 ) 00112 * .. 00113 * .. Local Scalars .. 00114 LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1 00115 INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT 00116 PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 ) 00117 REAL AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL, 00118 $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA, 00119 $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX, 00120 $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2 00121 * .. 00122 * .. External Functions .. 00123 LOGICAL SISNAN 00124 REAL SLAMCH 00125 EXTERNAL SISNAN, SLAMCH 00126 * .. 00127 * .. External Subroutines .. 00128 EXTERNAL SCOPY 00129 * .. 00130 * .. Intrinsic Functions .. 00131 INTRINSIC ABS 00132 * .. 00133 * .. Executable Statements .. 00134 * 00135 INFO = 0 00136 FACT = REAL(2**KTRYMAX) 00137 EPS = SLAMCH( 'Precision' ) 00138 SHIFT = 0 00139 FORCER = .FALSE. 00140 00141 00142 * Note that we cannot guarantee that for any of the shifts tried, 00143 * the factorization has a small or even moderate element growth. 00144 * There could be Ritz values at both ends of the cluster and despite 00145 * backing off, there are examples where all factorizations tried 00146 * (in IEEE mode, allowing zero pivots & infinities) have INFINITE 00147 * element growth. 00148 * For this reason, we should use PIVMIN in this subroutine so that at 00149 * least the L D L^T factorization exists. It can be checked afterwards 00150 * whether the element growth caused bad residuals/orthogonality. 00151 00152 * Decide whether the code should accept the best among all 00153 * representations despite large element growth or signal INFO=1 00154 NOFAIL = .TRUE. 00155 * 00156 00157 * Compute the average gap length of the cluster 00158 CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT) 00159 AVGAP = CLWDTH / REAL(CLEND-CLSTRT) 00160 MINGAP = MIN(CLGAPL, CLGAPR) 00161 * Initial values for shifts to both ends of cluster 00162 LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT ) 00163 RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND ) 00164 00165 * Use a small fudge to make sure that we really shift to the outside 00166 LSIGMA = LSIGMA - ABS(LSIGMA)* TWO * EPS 00167 RSIGMA = RSIGMA + ABS(RSIGMA)* TWO * EPS 00168 00169 * Compute upper bounds for how much to back off the initial shifts 00170 LDMAX = QUART * MINGAP + TWO * PIVMIN 00171 RDMAX = QUART * MINGAP + TWO * PIVMIN 00172 00173 LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT 00174 RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT 00175 * 00176 * Initialize the record of the best representation found 00177 * 00178 S = SLAMCH( 'S' ) 00179 SMLGROWTH = ONE / S 00180 FAIL = REAL(N-1)*MINGAP/(SPDIAM*EPS) 00181 FAIL2 = REAL(N-1)*MINGAP/(SPDIAM*SQRT(EPS)) 00182 BESTSHIFT = LSIGMA 00183 * 00184 * while (KTRY <= KTRYMAX) 00185 KTRY = 0 00186 GROWTHBOUND = MAXGROWTH1*SPDIAM 00187 00188 5 CONTINUE 00189 SAWNAN1 = .FALSE. 00190 SAWNAN2 = .FALSE. 00191 * Ensure that we do not back off too much of the initial shifts 00192 LDELTA = MIN(LDMAX,LDELTA) 00193 RDELTA = MIN(RDMAX,RDELTA) 00194 00195 * Compute the element growth when shifting to both ends of the cluster 00196 * accept the shift if there is no element growth at one of the two ends 00197 00198 * Left end 00199 S = -LSIGMA 00200 DPLUS( 1 ) = D( 1 ) + S 00201 IF(ABS(DPLUS(1)).LT.PIVMIN) THEN 00202 DPLUS(1) = -PIVMIN 00203 * Need to set SAWNAN1 because refined RRR test should not be used 00204 * in this case 00205 SAWNAN1 = .TRUE. 00206 ENDIF 00207 MAX1 = ABS( DPLUS( 1 ) ) 00208 DO 6 I = 1, N - 1 00209 LPLUS( I ) = LD( I ) / DPLUS( I ) 00210 S = S*LPLUS( I )*L( I ) - LSIGMA 00211 DPLUS( I+1 ) = D( I+1 ) + S 00212 IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN 00213 DPLUS(I+1) = -PIVMIN 00214 * Need to set SAWNAN1 because refined RRR test should not be used 00215 * in this case 00216 SAWNAN1 = .TRUE. 00217 ENDIF 00218 MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) ) 00219 6 CONTINUE 00220 SAWNAN1 = SAWNAN1 .OR. SISNAN( MAX1 ) 00221 00222 IF( FORCER .OR. 00223 $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN 00224 SIGMA = LSIGMA 00225 SHIFT = SLEFT 00226 GOTO 100 00227 ENDIF 00228 00229 * Right end 00230 S = -RSIGMA 00231 WORK( 1 ) = D( 1 ) + S 00232 IF(ABS(WORK(1)).LT.PIVMIN) THEN 00233 WORK(1) = -PIVMIN 00234 * Need to set SAWNAN2 because refined RRR test should not be used 00235 * in this case 00236 SAWNAN2 = .TRUE. 00237 ENDIF 00238 MAX2 = ABS( WORK( 1 ) ) 00239 DO 7 I = 1, N - 1 00240 WORK( N+I ) = LD( I ) / WORK( I ) 00241 S = S*WORK( N+I )*L( I ) - RSIGMA 00242 WORK( I+1 ) = D( I+1 ) + S 00243 IF(ABS(WORK(I+1)).LT.PIVMIN) THEN 00244 WORK(I+1) = -PIVMIN 00245 * Need to set SAWNAN2 because refined RRR test should not be used 00246 * in this case 00247 SAWNAN2 = .TRUE. 00248 ENDIF 00249 MAX2 = MAX( MAX2,ABS(WORK(I+1)) ) 00250 7 CONTINUE 00251 SAWNAN2 = SAWNAN2 .OR. SISNAN( MAX2 ) 00252 00253 IF( FORCER .OR. 00254 $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN 00255 SIGMA = RSIGMA 00256 SHIFT = SRIGHT 00257 GOTO 100 00258 ENDIF 00259 * If we are at this point, both shifts led to too much element growth 00260 00261 * Record the better of the two shifts (provided it didn't lead to NaN) 00262 IF(SAWNAN1.AND.SAWNAN2) THEN 00263 * both MAX1 and MAX2 are NaN 00264 GOTO 50 00265 ELSE 00266 IF( .NOT.SAWNAN1 ) THEN 00267 INDX = 1 00268 IF(MAX1.LE.SMLGROWTH) THEN 00269 SMLGROWTH = MAX1 00270 BESTSHIFT = LSIGMA 00271 ENDIF 00272 ENDIF 00273 IF( .NOT.SAWNAN2 ) THEN 00274 IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2 00275 IF(MAX2.LE.SMLGROWTH) THEN 00276 SMLGROWTH = MAX2 00277 BESTSHIFT = RSIGMA 00278 ENDIF 00279 ENDIF 00280 ENDIF 00281 00282 * If we are here, both the left and the right shift led to 00283 * element growth. If the element growth is moderate, then 00284 * we may still accept the representation, if it passes a 00285 * refined test for RRR. This test supposes that no NaN occurred. 00286 * Moreover, we use the refined RRR test only for isolated clusters. 00287 IF((CLWDTH.LT.MINGAP/REAL(128)) .AND. 00288 $ (MIN(MAX1,MAX2).LT.FAIL2) 00289 $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN 00290 DORRR1 = .TRUE. 00291 ELSE 00292 DORRR1 = .FALSE. 00293 ENDIF 00294 TRYRRR1 = .TRUE. 00295 IF( TRYRRR1 .AND. DORRR1 ) THEN 00296 IF(INDX.EQ.1) THEN 00297 TMP = ABS( DPLUS( N ) ) 00298 ZNM2 = ONE 00299 PROD = ONE 00300 OLDP = ONE 00301 DO 15 I = N-1, 1, -1 00302 IF( PROD .LE. EPS ) THEN 00303 PROD = 00304 $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP 00305 ELSE 00306 PROD = PROD*ABS(WORK(N+I)) 00307 END IF 00308 OLDP = PROD 00309 ZNM2 = ZNM2 + PROD**2 00310 TMP = MAX( TMP, ABS( DPLUS( I ) * PROD )) 00311 15 CONTINUE 00312 RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 00313 IF (RRR1.LE.MAXGROWTH2) THEN 00314 SIGMA = LSIGMA 00315 SHIFT = SLEFT 00316 GOTO 100 00317 ENDIF 00318 ELSE IF(INDX.EQ.2) THEN 00319 TMP = ABS( WORK( N ) ) 00320 ZNM2 = ONE 00321 PROD = ONE 00322 OLDP = ONE 00323 DO 16 I = N-1, 1, -1 00324 IF( PROD .LE. EPS ) THEN 00325 PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP 00326 ELSE 00327 PROD = PROD*ABS(LPLUS(I)) 00328 END IF 00329 OLDP = PROD 00330 ZNM2 = ZNM2 + PROD**2 00331 TMP = MAX( TMP, ABS( WORK( I ) * PROD )) 00332 16 CONTINUE 00333 RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 00334 IF (RRR2.LE.MAXGROWTH2) THEN 00335 SIGMA = RSIGMA 00336 SHIFT = SRIGHT 00337 GOTO 100 00338 ENDIF 00339 END IF 00340 ENDIF 00341 00342 50 CONTINUE 00343 00344 IF (KTRY.LT.KTRYMAX) THEN 00345 * If we are here, both shifts failed also the RRR test. 00346 * Back off to the outside 00347 LSIGMA = MAX( LSIGMA - LDELTA, 00348 $ LSIGMA - LDMAX) 00349 RSIGMA = MIN( RSIGMA + RDELTA, 00350 $ RSIGMA + RDMAX ) 00351 LDELTA = TWO * LDELTA 00352 RDELTA = TWO * RDELTA 00353 KTRY = KTRY + 1 00354 GOTO 5 00355 ELSE 00356 * None of the representations investigated satisfied our 00357 * criteria. Take the best one we found. 00358 IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN 00359 LSIGMA = BESTSHIFT 00360 RSIGMA = BESTSHIFT 00361 FORCER = .TRUE. 00362 GOTO 5 00363 ELSE 00364 INFO = 1 00365 RETURN 00366 ENDIF 00367 END IF 00368 00369 100 CONTINUE 00370 IF (SHIFT.EQ.SLEFT) THEN 00371 ELSEIF (SHIFT.EQ.SRIGHT) THEN 00372 * store new L and D back into DPLUS, LPLUS 00373 CALL SCOPY( N, WORK, 1, DPLUS, 1 ) 00374 CALL SCOPY( N-1, WORK(N+1), 1, LPLUS, 1 ) 00375 ENDIF 00376 00377 RETURN 00378 * 00379 * End of SLARRF 00380 * 00381 END
1.7.3 
