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