![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLAR1V 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLAR1V + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clar1v.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clar1v.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clar1v.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLAR1V( 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 * REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 00029 * $ RQCORR, ZTZ 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER ISUPPZ( * ) 00033 * REAL D( * ), L( * ), LD( * ), LLD( * ), 00034 * $ WORK( * ) 00035 * COMPLEX Z( * ) 00036 * .. 00037 * 00038 * 00039 *> \par Purpose: 00040 * ============= 00041 *> 00042 *> \verbatim 00043 *> 00044 *> CLAR1V 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 00118 *> The minimum pivot in the Sturm sequence. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] GAPTOL 00122 *> \verbatim 00123 *> GAPTOL is REAL 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 COMPLEX 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 REAL 00152 *> The square of the 2-norm of Z. 00153 *> \endverbatim 00154 *> 00155 *> \param[out] MINGMA 00156 *> \verbatim 00157 *> MINGMA is REAL 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 REAL 00185 *> NRMINV = 1/SQRT( ZTZ ) 00186 *> \endverbatim 00187 *> 00188 *> \param[out] RESID 00189 *> \verbatim 00190 *> RESID is REAL 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 REAL 00198 *> The Rayleigh Quotient correction to LAMBDA. 00199 *> RQCORR = MINGMA*TMP 00200 *> \endverbatim 00201 *> 00202 *> \param[out] WORK 00203 *> \verbatim 00204 *> WORK is REAL 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 complexOTHERauxiliary 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 CLAR1V( 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 REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 00242 $ RQCORR, ZTZ 00243 * .. 00244 * .. Array Arguments .. 00245 INTEGER ISUPPZ( * ) 00246 REAL D( * ), L( * ), LD( * ), LLD( * ), 00247 $ WORK( * ) 00248 COMPLEX Z( * ) 00249 * .. 00250 * 00251 * ===================================================================== 00252 * 00253 * .. Parameters .. 00254 REAL ZERO, ONE 00255 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00256 COMPLEX CONE 00257 PARAMETER ( CONE = ( 1.0E0, 0.0E0 ) ) 00258 00259 * .. 00260 * .. Local Scalars .. 00261 LOGICAL SAWNAN1, SAWNAN2 00262 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1, 00263 $ R2 00264 REAL DMINUS, DPLUS, EPS, S, TMP 00265 * .. 00266 * .. External Functions .. 00267 LOGICAL SISNAN 00268 REAL SLAMCH 00269 EXTERNAL SISNAN, SLAMCH 00270 * .. 00271 * .. Intrinsic Functions .. 00272 INTRINSIC ABS, REAL 00273 * .. 00274 * .. Executable Statements .. 00275 * 00276 EPS = SLAMCH( 'Precision' ) 00277 00278 00279 IF( R.EQ.0 ) THEN 00280 R1 = B1 00281 R2 = BN 00282 ELSE 00283 R1 = R 00284 R2 = R 00285 END IF 00286 00287 * Storage for LPLUS 00288 INDLPL = 0 00289 * Storage for UMINUS 00290 INDUMN = N 00291 INDS = 2*N + 1 00292 INDP = 3*N + 1 00293 00294 IF( B1.EQ.1 ) THEN 00295 WORK( INDS ) = ZERO 00296 ELSE 00297 WORK( INDS+B1-1 ) = LLD( B1-1 ) 00298 END IF 00299 00300 * 00301 * Compute the stationary transform (using the differential form) 00302 * until the index R2. 00303 * 00304 SAWNAN1 = .FALSE. 00305 NEG1 = 0 00306 S = WORK( INDS+B1-1 ) - LAMBDA 00307 DO 50 I = B1, R1 - 1 00308 DPLUS = D( I ) + S 00309 WORK( INDLPL+I ) = LD( I ) / DPLUS 00310 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00311 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00312 S = WORK( INDS+I ) - LAMBDA 00313 50 CONTINUE 00314 SAWNAN1 = SISNAN( S ) 00315 IF( SAWNAN1 ) GOTO 60 00316 DO 51 I = R1, R2 - 1 00317 DPLUS = D( I ) + S 00318 WORK( INDLPL+I ) = LD( I ) / DPLUS 00319 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00320 S = WORK( INDS+I ) - LAMBDA 00321 51 CONTINUE 00322 SAWNAN1 = SISNAN( S ) 00323 * 00324 60 CONTINUE 00325 IF( SAWNAN1 ) THEN 00326 * Runs a slower version of the above loop if a NaN is detected 00327 NEG1 = 0 00328 S = WORK( INDS+B1-1 ) - LAMBDA 00329 DO 70 I = B1, R1 - 1 00330 DPLUS = D( I ) + S 00331 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00332 WORK( INDLPL+I ) = LD( I ) / DPLUS 00333 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 00334 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00335 IF( WORK( INDLPL+I ).EQ.ZERO ) 00336 $ WORK( INDS+I ) = LLD( I ) 00337 S = WORK( INDS+I ) - LAMBDA 00338 70 CONTINUE 00339 DO 71 I = R1, R2 - 1 00340 DPLUS = D( I ) + S 00341 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 00342 WORK( INDLPL+I ) = LD( I ) / DPLUS 00343 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 00344 IF( WORK( INDLPL+I ).EQ.ZERO ) 00345 $ WORK( INDS+I ) = LLD( I ) 00346 S = WORK( INDS+I ) - LAMBDA 00347 71 CONTINUE 00348 END IF 00349 * 00350 * Compute the progressive transform (using the differential form) 00351 * until the index R1 00352 * 00353 SAWNAN2 = .FALSE. 00354 NEG2 = 0 00355 WORK( INDP+BN-1 ) = D( BN ) - LAMBDA 00356 DO 80 I = BN - 1, R1, -1 00357 DMINUS = LLD( I ) + WORK( INDP+I ) 00358 TMP = D( I ) / DMINUS 00359 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00360 WORK( INDUMN+I ) = L( I )*TMP 00361 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00362 80 CONTINUE 00363 TMP = WORK( INDP+R1-1 ) 00364 SAWNAN2 = SISNAN( TMP ) 00365 00366 IF( SAWNAN2 ) THEN 00367 * Runs a slower version of the above loop if a NaN is detected 00368 NEG2 = 0 00369 DO 100 I = BN-1, R1, -1 00370 DMINUS = LLD( I ) + WORK( INDP+I ) 00371 IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN 00372 TMP = D( I ) / DMINUS 00373 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 00374 WORK( INDUMN+I ) = L( I )*TMP 00375 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 00376 IF( TMP.EQ.ZERO ) 00377 $ WORK( INDP+I-1 ) = D( I ) - LAMBDA 00378 100 CONTINUE 00379 END IF 00380 * 00381 * Find the index (from R1 to R2) of the largest (in magnitude) 00382 * diagonal element of the inverse 00383 * 00384 MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 ) 00385 IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1 00386 IF( WANTNC ) THEN 00387 NEGCNT = NEG1 + NEG2 00388 ELSE 00389 NEGCNT = -1 00390 ENDIF 00391 IF( ABS(MINGMA).EQ.ZERO ) 00392 $ MINGMA = EPS*WORK( INDS+R1-1 ) 00393 R = R1 00394 DO 110 I = R1, R2 - 1 00395 TMP = WORK( INDS+I ) + WORK( INDP+I ) 00396 IF( TMP.EQ.ZERO ) 00397 $ TMP = EPS*WORK( INDS+I ) 00398 IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN 00399 MINGMA = TMP 00400 R = I + 1 00401 END IF 00402 110 CONTINUE 00403 * 00404 * Compute the FP vector: solve N^T v = e_r 00405 * 00406 ISUPPZ( 1 ) = B1 00407 ISUPPZ( 2 ) = BN 00408 Z( R ) = CONE 00409 ZTZ = ONE 00410 * 00411 * Compute the FP vector upwards from R 00412 * 00413 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00414 DO 210 I = R-1, B1, -1 00415 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00416 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00417 $ THEN 00418 Z( I ) = ZERO 00419 ISUPPZ( 1 ) = I + 1 00420 GOTO 220 00421 ENDIF 00422 ZTZ = ZTZ + REAL( Z( I )*Z( I ) ) 00423 210 CONTINUE 00424 220 CONTINUE 00425 ELSE 00426 * Run slower loop if NaN occurred. 00427 DO 230 I = R - 1, B1, -1 00428 IF( Z( I+1 ).EQ.ZERO ) THEN 00429 Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) 00430 ELSE 00431 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 00432 END IF 00433 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00434 $ THEN 00435 Z( I ) = ZERO 00436 ISUPPZ( 1 ) = I + 1 00437 GO TO 240 00438 END IF 00439 ZTZ = ZTZ + REAL( Z( I )*Z( I ) ) 00440 230 CONTINUE 00441 240 CONTINUE 00442 ENDIF 00443 00444 * Compute the FP vector downwards from R in blocks of size BLKSIZ 00445 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 00446 DO 250 I = R, BN-1 00447 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00448 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00449 $ THEN 00450 Z( I+1 ) = ZERO 00451 ISUPPZ( 2 ) = I 00452 GO TO 260 00453 END IF 00454 ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) ) 00455 250 CONTINUE 00456 260 CONTINUE 00457 ELSE 00458 * Run slower loop if NaN occurred. 00459 DO 270 I = R, BN - 1 00460 IF( Z( I ).EQ.ZERO ) THEN 00461 Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) 00462 ELSE 00463 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 00464 END IF 00465 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 00466 $ THEN 00467 Z( I+1 ) = ZERO 00468 ISUPPZ( 2 ) = I 00469 GO TO 280 00470 END IF 00471 ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) ) 00472 270 CONTINUE 00473 280 CONTINUE 00474 END IF 00475 * 00476 * Compute quantities for convergence test 00477 * 00478 TMP = ONE / ZTZ 00479 NRMINV = SQRT( TMP ) 00480 RESID = ABS( MINGMA )*NRMINV 00481 RQCORR = MINGMA*TMP 00482 * 00483 * 00484 RETURN 00485 * 00486 * End of CLAR1V 00487 * 00488 END