![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLARRV 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLARRV + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarrv.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarrv.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarrv.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLARRV( 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 * COMPLEX Z( LDZ, * ) 00037 * .. 00038 * 00039 * 00040 *> \par Purpose: 00041 * ============= 00042 *> 00043 *> \verbatim 00044 *> 00045 *> CLARRV 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 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 CLARRV. 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 complexOTHERauxiliary 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 CLARRV( 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 COMPLEX Z( LDZ, * ) 00301 * .. 00302 * 00303 * ===================================================================== 00304 * 00305 * .. Parameters .. 00306 INTEGER MAXITR 00307 PARAMETER ( MAXITR = 10 ) 00308 COMPLEX CZERO 00309 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ) ) 00310 REAL ZERO, ONE, TWO, THREE, FOUR, HALF 00311 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, 00312 $ TWO = 2.0E0, THREE = 3.0E0, 00313 $ FOUR = 4.0E0, HALF = 0.5E0) 00314 * .. 00315 * .. Local Scalars .. 00316 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ 00317 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1, 00318 $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG, 00319 $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER, 00320 $ ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS, 00321 $ NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST, 00322 $ NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST, 00323 $ OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX, 00324 $ WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU, 00325 $ ZUSEDW 00326 INTEGER INDIN1, INDIN2 00327 REAL BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU, 00328 $ LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID, 00329 $ RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF, 00330 $ SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ 00331 * .. 00332 * .. External Functions .. 00333 REAL SLAMCH 00334 EXTERNAL SLAMCH 00335 * .. 00336 * .. External Subroutines .. 00337 EXTERNAL CLAR1V, CLASET, CSSCAL, SCOPY, SLARRB, 00338 $ SLARRF 00339 * .. 00340 * .. Intrinsic Functions .. 00341 INTRINSIC ABS, REAL, MAX, MIN 00342 INTRINSIC CMPLX 00343 * .. 00344 * .. Executable Statements .. 00345 * .. 00346 00347 * The first N entries of WORK are reserved for the eigenvalues 00348 INDLD = N+1 00349 INDLLD= 2*N+1 00350 INDIN1 = 3*N + 1 00351 INDIN2 = 4*N + 1 00352 INDWRK = 5*N + 1 00353 MINWSIZE = 12 * N 00354 00355 DO 5 I= 1,MINWSIZE 00356 WORK( I ) = ZERO 00357 5 CONTINUE 00358 00359 * IWORK(IINDR+1:IINDR+N) hold the twist indices R for the 00360 * factorization used to compute the FP vector 00361 IINDR = 0 00362 * IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current 00363 * layer and the one above. 00364 IINDC1 = N 00365 IINDC2 = 2*N 00366 IINDWK = 3*N + 1 00367 00368 MINIWSIZE = 7 * N 00369 DO 10 I= 1,MINIWSIZE 00370 IWORK( I ) = 0 00371 10 CONTINUE 00372 00373 ZUSEDL = 1 00374 IF(DOL.GT.1) THEN 00375 * Set lower bound for use of Z 00376 ZUSEDL = DOL-1 00377 ENDIF 00378 ZUSEDU = M 00379 IF(DOU.LT.M) THEN 00380 * Set lower bound for use of Z 00381 ZUSEDU = DOU+1 00382 ENDIF 00383 * The width of the part of Z that is used 00384 ZUSEDW = ZUSEDU - ZUSEDL + 1 00385 00386 00387 CALL CLASET( 'Full', N, ZUSEDW, CZERO, CZERO, 00388 $ Z(1,ZUSEDL), LDZ ) 00389 00390 EPS = SLAMCH( 'Precision' ) 00391 RQTOL = TWO * EPS 00392 * 00393 * Set expert flags for standard code. 00394 TRYRQC = .TRUE. 00395 00396 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00397 ELSE 00398 * Only selected eigenpairs are computed. Since the other evalues 00399 * are not refined by RQ iteration, bisection has to compute to full 00400 * accuracy. 00401 RTOL1 = FOUR * EPS 00402 RTOL2 = FOUR * EPS 00403 ENDIF 00404 00405 * The entries WBEGIN:WEND in W, WERR, WGAP correspond to the 00406 * desired eigenvalues. The support of the nonzero eigenvector 00407 * entries is contained in the interval IBEGIN:IEND. 00408 * Remark that if k eigenpairs are desired, then the eigenvectors 00409 * are stored in k contiguous columns of Z. 00410 00411 * DONE is the number of eigenvectors already computed 00412 DONE = 0 00413 IBEGIN = 1 00414 WBEGIN = 1 00415 DO 170 JBLK = 1, IBLOCK( M ) 00416 IEND = ISPLIT( JBLK ) 00417 SIGMA = L( IEND ) 00418 * Find the eigenvectors of the submatrix indexed IBEGIN 00419 * through IEND. 00420 WEND = WBEGIN - 1 00421 15 CONTINUE 00422 IF( WEND.LT.M ) THEN 00423 IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN 00424 WEND = WEND + 1 00425 GO TO 15 00426 END IF 00427 END IF 00428 IF( WEND.LT.WBEGIN ) THEN 00429 IBEGIN = IEND + 1 00430 GO TO 170 00431 ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN 00432 IBEGIN = IEND + 1 00433 WBEGIN = WEND + 1 00434 GO TO 170 00435 END IF 00436 00437 * Find local spectral diameter of the block 00438 GL = GERS( 2*IBEGIN-1 ) 00439 GU = GERS( 2*IBEGIN ) 00440 DO 20 I = IBEGIN+1 , IEND 00441 GL = MIN( GERS( 2*I-1 ), GL ) 00442 GU = MAX( GERS( 2*I ), GU ) 00443 20 CONTINUE 00444 SPDIAM = GU - GL 00445 00446 * OLDIEN is the last index of the previous block 00447 OLDIEN = IBEGIN - 1 00448 * Calculate the size of the current block 00449 IN = IEND - IBEGIN + 1 00450 * The number of eigenvalues in the current block 00451 IM = WEND - WBEGIN + 1 00452 00453 * This is for a 1x1 block 00454 IF( IBEGIN.EQ.IEND ) THEN 00455 DONE = DONE+1 00456 Z( IBEGIN, WBEGIN ) = CMPLX( ONE, ZERO ) 00457 ISUPPZ( 2*WBEGIN-1 ) = IBEGIN 00458 ISUPPZ( 2*WBEGIN ) = IBEGIN 00459 W( WBEGIN ) = W( WBEGIN ) + SIGMA 00460 WORK( WBEGIN ) = W( WBEGIN ) 00461 IBEGIN = IEND + 1 00462 WBEGIN = WBEGIN + 1 00463 GO TO 170 00464 END IF 00465 00466 * The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) 00467 * Note that these can be approximations, in this case, the corresp. 00468 * entries of WERR give the size of the uncertainty interval. 00469 * The eigenvalue approximations will be refined when necessary as 00470 * high relative accuracy is required for the computation of the 00471 * corresponding eigenvectors. 00472 CALL SCOPY( IM, W( WBEGIN ), 1, 00473 $ WORK( WBEGIN ), 1 ) 00474 00475 * We store in W the eigenvalue approximations w.r.t. the original 00476 * matrix T. 00477 DO 30 I=1,IM 00478 W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA 00479 30 CONTINUE 00480 00481 00482 * NDEPTH is the current depth of the representation tree 00483 NDEPTH = 0 00484 * PARITY is either 1 or 0 00485 PARITY = 1 00486 * NCLUS is the number of clusters for the next level of the 00487 * representation tree, we start with NCLUS = 1 for the root 00488 NCLUS = 1 00489 IWORK( IINDC1+1 ) = 1 00490 IWORK( IINDC1+2 ) = IM 00491 00492 * IDONE is the number of eigenvectors already computed in the current 00493 * block 00494 IDONE = 0 00495 * loop while( IDONE.LT.IM ) 00496 * generate the representation tree for the current block and 00497 * compute the eigenvectors 00498 40 CONTINUE 00499 IF( IDONE.LT.IM ) THEN 00500 * This is a crude protection against infinitely deep trees 00501 IF( NDEPTH.GT.M ) THEN 00502 INFO = -2 00503 RETURN 00504 ENDIF 00505 * breadth first processing of the current level of the representation 00506 * tree: OLDNCL = number of clusters on current level 00507 OLDNCL = NCLUS 00508 * reset NCLUS to count the number of child clusters 00509 NCLUS = 0 00510 * 00511 PARITY = 1 - PARITY 00512 IF( PARITY.EQ.0 ) THEN 00513 OLDCLS = IINDC1 00514 NEWCLS = IINDC2 00515 ELSE 00516 OLDCLS = IINDC2 00517 NEWCLS = IINDC1 00518 END IF 00519 * Process the clusters on the current level 00520 DO 150 I = 1, OLDNCL 00521 J = OLDCLS + 2*I 00522 * OLDFST, OLDLST = first, last index of current cluster. 00523 * cluster indices start with 1 and are relative 00524 * to WBEGIN when accessing W, WGAP, WERR, Z 00525 OLDFST = IWORK( J-1 ) 00526 OLDLST = IWORK( J ) 00527 IF( NDEPTH.GT.0 ) THEN 00528 * Retrieve relatively robust representation (RRR) of cluster 00529 * that has been computed at the previous level 00530 * The RRR is stored in Z and overwritten once the eigenvectors 00531 * have been computed or when the cluster is refined 00532 00533 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00534 * Get representation from location of the leftmost evalue 00535 * of the cluster 00536 J = WBEGIN + OLDFST - 1 00537 ELSE 00538 IF(WBEGIN+OLDFST-1.LT.DOL) THEN 00539 * Get representation from the left end of Z array 00540 J = DOL - 1 00541 ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN 00542 * Get representation from the right end of Z array 00543 J = DOU 00544 ELSE 00545 J = WBEGIN + OLDFST - 1 00546 ENDIF 00547 ENDIF 00548 DO 45 K = 1, IN - 1 00549 D( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1, $ J ) ) 00550 L( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1, $ J+1 ) ) 00551 45 CONTINUE 00552 D( IEND ) = REAL( Z( IEND, J ) ) 00553 SIGMA = REAL( Z( IEND, J+1 ) ) 00554 00555 * Set the corresponding entries in Z to zero 00556 CALL CLASET( 'Full', IN, 2, CZERO, CZERO, 00557 $ Z( IBEGIN, J), LDZ ) 00558 END IF 00559 00560 * Compute DL and DLL of current RRR 00561 DO 50 J = IBEGIN, IEND-1 00562 TMP = D( J )*L( J ) 00563 WORK( INDLD-1+J ) = TMP 00564 WORK( INDLLD-1+J ) = TMP*L( J ) 00565 50 CONTINUE 00566 00567 IF( NDEPTH.GT.0 ) THEN 00568 * P and Q are index of the first and last eigenvalue to compute 00569 * within the current block 00570 P = INDEXW( WBEGIN-1+OLDFST ) 00571 Q = INDEXW( WBEGIN-1+OLDLST ) 00572 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET 00573 * through the Q-OFFSET elements of these arrays are to be used. 00574 * OFFSET = P-OLDFST 00575 OFFSET = INDEXW( WBEGIN ) - 1 00576 * perform limited bisection (if necessary) to get approximate 00577 * eigenvalues to the precision needed. 00578 CALL SLARRB( IN, D( IBEGIN ), 00579 $ WORK(INDLLD+IBEGIN-1), 00580 $ P, Q, RTOL1, RTOL2, OFFSET, 00581 $ WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN), 00582 $ WORK( INDWRK ), IWORK( IINDWK ), 00583 $ PIVMIN, SPDIAM, IN, IINFO ) 00584 IF( IINFO.NE.0 ) THEN 00585 INFO = -1 00586 RETURN 00587 ENDIF 00588 * We also recompute the extremal gaps. W holds all eigenvalues 00589 * of the unshifted matrix and must be used for computation 00590 * of WGAP, the entries of WORK might stem from RRRs with 00591 * different shifts. The gaps from WBEGIN-1+OLDFST to 00592 * WBEGIN-1+OLDLST are correctly computed in SLARRB. 00593 * However, we only allow the gaps to become greater since 00594 * this is what should happen when we decrease WERR 00595 IF( OLDFST.GT.1) THEN 00596 WGAP( WBEGIN+OLDFST-2 ) = 00597 $ MAX(WGAP(WBEGIN+OLDFST-2), 00598 $ W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1) 00599 $ - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) ) 00600 ENDIF 00601 IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN 00602 WGAP( WBEGIN+OLDLST-1 ) = 00603 $ MAX(WGAP(WBEGIN+OLDLST-1), 00604 $ W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST) 00605 $ - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) ) 00606 ENDIF 00607 * Each time the eigenvalues in WORK get refined, we store 00608 * the newly found approximation with all shifts applied in W 00609 DO 53 J=OLDFST,OLDLST 00610 W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA 00611 53 CONTINUE 00612 END IF 00613 00614 * Process the current node. 00615 NEWFST = OLDFST 00616 DO 140 J = OLDFST, OLDLST 00617 IF( J.EQ.OLDLST ) THEN 00618 * we are at the right end of the cluster, this is also the 00619 * boundary of the child cluster 00620 NEWLST = J 00621 ELSE IF ( WGAP( WBEGIN + J -1).GE. 00622 $ MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN 00623 * the right relative gap is big enough, the child cluster 00624 * (NEWFST,..,NEWLST) is well separated from the following 00625 NEWLST = J 00626 ELSE 00627 * inside a child cluster, the relative gap is not 00628 * big enough. 00629 GOTO 140 00630 END IF 00631 00632 * Compute size of child cluster found 00633 NEWSIZ = NEWLST - NEWFST + 1 00634 00635 * NEWFTT is the place in Z where the new RRR or the computed 00636 * eigenvector is to be stored 00637 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00638 * Store representation at location of the leftmost evalue 00639 * of the cluster 00640 NEWFTT = WBEGIN + NEWFST - 1 00641 ELSE 00642 IF(WBEGIN+NEWFST-1.LT.DOL) THEN 00643 * Store representation at the left end of Z array 00644 NEWFTT = DOL - 1 00645 ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN 00646 * Store representation at the right end of Z array 00647 NEWFTT = DOU 00648 ELSE 00649 NEWFTT = WBEGIN + NEWFST - 1 00650 ENDIF 00651 ENDIF 00652 00653 IF( NEWSIZ.GT.1) THEN 00654 * 00655 * Current child is not a singleton but a cluster. 00656 * Compute and store new representation of child. 00657 * 00658 * 00659 * Compute left and right cluster gap. 00660 * 00661 * LGAP and RGAP are not computed from WORK because 00662 * the eigenvalue approximations may stem from RRRs 00663 * different shifts. However, W hold all eigenvalues 00664 * of the unshifted matrix. Still, the entries in WGAP 00665 * have to be computed from WORK since the entries 00666 * in W might be of the same order so that gaps are not 00667 * exhibited correctly for very close eigenvalues. 00668 IF( NEWFST.EQ.1 ) THEN 00669 LGAP = MAX( ZERO, 00670 $ W(WBEGIN)-WERR(WBEGIN) - VL ) 00671 ELSE 00672 LGAP = WGAP( WBEGIN+NEWFST-2 ) 00673 ENDIF 00674 RGAP = WGAP( WBEGIN+NEWLST-1 ) 00675 * 00676 * Compute left- and rightmost eigenvalue of child 00677 * to high precision in order to shift as close 00678 * as possible and obtain as large relative gaps 00679 * as possible 00680 * 00681 DO 55 K =1,2 00682 IF(K.EQ.1) THEN 00683 P = INDEXW( WBEGIN-1+NEWFST ) 00684 ELSE 00685 P = INDEXW( WBEGIN-1+NEWLST ) 00686 ENDIF 00687 OFFSET = INDEXW( WBEGIN ) - 1 00688 CALL SLARRB( IN, D(IBEGIN), 00689 $ WORK( INDLLD+IBEGIN-1 ),P,P, 00690 $ RQTOL, RQTOL, OFFSET, 00691 $ WORK(WBEGIN),WGAP(WBEGIN), 00692 $ WERR(WBEGIN),WORK( INDWRK ), 00693 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 00694 $ IN, IINFO ) 00695 55 CONTINUE 00696 * 00697 IF((WBEGIN+NEWLST-1.LT.DOL).OR. 00698 $ (WBEGIN+NEWFST-1.GT.DOU)) THEN 00699 * if the cluster contains no desired eigenvalues 00700 * skip the computation of that branch of the rep. tree 00701 * 00702 * We could skip before the refinement of the extremal 00703 * eigenvalues of the child, but then the representation 00704 * tree could be different from the one when nothing is 00705 * skipped. For this reason we skip at this place. 00706 IDONE = IDONE + NEWLST - NEWFST + 1 00707 GOTO 139 00708 ENDIF 00709 * 00710 * Compute RRR of child cluster. 00711 * Note that the new RRR is stored in Z 00712 * 00713 * SLARRF needs LWORK = 2*N 00714 CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ), 00715 $ WORK(INDLD+IBEGIN-1), 00716 $ NEWFST, NEWLST, WORK(WBEGIN), 00717 $ WGAP(WBEGIN), WERR(WBEGIN), 00718 $ SPDIAM, LGAP, RGAP, PIVMIN, TAU, 00719 $ WORK( INDIN1 ), WORK( INDIN2 ), 00720 $ WORK( INDWRK ), IINFO ) 00721 * In the complex case, SLARRF cannot write 00722 * the new RRR directly into Z and needs an intermediate 00723 * workspace 00724 DO 56 K = 1, IN-1 00725 Z( IBEGIN+K-1, NEWFTT ) = 00726 $ CMPLX( WORK( INDIN1+K-1 ), ZERO ) 00727 Z( IBEGIN+K-1, NEWFTT+1 ) = 00728 $ CMPLX( WORK( INDIN2+K-1 ), ZERO ) 00729 56 CONTINUE 00730 Z( IEND, NEWFTT ) = 00731 $ CMPLX( WORK( INDIN1+IN-1 ), ZERO ) 00732 IF( IINFO.EQ.0 ) THEN 00733 * a new RRR for the cluster was found by SLARRF 00734 * update shift and store it 00735 SSIGMA = SIGMA + TAU 00736 Z( IEND, NEWFTT+1 ) = CMPLX( SSIGMA, ZERO ) 00737 * WORK() are the midpoints and WERR() the semi-width 00738 * Note that the entries in W are unchanged. 00739 DO 116 K = NEWFST, NEWLST 00740 FUDGE = 00741 $ THREE*EPS*ABS(WORK(WBEGIN+K-1)) 00742 WORK( WBEGIN + K - 1 ) = 00743 $ WORK( WBEGIN + K - 1) - TAU 00744 FUDGE = FUDGE + 00745 $ FOUR*EPS*ABS(WORK(WBEGIN+K-1)) 00746 * Fudge errors 00747 WERR( WBEGIN + K - 1 ) = 00748 $ WERR( WBEGIN + K - 1 ) + FUDGE 00749 * Gaps are not fudged. Provided that WERR is small 00750 * when eigenvalues are close, a zero gap indicates 00751 * that a new representation is needed for resolving 00752 * the cluster. A fudge could lead to a wrong decision 00753 * of judging eigenvalues 'separated' which in 00754 * reality are not. This could have a negative impact 00755 * on the orthogonality of the computed eigenvectors. 00756 116 CONTINUE 00757 00758 NCLUS = NCLUS + 1 00759 K = NEWCLS + 2*NCLUS 00760 IWORK( K-1 ) = NEWFST 00761 IWORK( K ) = NEWLST 00762 ELSE 00763 INFO = -2 00764 RETURN 00765 ENDIF 00766 ELSE 00767 * 00768 * Compute eigenvector of singleton 00769 * 00770 ITER = 0 00771 * 00772 TOL = FOUR * LOG(REAL(IN)) * EPS 00773 * 00774 K = NEWFST 00775 WINDEX = WBEGIN + K - 1 00776 WINDMN = MAX(WINDEX - 1,1) 00777 WINDPL = MIN(WINDEX + 1,M) 00778 LAMBDA = WORK( WINDEX ) 00779 DONE = DONE + 1 00780 * Check if eigenvector computation is to be skipped 00781 IF((WINDEX.LT.DOL).OR. 00782 $ (WINDEX.GT.DOU)) THEN 00783 ESKIP = .TRUE. 00784 GOTO 125 00785 ELSE 00786 ESKIP = .FALSE. 00787 ENDIF 00788 LEFT = WORK( WINDEX ) - WERR( WINDEX ) 00789 RIGHT = WORK( WINDEX ) + WERR( WINDEX ) 00790 INDEIG = INDEXW( WINDEX ) 00791 * Note that since we compute the eigenpairs for a child, 00792 * all eigenvalue approximations are w.r.t the same shift. 00793 * In this case, the entries in WORK should be used for 00794 * computing the gaps since they exhibit even very small 00795 * differences in the eigenvalues, as opposed to the 00796 * entries in W which might "look" the same. 00797 00798 IF( K .EQ. 1) THEN 00799 * In the case RANGE='I' and with not much initial 00800 * accuracy in LAMBDA and VL, the formula 00801 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) 00802 * can lead to an overestimation of the left gap and 00803 * thus to inadequately early RQI 'convergence'. 00804 * Prevent this by forcing a small left gap. 00805 LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 00806 ELSE 00807 LGAP = WGAP(WINDMN) 00808 ENDIF 00809 IF( K .EQ. IM) THEN 00810 * In the case RANGE='I' and with not much initial 00811 * accuracy in LAMBDA and VU, the formula 00812 * can lead to an overestimation of the right gap and 00813 * thus to inadequately early RQI 'convergence'. 00814 * Prevent this by forcing a small right gap. 00815 RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 00816 ELSE 00817 RGAP = WGAP(WINDEX) 00818 ENDIF 00819 GAP = MIN( LGAP, RGAP ) 00820 IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN 00821 * The eigenvector support can become wrong 00822 * because significant entries could be cut off due to a 00823 * large GAPTOL parameter in LAR1V. Prevent this. 00824 GAPTOL = ZERO 00825 ELSE 00826 GAPTOL = GAP * EPS 00827 ENDIF 00828 ISUPMN = IN 00829 ISUPMX = 1 00830 * Update WGAP so that it holds the minimum gap 00831 * to the left or the right. This is crucial in the 00832 * case where bisection is used to ensure that the 00833 * eigenvalue is refined up to the required precision. 00834 * The correct value is restored afterwards. 00835 SAVGAP = WGAP(WINDEX) 00836 WGAP(WINDEX) = GAP 00837 * We want to use the Rayleigh Quotient Correction 00838 * as often as possible since it converges quadratically 00839 * when we are close enough to the desired eigenvalue. 00840 * However, the Rayleigh Quotient can have the wrong sign 00841 * and lead us away from the desired eigenvalue. In this 00842 * case, the best we can do is to use bisection. 00843 USEDBS = .FALSE. 00844 USEDRQ = .FALSE. 00845 * Bisection is initially turned off unless it is forced 00846 NEEDBS = .NOT.TRYRQC 00847 120 CONTINUE 00848 * Check if bisection should be used to refine eigenvalue 00849 IF(NEEDBS) THEN 00850 * Take the bisection as new iterate 00851 USEDBS = .TRUE. 00852 ITMP1 = IWORK( IINDR+WINDEX ) 00853 OFFSET = INDEXW( WBEGIN ) - 1 00854 CALL SLARRB( IN, D(IBEGIN), 00855 $ WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG, 00856 $ ZERO, TWO*EPS, OFFSET, 00857 $ WORK(WBEGIN),WGAP(WBEGIN), 00858 $ WERR(WBEGIN),WORK( INDWRK ), 00859 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 00860 $ ITMP1, IINFO ) 00861 IF( IINFO.NE.0 ) THEN 00862 INFO = -3 00863 RETURN 00864 ENDIF 00865 LAMBDA = WORK( WINDEX ) 00866 * Reset twist index from inaccurate LAMBDA to 00867 * force computation of true MINGMA 00868 IWORK( IINDR+WINDEX ) = 0 00869 ENDIF 00870 * Given LAMBDA, compute the eigenvector. 00871 CALL CLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 00872 $ L( IBEGIN ), WORK(INDLD+IBEGIN-1), 00873 $ WORK(INDLLD+IBEGIN-1), 00874 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 00875 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 00876 $ IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ), 00877 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 00878 IF(ITER .EQ. 0) THEN 00879 BSTRES = RESID 00880 BSTW = LAMBDA 00881 ELSEIF(RESID.LT.BSTRES) THEN 00882 BSTRES = RESID 00883 BSTW = LAMBDA 00884 ENDIF 00885 ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 )) 00886 ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX )) 00887 ITER = ITER + 1 00888 00889 * sin alpha <= |resid|/gap 00890 * Note that both the residual and the gap are 00891 * proportional to the matrix, so ||T|| doesn't play 00892 * a role in the quotient 00893 00894 * 00895 * Convergence test for Rayleigh-Quotient iteration 00896 * (omitted when Bisection has been used) 00897 * 00898 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 00899 $ RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS) 00900 $ THEN 00901 * We need to check that the RQCORR update doesn't 00902 * move the eigenvalue away from the desired one and 00903 * towards a neighbor. -> protection with bisection 00904 IF(INDEIG.LE.NEGCNT) THEN 00905 * The wanted eigenvalue lies to the left 00906 SGNDEF = -ONE 00907 ELSE 00908 * The wanted eigenvalue lies to the right 00909 SGNDEF = ONE 00910 ENDIF 00911 * We only use the RQCORR if it improves the 00912 * the iterate reasonably. 00913 IF( ( RQCORR*SGNDEF.GE.ZERO ) 00914 $ .AND.( LAMBDA + RQCORR.LE. RIGHT) 00915 $ .AND.( LAMBDA + RQCORR.GE. LEFT) 00916 $ ) THEN 00917 USEDRQ = .TRUE. 00918 * Store new midpoint of bisection interval in WORK 00919 IF(SGNDEF.EQ.ONE) THEN 00920 * The current LAMBDA is on the left of the true 00921 * eigenvalue 00922 LEFT = LAMBDA 00923 * We prefer to assume that the error estimate 00924 * is correct. We could make the interval not 00925 * as a bracket but to be modified if the RQCORR 00926 * chooses to. In this case, the RIGHT side should 00927 * be modified as follows: 00928 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR) 00929 ELSE 00930 * The current LAMBDA is on the right of the true 00931 * eigenvalue 00932 RIGHT = LAMBDA 00933 * See comment about assuming the error estimate is 00934 * correct above. 00935 * LEFT = MIN(LEFT, LAMBDA + RQCORR) 00936 ENDIF 00937 WORK( WINDEX ) = 00938 $ HALF * (RIGHT + LEFT) 00939 * Take RQCORR since it has the correct sign and 00940 * improves the iterate reasonably 00941 LAMBDA = LAMBDA + RQCORR 00942 * Update width of error interval 00943 WERR( WINDEX ) = 00944 $ HALF * (RIGHT-LEFT) 00945 ELSE 00946 NEEDBS = .TRUE. 00947 ENDIF 00948 IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN 00949 * The eigenvalue is computed to bisection accuracy 00950 * compute eigenvector and stop 00951 USEDBS = .TRUE. 00952 GOTO 120 00953 ELSEIF( ITER.LT.MAXITR ) THEN 00954 GOTO 120 00955 ELSEIF( ITER.EQ.MAXITR ) THEN 00956 NEEDBS = .TRUE. 00957 GOTO 120 00958 ELSE 00959 INFO = 5 00960 RETURN 00961 END IF 00962 ELSE 00963 STP2II = .FALSE. 00964 IF(USEDRQ .AND. USEDBS .AND. 00965 $ BSTRES.LE.RESID) THEN 00966 LAMBDA = BSTW 00967 STP2II = .TRUE. 00968 ENDIF 00969 IF (STP2II) THEN 00970 * improve error angle by second step 00971 CALL CLAR1V( IN, 1, IN, LAMBDA, 00972 $ D( IBEGIN ), L( IBEGIN ), 00973 $ WORK(INDLD+IBEGIN-1), 00974 $ WORK(INDLLD+IBEGIN-1), 00975 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 00976 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 00977 $ IWORK( IINDR+WINDEX ), 00978 $ ISUPPZ( 2*WINDEX-1 ), 00979 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 00980 ENDIF 00981 WORK( WINDEX ) = LAMBDA 00982 END IF 00983 * 00984 * Compute FP-vector support w.r.t. whole matrix 00985 * 00986 ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN 00987 ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN 00988 ZFROM = ISUPPZ( 2*WINDEX-1 ) 00989 ZTO = ISUPPZ( 2*WINDEX ) 00990 ISUPMN = ISUPMN + OLDIEN 00991 ISUPMX = ISUPMX + OLDIEN 00992 * Ensure vector is ok if support in the RQI has changed 00993 IF(ISUPMN.LT.ZFROM) THEN 00994 DO 122 II = ISUPMN,ZFROM-1 00995 Z( II, WINDEX ) = ZERO 00996 122 CONTINUE 00997 ENDIF 00998 IF(ISUPMX.GT.ZTO) THEN 00999 DO 123 II = ZTO+1,ISUPMX 01000 Z( II, WINDEX ) = ZERO 01001 123 CONTINUE 01002 ENDIF 01003 CALL CSSCAL( ZTO-ZFROM+1, NRMINV, 01004 $ Z( ZFROM, WINDEX ), 1 ) 01005 125 CONTINUE 01006 * Update W 01007 W( WINDEX ) = LAMBDA+SIGMA 01008 * Recompute the gaps on the left and right 01009 * But only allow them to become larger and not 01010 * smaller (which can only happen through "bad" 01011 * cancellation and doesn't reflect the theory 01012 * where the initial gaps are underestimated due 01013 * to WERR being too crude.) 01014 IF(.NOT.ESKIP) THEN 01015 IF( K.GT.1) THEN 01016 WGAP( WINDMN ) = MAX( WGAP(WINDMN), 01017 $ W(WINDEX)-WERR(WINDEX) 01018 $ - W(WINDMN)-WERR(WINDMN) ) 01019 ENDIF 01020 IF( WINDEX.LT.WEND ) THEN 01021 WGAP( WINDEX ) = MAX( SAVGAP, 01022 $ W( WINDPL )-WERR( WINDPL ) 01023 $ - W( WINDEX )-WERR( WINDEX) ) 01024 ENDIF 01025 ENDIF 01026 IDONE = IDONE + 1 01027 ENDIF 01028 * here ends the code for the current child 01029 * 01030 139 CONTINUE 01031 * Proceed to any remaining child nodes 01032 NEWFST = J + 1 01033 140 CONTINUE 01034 150 CONTINUE 01035 NDEPTH = NDEPTH + 1 01036 GO TO 40 01037 END IF 01038 IBEGIN = IEND + 1 01039 WBEGIN = WEND + 1 01040 170 CONTINUE 01041 * 01042 01043 RETURN 01044 * 01045 * End of CLARRV 01046 * 01047 END 01048