![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLARRF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLARRF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND, 00022 * W, WGAP, WERR, 00023 * SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, 00024 * DPLUS, LPLUS, WORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * INTEGER CLSTRT, CLEND, INFO, N 00028 * DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 00029 * .. 00030 * .. Array Arguments .. 00031 * DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ), 00032 * $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> Given the initial representation L D L^T and its cluster of close 00042 *> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... 00043 *> W( CLEND ), DLARRF finds a new relatively robust representation 00044 *> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the 00045 *> eigenvalues of L(+) D(+) L(+)^T is relatively isolated. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] N 00052 *> \verbatim 00053 *> N is INTEGER 00054 *> The order of the matrix (subblock, if the matrix splitted). 00055 *> \endverbatim 00056 *> 00057 *> \param[in] D 00058 *> \verbatim 00059 *> D is DOUBLE PRECISION array, dimension (N) 00060 *> The N diagonal elements of the diagonal matrix D. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] L 00064 *> \verbatim 00065 *> L is DOUBLE PRECISION array, dimension (N-1) 00066 *> The (N-1) subdiagonal elements of the unit bidiagonal 00067 *> matrix L. 00068 *> \endverbatim 00069 *> 00070 *> \param[in] LD 00071 *> \verbatim 00072 *> LD is DOUBLE PRECISION array, dimension (N-1) 00073 *> The (N-1) elements L(i)*D(i). 00074 *> \endverbatim 00075 *> 00076 *> \param[in] CLSTRT 00077 *> \verbatim 00078 *> CLSTRT is INTEGER 00079 *> The index of the first eigenvalue in the cluster. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] CLEND 00083 *> \verbatim 00084 *> CLEND is INTEGER 00085 *> The index of the last eigenvalue in the cluster. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] W 00089 *> \verbatim 00090 *> W is DOUBLE PRECISION array, dimension 00091 *> dimension is >= (CLEND-CLSTRT+1) 00092 *> The eigenvalue APPROXIMATIONS of L D L^T in ascending order. 00093 *> W( CLSTRT ) through W( CLEND ) form the cluster of relatively 00094 *> close eigenalues. 00095 *> \endverbatim 00096 *> 00097 *> \param[in,out] WGAP 00098 *> \verbatim 00099 *> WGAP is DOUBLE PRECISION array, dimension 00100 *> dimension is >= (CLEND-CLSTRT+1) 00101 *> The separation from the right neighbor eigenvalue in W. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] WERR 00105 *> \verbatim 00106 *> WERR is DOUBLE PRECISION array, dimension 00107 *> dimension is >= (CLEND-CLSTRT+1) 00108 *> WERR contain the semiwidth of the uncertainty 00109 *> interval of the corresponding eigenvalue APPROXIMATION in W 00110 *> \endverbatim 00111 *> 00112 *> \param[in] SPDIAM 00113 *> \verbatim 00114 *> SPDIAM is DOUBLE PRECISION 00115 *> estimate of the spectral diameter obtained from the 00116 *> Gerschgorin intervals 00117 *> \endverbatim 00118 *> 00119 *> \param[in] CLGAPL 00120 *> \verbatim 00121 *> CLGAPL is DOUBLE PRECISION 00122 *> \endverbatim 00123 *> 00124 *> \param[in] CLGAPR 00125 *> \verbatim 00126 *> CLGAPR is DOUBLE PRECISION 00127 *> absolute gap on each end of the cluster. 00128 *> Set by the calling routine to protect against shifts too close 00129 *> to eigenvalues outside the cluster. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] PIVMIN 00133 *> \verbatim 00134 *> PIVMIN is DOUBLE PRECISION 00135 *> The minimum pivot allowed in the Sturm sequence. 00136 *> \endverbatim 00137 *> 00138 *> \param[out] SIGMA 00139 *> \verbatim 00140 *> SIGMA is DOUBLE PRECISION 00141 *> The shift used to form L(+) D(+) L(+)^T. 00142 *> \endverbatim 00143 *> 00144 *> \param[out] DPLUS 00145 *> \verbatim 00146 *> DPLUS is DOUBLE PRECISION array, dimension (N) 00147 *> The N diagonal elements of the diagonal matrix D(+). 00148 *> \endverbatim 00149 *> 00150 *> \param[out] LPLUS 00151 *> \verbatim 00152 *> LPLUS is DOUBLE PRECISION array, dimension (N-1) 00153 *> The first (N-1) elements of LPLUS contain the subdiagonal 00154 *> elements of the unit bidiagonal matrix L(+). 00155 *> \endverbatim 00156 *> 00157 *> \param[out] WORK 00158 *> \verbatim 00159 *> WORK is DOUBLE PRECISION array, dimension (2*N) 00160 *> Workspace. 00161 *> \endverbatim 00162 *> 00163 *> \param[out] INFO 00164 *> \verbatim 00165 *> INFO is INTEGER 00166 *> Signals processing OK (=0) or failure (=1) 00167 *> \endverbatim 00168 * 00169 * Authors: 00170 * ======== 00171 * 00172 *> \author Univ. of Tennessee 00173 *> \author Univ. of California Berkeley 00174 *> \author Univ. of Colorado Denver 00175 *> \author NAG Ltd. 00176 * 00177 *> \date November 2011 00178 * 00179 *> \ingroup auxOTHERauxiliary 00180 * 00181 *> \par Contributors: 00182 * ================== 00183 *> 00184 *> Beresford Parlett, University of California, Berkeley, USA \n 00185 *> Jim Demmel, University of California, Berkeley, USA \n 00186 *> Inderjit Dhillon, University of Texas, Austin, USA \n 00187 *> Osni Marques, LBNL/NERSC, USA \n 00188 *> Christof Voemel, University of California, Berkeley, USA 00189 * 00190 * ===================================================================== 00191 SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND, 00192 $ W, WGAP, WERR, 00193 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, 00194 $ DPLUS, LPLUS, WORK, INFO ) 00195 * 00196 * -- LAPACK auxiliary routine (version 3.4.0) -- 00197 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00198 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00199 * November 2011 00200 * 00201 * .. Scalar Arguments .. 00202 INTEGER CLSTRT, CLEND, INFO, N 00203 DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM 00204 * .. 00205 * .. Array Arguments .. 00206 DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ), 00207 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * ) 00208 * .. 00209 * 00210 * ===================================================================== 00211 * 00212 * .. Parameters .. 00213 DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO 00214 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0, FOUR = 4.0D0, 00215 $ QUART = 0.25D0, 00216 $ MAXGROWTH1 = 8.D0, 00217 $ MAXGROWTH2 = 8.D0 ) 00218 * .. 00219 * .. Local Scalars .. 00220 LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1 00221 INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT 00222 PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 ) 00223 DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL, 00224 $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA, 00225 $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX, 00226 $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2 00227 * .. 00228 * .. External Functions .. 00229 LOGICAL DISNAN 00230 DOUBLE PRECISION DLAMCH 00231 EXTERNAL DISNAN, DLAMCH 00232 * .. 00233 * .. External Subroutines .. 00234 EXTERNAL DCOPY 00235 * .. 00236 * .. Intrinsic Functions .. 00237 INTRINSIC ABS 00238 * .. 00239 * .. Executable Statements .. 00240 * 00241 INFO = 0 00242 FACT = DBLE(2**KTRYMAX) 00243 EPS = DLAMCH( 'Precision' ) 00244 SHIFT = 0 00245 FORCER = .FALSE. 00246 00247 00248 * Note that we cannot guarantee that for any of the shifts tried, 00249 * the factorization has a small or even moderate element growth. 00250 * There could be Ritz values at both ends of the cluster and despite 00251 * backing off, there are examples where all factorizations tried 00252 * (in IEEE mode, allowing zero pivots & infinities) have INFINITE 00253 * element growth. 00254 * For this reason, we should use PIVMIN in this subroutine so that at 00255 * least the L D L^T factorization exists. It can be checked afterwards 00256 * whether the element growth caused bad residuals/orthogonality. 00257 00258 * Decide whether the code should accept the best among all 00259 * representations despite large element growth or signal INFO=1 00260 NOFAIL = .TRUE. 00261 * 00262 00263 * Compute the average gap length of the cluster 00264 CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT) 00265 AVGAP = CLWDTH / DBLE(CLEND-CLSTRT) 00266 MINGAP = MIN(CLGAPL, CLGAPR) 00267 * Initial values for shifts to both ends of cluster 00268 LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT ) 00269 RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND ) 00270 00271 * Use a small fudge to make sure that we really shift to the outside 00272 LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS 00273 RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS 00274 00275 * Compute upper bounds for how much to back off the initial shifts 00276 LDMAX = QUART * MINGAP + TWO * PIVMIN 00277 RDMAX = QUART * MINGAP + TWO * PIVMIN 00278 00279 LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT 00280 RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT 00281 * 00282 * Initialize the record of the best representation found 00283 * 00284 S = DLAMCH( 'S' ) 00285 SMLGROWTH = ONE / S 00286 FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS) 00287 FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS)) 00288 BESTSHIFT = LSIGMA 00289 * 00290 * while (KTRY <= KTRYMAX) 00291 KTRY = 0 00292 GROWTHBOUND = MAXGROWTH1*SPDIAM 00293 00294 5 CONTINUE 00295 SAWNAN1 = .FALSE. 00296 SAWNAN2 = .FALSE. 00297 * Ensure that we do not back off too much of the initial shifts 00298 LDELTA = MIN(LDMAX,LDELTA) 00299 RDELTA = MIN(RDMAX,RDELTA) 00300 00301 * Compute the element growth when shifting to both ends of the cluster 00302 * accept the shift if there is no element growth at one of the two ends 00303 00304 * Left end 00305 S = -LSIGMA 00306 DPLUS( 1 ) = D( 1 ) + S 00307 IF(ABS(DPLUS(1)).LT.PIVMIN) THEN 00308 DPLUS(1) = -PIVMIN 00309 * Need to set SAWNAN1 because refined RRR test should not be used 00310 * in this case 00311 SAWNAN1 = .TRUE. 00312 ENDIF 00313 MAX1 = ABS( DPLUS( 1 ) ) 00314 DO 6 I = 1, N - 1 00315 LPLUS( I ) = LD( I ) / DPLUS( I ) 00316 S = S*LPLUS( I )*L( I ) - LSIGMA 00317 DPLUS( I+1 ) = D( I+1 ) + S 00318 IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN 00319 DPLUS(I+1) = -PIVMIN 00320 * Need to set SAWNAN1 because refined RRR test should not be used 00321 * in this case 00322 SAWNAN1 = .TRUE. 00323 ENDIF 00324 MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) ) 00325 6 CONTINUE 00326 SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 ) 00327 00328 IF( FORCER .OR. 00329 $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN 00330 SIGMA = LSIGMA 00331 SHIFT = SLEFT 00332 GOTO 100 00333 ENDIF 00334 00335 * Right end 00336 S = -RSIGMA 00337 WORK( 1 ) = D( 1 ) + S 00338 IF(ABS(WORK(1)).LT.PIVMIN) THEN 00339 WORK(1) = -PIVMIN 00340 * Need to set SAWNAN2 because refined RRR test should not be used 00341 * in this case 00342 SAWNAN2 = .TRUE. 00343 ENDIF 00344 MAX2 = ABS( WORK( 1 ) ) 00345 DO 7 I = 1, N - 1 00346 WORK( N+I ) = LD( I ) / WORK( I ) 00347 S = S*WORK( N+I )*L( I ) - RSIGMA 00348 WORK( I+1 ) = D( I+1 ) + S 00349 IF(ABS(WORK(I+1)).LT.PIVMIN) THEN 00350 WORK(I+1) = -PIVMIN 00351 * Need to set SAWNAN2 because refined RRR test should not be used 00352 * in this case 00353 SAWNAN2 = .TRUE. 00354 ENDIF 00355 MAX2 = MAX( MAX2,ABS(WORK(I+1)) ) 00356 7 CONTINUE 00357 SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 ) 00358 00359 IF( FORCER .OR. 00360 $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN 00361 SIGMA = RSIGMA 00362 SHIFT = SRIGHT 00363 GOTO 100 00364 ENDIF 00365 * If we are at this point, both shifts led to too much element growth 00366 00367 * Record the better of the two shifts (provided it didn't lead to NaN) 00368 IF(SAWNAN1.AND.SAWNAN2) THEN 00369 * both MAX1 and MAX2 are NaN 00370 GOTO 50 00371 ELSE 00372 IF( .NOT.SAWNAN1 ) THEN 00373 INDX = 1 00374 IF(MAX1.LE.SMLGROWTH) THEN 00375 SMLGROWTH = MAX1 00376 BESTSHIFT = LSIGMA 00377 ENDIF 00378 ENDIF 00379 IF( .NOT.SAWNAN2 ) THEN 00380 IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2 00381 IF(MAX2.LE.SMLGROWTH) THEN 00382 SMLGROWTH = MAX2 00383 BESTSHIFT = RSIGMA 00384 ENDIF 00385 ENDIF 00386 ENDIF 00387 00388 * If we are here, both the left and the right shift led to 00389 * element growth. If the element growth is moderate, then 00390 * we may still accept the representation, if it passes a 00391 * refined test for RRR. This test supposes that no NaN occurred. 00392 * Moreover, we use the refined RRR test only for isolated clusters. 00393 IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND. 00394 $ (MIN(MAX1,MAX2).LT.FAIL2) 00395 $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN 00396 DORRR1 = .TRUE. 00397 ELSE 00398 DORRR1 = .FALSE. 00399 ENDIF 00400 TRYRRR1 = .TRUE. 00401 IF( TRYRRR1 .AND. DORRR1 ) THEN 00402 IF(INDX.EQ.1) THEN 00403 TMP = ABS( DPLUS( N ) ) 00404 ZNM2 = ONE 00405 PROD = ONE 00406 OLDP = ONE 00407 DO 15 I = N-1, 1, -1 00408 IF( PROD .LE. EPS ) THEN 00409 PROD = 00410 $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP 00411 ELSE 00412 PROD = PROD*ABS(WORK(N+I)) 00413 END IF 00414 OLDP = PROD 00415 ZNM2 = ZNM2 + PROD**2 00416 TMP = MAX( TMP, ABS( DPLUS( I ) * PROD )) 00417 15 CONTINUE 00418 RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 00419 IF (RRR1.LE.MAXGROWTH2) THEN 00420 SIGMA = LSIGMA 00421 SHIFT = SLEFT 00422 GOTO 100 00423 ENDIF 00424 ELSE IF(INDX.EQ.2) THEN 00425 TMP = ABS( WORK( N ) ) 00426 ZNM2 = ONE 00427 PROD = ONE 00428 OLDP = ONE 00429 DO 16 I = N-1, 1, -1 00430 IF( PROD .LE. EPS ) THEN 00431 PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP 00432 ELSE 00433 PROD = PROD*ABS(LPLUS(I)) 00434 END IF 00435 OLDP = PROD 00436 ZNM2 = ZNM2 + PROD**2 00437 TMP = MAX( TMP, ABS( WORK( I ) * PROD )) 00438 16 CONTINUE 00439 RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) ) 00440 IF (RRR2.LE.MAXGROWTH2) THEN 00441 SIGMA = RSIGMA 00442 SHIFT = SRIGHT 00443 GOTO 100 00444 ENDIF 00445 END IF 00446 ENDIF 00447 00448 50 CONTINUE 00449 00450 IF (KTRY.LT.KTRYMAX) THEN 00451 * If we are here, both shifts failed also the RRR test. 00452 * Back off to the outside 00453 LSIGMA = MAX( LSIGMA - LDELTA, 00454 $ LSIGMA - LDMAX) 00455 RSIGMA = MIN( RSIGMA + RDELTA, 00456 $ RSIGMA + RDMAX ) 00457 LDELTA = TWO * LDELTA 00458 RDELTA = TWO * RDELTA 00459 KTRY = KTRY + 1 00460 GOTO 5 00461 ELSE 00462 * None of the representations investigated satisfied our 00463 * criteria. Take the best one we found. 00464 IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN 00465 LSIGMA = BESTSHIFT 00466 RSIGMA = BESTSHIFT 00467 FORCER = .TRUE. 00468 GOTO 5 00469 ELSE 00470 INFO = 1 00471 RETURN 00472 ENDIF 00473 END IF 00474 00475 100 CONTINUE 00476 IF (SHIFT.EQ.SLEFT) THEN 00477 ELSEIF (SHIFT.EQ.SRIGHT) THEN 00478 * store new L and D back into DPLUS, LPLUS 00479 CALL DCOPY( N, WORK, 1, DPLUS, 1 ) 00480 CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 ) 00481 ENDIF 00482 00483 RETURN 00484 * 00485 * End of DLARRF 00486 * 00487 END