![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAR1V 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAR1V + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlar1v.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlar1v.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlar1v.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD, 00022 * PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, 00023 * R, ISUPPZ, NRMINV, RESID, RQCORR, WORK ) 00024 * 00025 * .. Scalar Arguments .. 00026 * LOGICAL WANTNC 00027 * INTEGER B1, BN, N, NEGCNT, R 00028 * DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 00029 * $ RQCORR, ZTZ 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER ISUPPZ( * ) 00033 * DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ), 00034 * $ WORK( * ) 00035 * DOUBLE PRECISION Z( * ) 00036 * .. 00037 * 00038 * 00039 *> \par Purpose: 00040 * ============= 00041 *> 00042 *> \verbatim 00043 *> 00044 *> DLAR1V computes the (scaled) r-th column of the inverse of 00045 *> the sumbmatrix in rows B1 through BN of the tridiagonal matrix 00046 *> L D L**T - sigma I. When sigma is close to an eigenvalue, the 00047 *> computed vector is an accurate eigenvector. Usually, r corresponds 00048 *> to the index where the eigenvector is largest in magnitude. 00049 *> The following steps accomplish this computation : 00050 *> (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T, 00051 *> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T, 00052 *> (c) Computation of the diagonal elements of the inverse of 00053 *> L D L**T - sigma I by combining the above transforms, and choosing 00054 *> r as the index where the diagonal of the inverse is (one of the) 00055 *> largest in magnitude. 00056 *> (d) Computation of the (scaled) r-th column of the inverse using the 00057 *> twisted factorization obtained by combining the top part of the 00058 *> the stationary and the bottom part of the progressive transform. 00059 *> \endverbatim 00060 * 00061 * Arguments: 00062 * ========== 00063 * 00064 *> \param[in] N 00065 *> \verbatim 00066 *> N is INTEGER 00067 *> The order of the matrix L D L**T. 00068 *> \endverbatim 00069 *> 00070 *> \param[in] B1 00071 *> \verbatim 00072 *> B1 is INTEGER 00073 *> First index of the submatrix of L D L**T. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] BN 00077 *> \verbatim 00078 *> BN is INTEGER 00079 *> Last index of the submatrix of L D L**T. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LAMBDA 00083 *> \verbatim 00084 *> LAMBDA is DOUBLE PRECISION 00085 *> The shift. In order to compute an accurate eigenvector, 00086 *> LAMBDA should be a good approximation to an eigenvalue 00087 *> of L D L**T. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] L 00091 *> \verbatim 00092 *> L is DOUBLE PRECISION array, dimension (N-1) 00093 *> The (n-1) subdiagonal elements of the unit bidiagonal matrix 00094 *> L, in elements 1 to N-1. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] D 00098 *> \verbatim 00099 *> D is DOUBLE PRECISION array, dimension (N) 00100 *> The n diagonal elements of the diagonal matrix D. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] LD 00104 *> \verbatim 00105 *> LD is DOUBLE PRECISION array, dimension (N-1) 00106 *> The n-1 elements L(i)*D(i). 00107 *> \endverbatim 00108 *> 00109 *> \param[in] LLD 00110 *> \verbatim 00111 *> LLD is DOUBLE PRECISION array, dimension (N-1) 00112 *> The n-1 elements L(i)*L(i)*D(i). 00113 *> \endverbatim 00114 *> 00115 *> \param[in] PIVMIN 00116 *> \verbatim 00117 *> PIVMIN is DOUBLE PRECISION 00118 *> The minimum pivot in the Sturm sequence. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] GAPTOL 00122 *> \verbatim 00123 *> GAPTOL is DOUBLE PRECISION 00124 *> Tolerance that indicates when eigenvector entries are negligible 00125 *> w.r.t. their contribution to the residual. 00126 *> \endverbatim 00127 *> 00128 *> \param[in,out] Z 00129 *> \verbatim 00130 *> Z is DOUBLE PRECISION array, dimension (N) 00131 *> On input, all entries of Z must be set to 0. 00132 *> On output, Z contains the (scaled) r-th column of the 00133 *> inverse. The scaling is such that Z(R) equals 1. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] WANTNC 00137 *> \verbatim 00138 *> WANTNC is LOGICAL 00139 *> Specifies whether NEGCNT has to be computed. 00140 *> \endverbatim 00141 *> 00142 *> \param[out] NEGCNT 00143 *> \verbatim 00144 *> NEGCNT is INTEGER 00145 *> If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin 00146 *> in the matrix factorization L D L**T, and NEGCNT = -1 otherwise. 00147 *> \endverbatim 00148 *> 00149 *> \param[out] ZTZ 00150 *> \verbatim 00151 *> ZTZ is DOUBLE PRECISION 00152 *> The square of the 2-norm of Z. 00153 *> \endverbatim 00154 *> 00155 *> \param[out] MINGMA 00156 *> \verbatim 00157 *> MINGMA is DOUBLE PRECISION 00158 *> The reciprocal of the largest (in magnitude) diagonal 00159 *> element of the inverse of L D L**T - sigma I. 00160 *> \endverbatim 00161 *> 00162 *> \param[in,out] R 00163 *> \verbatim 00164 *> R is INTEGER 00165 *> The twist index for the twisted factorization used to 00166 *> compute Z. 00167 *> On input, 0 <= R <= N. If R is input as 0, R is set to 00168 *> the index where (L D L**T - sigma I)^{-1} is largest 00169 *> in magnitude. If 1 <= R <= N, R is unchanged. 00170 *> On output, R contains the twist index used to compute Z. 00171 *> Ideally, R designates the position of the maximum entry in the 00172 *> eigenvector. 00173 *> \endverbatim 00174 *> 00175 *> \param[out] ISUPPZ 00176 *> \verbatim 00177 *> ISUPPZ is INTEGER array, dimension (2) 00178 *> The support of the vector in Z, i.e., the vector Z is 00179 *> nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). 00180 *> \endverbatim 00181 *> 00182 *> \param[out] NRMINV 00183 *> \verbatim 00184 *> NRMINV is DOUBLE PRECISION 00185 *> NRMINV = 1/SQRT( ZTZ ) 00186 *> \endverbatim 00187 *> 00188 *> \param[out] RESID 00189 *> \verbatim 00190 *> RESID is DOUBLE PRECISION 00191 *> The residual of the FP vector. 00192 *> RESID = ABS( MINGMA )/SQRT( ZTZ ) 00193 *> \endverbatim 00194 *> 00195 *> \param[out] RQCORR 00196 *> \verbatim 00197 *> RQCORR is DOUBLE PRECISION 00198 *> The Rayleigh Quotient correction to LAMBDA. 00199 *> RQCORR = MINGMA*TMP 00200 *> \endverbatim 00201 *> 00202 *> \param[out] WORK 00203 *> \verbatim 00204 *> WORK is DOUBLE PRECISION array, dimension (4*N) 00205 *> \endverbatim 00206 * 00207 * Authors: 00208 * ======== 00209 * 00210 *> \author Univ. of Tennessee 00211 *> \author Univ. of California Berkeley 00212 *> \author Univ. of Colorado Denver 00213 *> \author NAG Ltd. 00214 * 00215 *> \date November 2011 00216 * 00217 *> \ingroup doubleOTHERauxiliary 00218 * 00219 *> \par Contributors: 00220 * ================== 00221 *> 00222 *> Beresford Parlett, University of California, Berkeley, USA \n 00223 *> Jim Demmel, University of California, Berkeley, USA \n 00224 *> Inderjit Dhillon, University of Texas, Austin, USA \n 00225 *> Osni Marques, LBNL/NERSC, USA \n 00226 *> Christof Voemel, University of California, Berkeley, USA 00227 * 00228 * ===================================================================== 00229 SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD, 00230 $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, 00231 $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK ) 00232 * 00233 * -- LAPACK auxiliary routine (version 3.4.0) -- 00234 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00235 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00236 * November 2011 00237 * 00238 * .. Scalar Arguments .. 00239 LOGICAL WANTNC 00240 INTEGER B1, BN, N, NEGCNT, R 00241 DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 00242 $ RQCORR, ZTZ 00243 * .. 00244 * .. Array Arguments .. 00245 INTEGER ISUPPZ( * ) 00246 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ), 00247 $ WORK( * ) 00248 DOUBLE PRECISION Z( * ) 00249 * .. 00250 * 00251 * ===================================================================== 00252 * 00253 * .. Parameters .. 00254 DOUBLE PRECISION ZERO, ONE 00255 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00256 00257 * .. 00258 * .. Local Scalars .. 00259 LOGICAL SAWNAN1, SAWNAN2 00260 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1, 00261 $ R2 00262 DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP 00263 * .. 00264 * .. External Functions .. 00265 LOGICAL DISNAN 00266 DOUBLE PRECISION DLAMCH 00267 EXTERNAL DISNAN, DLAMCH 00268 * .. 00269 * .. Intrinsic Functions .. 00270 INTRINSIC ABS 00271 * .. 00272 * .. Executable Statements .. 00273 * 00274 EPS = DLAMCH( 'Precision' ) 00275 00276 00277 IF( R.EQ.0 ) THEN 00278 R1 = B1 00279 R2 = BN 00280 ELSE 00281 R1 = R 00282 R2 = R 00283 END IF 00284 00285 * Storage for LPLUS 00286 INDLPL = 0 00287 * Storage for UMINUS 00288 INDUMN = N 00289 INDS = 2*N + 1 00290 INDP = 3*N + 1 00291 00292 IF( B1.EQ.1 ) THEN 00293 WORK( INDS ) = ZERO 00294 ELSE 00295 WORK( INDS+B1-1 ) = LLD( B1-1 ) 00296 END IF 00297 00298 * 00299 * Compute the stationary transform (using the differential form) 00300 * until the index R2. 00301 * 00302 SAWNAN1 = .FALSE. 00303 NEG1 = 0 00304 S = WORK( INDS+B1-1 ) - LAMBDA 00305 DO 50 I = B1, R1 - 1 00306 DPLUS = D( I ) + S 00307 WORK( INDLPL+I ) = LD( I ) / DPLUS 00308 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00309 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00310 S = WORK( INDS+I ) - LAMBDA 00311 50 CONTINUE 00312 SAWNAN1 = DISNAN( S ) 00313 IF( SAWNAN1 ) GOTO 60 00314 DO 51 I = R1, R2 - 1 00315 DPLUS = D( I ) + S 00316 WORK( INDLPL+I ) = LD( I ) / DPLUS 00317 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00318 S = WORK( INDS+I ) - LAMBDA 00319 51 CONTINUE 00320 SAWNAN1 = DISNAN( S ) 00321 * 00322 60 CONTINUE 00323 IF( SAWNAN1 ) THEN 00324 * Runs a slower version of the above loop if a NaN is detected 00325 NEG1 = 0 00326 S = WORK( INDS+B1-1 ) - LAMBDA 00327 DO 70 I = B1, R1 - 1 00328 DPLUS = D( I ) + S 00329 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00330 WORK( INDLPL+I ) = LD( I ) / DPLUS 00331 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00332 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00333 IF( WORK( INDLPL+I ).EQ.ZERO ) 00334 $ WORK( INDS+I ) = LLD( I ) 00335 S = WORK( INDS+I ) - LAMBDA 00336 70 CONTINUE 00337 DO 71 I = R1, R2 - 1 00338 DPLUS = D( I ) + S 00339 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00340 WORK( INDLPL+I ) = LD( I ) / DPLUS 00341 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00342 IF( WORK( INDLPL+I ).EQ.ZERO ) 00343 $ WORK( INDS+I ) = LLD( I ) 00344 S = WORK( INDS+I ) - LAMBDA 00345 71 CONTINUE 00346 END IF 00347 * 00348 * Compute the progressive transform (using the differential form) 00349 * until the index R1 00350 * 00351 SAWNAN2 = .FALSE. 00352 NEG2 = 0 00353 WORK( INDP+BN-1 ) = D( BN ) - LAMBDA 00354 DO 80 I = BN - 1, R1, -1 00355 DMINUS = LLD( I ) + WORK( INDP+I ) 00356 TMP = D( I ) / DMINUS 00357 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00358 WORK( INDUMN+I ) = L( I )*TMP 00359 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00360 80 CONTINUE 00361 TMP = WORK( INDP+R1-1 ) 00362 SAWNAN2 = DISNAN( TMP ) 00363 00364 IF( SAWNAN2 ) THEN 00365 * Runs a slower version of the above loop if a NaN is detected 00366 NEG2 = 0 00367 DO 100 I = BN-1, R1, -1 00368 DMINUS = LLD( I ) + WORK( INDP+I ) 00369 IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN 00370 TMP = D( I ) / DMINUS 00371 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00372 WORK( INDUMN+I ) = L( I )*TMP 00373 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00374 IF( TMP.EQ.ZERO ) 00375 $ WORK( INDP+I-1 ) = D( I ) - LAMBDA 00376 100 CONTINUE 00377 END IF 00378 * 00379 * Find the index (from R1 to R2) of the largest (in magnitude) 00380 * diagonal element of the inverse 00381 * 00382 MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 ) 00383 IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1 00384 IF( WANTNC ) THEN 00385 NEGCNT = NEG1 + NEG2 00386 ELSE 00387 NEGCNT = -1 00388 ENDIF 00389 IF( ABS(MINGMA).EQ.ZERO ) 00390 $ MINGMA = EPS*WORK( INDS+R1-1 ) 00391 R = R1 00392 DO 110 I = R1, R2 - 1 00393 TMP = WORK( INDS+I ) + WORK( INDP+I ) 00394 IF( TMP.EQ.ZERO ) 00395 $ TMP = EPS*WORK( INDS+I ) 00396 IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN 00397 MINGMA = TMP 00398 R = I + 1 00399 END IF 00400 110 CONTINUE 00401 * 00402 * Compute the FP vector: solve N^T v = e_r 00403 * 00404 ISUPPZ( 1 ) = B1 00405 ISUPPZ( 2 ) = BN 00406 Z( R ) = ONE 00407 ZTZ = ONE 00408 * 00409 * Compute the FP vector upwards from R 00410 * 00411 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00412 DO 210 I = R-1, B1, -1 00413 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00414 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00415 $ THEN 00416 Z( I ) = ZERO 00417 ISUPPZ( 1 ) = I + 1 00418 GOTO 220 00419 ENDIF 00420 ZTZ = ZTZ + Z( I )*Z( I ) 00421 210 CONTINUE 00422 220 CONTINUE 00423 ELSE 00424 * Run slower loop if NaN occurred. 00425 DO 230 I = R - 1, B1, -1 00426 IF( Z( I+1 ).EQ.ZERO ) THEN 00427 Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) 00428 ELSE 00429 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00430 END IF 00431 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00432 $ THEN 00433 Z( I ) = ZERO 00434 ISUPPZ( 1 ) = I + 1 00435 GO TO 240 00436 END IF 00437 ZTZ = ZTZ + Z( I )*Z( I ) 00438 230 CONTINUE 00439 240 CONTINUE 00440 ENDIF 00441 00442 * Compute the FP vector downwards from R in blocks of size BLKSIZ 00443 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00444 DO 250 I = R, BN-1 00445 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00446 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00447 $ THEN 00448 Z( I+1 ) = ZERO 00449 ISUPPZ( 2 ) = I 00450 GO TO 260 00451 END IF 00452 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 00453 250 CONTINUE 00454 260 CONTINUE 00455 ELSE 00456 * Run slower loop if NaN occurred. 00457 DO 270 I = R, BN - 1 00458 IF( Z( I ).EQ.ZERO ) THEN 00459 Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) 00460 ELSE 00461 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00462 END IF 00463 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00464 $ THEN 00465 Z( I+1 ) = ZERO 00466 ISUPPZ( 2 ) = I 00467 GO TO 280 00468 END IF 00469 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 00470 270 CONTINUE 00471 280 CONTINUE 00472 END IF 00473 * 00474 * Compute quantities for convergence test 00475 * 00476 TMP = ONE / ZTZ 00477 NRMINV = SQRT( TMP ) 00478 RESID = ABS( MINGMA )*NRMINV 00479 RQCORR = MINGMA*TMP 00480 * 00481 * 00482 RETURN 00483 * 00484 * End of DLAR1V 00485 * 00486 END