![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLARRV 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLARRV + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrv.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrv.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrv.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLARRV( N, VL, VU, D, L, PIVMIN, 00022 * ISPLIT, M, DOL, DOU, MINRGP, 00023 * RTOL1, RTOL2, W, WERR, WGAP, 00024 * IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, 00025 * WORK, IWORK, INFO ) 00026 * 00027 * .. Scalar Arguments .. 00028 * INTEGER DOL, DOU, INFO, LDZ, M, N 00029 * REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 00033 * $ ISUPPZ( * ), IWORK( * ) 00034 * REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 00035 * $ WGAP( * ), WORK( * ) 00036 * REAL Z( LDZ, * ) 00037 * .. 00038 * 00039 * 00040 *> \par Purpose: 00041 * ============= 00042 *> 00043 *> \verbatim 00044 *> 00045 *> SLARRV computes the eigenvectors of the tridiagonal matrix 00046 *> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. 00047 *> The input eigenvalues should have been computed by SLARRE. 00048 *> \endverbatim 00049 * 00050 * Arguments: 00051 * ========== 00052 * 00053 *> \param[in] N 00054 *> \verbatim 00055 *> N is INTEGER 00056 *> The order of the matrix. N >= 0. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] VL 00060 *> \verbatim 00061 *> VL is REAL 00062 *> \endverbatim 00063 *> 00064 *> \param[in] VU 00065 *> \verbatim 00066 *> VU is REAL 00067 *> Lower and upper bounds of the interval that contains the desired 00068 *> eigenvalues. VL < VU. Needed to compute gaps on the left or right 00069 *> end of the extremal eigenvalues in the desired RANGE. 00070 *> \endverbatim 00071 *> 00072 *> \param[in,out] D 00073 *> \verbatim 00074 *> D is REAL array, dimension (N) 00075 *> On entry, the N diagonal elements of the diagonal matrix D. 00076 *> On exit, D may be overwritten. 00077 *> \endverbatim 00078 *> 00079 *> \param[in,out] L 00080 *> \verbatim 00081 *> L is REAL array, dimension (N) 00082 *> On entry, the (N-1) subdiagonal elements of the unit 00083 *> bidiagonal matrix L are in elements 1 to N-1 of L 00084 *> (if the matrix is not splitted.) At the end of each block 00085 *> is stored the corresponding shift as given by SLARRE. 00086 *> On exit, L is overwritten. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] PIVMIN 00090 *> \verbatim 00091 *> PIVMIN is REAL 00092 *> The minimum pivot allowed in the Sturm sequence. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] ISPLIT 00096 *> \verbatim 00097 *> ISPLIT is INTEGER array, dimension (N) 00098 *> The splitting points, at which T breaks up into blocks. 00099 *> The first block consists of rows/columns 1 to 00100 *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 00101 *> through ISPLIT( 2 ), etc. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] M 00105 *> \verbatim 00106 *> M is INTEGER 00107 *> The total number of input eigenvalues. 0 <= M <= N. 00108 *> \endverbatim 00109 *> 00110 *> \param[in] DOL 00111 *> \verbatim 00112 *> DOL is INTEGER 00113 *> \endverbatim 00114 *> 00115 *> \param[in] DOU 00116 *> \verbatim 00117 *> DOU is INTEGER 00118 *> If the user wants to compute only selected eigenvectors from all 00119 *> the eigenvalues supplied, he can specify an index range DOL:DOU. 00120 *> Or else the setting DOL=1, DOU=M should be applied. 00121 *> Note that DOL and DOU refer to the order in which the eigenvalues 00122 *> are stored in W. 00123 *> If the user wants to compute only selected eigenpairs, then 00124 *> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the 00125 *> computed eigenvectors. All other columns of Z are set to zero. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] MINRGP 00129 *> \verbatim 00130 *> MINRGP is REAL 00131 *> \endverbatim 00132 *> 00133 *> \param[in] RTOL1 00134 *> \verbatim 00135 *> RTOL1 is REAL 00136 *> \endverbatim 00137 *> 00138 *> \param[in] RTOL2 00139 *> \verbatim 00140 *> RTOL2 is REAL 00141 *> Parameters for bisection. 00142 *> An interval [LEFT,RIGHT] has converged if 00143 *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) 00144 *> \endverbatim 00145 *> 00146 *> \param[in,out] W 00147 *> \verbatim 00148 *> W is REAL array, dimension (N) 00149 *> The first M elements of W contain the APPROXIMATE eigenvalues for 00150 *> which eigenvectors are to be computed. The eigenvalues 00151 *> should be grouped by split-off block and ordered from 00152 *> smallest to largest within the block ( The output array 00153 *> W from SLARRE is expected here ). Furthermore, they are with 00154 *> respect to the shift of the corresponding root representation 00155 *> for their block. On exit, W holds the eigenvalues of the 00156 *> UNshifted matrix. 00157 *> \endverbatim 00158 *> 00159 *> \param[in,out] WERR 00160 *> \verbatim 00161 *> WERR is REAL array, dimension (N) 00162 *> The first M elements contain the semiwidth of the uncertainty 00163 *> interval of the corresponding eigenvalue in W 00164 *> \endverbatim 00165 *> 00166 *> \param[in,out] WGAP 00167 *> \verbatim 00168 *> WGAP is REAL array, dimension (N) 00169 *> The separation from the right neighbor eigenvalue in W. 00170 *> \endverbatim 00171 *> 00172 *> \param[in] IBLOCK 00173 *> \verbatim 00174 *> IBLOCK is INTEGER array, dimension (N) 00175 *> The indices of the blocks (submatrices) associated with the 00176 *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue 00177 *> W(i) belongs to the first block from the top, =2 if W(i) 00178 *> belongs to the second block, etc. 00179 *> \endverbatim 00180 *> 00181 *> \param[in] INDEXW 00182 *> \verbatim 00183 *> INDEXW is INTEGER array, dimension (N) 00184 *> The indices of the eigenvalues within each block (submatrix); 00185 *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the 00186 *> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. 00187 *> \endverbatim 00188 *> 00189 *> \param[in] GERS 00190 *> \verbatim 00191 *> GERS is REAL array, dimension (2*N) 00192 *> The N Gerschgorin intervals (the i-th Gerschgorin interval 00193 *> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should 00194 *> be computed from the original UNshifted matrix. 00195 *> \endverbatim 00196 *> 00197 *> \param[out] Z 00198 *> \verbatim 00199 *> Z is REAL array, dimension (LDZ, max(1,M) ) 00200 *> If INFO = 0, the first M columns of Z contain the 00201 *> orthonormal eigenvectors of the matrix T 00202 *> corresponding to the input eigenvalues, with the i-th 00203 *> column of Z holding the eigenvector associated with W(i). 00204 *> Note: the user must ensure that at least max(1,M) columns are 00205 *> supplied in the array Z. 00206 *> \endverbatim 00207 *> 00208 *> \param[in] LDZ 00209 *> \verbatim 00210 *> LDZ is INTEGER 00211 *> The leading dimension of the array Z. LDZ >= 1, and if 00212 *> JOBZ = 'V', LDZ >= max(1,N). 00213 *> \endverbatim 00214 *> 00215 *> \param[out] ISUPPZ 00216 *> \verbatim 00217 *> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) 00218 *> The support of the eigenvectors in Z, i.e., the indices 00219 *> indicating the nonzero elements in Z. The I-th eigenvector 00220 *> is nonzero only in elements ISUPPZ( 2*I-1 ) through 00221 *> ISUPPZ( 2*I ). 00222 *> \endverbatim 00223 *> 00224 *> \param[out] WORK 00225 *> \verbatim 00226 *> WORK is REAL array, dimension (12*N) 00227 *> \endverbatim 00228 *> 00229 *> \param[out] IWORK 00230 *> \verbatim 00231 *> IWORK is INTEGER array, dimension (7*N) 00232 *> \endverbatim 00233 *> 00234 *> \param[out] INFO 00235 *> \verbatim 00236 *> INFO is INTEGER 00237 *> = 0: successful exit 00238 *> 00239 *> > 0: A problem occured in SLARRV. 00240 *> < 0: One of the called subroutines signaled an internal problem. 00241 *> Needs inspection of the corresponding parameter IINFO 00242 *> for further information. 00243 *> 00244 *> =-1: Problem in SLARRB when refining a child's eigenvalues. 00245 *> =-2: Problem in SLARRF when computing the RRR of a child. 00246 *> When a child is inside a tight cluster, it can be difficult 00247 *> to find an RRR. A partial remedy from the user's point of 00248 *> view is to make the parameter MINRGP smaller and recompile. 00249 *> However, as the orthogonality of the computed vectors is 00250 *> proportional to 1/MINRGP, the user should be aware that 00251 *> he might be trading in precision when he decreases MINRGP. 00252 *> =-3: Problem in SLARRB when refining a single eigenvalue 00253 *> after the Rayleigh correction was rejected. 00254 *> = 5: The Rayleigh Quotient Iteration failed to converge to 00255 *> full accuracy in MAXITR steps. 00256 *> \endverbatim 00257 * 00258 * Authors: 00259 * ======== 00260 * 00261 *> \author Univ. of Tennessee 00262 *> \author Univ. of California Berkeley 00263 *> \author Univ. of Colorado Denver 00264 *> \author NAG Ltd. 00265 * 00266 *> \date November 2011 00267 * 00268 *> \ingroup realOTHERauxiliary 00269 * 00270 *> \par Contributors: 00271 * ================== 00272 *> 00273 *> Beresford Parlett, University of California, Berkeley, USA \n 00274 *> Jim Demmel, University of California, Berkeley, USA \n 00275 *> Inderjit Dhillon, University of Texas, Austin, USA \n 00276 *> Osni Marques, LBNL/NERSC, USA \n 00277 *> Christof Voemel, University of California, Berkeley, USA 00278 * 00279 * ===================================================================== 00280 SUBROUTINE SLARRV( N, VL, VU, D, L, PIVMIN, 00281 $ ISPLIT, M, DOL, DOU, MINRGP, 00282 $ RTOL1, RTOL2, W, WERR, WGAP, 00283 $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, 00284 $ WORK, IWORK, INFO ) 00285 * 00286 * -- LAPACK auxiliary routine (version 3.4.0) -- 00287 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00288 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00289 * November 2011 00290 * 00291 * .. Scalar Arguments .. 00292 INTEGER DOL, DOU, INFO, LDZ, M, N 00293 REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 00294 * .. 00295 * .. Array Arguments .. 00296 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 00297 $ ISUPPZ( * ), IWORK( * ) 00298 REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 00299 $ WGAP( * ), WORK( * ) 00300 REAL Z( LDZ, * ) 00301 * .. 00302 * 00303 * ===================================================================== 00304 * 00305 * .. Parameters .. 00306 INTEGER MAXITR 00307 PARAMETER ( MAXITR = 10 ) 00308 REAL ZERO, ONE, TWO, THREE, FOUR, HALF 00309 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, 00310 $ TWO = 2.0E0, THREE = 3.0E0, 00311 $ FOUR = 4.0E0, HALF = 0.5E0) 00312 * .. 00313 * .. Local Scalars .. 00314 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ 00315 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1, 00316 $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG, 00317 $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER, 00318 $ ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS, 00319 $ NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST, 00320 $ NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST, 00321 $ OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX, 00322 $ WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU, 00323 $ ZUSEDW 00324 REAL BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU, 00325 $ LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID, 00326 $ RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF, 00327 $ SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ 00328 * .. 00329 * .. External Functions .. 00330 REAL SLAMCH 00331 EXTERNAL SLAMCH 00332 * .. 00333 * .. External Subroutines .. 00334 EXTERNAL SCOPY, SLAR1V, SLARRB, SLARRF, SLASET, 00335 $ SSCAL 00336 * .. 00337 * .. Intrinsic Functions .. 00338 INTRINSIC ABS, REAL, MAX, MIN 00339 * .. 00340 * .. Executable Statements .. 00341 * .. 00342 00343 * The first N entries of WORK are reserved for the eigenvalues 00344 INDLD = N+1 00345 INDLLD= 2*N+1 00346 INDWRK= 3*N+1 00347 MINWSIZE = 12 * N 00348 00349 DO 5 I= 1,MINWSIZE 00350 WORK( I ) = ZERO 00351 5 CONTINUE 00352 00353 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the 00354 * factorization used to compute the FP vector 00355 IINDR = 0 00356 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current 00357 * layer and the one above. 00358 IINDC1 = N 00359 IINDC2 = 2*N 00360 IINDWK = 3*N + 1 00361 00362 MINIWSIZE = 7 * N 00363 DO 10 I= 1,MINIWSIZE 00364 IWORK( I ) = 0 00365 10 CONTINUE 00366 00367 ZUSEDL = 1 00368 IF(DOL.GT.1) THEN 00369 * Set lower bound for use of Z 00370 ZUSEDL = DOL-1 00371 ENDIF 00372 ZUSEDU = M 00373 IF(DOU.LT.M) THEN 00374 * Set lower bound for use of Z 00375 ZUSEDU = DOU+1 00376 ENDIF 00377 * The width of the part of Z that is used 00378 ZUSEDW = ZUSEDU - ZUSEDL + 1 00379 00380 00381 CALL SLASET( 'Full', N, ZUSEDW, ZERO, ZERO, 00382 $ Z(1,ZUSEDL), LDZ ) 00383 00384 EPS = SLAMCH( 'Precision' ) 00385 RQTOL = TWO * EPS 00386 * 00387 * Set expert flags for standard code. 00388 TRYRQC = .TRUE. 00389 00390 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00391 ELSE 00392 * Only selected eigenpairs are computed. Since the other evalues 00393 * are not refined by RQ iteration, bisection has to compute to full 00394 * accuracy. 00395 RTOL1 = FOUR * EPS 00396 RTOL2 = FOUR * EPS 00397 ENDIF 00398 00399 * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the 00400 * desired eigenvalues. The support of the nonzero eigenvector 00401 * entries is contained in the interval IBEGIN:IEND. 00402 * Remark that if k eigenpairs are desired, then the eigenvectors 00403 * are stored in k contiguous columns of Z. 00404 00405 * DONE is the number of eigenvectors already computed 00406 DONE = 0 00407 IBEGIN = 1 00408 WBEGIN = 1 00409 DO 170 JBLK = 1, IBLOCK( M ) 00410 IEND = ISPLIT( JBLK ) 00411 SIGMA = L( IEND ) 00412 * Find the eigenvectors of the submatrix indexed IBEGIN 00413 * through IEND. 00414 WEND = WBEGIN - 1 00415 15 CONTINUE 00416 IF( WEND.LT.M ) THEN 00417 IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN 00418 WEND = WEND + 1 00419 GO TO 15 00420 END IF 00421 END IF 00422 IF( WEND.LT.WBEGIN ) THEN 00423 IBEGIN = IEND + 1 00424 GO TO 170 00425 ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN 00426 IBEGIN = IEND + 1 00427 WBEGIN = WEND + 1 00428 GO TO 170 00429 END IF 00430 00431 * Find local spectral diameter of the block 00432 GL = GERS( 2*IBEGIN-1 ) 00433 GU = GERS( 2*IBEGIN ) 00434 DO 20 I = IBEGIN+1 , IEND 00435 GL = MIN( GERS( 2*I-1 ), GL ) 00436 GU = MAX( GERS( 2*I ), GU ) 00437 20 CONTINUE 00438 SPDIAM = GU - GL 00439 00440 * OLDIEN is the last index of the previous block 00441 OLDIEN = IBEGIN - 1 00442 * Calculate the size of the current block 00443 IN = IEND - IBEGIN + 1 00444 * The number of eigenvalues in the current block 00445 IM = WEND - WBEGIN + 1 00446 00447 * This is for a 1x1 block 00448 IF( IBEGIN.EQ.IEND ) THEN 00449 DONE = DONE+1 00450 Z( IBEGIN, WBEGIN ) = ONE 00451 ISUPPZ( 2*WBEGIN-1 ) = IBEGIN 00452 ISUPPZ( 2*WBEGIN ) = IBEGIN 00453 W( WBEGIN ) = W( WBEGIN ) + SIGMA 00454 WORK( WBEGIN ) = W( WBEGIN ) 00455 IBEGIN = IEND + 1 00456 WBEGIN = WBEGIN + 1 00457 GO TO 170 00458 END IF 00459 00460 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) 00461 * Note that these can be approximations, in this case, the corresp. 00462 * entries of WERR give the size of the uncertainty interval. 00463 * The eigenvalue approximations will be refined when necessary as 00464 * high relative accuracy is required for the computation of the 00465 * corresponding eigenvectors. 00466 CALL SCOPY( IM, W( WBEGIN ), 1, 00467 $ WORK( WBEGIN ), 1 ) 00468 00469 * We store in W the eigenvalue approximations w.r.t. the original 00470 * matrix T. 00471 DO 30 I=1,IM 00472 W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA 00473 30 CONTINUE 00474 00475 00476 * NDEPTH is the current depth of the representation tree 00477 NDEPTH = 0 00478 * PARITY is either 1 or 0 00479 PARITY = 1 00480 * NCLUS is the number of clusters for the next level of the 00481 * representation tree, we start with NCLUS = 1 for the root 00482 NCLUS = 1 00483 IWORK( IINDC1+1 ) = 1 00484 IWORK( IINDC1+2 ) = IM 00485 00486 * IDONE is the number of eigenvectors already computed in the current 00487 * block 00488 IDONE = 0 00489 * loop while( IDONE.LT.IM ) 00490 * generate the representation tree for the current block and 00491 * compute the eigenvectors 00492 40 CONTINUE 00493 IF( IDONE.LT.IM ) THEN 00494 * This is a crude protection against infinitely deep trees 00495 IF( NDEPTH.GT.M ) THEN 00496 INFO = -2 00497 RETURN 00498 ENDIF 00499 * breadth first processing of the current level of the representation 00500 * tree: OLDNCL = number of clusters on current level 00501 OLDNCL = NCLUS 00502 * reset NCLUS to count the number of child clusters 00503 NCLUS = 0 00504 * 00505 PARITY = 1 - PARITY 00506 IF( PARITY.EQ.0 ) THEN 00507 OLDCLS = IINDC1 00508 NEWCLS = IINDC2 00509 ELSE 00510 OLDCLS = IINDC2 00511 NEWCLS = IINDC1 00512 END IF 00513 * Process the clusters on the current level 00514 DO 150 I = 1, OLDNCL 00515 J = OLDCLS + 2*I 00516 * OLDFST, OLDLST = first, last index of current cluster. 00517 * cluster indices start with 1 and are relative 00518 * to WBEGIN when accessing W, WGAP, WERR, Z 00519 OLDFST = IWORK( J-1 ) 00520 OLDLST = IWORK( J ) 00521 IF( NDEPTH.GT.0 ) THEN 00522 * Retrieve relatively robust representation (RRR) of cluster 00523 * that has been computed at the previous level 00524 * The RRR is stored in Z and overwritten once the eigenvectors 00525 * have been computed or when the cluster is refined 00526 00527 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00528 * Get representation from location of the leftmost evalue 00529 * of the cluster 00530 J = WBEGIN + OLDFST - 1 00531 ELSE 00532 IF(WBEGIN+OLDFST-1.LT.DOL) THEN 00533 * Get representation from the left end of Z array 00534 J = DOL - 1 00535 ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN 00536 * Get representation from the right end of Z array 00537 J = DOU 00538 ELSE 00539 J = WBEGIN + OLDFST - 1 00540 ENDIF 00541 ENDIF 00542 CALL SCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 ) 00543 CALL SCOPY( IN-1, Z( IBEGIN, J+1 ), 1, L( IBEGIN ), 00544 $ 1 ) 00545 SIGMA = Z( IEND, J+1 ) 00546 00547 * Set the corresponding entries in Z to zero 00548 CALL SLASET( 'Full', IN, 2, ZERO, ZERO, 00549 $ Z( IBEGIN, J), LDZ ) 00550 END IF 00551 00552 * Compute DL and DLL of current RRR 00553 DO 50 J = IBEGIN, IEND-1 00554 TMP = D( J )*L( J ) 00555 WORK( INDLD-1+J ) = TMP 00556 WORK( INDLLD-1+J ) = TMP*L( J ) 00557 50 CONTINUE 00558 00559 IF( NDEPTH.GT.0 ) THEN 00560 * P and Q are index of the first and last eigenvalue to compute 00561 * within the current block 00562 P = INDEXW( WBEGIN-1+OLDFST ) 00563 Q = INDEXW( WBEGIN-1+OLDLST ) 00564 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET 00565 * through the Q-OFFSET elements of these arrays are to be used. 00566 * OFFSET = P-OLDFST 00567 OFFSET = INDEXW( WBEGIN ) - 1 00568 * perform limited bisection (if necessary) to get approximate 00569 * eigenvalues to the precision needed. 00570 CALL SLARRB( IN, D( IBEGIN ), 00571 $ WORK(INDLLD+IBEGIN-1), 00572 $ P, Q, RTOL1, RTOL2, OFFSET, 00573 $ WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN), 00574 $ WORK( INDWRK ), IWORK( IINDWK ), 00575 $ PIVMIN, SPDIAM, IN, IINFO ) 00576 IF( IINFO.NE.0 ) THEN 00577 INFO = -1 00578 RETURN 00579 ENDIF 00580 * We also recompute the extremal gaps. W holds all eigenvalues 00581 * of the unshifted matrix and must be used for computation 00582 * of WGAP, the entries of WORK might stem from RRRs with 00583 * different shifts. The gaps from WBEGIN-1+OLDFST to 00584 * WBEGIN-1+OLDLST are correctly computed in SLARRB. 00585 * However, we only allow the gaps to become greater since 00586 * this is what should happen when we decrease WERR 00587 IF( OLDFST.GT.1) THEN 00588 WGAP( WBEGIN+OLDFST-2 ) = 00589 $ MAX(WGAP(WBEGIN+OLDFST-2), 00590 $ W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1) 00591 $ - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) ) 00592 ENDIF 00593 IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN 00594 WGAP( WBEGIN+OLDLST-1 ) = 00595 $ MAX(WGAP(WBEGIN+OLDLST-1), 00596 $ W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST) 00597 $ - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) ) 00598 ENDIF 00599 * Each time the eigenvalues in WORK get refined, we store 00600 * the newly found approximation with all shifts applied in W 00601 DO 53 J=OLDFST,OLDLST 00602 W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA 00603 53 CONTINUE 00604 END IF 00605 00606 * Process the current node. 00607 NEWFST = OLDFST 00608 DO 140 J = OLDFST, OLDLST 00609 IF( J.EQ.OLDLST ) THEN 00610 * we are at the right end of the cluster, this is also the 00611 * boundary of the child cluster 00612 NEWLST = J 00613 ELSE IF ( WGAP( WBEGIN + J -1).GE. 00614 $ MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN 00615 * the right relative gap is big enough, the child cluster 00616 * (NEWFST,..,NEWLST) is well separated from the following 00617 NEWLST = J 00618 ELSE 00619 * inside a child cluster, the relative gap is not 00620 * big enough. 00621 GOTO 140 00622 END IF 00623 00624 * Compute size of child cluster found 00625 NEWSIZ = NEWLST - NEWFST + 1 00626 00627 * NEWFTT is the place in Z where the new RRR or the computed 00628 * eigenvector is to be stored 00629 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00630 * Store representation at location of the leftmost evalue 00631 * of the cluster 00632 NEWFTT = WBEGIN + NEWFST - 1 00633 ELSE 00634 IF(WBEGIN+NEWFST-1.LT.DOL) THEN 00635 * Store representation at the left end of Z array 00636 NEWFTT = DOL - 1 00637 ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN 00638 * Store representation at the right end of Z array 00639 NEWFTT = DOU 00640 ELSE 00641 NEWFTT = WBEGIN + NEWFST - 1 00642 ENDIF 00643 ENDIF 00644 00645 IF( NEWSIZ.GT.1) THEN 00646 * 00647 * Current child is not a singleton but a cluster. 00648 * Compute and store new representation of child. 00649 * 00650 * 00651 * Compute left and right cluster gap. 00652 * 00653 * LGAP and RGAP are not computed from WORK because 00654 * the eigenvalue approximations may stem from RRRs 00655 * different shifts. However, W hold all eigenvalues 00656 * of the unshifted matrix. Still, the entries in WGAP 00657 * have to be computed from WORK since the entries 00658 * in W might be of the same order so that gaps are not 00659 * exhibited correctly for very close eigenvalues. 00660 IF( NEWFST.EQ.1 ) THEN 00661 LGAP = MAX( ZERO, 00662 $ W(WBEGIN)-WERR(WBEGIN) - VL ) 00663 ELSE 00664 LGAP = WGAP( WBEGIN+NEWFST-2 ) 00665 ENDIF 00666 RGAP = WGAP( WBEGIN+NEWLST-1 ) 00667 * 00668 * Compute left- and rightmost eigenvalue of child 00669 * to high precision in order to shift as close 00670 * as possible and obtain as large relative gaps 00671 * as possible 00672 * 00673 DO 55 K =1,2 00674 IF(K.EQ.1) THEN 00675 P = INDEXW( WBEGIN-1+NEWFST ) 00676 ELSE 00677 P = INDEXW( WBEGIN-1+NEWLST ) 00678 ENDIF 00679 OFFSET = INDEXW( WBEGIN ) - 1 00680 CALL SLARRB( IN, D(IBEGIN), 00681 $ WORK( INDLLD+IBEGIN-1 ),P,P, 00682 $ RQTOL, RQTOL, OFFSET, 00683 $ WORK(WBEGIN),WGAP(WBEGIN), 00684 $ WERR(WBEGIN),WORK( INDWRK ), 00685 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 00686 $ IN, IINFO ) 00687 55 CONTINUE 00688 * 00689 IF((WBEGIN+NEWLST-1.LT.DOL).OR. 00690 $ (WBEGIN+NEWFST-1.GT.DOU)) THEN 00691 * if the cluster contains no desired eigenvalues 00692 * skip the computation of that branch of the rep. tree 00693 * 00694 * We could skip before the refinement of the extremal 00695 * eigenvalues of the child, but then the representation 00696 * tree could be different from the one when nothing is 00697 * skipped. For this reason we skip at this place. 00698 IDONE = IDONE + NEWLST - NEWFST + 1 00699 GOTO 139 00700 ENDIF 00701 * 00702 * Compute RRR of child cluster. 00703 * Note that the new RRR is stored in Z 00704 * 00705 * SLARRF needs LWORK = 2*N 00706 CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ), 00707 $ WORK(INDLD+IBEGIN-1), 00708 $ NEWFST, NEWLST, WORK(WBEGIN), 00709 $ WGAP(WBEGIN), WERR(WBEGIN), 00710 $ SPDIAM, LGAP, RGAP, PIVMIN, TAU, 00711 $ Z(IBEGIN, NEWFTT),Z(IBEGIN, NEWFTT+1), 00712 $ WORK( INDWRK ), IINFO ) 00713 IF( IINFO.EQ.0 ) THEN 00714 * a new RRR for the cluster was found by SLARRF 00715 * update shift and store it 00716 SSIGMA = SIGMA + TAU 00717 Z( IEND, NEWFTT+1 ) = SSIGMA 00718 * WORK() are the midpoints and WERR() the semi-width 00719 * Note that the entries in W are unchanged. 00720 DO 116 K = NEWFST, NEWLST 00721 FUDGE = 00722 $ THREE*EPS*ABS(WORK(WBEGIN+K-1)) 00723 WORK( WBEGIN + K - 1 ) = 00724 $ WORK( WBEGIN + K - 1) - TAU 00725 FUDGE = FUDGE + 00726 $ FOUR*EPS*ABS(WORK(WBEGIN+K-1)) 00727 * Fudge errors 00728 WERR( WBEGIN + K - 1 ) = 00729 $ WERR( WBEGIN + K - 1 ) + FUDGE 00730 * Gaps are not fudged. Provided that WERR is small 00731 * when eigenvalues are close, a zero gap indicates 00732 * that a new representation is needed for resolving 00733 * the cluster. A fudge could lead to a wrong decision 00734 * of judging eigenvalues 'separated' which in 00735 * reality are not. This could have a negative impact 00736 * on the orthogonality of the computed eigenvectors. 00737 116 CONTINUE 00738 00739 NCLUS = NCLUS + 1 00740 K = NEWCLS + 2*NCLUS 00741 IWORK( K-1 ) = NEWFST 00742 IWORK( K ) = NEWLST 00743 ELSE 00744 INFO = -2 00745 RETURN 00746 ENDIF 00747 ELSE 00748 * 00749 * Compute eigenvector of singleton 00750 * 00751 ITER = 0 00752 * 00753 TOL = FOUR * LOG(REAL(IN)) * EPS 00754 * 00755 K = NEWFST 00756 WINDEX = WBEGIN + K - 1 00757 WINDMN = MAX(WINDEX - 1,1) 00758 WINDPL = MIN(WINDEX + 1,M) 00759 LAMBDA = WORK( WINDEX ) 00760 DONE = DONE + 1 00761 * Check if eigenvector computation is to be skipped 00762 IF((WINDEX.LT.DOL).OR. 00763 $ (WINDEX.GT.DOU)) THEN 00764 ESKIP = .TRUE. 00765 GOTO 125 00766 ELSE 00767 ESKIP = .FALSE. 00768 ENDIF 00769 LEFT = WORK( WINDEX ) - WERR( WINDEX ) 00770 RIGHT = WORK( WINDEX ) + WERR( WINDEX ) 00771 INDEIG = INDEXW( WINDEX ) 00772 * Note that since we compute the eigenpairs for a child, 00773 * all eigenvalue approximations are w.r.t the same shift. 00774 * In this case, the entries in WORK should be used for 00775 * computing the gaps since they exhibit even very small 00776 * differences in the eigenvalues, as opposed to the 00777 * entries in W which might "look" the same. 00778 00779 IF( K .EQ. 1) THEN 00780 * In the case RANGE='I' and with not much initial 00781 * accuracy in LAMBDA and VL, the formula 00782 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) 00783 * can lead to an overestimation of the left gap and 00784 * thus to inadequately early RQI 'convergence'. 00785 * Prevent this by forcing a small left gap. 00786 LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 00787 ELSE 00788 LGAP = WGAP(WINDMN) 00789 ENDIF 00790 IF( K .EQ. IM) THEN 00791 * In the case RANGE='I' and with not much initial 00792 * accuracy in LAMBDA and VU, the formula 00793 * can lead to an overestimation of the right gap and 00794 * thus to inadequately early RQI 'convergence'. 00795 * Prevent this by forcing a small right gap. 00796 RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 00797 ELSE 00798 RGAP = WGAP(WINDEX) 00799 ENDIF 00800 GAP = MIN( LGAP, RGAP ) 00801 IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN 00802 * The eigenvector support can become wrong 00803 * because significant entries could be cut off due to a 00804 * large GAPTOL parameter in LAR1V. Prevent this. 00805 GAPTOL = ZERO 00806 ELSE 00807 GAPTOL = GAP * EPS 00808 ENDIF 00809 ISUPMN = IN 00810 ISUPMX = 1 00811 * Update WGAP so that it holds the minimum gap 00812 * to the left or the right. This is crucial in the 00813 * case where bisection is used to ensure that the 00814 * eigenvalue is refined up to the required precision. 00815 * The correct value is restored afterwards. 00816 SAVGAP = WGAP(WINDEX) 00817 WGAP(WINDEX) = GAP 00818 * We want to use the Rayleigh Quotient Correction 00819 * as often as possible since it converges quadratically 00820 * when we are close enough to the desired eigenvalue. 00821 * However, the Rayleigh Quotient can have the wrong sign 00822 * and lead us away from the desired eigenvalue. In this 00823 * case, the best we can do is to use bisection. 00824 USEDBS = .FALSE. 00825 USEDRQ = .FALSE. 00826 * Bisection is initially turned off unless it is forced 00827 NEEDBS = .NOT.TRYRQC 00828 120 CONTINUE 00829 * Check if bisection should be used to refine eigenvalue 00830 IF(NEEDBS) THEN 00831 * Take the bisection as new iterate 00832 USEDBS = .TRUE. 00833 ITMP1 = IWORK( IINDR+WINDEX ) 00834 OFFSET = INDEXW( WBEGIN ) - 1 00835 CALL SLARRB( IN, D(IBEGIN), 00836 $ WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG, 00837 $ ZERO, TWO*EPS, OFFSET, 00838 $ WORK(WBEGIN),WGAP(WBEGIN), 00839 $ WERR(WBEGIN),WORK( INDWRK ), 00840 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 00841 $ ITMP1, IINFO ) 00842 IF( IINFO.NE.0 ) THEN 00843 INFO = -3 00844 RETURN 00845 ENDIF 00846 LAMBDA = WORK( WINDEX ) 00847 * Reset twist index from inaccurate LAMBDA to 00848 * force computation of true MINGMA 00849 IWORK( IINDR+WINDEX ) = 0 00850 ENDIF 00851 * Given LAMBDA, compute the eigenvector. 00852 CALL SLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 00853 $ L( IBEGIN ), WORK(INDLD+IBEGIN-1), 00854 $ WORK(INDLLD+IBEGIN-1), 00855 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 00856 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 00857 $ IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ), 00858 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 00859 IF(ITER .EQ. 0) THEN 00860 BSTRES = RESID 00861 BSTW = LAMBDA 00862 ELSEIF(RESID.LT.BSTRES) THEN 00863 BSTRES = RESID 00864 BSTW = LAMBDA 00865 ENDIF 00866 ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 )) 00867 ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX )) 00868 ITER = ITER + 1 00869 00870 * sin alpha <= |resid|/gap 00871 * Note that both the residual and the gap are 00872 * proportional to the matrix, so ||T|| doesn't play 00873 * a role in the quotient 00874 00875 * 00876 * Convergence test for Rayleigh-Quotient iteration 00877 * (omitted when Bisection has been used) 00878 * 00879 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 00880 $ RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS) 00881 $ THEN 00882 * We need to check that the RQCORR update doesn't 00883 * move the eigenvalue away from the desired one and 00884 * towards a neighbor. -> protection with bisection 00885 IF(INDEIG.LE.NEGCNT) THEN 00886 * The wanted eigenvalue lies to the left 00887 SGNDEF = -ONE 00888 ELSE 00889 * The wanted eigenvalue lies to the right 00890 SGNDEF = ONE 00891 ENDIF 00892 * We only use the RQCORR if it improves the 00893 * the iterate reasonably. 00894 IF( ( RQCORR*SGNDEF.GE.ZERO ) 00895 $ .AND.( LAMBDA + RQCORR.LE. RIGHT) 00896 $ .AND.( LAMBDA + RQCORR.GE. LEFT) 00897 $ ) THEN 00898 USEDRQ = .TRUE. 00899 * Store new midpoint of bisection interval in WORK 00900 IF(SGNDEF.EQ.ONE) THEN 00901 * The current LAMBDA is on the left of the true 00902 * eigenvalue 00903 LEFT = LAMBDA 00904 * We prefer to assume that the error estimate 00905 * is correct. We could make the interval not 00906 * as a bracket but to be modified if the RQCORR 00907 * chooses to. In this case, the RIGHT side should 00908 * be modified as follows: 00909 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR) 00910 ELSE 00911 * The current LAMBDA is on the right of the true 00912 * eigenvalue 00913 RIGHT = LAMBDA 00914 * See comment about assuming the error estimate is 00915 * correct above. 00916 * LEFT = MIN(LEFT, LAMBDA + RQCORR) 00917 ENDIF 00918 WORK( WINDEX ) = 00919 $ HALF * (RIGHT + LEFT) 00920 * Take RQCORR since it has the correct sign and 00921 * improves the iterate reasonably 00922 LAMBDA = LAMBDA + RQCORR 00923 * Update width of error interval 00924 WERR( WINDEX ) = 00925 $ HALF * (RIGHT-LEFT) 00926 ELSE 00927 NEEDBS = .TRUE. 00928 ENDIF 00929 IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN 00930 * The eigenvalue is computed to bisection accuracy 00931 * compute eigenvector and stop 00932 USEDBS = .TRUE. 00933 GOTO 120 00934 ELSEIF( ITER.LT.MAXITR ) THEN 00935 GOTO 120 00936 ELSEIF( ITER.EQ.MAXITR ) THEN 00937 NEEDBS = .TRUE. 00938 GOTO 120 00939 ELSE 00940 INFO = 5 00941 RETURN 00942 END IF 00943 ELSE 00944 STP2II = .FALSE. 00945 IF(USEDRQ .AND. USEDBS .AND. 00946 $ BSTRES.LE.RESID) THEN 00947 LAMBDA = BSTW 00948 STP2II = .TRUE. 00949 ENDIF 00950 IF (STP2II) THEN 00951 * improve error angle by second step 00952 CALL SLAR1V( IN, 1, IN, LAMBDA, 00953 $ D( IBEGIN ), L( IBEGIN ), 00954 $ WORK(INDLD+IBEGIN-1), 00955 $ WORK(INDLLD+IBEGIN-1), 00956 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 00957 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 00958 $ IWORK( IINDR+WINDEX ), 00959 $ ISUPPZ( 2*WINDEX-1 ), 00960 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 00961 ENDIF 00962 WORK( WINDEX ) = LAMBDA 00963 END IF 00964 * 00965 * Compute FP-vector support w.r.t. whole matrix 00966 * 00967 ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN 00968 ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN 00969 ZFROM = ISUPPZ( 2*WINDEX-1 ) 00970 ZTO = ISUPPZ( 2*WINDEX ) 00971 ISUPMN = ISUPMN + OLDIEN 00972 ISUPMX = ISUPMX + OLDIEN 00973 * Ensure vector is ok if support in the RQI has changed 00974 IF(ISUPMN.LT.ZFROM) THEN 00975 DO 122 II = ISUPMN,ZFROM-1 00976 Z( II, WINDEX ) = ZERO 00977 122 CONTINUE 00978 ENDIF 00979 IF(ISUPMX.GT.ZTO) THEN 00980 DO 123 II = ZTO+1,ISUPMX 00981 Z( II, WINDEX ) = ZERO 00982 123 CONTINUE 00983 ENDIF 00984 CALL SSCAL( ZTO-ZFROM+1, NRMINV, 00985 $ Z( ZFROM, WINDEX ), 1 ) 00986 125 CONTINUE 00987 * Update W 00988 W( WINDEX ) = LAMBDA+SIGMA 00989 * Recompute the gaps on the left and right 00990 * But only allow them to become larger and not 00991 * smaller (which can only happen through "bad" 00992 * cancellation and doesn't reflect the theory 00993 * where the initial gaps are underestimated due 00994 * to WERR being too crude.) 00995 IF(.NOT.ESKIP) THEN 00996 IF( K.GT.1) THEN 00997 WGAP( WINDMN ) = MAX( WGAP(WINDMN), 00998 $ W(WINDEX)-WERR(WINDEX) 00999 $ - W(WINDMN)-WERR(WINDMN) ) 01000 ENDIF 01001 IF( WINDEX.LT.WEND ) THEN 01002 WGAP( WINDEX ) = MAX( SAVGAP, 01003 $ W( WINDPL )-WERR( WINDPL ) 01004 $ - W( WINDEX )-WERR( WINDEX) ) 01005 ENDIF 01006 ENDIF 01007 IDONE = IDONE + 1 01008 ENDIF 01009 * here ends the code for the current child 01010 * 01011 139 CONTINUE 01012 * Proceed to any remaining child nodes 01013 NEWFST = J + 1 01014 140 CONTINUE 01015 150 CONTINUE 01016 NDEPTH = NDEPTH + 1 01017 GO TO 40 01018 END IF 01019 IBEGIN = IEND + 1 01020 WBEGIN = WEND + 1 01021 170 CONTINUE 01022 * 01023 01024 RETURN 01025 * 01026 * End of SLARRV 01027 * 01028 END