 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CLAR1V( 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.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 LOGICAL WANTNC 00012 INTEGER B1, BN, N, NEGCNT, R 00013 REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 00014 $ RQCORR, ZTZ 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER ISUPPZ( * ) 00018 REAL D( * ), L( * ), LD( * ), LLD( * ), 00019 $ WORK( * ) 00020 COMPLEX Z( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * CLAR1V 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) REAL 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) REAL 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) REAL array, dimension (N) 00064 * The n diagonal elements of the diagonal matrix D. 00065 * 00066 * LD (input) REAL array, dimension (N-1) 00067 * The n-1 elements L(i)*D(i). 00068 * 00069 * LLD (input) REAL array, dimension (N-1) 00070 * The n-1 elements L(i)*L(i)*D(i). 00071 * 00072 * PIVMIN (input) REAL 00073 * The minimum pivot in the Sturm sequence. 00074 * 00075 * GAPTOL (input) REAL 00076 * Tolerance that indicates when eigenvector entries are negligible 00077 * w.r.t. their contribution to the residual. 00078 * 00079 * Z (input/output) COMPLEX 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) REAL 00092 * The square of the 2-norm of Z. 00093 * 00094 * MINGMA (output) REAL 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) REAL 00113 * NRMINV = 1/SQRT( ZTZ ) 00114 * 00115 * RESID (output) REAL 00116 * The residual of the FP vector. 00117 * RESID = ABS( MINGMA )/SQRT( ZTZ ) 00118 * 00119 * RQCORR (output) REAL 00120 * The Rayleigh Quotient correction to LAMBDA. 00121 * RQCORR = MINGMA*TMP 00122 * 00123 * WORK (workspace) REAL 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 REAL ZERO, ONE 00139 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00140 COMPLEX CONE 00141 PARAMETER ( CONE = ( 1.0E0, 0.0E0 ) ) 00142 00143 * .. 00144 * .. Local Scalars .. 00145 LOGICAL SAWNAN1, SAWNAN2 00146 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1, 00147 $ R2 00148 REAL DMINUS, DPLUS, EPS, S, TMP 00149 * .. 00150 * .. External Functions .. 00151 LOGICAL SISNAN 00152 REAL SLAMCH 00153 EXTERNAL SISNAN, SLAMCH 00154 * .. 00155 * .. Intrinsic Functions .. 00156 INTRINSIC ABS, REAL 00157 * .. 00158 * .. Executable Statements .. 00159 * 00160 EPS = SLAMCH( 'Precision' ) 00161 00162 00163 IF( R.EQ.0 ) THEN 00164 R1 = B1 00165 R2 = BN 00166 ELSE 00167 R1 = R 00168 R2 = R 00169 END IF 00170 00171 * Storage for LPLUS 00172 INDLPL = 0 00173 * Storage for UMINUS 00174 INDUMN = N 00175 INDS = 2*N + 1 00176 INDP = 3*N + 1 00177 00178 IF( B1.EQ.1 ) THEN 00179 WORK( INDS ) = ZERO 00180 ELSE 00181 WORK( INDS+B1-1 ) = LLD( B1-1 ) 00182 END IF 00183 00184 * 00185 * Compute the stationary transform (using the differential form) 00186 * until the index R2. 00187 * 00188 SAWNAN1 = .FALSE. 00189 NEG1 = 0 00190 S = WORK( INDS+B1-1 ) - LAMBDA 00191 DO 50 I = B1, R1 - 1 00192 DPLUS = D( I ) + S 00193 WORK( INDLPL+I ) = LD( I ) / DPLUS 00194 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00195 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00196 S = WORK( INDS+I ) - LAMBDA 00197 50 CONTINUE 00198 SAWNAN1 = SISNAN( S ) 00199 IF( SAWNAN1 ) GOTO 60 00200 DO 51 I = R1, R2 - 1 00201 DPLUS = D( I ) + S 00202 WORK( INDLPL+I ) = LD( I ) / DPLUS 00203 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00204 S = WORK( INDS+I ) - LAMBDA 00205 51 CONTINUE 00206 SAWNAN1 = SISNAN( S ) 00207 * 00208 60 CONTINUE 00209 IF( SAWNAN1 ) THEN 00210 * Runs a slower version of the above loop if a NaN is detected 00211 NEG1 = 0 00212 S = WORK( INDS+B1-1 ) - LAMBDA 00213 DO 70 I = B1, R1 - 1 00214 DPLUS = D( I ) + S 00215 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00216 WORK( INDLPL+I ) = LD( I ) / DPLUS 00217 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00218 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00219 IF( WORK( INDLPL+I ).EQ.ZERO ) 00220 $ WORK( INDS+I ) = LLD( I ) 00221 S = WORK( INDS+I ) - LAMBDA 00222 70 CONTINUE 00223 DO 71 I = R1, R2 - 1 00224 DPLUS = D( I ) + S 00225 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00226 WORK( INDLPL+I ) = LD( I ) / DPLUS 00227 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00228 IF( WORK( INDLPL+I ).EQ.ZERO ) 00229 $ WORK( INDS+I ) = LLD( I ) 00230 S = WORK( INDS+I ) - LAMBDA 00231 71 CONTINUE 00232 END IF 00233 * 00234 * Compute the progressive transform (using the differential form) 00235 * until the index R1 00236 * 00237 SAWNAN2 = .FALSE. 00238 NEG2 = 0 00239 WORK( INDP+BN-1 ) = D( BN ) - LAMBDA 00240 DO 80 I = BN - 1, R1, -1 00241 DMINUS = LLD( I ) + WORK( INDP+I ) 00242 TMP = D( I ) / DMINUS 00243 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00244 WORK( INDUMN+I ) = L( I )*TMP 00245 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00246 80 CONTINUE 00247 TMP = WORK( INDP+R1-1 ) 00248 SAWNAN2 = SISNAN( TMP ) 00249 00250 IF( SAWNAN2 ) THEN 00251 * Runs a slower version of the above loop if a NaN is detected 00252 NEG2 = 0 00253 DO 100 I = BN-1, R1, -1 00254 DMINUS = LLD( I ) + WORK( INDP+I ) 00255 IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN 00256 TMP = D( I ) / DMINUS 00257 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00258 WORK( INDUMN+I ) = L( I )*TMP 00259 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00260 IF( TMP.EQ.ZERO ) 00261 $ WORK( INDP+I-1 ) = D( I ) - LAMBDA 00262 100 CONTINUE 00263 END IF 00264 * 00265 * Find the index (from R1 to R2) of the largest (in magnitude) 00266 * diagonal element of the inverse 00267 * 00268 MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 ) 00269 IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1 00270 IF( WANTNC ) THEN 00271 NEGCNT = NEG1 + NEG2 00272 ELSE 00273 NEGCNT = -1 00274 ENDIF 00275 IF( ABS(MINGMA).EQ.ZERO ) 00276 $ MINGMA = EPS*WORK( INDS+R1-1 ) 00277 R = R1 00278 DO 110 I = R1, R2 - 1 00279 TMP = WORK( INDS+I ) + WORK( INDP+I ) 00280 IF( TMP.EQ.ZERO ) 00281 $ TMP = EPS*WORK( INDS+I ) 00282 IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN 00283 MINGMA = TMP 00284 R = I + 1 00285 END IF 00286 110 CONTINUE 00287 * 00288 * Compute the FP vector: solve N^T v = e_r 00289 * 00290 ISUPPZ( 1 ) = B1 00291 ISUPPZ( 2 ) = BN 00292 Z( R ) = CONE 00293 ZTZ = ONE 00294 * 00295 * Compute the FP vector upwards from R 00296 * 00297 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00298 DO 210 I = R-1, B1, -1 00299 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00300 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00301 $ THEN 00302 Z( I ) = ZERO 00303 ISUPPZ( 1 ) = I + 1 00304 GOTO 220 00305 ENDIF 00306 ZTZ = ZTZ + REAL( Z( I )*Z( I ) ) 00307 210 CONTINUE 00308 220 CONTINUE 00309 ELSE 00310 * Run slower loop if NaN occurred. 00311 DO 230 I = R - 1, B1, -1 00312 IF( Z( I+1 ).EQ.ZERO ) THEN 00313 Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) 00314 ELSE 00315 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00316 END IF 00317 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00318 $ THEN 00319 Z( I ) = ZERO 00320 ISUPPZ( 1 ) = I + 1 00321 GO TO 240 00322 END IF 00323 ZTZ = ZTZ + REAL( Z( I )*Z( I ) ) 00324 230 CONTINUE 00325 240 CONTINUE 00326 ENDIF 00327 00328 * Compute the FP vector downwards from R in blocks of size BLKSIZ 00329 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00330 DO 250 I = R, BN-1 00331 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00332 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00333 $ THEN 00334 Z( I+1 ) = ZERO 00335 ISUPPZ( 2 ) = I 00336 GO TO 260 00337 END IF 00338 ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) ) 00339 250 CONTINUE 00340 260 CONTINUE 00341 ELSE 00342 * Run slower loop if NaN occurred. 00343 DO 270 I = R, BN - 1 00344 IF( Z( I ).EQ.ZERO ) THEN 00345 Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) 00346 ELSE 00347 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00348 END IF 00349 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00350 $ THEN 00351 Z( I+1 ) = ZERO 00352 ISUPPZ( 2 ) = I 00353 GO TO 280 00354 END IF 00355 ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) ) 00356 270 CONTINUE 00357 280 CONTINUE 00358 END IF 00359 * 00360 * Compute quantities for convergence test 00361 * 00362 TMP = ONE / ZTZ 00363 NRMINV = SQRT( TMP ) 00364 RESID = ABS( MINGMA )*NRMINV 00365 RQCORR = MINGMA*TMP 00366 * 00367 * 00368 RETURN 00369 * 00370 * End of CLAR1V 00371 * 00372 END
1.7.3 
