![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLARRV 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLARRV + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarrv.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarrv.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarrv.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLARRV( 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 * DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 00033 * $ ISUPPZ( * ), IWORK( * ) 00034 * DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 00035 * $ WGAP( * ), WORK( * ) 00036 * COMPLEX*16 Z( LDZ, * ) 00037 * .. 00038 * 00039 * 00040 *> \par Purpose: 00041 * ============= 00042 *> 00043 *> \verbatim 00044 *> 00045 *> ZLARRV 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 DLARRE. 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 DOUBLE PRECISION 00062 *> \endverbatim 00063 *> 00064 *> \param[in] VU 00065 *> \verbatim 00066 *> VU is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLARRE. 00086 *> On exit, L is overwritten. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] PIVMIN 00090 *> \verbatim 00091 *> PIVMIN is DOUBLE PRECISION 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 DOUBLE PRECISION 00131 *> \endverbatim 00132 *> 00133 *> \param[in] RTOL1 00134 *> \verbatim 00135 *> RTOL1 is DOUBLE PRECISION 00136 *> \endverbatim 00137 *> 00138 *> \param[in] RTOL2 00139 *> \verbatim 00140 *> RTOL2 is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLARRE 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 COMPLEX*16 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 DOUBLE PRECISION 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 ZLARRV. 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 DLARRB when refining a child's eigenvalues. 00245 *> =-2: Problem in DLARRF 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 DLARRB 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 complex16OTHERauxiliary 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 ZLARRV( 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 DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 00294 * .. 00295 * .. Array Arguments .. 00296 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 00297 $ ISUPPZ( * ), IWORK( * ) 00298 DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 00299 $ WGAP( * ), WORK( * ) 00300 COMPLEX*16 Z( LDZ, * ) 00301 * .. 00302 * 00303 * ===================================================================== 00304 * 00305 * .. Parameters .. 00306 INTEGER MAXITR 00307 PARAMETER ( MAXITR = 10 ) 00308 COMPLEX*16 CZERO 00309 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) ) 00310 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF 00311 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 00312 $ TWO = 2.0D0, THREE = 3.0D0, 00313 $ FOUR = 4.0D0, HALF = 0.5D0) 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00334 EXTERNAL DLAMCH 00335 * .. 00336 * .. External Subroutines .. 00337 EXTERNAL DCOPY, DLARRB, DLARRF, ZDSCAL, ZLAR1V, 00338 $ ZLASET 00339 * .. 00340 * .. Intrinsic Functions .. 00341 INTRINSIC ABS, DBLE, MAX, MIN 00342 INTRINSIC DCMPLX 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 ZLASET( 'Full', N, ZUSEDW, CZERO, CZERO, 00388 $ Z(1,ZUSEDL), LDZ ) 00389 00390 EPS = DLAMCH( '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 ) = DCMPLX( 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 DCOPY( 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 ) = DBLE( Z( IBEGIN+K-1, 00550 $ J ) ) 00551 L( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1, 00552 $ J+1 ) ) 00553 45 CONTINUE 00554 D( IEND ) = DBLE( Z( IEND, J ) ) 00555 SIGMA = DBLE( Z( IEND, J+1 ) ) 00556 00557 * Set the corresponding entries in Z to zero 00558 CALL ZLASET( 'Full', IN, 2, CZERO, CZERO, 00559 $ Z( IBEGIN, J), LDZ ) 00560 END IF 00561 00562 * Compute DL and DLL of current RRR 00563 DO 50 J = IBEGIN, IEND-1 00564 TMP = D( J )*L( J ) 00565 WORK( INDLD-1+J ) = TMP 00566 WORK( INDLLD-1+J ) = TMP*L( J ) 00567 50 CONTINUE 00568 00569 IF( NDEPTH.GT.0 ) THEN 00570 * P and Q are index of the first and last eigenvalue to compute 00571 * within the current block 00572 P = INDEXW( WBEGIN-1+OLDFST ) 00573 Q = INDEXW( WBEGIN-1+OLDLST ) 00574 * Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET 00575 * through the Q-OFFSET elements of these arrays are to be used. 00576 * OFFSET = P-OLDFST 00577 OFFSET = INDEXW( WBEGIN ) - 1 00578 * perform limited bisection (if necessary) to get approximate 00579 * eigenvalues to the precision needed. 00580 CALL DLARRB( IN, D( IBEGIN ), 00581 $ WORK(INDLLD+IBEGIN-1), 00582 $ P, Q, RTOL1, RTOL2, OFFSET, 00583 $ WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN), 00584 $ WORK( INDWRK ), IWORK( IINDWK ), 00585 $ PIVMIN, SPDIAM, IN, IINFO ) 00586 IF( IINFO.NE.0 ) THEN 00587 INFO = -1 00588 RETURN 00589 ENDIF 00590 * We also recompute the extremal gaps. W holds all eigenvalues 00591 * of the unshifted matrix and must be used for computation 00592 * of WGAP, the entries of WORK might stem from RRRs with 00593 * different shifts. The gaps from WBEGIN-1+OLDFST to 00594 * WBEGIN-1+OLDLST are correctly computed in DLARRB. 00595 * However, we only allow the gaps to become greater since 00596 * this is what should happen when we decrease WERR 00597 IF( OLDFST.GT.1) THEN 00598 WGAP( WBEGIN+OLDFST-2 ) = 00599 $ MAX(WGAP(WBEGIN+OLDFST-2), 00600 $ W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1) 00601 $ - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) ) 00602 ENDIF 00603 IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN 00604 WGAP( WBEGIN+OLDLST-1 ) = 00605 $ MAX(WGAP(WBEGIN+OLDLST-1), 00606 $ W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST) 00607 $ - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) ) 00608 ENDIF 00609 * Each time the eigenvalues in WORK get refined, we store 00610 * the newly found approximation with all shifts applied in W 00611 DO 53 J=OLDFST,OLDLST 00612 W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA 00613 53 CONTINUE 00614 END IF 00615 00616 * Process the current node. 00617 NEWFST = OLDFST 00618 DO 140 J = OLDFST, OLDLST 00619 IF( J.EQ.OLDLST ) THEN 00620 * we are at the right end of the cluster, this is also the 00621 * boundary of the child cluster 00622 NEWLST = J 00623 ELSE IF ( WGAP( WBEGIN + J -1).GE. 00624 $ MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN 00625 * the right relative gap is big enough, the child cluster 00626 * (NEWFST,..,NEWLST) is well separated from the following 00627 NEWLST = J 00628 ELSE 00629 * inside a child cluster, the relative gap is not 00630 * big enough. 00631 GOTO 140 00632 END IF 00633 00634 * Compute size of child cluster found 00635 NEWSIZ = NEWLST - NEWFST + 1 00636 00637 * NEWFTT is the place in Z where the new RRR or the computed 00638 * eigenvector is to be stored 00639 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 00640 * Store representation at location of the leftmost evalue 00641 * of the cluster 00642 NEWFTT = WBEGIN + NEWFST - 1 00643 ELSE 00644 IF(WBEGIN+NEWFST-1.LT.DOL) THEN 00645 * Store representation at the left end of Z array 00646 NEWFTT = DOL - 1 00647 ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN 00648 * Store representation at the right end of Z array 00649 NEWFTT = DOU 00650 ELSE 00651 NEWFTT = WBEGIN + NEWFST - 1 00652 ENDIF 00653 ENDIF 00654 00655 IF( NEWSIZ.GT.1) THEN 00656 * 00657 * Current child is not a singleton but a cluster. 00658 * Compute and store new representation of child. 00659 * 00660 * 00661 * Compute left and right cluster gap. 00662 * 00663 * LGAP and RGAP are not computed from WORK because 00664 * the eigenvalue approximations may stem from RRRs 00665 * different shifts. However, W hold all eigenvalues 00666 * of the unshifted matrix. Still, the entries in WGAP 00667 * have to be computed from WORK since the entries 00668 * in W might be of the same order so that gaps are not 00669 * exhibited correctly for very close eigenvalues. 00670 IF( NEWFST.EQ.1 ) THEN 00671 LGAP = MAX( ZERO, 00672 $ W(WBEGIN)-WERR(WBEGIN) - VL ) 00673 ELSE 00674 LGAP = WGAP( WBEGIN+NEWFST-2 ) 00675 ENDIF 00676 RGAP = WGAP( WBEGIN+NEWLST-1 ) 00677 * 00678 * Compute left- and rightmost eigenvalue of child 00679 * to high precision in order to shift as close 00680 * as possible and obtain as large relative gaps 00681 * as possible 00682 * 00683 DO 55 K =1,2 00684 IF(K.EQ.1) THEN 00685 P = INDEXW( WBEGIN-1+NEWFST ) 00686 ELSE 00687 P = INDEXW( WBEGIN-1+NEWLST ) 00688 ENDIF 00689 OFFSET = INDEXW( WBEGIN ) - 1 00690 CALL DLARRB( IN, D(IBEGIN), 00691 $ WORK( INDLLD+IBEGIN-1 ),P,P, 00692 $ RQTOL, RQTOL, OFFSET, 00693 $ WORK(WBEGIN),WGAP(WBEGIN), 00694 $ WERR(WBEGIN),WORK( INDWRK ), 00695 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 00696 $ IN, IINFO ) 00697 55 CONTINUE 00698 * 00699 IF((WBEGIN+NEWLST-1.LT.DOL).OR. 00700 $ (WBEGIN+NEWFST-1.GT.DOU)) THEN 00701 * if the cluster contains no desired eigenvalues 00702 * skip the computation of that branch of the rep. tree 00703 * 00704 * We could skip before the refinement of the extremal 00705 * eigenvalues of the child, but then the representation 00706 * tree could be different from the one when nothing is 00707 * skipped. For this reason we skip at this place. 00708 IDONE = IDONE + NEWLST - NEWFST + 1 00709 GOTO 139 00710 ENDIF 00711 * 00712 * Compute RRR of child cluster. 00713 * Note that the new RRR is stored in Z 00714 * 00715 * DLARRF needs LWORK = 2*N 00716 CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ), 00717 $ WORK(INDLD+IBEGIN-1), 00718 $ NEWFST, NEWLST, WORK(WBEGIN), 00719 $ WGAP(WBEGIN), WERR(WBEGIN), 00720 $ SPDIAM, LGAP, RGAP, PIVMIN, TAU, 00721 $ WORK( INDIN1 ), WORK( INDIN2 ), 00722 $ WORK( INDWRK ), IINFO ) 00723 * In the complex case, DLARRF cannot write 00724 * the new RRR directly into Z and needs an intermediate 00725 * workspace 00726 DO 56 K = 1, IN-1 00727 Z( IBEGIN+K-1, NEWFTT ) = 00728 $ DCMPLX( WORK( INDIN1+K-1 ), ZERO ) 00729 Z( IBEGIN+K-1, NEWFTT+1 ) = 00730 $ DCMPLX( WORK( INDIN2+K-1 ), ZERO ) 00731 56 CONTINUE 00732 Z( IEND, NEWFTT ) = 00733 $ DCMPLX( WORK( INDIN1+IN-1 ), ZERO ) 00734 IF( IINFO.EQ.0 ) THEN 00735 * a new RRR for the cluster was found by DLARRF 00736 * update shift and store it 00737 SSIGMA = SIGMA + TAU 00738 Z( IEND, NEWFTT+1 ) = DCMPLX( SSIGMA, ZERO ) 00739 * WORK() are the midpoints and WERR() the semi-width 00740 * Note that the entries in W are unchanged. 00741 DO 116 K = NEWFST, NEWLST 00742 FUDGE = 00743 $ THREE*EPS*ABS(WORK(WBEGIN+K-1)) 00744 WORK( WBEGIN + K - 1 ) = 00745 $ WORK( WBEGIN + K - 1) - TAU 00746 FUDGE = FUDGE + 00747 $ FOUR*EPS*ABS(WORK(WBEGIN+K-1)) 00748 * Fudge errors 00749 WERR( WBEGIN + K - 1 ) = 00750 $ WERR( WBEGIN + K - 1 ) + FUDGE 00751 * Gaps are not fudged. Provided that WERR is small 00752 * when eigenvalues are close, a zero gap indicates 00753 * that a new representation is needed for resolving 00754 * the cluster. A fudge could lead to a wrong decision 00755 * of judging eigenvalues 'separated' which in 00756 * reality are not. This could have a negative impact 00757 * on the orthogonality of the computed eigenvectors. 00758 116 CONTINUE 00759 00760 NCLUS = NCLUS + 1 00761 K = NEWCLS + 2*NCLUS 00762 IWORK( K-1 ) = NEWFST 00763 IWORK( K ) = NEWLST 00764 ELSE 00765 INFO = -2 00766 RETURN 00767 ENDIF 00768 ELSE 00769 * 00770 * Compute eigenvector of singleton 00771 * 00772 ITER = 0 00773 * 00774 TOL = FOUR * LOG(DBLE(IN)) * EPS 00775 * 00776 K = NEWFST 00777 WINDEX = WBEGIN + K - 1 00778 WINDMN = MAX(WINDEX - 1,1) 00779 WINDPL = MIN(WINDEX + 1,M) 00780 LAMBDA = WORK( WINDEX ) 00781 DONE = DONE + 1 00782 * Check if eigenvector computation is to be skipped 00783 IF((WINDEX.LT.DOL).OR. 00784 $ (WINDEX.GT.DOU)) THEN 00785 ESKIP = .TRUE. 00786 GOTO 125 00787 ELSE 00788 ESKIP = .FALSE. 00789 ENDIF 00790 LEFT = WORK( WINDEX ) - WERR( WINDEX ) 00791 RIGHT = WORK( WINDEX ) + WERR( WINDEX ) 00792 INDEIG = INDEXW( WINDEX ) 00793 * Note that since we compute the eigenpairs for a child, 00794 * all eigenvalue approximations are w.r.t the same shift. 00795 * In this case, the entries in WORK should be used for 00796 * computing the gaps since they exhibit even very small 00797 * differences in the eigenvalues, as opposed to the 00798 * entries in W which might "look" the same. 00799 00800 IF( K .EQ. 1) THEN 00801 * In the case RANGE='I' and with not much initial 00802 * accuracy in LAMBDA and VL, the formula 00803 * LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) 00804 * can lead to an overestimation of the left gap and 00805 * thus to inadequately early RQI 'convergence'. 00806 * Prevent this by forcing a small left gap. 00807 LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 00808 ELSE 00809 LGAP = WGAP(WINDMN) 00810 ENDIF 00811 IF( K .EQ. IM) THEN 00812 * In the case RANGE='I' and with not much initial 00813 * accuracy in LAMBDA and VU, the formula 00814 * can lead to an overestimation of the right gap and 00815 * thus to inadequately early RQI 'convergence'. 00816 * Prevent this by forcing a small right gap. 00817 RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 00818 ELSE 00819 RGAP = WGAP(WINDEX) 00820 ENDIF 00821 GAP = MIN( LGAP, RGAP ) 00822 IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN 00823 * The eigenvector support can become wrong 00824 * because significant entries could be cut off due to a 00825 * large GAPTOL parameter in LAR1V. Prevent this. 00826 GAPTOL = ZERO 00827 ELSE 00828 GAPTOL = GAP * EPS 00829 ENDIF 00830 ISUPMN = IN 00831 ISUPMX = 1 00832 * Update WGAP so that it holds the minimum gap 00833 * to the left or the right. This is crucial in the 00834 * case where bisection is used to ensure that the 00835 * eigenvalue is refined up to the required precision. 00836 * The correct value is restored afterwards. 00837 SAVGAP = WGAP(WINDEX) 00838 WGAP(WINDEX) = GAP 00839 * We want to use the Rayleigh Quotient Correction 00840 * as often as possible since it converges quadratically 00841 * when we are close enough to the desired eigenvalue. 00842 * However, the Rayleigh Quotient can have the wrong sign 00843 * and lead us away from the desired eigenvalue. In this 00844 * case, the best we can do is to use bisection. 00845 USEDBS = .FALSE. 00846 USEDRQ = .FALSE. 00847 * Bisection is initially turned off unless it is forced 00848 NEEDBS = .NOT.TRYRQC 00849 120 CONTINUE 00850 * Check if bisection should be used to refine eigenvalue 00851 IF(NEEDBS) THEN 00852 * Take the bisection as new iterate 00853 USEDBS = .TRUE. 00854 ITMP1 = IWORK( IINDR+WINDEX ) 00855 OFFSET = INDEXW( WBEGIN ) - 1 00856 CALL DLARRB( IN, D(IBEGIN), 00857 $ WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG, 00858 $ ZERO, TWO*EPS, OFFSET, 00859 $ WORK(WBEGIN),WGAP(WBEGIN), 00860 $ WERR(WBEGIN),WORK( INDWRK ), 00861 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 00862 $ ITMP1, IINFO ) 00863 IF( IINFO.NE.0 ) THEN 00864 INFO = -3 00865 RETURN 00866 ENDIF 00867 LAMBDA = WORK( WINDEX ) 00868 * Reset twist index from inaccurate LAMBDA to 00869 * force computation of true MINGMA 00870 IWORK( IINDR+WINDEX ) = 0 00871 ENDIF 00872 * Given LAMBDA, compute the eigenvector. 00873 CALL ZLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 00874 $ L( IBEGIN ), WORK(INDLD+IBEGIN-1), 00875 $ WORK(INDLLD+IBEGIN-1), 00876 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 00877 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 00878 $ IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ), 00879 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 00880 IF(ITER .EQ. 0) THEN 00881 BSTRES = RESID 00882 BSTW = LAMBDA 00883 ELSEIF(RESID.LT.BSTRES) THEN 00884 BSTRES = RESID 00885 BSTW = LAMBDA 00886 ENDIF 00887 ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 )) 00888 ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX )) 00889 ITER = ITER + 1 00890 00891 * sin alpha <= |resid|/gap 00892 * Note that both the residual and the gap are 00893 * proportional to the matrix, so ||T|| doesn't play 00894 * a role in the quotient 00895 00896 * 00897 * Convergence test for Rayleigh-Quotient iteration 00898 * (omitted when Bisection has been used) 00899 * 00900 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 00901 $ RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS) 00902 $ THEN 00903 * We need to check that the RQCORR update doesn't 00904 * move the eigenvalue away from the desired one and 00905 * towards a neighbor. -> protection with bisection 00906 IF(INDEIG.LE.NEGCNT) THEN 00907 * The wanted eigenvalue lies to the left 00908 SGNDEF = -ONE 00909 ELSE 00910 * The wanted eigenvalue lies to the right 00911 SGNDEF = ONE 00912 ENDIF 00913 * We only use the RQCORR if it improves the 00914 * the iterate reasonably. 00915 IF( ( RQCORR*SGNDEF.GE.ZERO ) 00916 $ .AND.( LAMBDA + RQCORR.LE. RIGHT) 00917 $ .AND.( LAMBDA + RQCORR.GE. LEFT) 00918 $ ) THEN 00919 USEDRQ = .TRUE. 00920 * Store new midpoint of bisection interval in WORK 00921 IF(SGNDEF.EQ.ONE) THEN 00922 * The current LAMBDA is on the left of the true 00923 * eigenvalue 00924 LEFT = LAMBDA 00925 * We prefer to assume that the error estimate 00926 * is correct. We could make the interval not 00927 * as a bracket but to be modified if the RQCORR 00928 * chooses to. In this case, the RIGHT side should 00929 * be modified as follows: 00930 * RIGHT = MAX(RIGHT, LAMBDA + RQCORR) 00931 ELSE 00932 * The current LAMBDA is on the right of the true 00933 * eigenvalue 00934 RIGHT = LAMBDA 00935 * See comment about assuming the error estimate is 00936 * correct above. 00937 * LEFT = MIN(LEFT, LAMBDA + RQCORR) 00938 ENDIF 00939 WORK( WINDEX ) = 00940 $ HALF * (RIGHT + LEFT) 00941 * Take RQCORR since it has the correct sign and 00942 * improves the iterate reasonably 00943 LAMBDA = LAMBDA + RQCORR 00944 * Update width of error interval 00945 WERR( WINDEX ) = 00946 $ HALF * (RIGHT-LEFT) 00947 ELSE 00948 NEEDBS = .TRUE. 00949 ENDIF 00950 IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN 00951 * The eigenvalue is computed to bisection accuracy 00952 * compute eigenvector and stop 00953 USEDBS = .TRUE. 00954 GOTO 120 00955 ELSEIF( ITER.LT.MAXITR ) THEN 00956 GOTO 120 00957 ELSEIF( ITER.EQ.MAXITR ) THEN 00958 NEEDBS = .TRUE. 00959 GOTO 120 00960 ELSE 00961 INFO = 5 00962 RETURN 00963 END IF 00964 ELSE 00965 STP2II = .FALSE. 00966 IF(USEDRQ .AND. USEDBS .AND. 00967 $ BSTRES.LE.RESID) THEN 00968 LAMBDA = BSTW 00969 STP2II = .TRUE. 00970 ENDIF 00971 IF (STP2II) THEN 00972 * improve error angle by second step 00973 CALL ZLAR1V( IN, 1, IN, LAMBDA, 00974 $ D( IBEGIN ), L( IBEGIN ), 00975 $ WORK(INDLD+IBEGIN-1), 00976 $ WORK(INDLLD+IBEGIN-1), 00977 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 00978 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 00979 $ IWORK( IINDR+WINDEX ), 00980 $ ISUPPZ( 2*WINDEX-1 ), 00981 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 00982 ENDIF 00983 WORK( WINDEX ) = LAMBDA 00984 END IF 00985 * 00986 * Compute FP-vector support w.r.t. whole matrix 00987 * 00988 ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN 00989 ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN 00990 ZFROM = ISUPPZ( 2*WINDEX-1 ) 00991 ZTO = ISUPPZ( 2*WINDEX ) 00992 ISUPMN = ISUPMN + OLDIEN 00993 ISUPMX = ISUPMX + OLDIEN 00994 * Ensure vector is ok if support in the RQI has changed 00995 IF(ISUPMN.LT.ZFROM) THEN 00996 DO 122 II = ISUPMN,ZFROM-1 00997 Z( II, WINDEX ) = ZERO 00998 122 CONTINUE 00999 ENDIF 01000 IF(ISUPMX.GT.ZTO) THEN 01001 DO 123 II = ZTO+1,ISUPMX 01002 Z( II, WINDEX ) = ZERO 01003 123 CONTINUE 01004 ENDIF 01005 CALL ZDSCAL( ZTO-ZFROM+1, NRMINV, 01006 $ Z( ZFROM, WINDEX ), 1 ) 01007 125 CONTINUE 01008 * Update W 01009 W( WINDEX ) = LAMBDA+SIGMA 01010 * Recompute the gaps on the left and right 01011 * But only allow them to become larger and not 01012 * smaller (which can only happen through "bad" 01013 * cancellation and doesn't reflect the theory 01014 * where the initial gaps are underestimated due 01015 * to WERR being too crude.) 01016 IF(.NOT.ESKIP) THEN 01017 IF( K.GT.1) THEN 01018 WGAP( WINDMN ) = MAX( WGAP(WINDMN), 01019 $ W(WINDEX)-WERR(WINDEX) 01020 $ - W(WINDMN)-WERR(WINDMN) ) 01021 ENDIF 01022 IF( WINDEX.LT.WEND ) THEN 01023 WGAP( WINDEX ) = MAX( SAVGAP, 01024 $ W( WINDPL )-WERR( WINDPL ) 01025 $ - W( WINDEX )-WERR( WINDEX) ) 01026 ENDIF 01027 ENDIF 01028 IDONE = IDONE + 1 01029 ENDIF 01030 * here ends the code for the current child 01031 * 01032 139 CONTINUE 01033 * Proceed to any remaining child nodes 01034 NEWFST = J + 1 01035 140 CONTINUE 01036 150 CONTINUE 01037 NDEPTH = NDEPTH + 1 01038 GO TO 40 01039 END IF 01040 IBEGIN = IEND + 1 01041 WBEGIN = WEND + 1 01042 170 CONTINUE 01043 * 01044 01045 RETURN 01046 * 01047 * End of ZLARRV 01048 * 01049 END