|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD, 00002 $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, 00003 $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK ) 00004 * 00005 * -- LAPACK auxiliary routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 LOGICAL WANTNC 00012 INTEGER B1, BN, N, NEGCNT, R 00013 DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 00014 $ RQCORR, ZTZ 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER ISUPPZ( * ) 00018 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ), 00019 $ WORK( * ) 00020 DOUBLE PRECISION Z( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * DLAR1V computes the (scaled) r-th column of the inverse of 00027 * the sumbmatrix in rows B1 through BN of the tridiagonal matrix 00028 * L D L^T - sigma I. When sigma is close to an eigenvalue, the 00029 * computed vector is an accurate eigenvector. Usually, r corresponds 00030 * to the index where the eigenvector is largest in magnitude. 00031 * The following steps accomplish this computation : 00032 * (a) Stationary qd transform, L D L^T - sigma I = L(+) D(+) L(+)^T, 00033 * (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T, 00034 * (c) Computation of the diagonal elements of the inverse of 00035 * L D L^T - sigma I by combining the above transforms, and choosing 00036 * r as the index where the diagonal of the inverse is (one of the) 00037 * largest in magnitude. 00038 * (d) Computation of the (scaled) r-th column of the inverse using the 00039 * twisted factorization obtained by combining the top part of the 00040 * the stationary and the bottom part of the progressive transform. 00041 * 00042 * Arguments 00043 * ========= 00044 * 00045 * N (input) INTEGER 00046 * The order of the matrix L D L^T. 00047 * 00048 * B1 (input) INTEGER 00049 * First index of the submatrix of L D L^T. 00050 * 00051 * BN (input) INTEGER 00052 * Last index of the submatrix of L D L^T. 00053 * 00054 * LAMBDA (input) DOUBLE PRECISION 00055 * The shift. In order to compute an accurate eigenvector, 00056 * LAMBDA should be a good approximation to an eigenvalue 00057 * of L D L^T. 00058 * 00059 * L (input) DOUBLE PRECISION array, dimension (N-1) 00060 * The (n-1) subdiagonal elements of the unit bidiagonal matrix 00061 * L, in elements 1 to N-1. 00062 * 00063 * D (input) DOUBLE PRECISION array, dimension (N) 00064 * The n diagonal elements of the diagonal matrix D. 00065 * 00066 * LD (input) DOUBLE PRECISION array, dimension (N-1) 00067 * The n-1 elements L(i)*D(i). 00068 * 00069 * LLD (input) DOUBLE PRECISION array, dimension (N-1) 00070 * The n-1 elements L(i)*L(i)*D(i). 00071 * 00072 * PIVMIN (input) DOUBLE PRECISION 00073 * The minimum pivot in the Sturm sequence. 00074 * 00075 * GAPTOL (input) DOUBLE PRECISION 00076 * Tolerance that indicates when eigenvector entries are negligible 00077 * w.r.t. their contribution to the residual. 00078 * 00079 * Z (input/output) DOUBLE PRECISION array, dimension (N) 00080 * On input, all entries of Z must be set to 0. 00081 * On output, Z contains the (scaled) r-th column of the 00082 * inverse. The scaling is such that Z(R) equals 1. 00083 * 00084 * WANTNC (input) LOGICAL 00085 * Specifies whether NEGCNT has to be computed. 00086 * 00087 * NEGCNT (output) INTEGER 00088 * If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin 00089 * in the matrix factorization L D L^T, and NEGCNT = -1 otherwise. 00090 * 00091 * ZTZ (output) DOUBLE PRECISION 00092 * The square of the 2-norm of Z. 00093 * 00094 * MINGMA (output) DOUBLE PRECISION 00095 * The reciprocal of the largest (in magnitude) diagonal 00096 * element of the inverse of L D L^T - sigma I. 00097 * 00098 * R (input/output) INTEGER 00099 * The twist index for the twisted factorization used to 00100 * compute Z. 00101 * On input, 0 <= R <= N. If R is input as 0, R is set to 00102 * the index where (L D L^T - sigma I)^{-1} is largest 00103 * in magnitude. If 1 <= R <= N, R is unchanged. 00104 * On output, R contains the twist index used to compute Z. 00105 * Ideally, R designates the position of the maximum entry in the 00106 * eigenvector. 00107 * 00108 * ISUPPZ (output) INTEGER array, dimension (2) 00109 * The support of the vector in Z, i.e., the vector Z is 00110 * nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). 00111 * 00112 * NRMINV (output) DOUBLE PRECISION 00113 * NRMINV = 1/SQRT( ZTZ ) 00114 * 00115 * RESID (output) DOUBLE PRECISION 00116 * The residual of the FP vector. 00117 * RESID = ABS( MINGMA )/SQRT( ZTZ ) 00118 * 00119 * RQCORR (output) DOUBLE PRECISION 00120 * The Rayleigh Quotient correction to LAMBDA. 00121 * RQCORR = MINGMA*TMP 00122 * 00123 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 00124 * 00125 * Further Details 00126 * =============== 00127 * 00128 * Based on contributions by 00129 * Beresford Parlett, University of California, Berkeley, USA 00130 * Jim Demmel, University of California, Berkeley, USA 00131 * Inderjit Dhillon, University of Texas, Austin, USA 00132 * Osni Marques, LBNL/NERSC, USA 00133 * Christof Voemel, University of California, Berkeley, USA 00134 * 00135 * ===================================================================== 00136 * 00137 * .. Parameters .. 00138 DOUBLE PRECISION ZERO, ONE 00139 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00140 00141 * .. 00142 * .. Local Scalars .. 00143 LOGICAL SAWNAN1, SAWNAN2 00144 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1, 00145 $ R2 00146 DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP 00147 * .. 00148 * .. External Functions .. 00149 LOGICAL DISNAN 00150 DOUBLE PRECISION DLAMCH 00151 EXTERNAL DISNAN, DLAMCH 00152 * .. 00153 * .. Intrinsic Functions .. 00154 INTRINSIC ABS 00155 * .. 00156 * .. Executable Statements .. 00157 * 00158 EPS = DLAMCH( 'Precision' ) 00159 00160 00161 IF( R.EQ.0 ) THEN 00162 R1 = B1 00163 R2 = BN 00164 ELSE 00165 R1 = R 00166 R2 = R 00167 END IF 00168 00169 * Storage for LPLUS 00170 INDLPL = 0 00171 * Storage for UMINUS 00172 INDUMN = N 00173 INDS = 2*N + 1 00174 INDP = 3*N + 1 00175 00176 IF( B1.EQ.1 ) THEN 00177 WORK( INDS ) = ZERO 00178 ELSE 00179 WORK( INDS+B1-1 ) = LLD( B1-1 ) 00180 END IF 00181 00182 * 00183 * Compute the stationary transform (using the differential form) 00184 * until the index R2. 00185 * 00186 SAWNAN1 = .FALSE. 00187 NEG1 = 0 00188 S = WORK( INDS+B1-1 ) - LAMBDA 00189 DO 50 I = B1, R1 - 1 00190 DPLUS = D( I ) + S 00191 WORK( INDLPL+I ) = LD( I ) / DPLUS 00192 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00193 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00194 S = WORK( INDS+I ) - LAMBDA 00195 50 CONTINUE 00196 SAWNAN1 = DISNAN( S ) 00197 IF( SAWNAN1 ) GOTO 60 00198 DO 51 I = R1, R2 - 1 00199 DPLUS = D( I ) + S 00200 WORK( INDLPL+I ) = LD( I ) / DPLUS 00201 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00202 S = WORK( INDS+I ) - LAMBDA 00203 51 CONTINUE 00204 SAWNAN1 = DISNAN( S ) 00205 * 00206 60 CONTINUE 00207 IF( SAWNAN1 ) THEN 00208 * Runs a slower version of the above loop if a NaN is detected 00209 NEG1 = 0 00210 S = WORK( INDS+B1-1 ) - LAMBDA 00211 DO 70 I = B1, R1 - 1 00212 DPLUS = D( I ) + S 00213 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00214 WORK( INDLPL+I ) = LD( I ) / DPLUS 00215 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00216 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00217 IF( WORK( INDLPL+I ).EQ.ZERO ) 00218 $ WORK( INDS+I ) = LLD( I ) 00219 S = WORK( INDS+I ) - LAMBDA 00220 70 CONTINUE 00221 DO 71 I = R1, R2 - 1 00222 DPLUS = D( I ) + S 00223 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00224 WORK( INDLPL+I ) = LD( I ) / DPLUS 00225 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00226 IF( WORK( INDLPL+I ).EQ.ZERO ) 00227 $ WORK( INDS+I ) = LLD( I ) 00228 S = WORK( INDS+I ) - LAMBDA 00229 71 CONTINUE 00230 END IF 00231 * 00232 * Compute the progressive transform (using the differential form) 00233 * until the index R1 00234 * 00235 SAWNAN2 = .FALSE. 00236 NEG2 = 0 00237 WORK( INDP+BN-1 ) = D( BN ) - LAMBDA 00238 DO 80 I = BN - 1, R1, -1 00239 DMINUS = LLD( I ) + WORK( INDP+I ) 00240 TMP = D( I ) / DMINUS 00241 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00242 WORK( INDUMN+I ) = L( I )*TMP 00243 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00244 80 CONTINUE 00245 TMP = WORK( INDP+R1-1 ) 00246 SAWNAN2 = DISNAN( TMP ) 00247 00248 IF( SAWNAN2 ) THEN 00249 * Runs a slower version of the above loop if a NaN is detected 00250 NEG2 = 0 00251 DO 100 I = BN-1, R1, -1 00252 DMINUS = LLD( I ) + WORK( INDP+I ) 00253 IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN 00254 TMP = D( I ) / DMINUS 00255 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00256 WORK( INDUMN+I ) = L( I )*TMP 00257 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00258 IF( TMP.EQ.ZERO ) 00259 $ WORK( INDP+I-1 ) = D( I ) - LAMBDA 00260 100 CONTINUE 00261 END IF 00262 * 00263 * Find the index (from R1 to R2) of the largest (in magnitude) 00264 * diagonal element of the inverse 00265 * 00266 MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 ) 00267 IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1 00268 IF( WANTNC ) THEN 00269 NEGCNT = NEG1 + NEG2 00270 ELSE 00271 NEGCNT = -1 00272 ENDIF 00273 IF( ABS(MINGMA).EQ.ZERO ) 00274 $ MINGMA = EPS*WORK( INDS+R1-1 ) 00275 R = R1 00276 DO 110 I = R1, R2 - 1 00277 TMP = WORK( INDS+I ) + WORK( INDP+I ) 00278 IF( TMP.EQ.ZERO ) 00279 $ TMP = EPS*WORK( INDS+I ) 00280 IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN 00281 MINGMA = TMP 00282 R = I + 1 00283 END IF 00284 110 CONTINUE 00285 * 00286 * Compute the FP vector: solve N^T v = e_r 00287 * 00288 ISUPPZ( 1 ) = B1 00289 ISUPPZ( 2 ) = BN 00290 Z( R ) = ONE 00291 ZTZ = ONE 00292 * 00293 * Compute the FP vector upwards from R 00294 * 00295 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00296 DO 210 I = R-1, B1, -1 00297 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00298 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00299 $ THEN 00300 Z( I ) = ZERO 00301 ISUPPZ( 1 ) = I + 1 00302 GOTO 220 00303 ENDIF 00304 ZTZ = ZTZ + Z( I )*Z( I ) 00305 210 CONTINUE 00306 220 CONTINUE 00307 ELSE 00308 * Run slower loop if NaN occurred. 00309 DO 230 I = R - 1, B1, -1 00310 IF( Z( I+1 ).EQ.ZERO ) THEN 00311 Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) 00312 ELSE 00313 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00314 END IF 00315 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00316 $ THEN 00317 Z( I ) = ZERO 00318 ISUPPZ( 1 ) = I + 1 00319 GO TO 240 00320 END IF 00321 ZTZ = ZTZ + Z( I )*Z( I ) 00322 230 CONTINUE 00323 240 CONTINUE 00324 ENDIF 00325 00326 * Compute the FP vector downwards from R in blocks of size BLKSIZ 00327 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00328 DO 250 I = R, BN-1 00329 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00330 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00331 $ THEN 00332 Z( I+1 ) = ZERO 00333 ISUPPZ( 2 ) = I 00334 GO TO 260 00335 END IF 00336 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 00337 250 CONTINUE 00338 260 CONTINUE 00339 ELSE 00340 * Run slower loop if NaN occurred. 00341 DO 270 I = R, BN - 1 00342 IF( Z( I ).EQ.ZERO ) THEN 00343 Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) 00344 ELSE 00345 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00346 END IF 00347 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00348 $ THEN 00349 Z( I+1 ) = ZERO 00350 ISUPPZ( 2 ) = I 00351 GO TO 280 00352 END IF 00353 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 00354 270 CONTINUE 00355 280 CONTINUE 00356 END IF 00357 * 00358 * Compute quantities for convergence test 00359 * 00360 TMP = ONE / ZTZ 00361 NRMINV = SQRT( TMP ) 00362 RESID = ABS( MINGMA )*NRMINV 00363 RQCORR = MINGMA*TMP 00364 * 00365 * 00366 RETURN 00367 * 00368 * End of DLAR1V 00369 * 00370 END
1.7.3 
