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