00001 /* slarrv.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static real c_b5 = 0.f; 00019 static integer c__1 = 1; 00020 static integer c__2 = 2; 00021 00022 /* Subroutine */ int slarrv_(integer *n, real *vl, real *vu, real *d__, real * 00023 l, real *pivmin, integer *isplit, integer *m, integer *dol, integer * 00024 dou, real *minrgp, real *rtol1, real *rtol2, real *w, real *werr, 00025 real *wgap, integer *iblock, integer *indexw, real *gers, real *z__, 00026 integer *ldz, integer *isuppz, real *work, integer *iwork, integer * 00027 info) 00028 { 00029 /* System generated locals */ 00030 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5; 00031 real r__1, r__2; 00032 logical L__1; 00033 00034 /* Builtin functions */ 00035 double log(doublereal); 00036 00037 /* Local variables */ 00038 integer minwsize, i__, j, k, p, q, miniwsize, ii; 00039 real gl; 00040 integer im, in; 00041 real gu, gap, eps, tau, tol, tmp; 00042 integer zto; 00043 real ztz; 00044 integer iend, jblk; 00045 real lgap; 00046 integer done; 00047 real rgap, left; 00048 integer wend, iter; 00049 real bstw; 00050 integer itmp1, indld; 00051 real fudge; 00052 integer idone; 00053 real sigma; 00054 integer iinfo, iindr; 00055 real resid; 00056 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *); 00057 logical eskip; 00058 real right; 00059 integer nclus, zfrom; 00060 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 00061 integer *); 00062 real rqtol; 00063 integer iindc1, iindc2; 00064 extern /* Subroutine */ int slar1v_(integer *, integer *, integer *, real 00065 *, real *, real *, real *, real *, real *, real *, real *, 00066 logical *, integer *, real *, real *, integer *, integer *, real * 00067 , real *, real *, real *); 00068 logical stp2ii; 00069 real lambda; 00070 integer ibegin, indeig; 00071 logical needbs; 00072 integer indlld; 00073 real sgndef, mingma; 00074 extern doublereal slamch_(char *); 00075 integer oldien, oldncl, wbegin; 00076 real spdiam; 00077 integer negcnt, oldcls; 00078 real savgap; 00079 integer ndepth; 00080 real ssigma; 00081 logical usedbs; 00082 integer iindwk, offset; 00083 real gaptol; 00084 extern /* Subroutine */ int slarrb_(integer *, real *, real *, integer *, 00085 integer *, real *, real *, integer *, real *, real *, real *, 00086 real *, integer *, real *, real *, integer *, integer *), slarrf_( 00087 integer *, real *, real *, real *, integer *, integer *, real *, 00088 real *, real *, real *, real *, real *, real *, real *, real *, 00089 real *, real *, integer *); 00090 integer newcls, oldfst, indwrk, windex, oldlst; 00091 logical usedrq; 00092 integer newfst, newftt, parity, windmn, isupmn, newlst, windpl, zusedl, 00093 newsiz, zusedu, zusedw; 00094 real bstres, nrminv; 00095 logical tryrqc; 00096 integer isupmx; 00097 real rqcorr; 00098 extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *, 00099 real *, real *, integer *); 00100 00101 00102 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00103 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00104 /* November 2006 */ 00105 00106 /* .. Scalar Arguments .. */ 00107 /* .. */ 00108 /* .. Array Arguments .. */ 00109 /* .. */ 00110 00111 /* Purpose */ 00112 /* ======= */ 00113 00114 /* SLARRV computes the eigenvectors of the tridiagonal matrix */ 00115 /* T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T. */ 00116 /* The input eigenvalues should have been computed by SLARRE. */ 00117 00118 /* Arguments */ 00119 /* ========= */ 00120 00121 /* N (input) INTEGER */ 00122 /* The order of the matrix. N >= 0. */ 00123 00124 /* VL (input) REAL */ 00125 /* VU (input) REAL */ 00126 /* Lower and upper bounds of the interval that contains the desired */ 00127 /* eigenvalues. VL < VU. Needed to compute gaps on the left or right */ 00128 /* end of the extremal eigenvalues in the desired RANGE. */ 00129 00130 /* D (input/output) REAL array, dimension (N) */ 00131 /* On entry, the N diagonal elements of the diagonal matrix D. */ 00132 /* On exit, D may be overwritten. */ 00133 00134 /* L (input/output) REAL array, dimension (N) */ 00135 /* On entry, the (N-1) subdiagonal elements of the unit */ 00136 /* bidiagonal matrix L are in elements 1 to N-1 of L */ 00137 /* (if the matrix is not splitted.) At the end of each block */ 00138 /* is stored the corresponding shift as given by SLARRE. */ 00139 /* On exit, L is overwritten. */ 00140 00141 /* PIVMIN (in) DOUBLE PRECISION */ 00142 /* The minimum pivot allowed in the Sturm sequence. */ 00143 00144 /* ISPLIT (input) INTEGER array, dimension (N) */ 00145 /* The splitting points, at which T breaks up into blocks. */ 00146 /* The first block consists of rows/columns 1 to */ 00147 /* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */ 00148 /* through ISPLIT( 2 ), etc. */ 00149 00150 /* M (input) INTEGER */ 00151 /* The total number of input eigenvalues. 0 <= M <= N. */ 00152 00153 /* DOL (input) INTEGER */ 00154 /* DOU (input) INTEGER */ 00155 /* If the user wants to compute only selected eigenvectors from all */ 00156 /* the eigenvalues supplied, he can specify an index range DOL:DOU. */ 00157 /* Or else the setting DOL=1, DOU=M should be applied. */ 00158 /* Note that DOL and DOU refer to the order in which the eigenvalues */ 00159 /* are stored in W. */ 00160 /* If the user wants to compute only selected eigenpairs, then */ 00161 /* the columns DOL-1 to DOU+1 of the eigenvector space Z contain the */ 00162 /* computed eigenvectors. All other columns of Z are set to zero. */ 00163 00164 /* MINRGP (input) REAL */ 00165 00166 /* RTOL1 (input) REAL */ 00167 /* RTOL2 (input) REAL */ 00168 /* Parameters for bisection. */ 00169 /* An interval [LEFT,RIGHT] has converged if */ 00170 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */ 00171 00172 /* W (input/output) REAL array, dimension (N) */ 00173 /* The first M elements of W contain the APPROXIMATE eigenvalues for */ 00174 /* which eigenvectors are to be computed. The eigenvalues */ 00175 /* should be grouped by split-off block and ordered from */ 00176 /* smallest to largest within the block ( The output array */ 00177 /* W from SLARRE is expected here ). Furthermore, they are with */ 00178 /* respect to the shift of the corresponding root representation */ 00179 /* for their block. On exit, W holds the eigenvalues of the */ 00180 /* UNshifted matrix. */ 00181 00182 /* WERR (input/output) REAL array, dimension (N) */ 00183 /* The first M elements contain the semiwidth of the uncertainty */ 00184 /* interval of the corresponding eigenvalue in W */ 00185 00186 /* WGAP (input/output) REAL array, dimension (N) */ 00187 /* The separation from the right neighbor eigenvalue in W. */ 00188 00189 /* IBLOCK (input) INTEGER array, dimension (N) */ 00190 /* The indices of the blocks (submatrices) associated with the */ 00191 /* corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue */ 00192 /* W(i) belongs to the first block from the top, =2 if W(i) */ 00193 /* belongs to the second block, etc. */ 00194 00195 /* INDEXW (input) INTEGER array, dimension (N) */ 00196 /* The indices of the eigenvalues within each block (submatrix); */ 00197 /* for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the */ 00198 /* i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. */ 00199 00200 /* GERS (input) REAL array, dimension (2*N) */ 00201 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */ 00202 /* is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should */ 00203 /* be computed from the original UNshifted matrix. */ 00204 00205 /* Z (output) REAL array, dimension (LDZ, max(1,M) ) */ 00206 /* If INFO = 0, the first M columns of Z contain the */ 00207 /* orthonormal eigenvectors of the matrix T */ 00208 /* corresponding to the input eigenvalues, with the i-th */ 00209 /* column of Z holding the eigenvector associated with W(i). */ 00210 /* Note: the user must ensure that at least max(1,M) columns are */ 00211 /* supplied in the array Z. */ 00212 00213 /* LDZ (input) INTEGER */ 00214 /* The leading dimension of the array Z. LDZ >= 1, and if */ 00215 /* JOBZ = 'V', LDZ >= max(1,N). */ 00216 00217 /* ISUPPZ (output) 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 00223 /* WORK (workspace) REAL array, dimension (12*N) */ 00224 00225 /* IWORK (workspace) INTEGER array, dimension (7*N) */ 00226 00227 /* INFO (output) INTEGER */ 00228 /* = 0: successful exit */ 00229 00230 /* > 0: A problem occured in SLARRV. */ 00231 /* < 0: One of the called subroutines signaled an internal problem. */ 00232 /* Needs inspection of the corresponding parameter IINFO */ 00233 /* for further information. */ 00234 00235 /* =-1: Problem in SLARRB when refining a child's eigenvalues. */ 00236 /* =-2: Problem in SLARRF when computing the RRR of a child. */ 00237 /* When a child is inside a tight cluster, it can be difficult */ 00238 /* to find an RRR. A partial remedy from the user's point of */ 00239 /* view is to make the parameter MINRGP smaller and recompile. */ 00240 /* However, as the orthogonality of the computed vectors is */ 00241 /* proportional to 1/MINRGP, the user should be aware that */ 00242 /* he might be trading in precision when he decreases MINRGP. */ 00243 /* =-3: Problem in SLARRB when refining a single eigenvalue */ 00244 /* after the Rayleigh correction was rejected. */ 00245 /* = 5: The Rayleigh Quotient Iteration failed to converge to */ 00246 /* full accuracy in MAXITR steps. */ 00247 00248 /* Further Details */ 00249 /* =============== */ 00250 00251 /* Based on contributions by */ 00252 /* Beresford Parlett, University of California, Berkeley, USA */ 00253 /* Jim Demmel, University of California, Berkeley, USA */ 00254 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00255 /* Osni Marques, LBNL/NERSC, USA */ 00256 /* Christof Voemel, University of California, Berkeley, USA */ 00257 00258 /* ===================================================================== */ 00259 00260 /* .. Parameters .. */ 00261 /* .. */ 00262 /* .. Local Scalars .. */ 00263 /* .. */ 00264 /* .. External Functions .. */ 00265 /* .. */ 00266 /* .. External Subroutines .. */ 00267 /* .. */ 00268 /* .. Intrinsic Functions .. */ 00269 /* .. */ 00270 /* .. Executable Statements .. */ 00271 /* .. */ 00272 /* The first N entries of WORK are reserved for the eigenvalues */ 00273 /* Parameter adjustments */ 00274 --d__; 00275 --l; 00276 --isplit; 00277 --w; 00278 --werr; 00279 --wgap; 00280 --iblock; 00281 --indexw; 00282 --gers; 00283 z_dim1 = *ldz; 00284 z_offset = 1 + z_dim1; 00285 z__ -= z_offset; 00286 --isuppz; 00287 --work; 00288 --iwork; 00289 00290 /* Function Body */ 00291 indld = *n + 1; 00292 indlld = (*n << 1) + 1; 00293 indwrk = *n * 3 + 1; 00294 minwsize = *n * 12; 00295 i__1 = minwsize; 00296 for (i__ = 1; i__ <= i__1; ++i__) { 00297 work[i__] = 0.f; 00298 /* L5: */ 00299 } 00300 /* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the */ 00301 /* factorization used to compute the FP vector */ 00302 iindr = 0; 00303 /* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current */ 00304 /* layer and the one above. */ 00305 iindc1 = *n; 00306 iindc2 = *n << 1; 00307 iindwk = *n * 3 + 1; 00308 miniwsize = *n * 7; 00309 i__1 = miniwsize; 00310 for (i__ = 1; i__ <= i__1; ++i__) { 00311 iwork[i__] = 0; 00312 /* L10: */ 00313 } 00314 zusedl = 1; 00315 if (*dol > 1) { 00316 /* Set lower bound for use of Z */ 00317 zusedl = *dol - 1; 00318 } 00319 zusedu = *m; 00320 if (*dou < *m) { 00321 /* Set lower bound for use of Z */ 00322 zusedu = *dou + 1; 00323 } 00324 /* The width of the part of Z that is used */ 00325 zusedw = zusedu - zusedl + 1; 00326 slaset_("Full", n, &zusedw, &c_b5, &c_b5, &z__[zusedl * z_dim1 + 1], ldz); 00327 eps = slamch_("Precision"); 00328 rqtol = eps * 2.f; 00329 00330 /* Set expert flags for standard code. */ 00331 tryrqc = TRUE_; 00332 if (*dol == 1 && *dou == *m) { 00333 } else { 00334 /* Only selected eigenpairs are computed. Since the other evalues */ 00335 /* are not refined by RQ iteration, bisection has to compute to full */ 00336 /* accuracy. */ 00337 *rtol1 = eps * 4.f; 00338 *rtol2 = eps * 4.f; 00339 } 00340 /* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the */ 00341 /* desired eigenvalues. The support of the nonzero eigenvector */ 00342 /* entries is contained in the interval IBEGIN:IEND. */ 00343 /* Remark that if k eigenpairs are desired, then the eigenvectors */ 00344 /* are stored in k contiguous columns of Z. */ 00345 /* DONE is the number of eigenvectors already computed */ 00346 done = 0; 00347 ibegin = 1; 00348 wbegin = 1; 00349 i__1 = iblock[*m]; 00350 for (jblk = 1; jblk <= i__1; ++jblk) { 00351 iend = isplit[jblk]; 00352 sigma = l[iend]; 00353 /* Find the eigenvectors of the submatrix indexed IBEGIN */ 00354 /* through IEND. */ 00355 wend = wbegin - 1; 00356 L15: 00357 if (wend < *m) { 00358 if (iblock[wend + 1] == jblk) { 00359 ++wend; 00360 goto L15; 00361 } 00362 } 00363 if (wend < wbegin) { 00364 ibegin = iend + 1; 00365 goto L170; 00366 } else if (wend < *dol || wbegin > *dou) { 00367 ibegin = iend + 1; 00368 wbegin = wend + 1; 00369 goto L170; 00370 } 00371 /* Find local spectral diameter of the block */ 00372 gl = gers[(ibegin << 1) - 1]; 00373 gu = gers[ibegin * 2]; 00374 i__2 = iend; 00375 for (i__ = ibegin + 1; i__ <= i__2; ++i__) { 00376 /* Computing MIN */ 00377 r__1 = gers[(i__ << 1) - 1]; 00378 gl = dmin(r__1,gl); 00379 /* Computing MAX */ 00380 r__1 = gers[i__ * 2]; 00381 gu = dmax(r__1,gu); 00382 /* L20: */ 00383 } 00384 spdiam = gu - gl; 00385 /* OLDIEN is the last index of the previous block */ 00386 oldien = ibegin - 1; 00387 /* Calculate the size of the current block */ 00388 in = iend - ibegin + 1; 00389 /* The number of eigenvalues in the current block */ 00390 im = wend - wbegin + 1; 00391 /* This is for a 1x1 block */ 00392 if (ibegin == iend) { 00393 ++done; 00394 z__[ibegin + wbegin * z_dim1] = 1.f; 00395 isuppz[(wbegin << 1) - 1] = ibegin; 00396 isuppz[wbegin * 2] = ibegin; 00397 w[wbegin] += sigma; 00398 work[wbegin] = w[wbegin]; 00399 ibegin = iend + 1; 00400 ++wbegin; 00401 goto L170; 00402 } 00403 /* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) */ 00404 /* Note that these can be approximations, in this case, the corresp. */ 00405 /* entries of WERR give the size of the uncertainty interval. */ 00406 /* The eigenvalue approximations will be refined when necessary as */ 00407 /* high relative accuracy is required for the computation of the */ 00408 /* corresponding eigenvectors. */ 00409 scopy_(&im, &w[wbegin], &c__1, &work[wbegin], &c__1); 00410 /* We store in W the eigenvalue approximations w.r.t. the original */ 00411 /* matrix T. */ 00412 i__2 = im; 00413 for (i__ = 1; i__ <= i__2; ++i__) { 00414 w[wbegin + i__ - 1] += sigma; 00415 /* L30: */ 00416 } 00417 /* NDEPTH is the current depth of the representation tree */ 00418 ndepth = 0; 00419 /* PARITY is either 1 or 0 */ 00420 parity = 1; 00421 /* NCLUS is the number of clusters for the next level of the */ 00422 /* representation tree, we start with NCLUS = 1 for the root */ 00423 nclus = 1; 00424 iwork[iindc1 + 1] = 1; 00425 iwork[iindc1 + 2] = im; 00426 /* IDONE is the number of eigenvectors already computed in the current */ 00427 /* block */ 00428 idone = 0; 00429 /* loop while( IDONE.LT.IM ) */ 00430 /* generate the representation tree for the current block and */ 00431 /* compute the eigenvectors */ 00432 L40: 00433 if (idone < im) { 00434 /* This is a crude protection against infinitely deep trees */ 00435 if (ndepth > *m) { 00436 *info = -2; 00437 return 0; 00438 } 00439 /* breadth first processing of the current level of the representation */ 00440 /* tree: OLDNCL = number of clusters on current level */ 00441 oldncl = nclus; 00442 /* reset NCLUS to count the number of child clusters */ 00443 nclus = 0; 00444 00445 parity = 1 - parity; 00446 if (parity == 0) { 00447 oldcls = iindc1; 00448 newcls = iindc2; 00449 } else { 00450 oldcls = iindc2; 00451 newcls = iindc1; 00452 } 00453 /* Process the clusters on the current level */ 00454 i__2 = oldncl; 00455 for (i__ = 1; i__ <= i__2; ++i__) { 00456 j = oldcls + (i__ << 1); 00457 /* OLDFST, OLDLST = first, last index of current cluster. */ 00458 /* cluster indices start with 1 and are relative */ 00459 /* to WBEGIN when accessing W, WGAP, WERR, Z */ 00460 oldfst = iwork[j - 1]; 00461 oldlst = iwork[j]; 00462 if (ndepth > 0) { 00463 /* Retrieve relatively robust representation (RRR) of cluster */ 00464 /* that has been computed at the previous level */ 00465 /* The RRR is stored in Z and overwritten once the eigenvectors */ 00466 /* have been computed or when the cluster is refined */ 00467 if (*dol == 1 && *dou == *m) { 00468 /* Get representation from location of the leftmost evalue */ 00469 /* of the cluster */ 00470 j = wbegin + oldfst - 1; 00471 } else { 00472 if (wbegin + oldfst - 1 < *dol) { 00473 /* Get representation from the left end of Z array */ 00474 j = *dol - 1; 00475 } else if (wbegin + oldfst - 1 > *dou) { 00476 /* Get representation from the right end of Z array */ 00477 j = *dou; 00478 } else { 00479 j = wbegin + oldfst - 1; 00480 } 00481 } 00482 scopy_(&in, &z__[ibegin + j * z_dim1], &c__1, &d__[ibegin] 00483 , &c__1); 00484 i__3 = in - 1; 00485 scopy_(&i__3, &z__[ibegin + (j + 1) * z_dim1], &c__1, &l[ 00486 ibegin], &c__1); 00487 sigma = z__[iend + (j + 1) * z_dim1]; 00488 /* Set the corresponding entries in Z to zero */ 00489 slaset_("Full", &in, &c__2, &c_b5, &c_b5, &z__[ibegin + j 00490 * z_dim1], ldz); 00491 } 00492 /* Compute DL and DLL of current RRR */ 00493 i__3 = iend - 1; 00494 for (j = ibegin; j <= i__3; ++j) { 00495 tmp = d__[j] * l[j]; 00496 work[indld - 1 + j] = tmp; 00497 work[indlld - 1 + j] = tmp * l[j]; 00498 /* L50: */ 00499 } 00500 if (ndepth > 0) { 00501 /* P and Q are index of the first and last eigenvalue to compute */ 00502 /* within the current block */ 00503 p = indexw[wbegin - 1 + oldfst]; 00504 q = indexw[wbegin - 1 + oldlst]; 00505 /* Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET */ 00506 /* thru' Q-OFFSET elements of these arrays are to be used. */ 00507 /* OFFSET = P-OLDFST */ 00508 offset = indexw[wbegin] - 1; 00509 /* perform limited bisection (if necessary) to get approximate */ 00510 /* eigenvalues to the precision needed. */ 00511 slarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &p, 00512 &q, rtol1, rtol2, &offset, &work[wbegin], &wgap[ 00513 wbegin], &werr[wbegin], &work[indwrk], &iwork[ 00514 iindwk], pivmin, &spdiam, &in, &iinfo); 00515 if (iinfo != 0) { 00516 *info = -1; 00517 return 0; 00518 } 00519 /* We also recompute the extremal gaps. W holds all eigenvalues */ 00520 /* of the unshifted matrix and must be used for computation */ 00521 /* of WGAP, the entries of WORK might stem from RRRs with */ 00522 /* different shifts. The gaps from WBEGIN-1+OLDFST to */ 00523 /* WBEGIN-1+OLDLST are correctly computed in SLARRB. */ 00524 /* However, we only allow the gaps to become greater since */ 00525 /* this is what should happen when we decrease WERR */ 00526 if (oldfst > 1) { 00527 /* Computing MAX */ 00528 r__1 = wgap[wbegin + oldfst - 2], r__2 = w[wbegin + 00529 oldfst - 1] - werr[wbegin + oldfst - 1] - w[ 00530 wbegin + oldfst - 2] - werr[wbegin + oldfst - 00531 2]; 00532 wgap[wbegin + oldfst - 2] = dmax(r__1,r__2); 00533 } 00534 if (wbegin + oldlst - 1 < wend) { 00535 /* Computing MAX */ 00536 r__1 = wgap[wbegin + oldlst - 1], r__2 = w[wbegin + 00537 oldlst] - werr[wbegin + oldlst] - w[wbegin + 00538 oldlst - 1] - werr[wbegin + oldlst - 1]; 00539 wgap[wbegin + oldlst - 1] = dmax(r__1,r__2); 00540 } 00541 /* Each time the eigenvalues in WORK get refined, we store */ 00542 /* the newly found approximation with all shifts applied in W */ 00543 i__3 = oldlst; 00544 for (j = oldfst; j <= i__3; ++j) { 00545 w[wbegin + j - 1] = work[wbegin + j - 1] + sigma; 00546 /* L53: */ 00547 } 00548 } 00549 /* Process the current node. */ 00550 newfst = oldfst; 00551 i__3 = oldlst; 00552 for (j = oldfst; j <= i__3; ++j) { 00553 if (j == oldlst) { 00554 /* we are at the right end of the cluster, this is also the */ 00555 /* boundary of the child cluster */ 00556 newlst = j; 00557 } else if (wgap[wbegin + j - 1] >= *minrgp * (r__1 = work[ 00558 wbegin + j - 1], dabs(r__1))) { 00559 /* the right relative gap is big enough, the child cluster */ 00560 /* (NEWFST,..,NEWLST) is well separated from the following */ 00561 newlst = j; 00562 } else { 00563 /* inside a child cluster, the relative gap is not */ 00564 /* big enough. */ 00565 goto L140; 00566 } 00567 /* Compute size of child cluster found */ 00568 newsiz = newlst - newfst + 1; 00569 /* NEWFTT is the place in Z where the new RRR or the computed */ 00570 /* eigenvector is to be stored */ 00571 if (*dol == 1 && *dou == *m) { 00572 /* Store representation at location of the leftmost evalue */ 00573 /* of the cluster */ 00574 newftt = wbegin + newfst - 1; 00575 } else { 00576 if (wbegin + newfst - 1 < *dol) { 00577 /* Store representation at the left end of Z array */ 00578 newftt = *dol - 1; 00579 } else if (wbegin + newfst - 1 > *dou) { 00580 /* Store representation at the right end of Z array */ 00581 newftt = *dou; 00582 } else { 00583 newftt = wbegin + newfst - 1; 00584 } 00585 } 00586 if (newsiz > 1) { 00587 00588 /* Current child is not a singleton but a cluster. */ 00589 /* Compute and store new representation of child. */ 00590 00591 00592 /* Compute left and right cluster gap. */ 00593 00594 /* LGAP and RGAP are not computed from WORK because */ 00595 /* the eigenvalue approximations may stem from RRRs */ 00596 /* different shifts. However, W hold all eigenvalues */ 00597 /* of the unshifted matrix. Still, the entries in WGAP */ 00598 /* have to be computed from WORK since the entries */ 00599 /* in W might be of the same order so that gaps are not */ 00600 /* exhibited correctly for very close eigenvalues. */ 00601 if (newfst == 1) { 00602 /* Computing MAX */ 00603 r__1 = 0.f, r__2 = w[wbegin] - werr[wbegin] - *vl; 00604 lgap = dmax(r__1,r__2); 00605 } else { 00606 lgap = wgap[wbegin + newfst - 2]; 00607 } 00608 rgap = wgap[wbegin + newlst - 1]; 00609 00610 /* Compute left- and rightmost eigenvalue of child */ 00611 /* to high precision in order to shift as close */ 00612 /* as possible and obtain as large relative gaps */ 00613 /* as possible */ 00614 00615 for (k = 1; k <= 2; ++k) { 00616 if (k == 1) { 00617 p = indexw[wbegin - 1 + newfst]; 00618 } else { 00619 p = indexw[wbegin - 1 + newlst]; 00620 } 00621 offset = indexw[wbegin] - 1; 00622 slarrb_(&in, &d__[ibegin], &work[indlld + ibegin 00623 - 1], &p, &p, &rqtol, &rqtol, &offset, & 00624 work[wbegin], &wgap[wbegin], &werr[wbegin] 00625 , &work[indwrk], &iwork[iindwk], pivmin, & 00626 spdiam, &in, &iinfo); 00627 /* L55: */ 00628 } 00629 00630 if (wbegin + newlst - 1 < *dol || wbegin + newfst - 1 00631 > *dou) { 00632 /* if the cluster contains no desired eigenvalues */ 00633 /* skip the computation of that branch of the rep. tree */ 00634 00635 /* We could skip before the refinement of the extremal */ 00636 /* eigenvalues of the child, but then the representation */ 00637 /* tree could be different from the one when nothing is */ 00638 /* skipped. For this reason we skip at this place. */ 00639 idone = idone + newlst - newfst + 1; 00640 goto L139; 00641 } 00642 00643 /* Compute RRR of child cluster. */ 00644 /* Note that the new RRR is stored in Z */ 00645 00646 /* SLARRF needs LWORK = 2*N */ 00647 slarrf_(&in, &d__[ibegin], &l[ibegin], &work[indld + 00648 ibegin - 1], &newfst, &newlst, &work[wbegin], 00649 &wgap[wbegin], &werr[wbegin], &spdiam, &lgap, 00650 &rgap, pivmin, &tau, &z__[ibegin + newftt * 00651 z_dim1], &z__[ibegin + (newftt + 1) * z_dim1], 00652 &work[indwrk], &iinfo); 00653 if (iinfo == 0) { 00654 /* a new RRR for the cluster was found by SLARRF */ 00655 /* update shift and store it */ 00656 ssigma = sigma + tau; 00657 z__[iend + (newftt + 1) * z_dim1] = ssigma; 00658 /* WORK() are the midpoints and WERR() the semi-width */ 00659 /* Note that the entries in W are unchanged. */ 00660 i__4 = newlst; 00661 for (k = newfst; k <= i__4; ++k) { 00662 fudge = eps * 3.f * (r__1 = work[wbegin + k - 00663 1], dabs(r__1)); 00664 work[wbegin + k - 1] -= tau; 00665 fudge += eps * 4.f * (r__1 = work[wbegin + k 00666 - 1], dabs(r__1)); 00667 /* Fudge errors */ 00668 werr[wbegin + k - 1] += fudge; 00669 /* Gaps are not fudged. Provided that WERR is small */ 00670 /* when eigenvalues are close, a zero gap indicates */ 00671 /* that a new representation is needed for resolving */ 00672 /* the cluster. A fudge could lead to a wrong decision */ 00673 /* of judging eigenvalues 'separated' which in */ 00674 /* reality are not. This could have a negative impact */ 00675 /* on the orthogonality of the computed eigenvectors. */ 00676 /* L116: */ 00677 } 00678 ++nclus; 00679 k = newcls + (nclus << 1); 00680 iwork[k - 1] = newfst; 00681 iwork[k] = newlst; 00682 } else { 00683 *info = -2; 00684 return 0; 00685 } 00686 } else { 00687 00688 /* Compute eigenvector of singleton */ 00689 00690 iter = 0; 00691 00692 tol = log((real) in) * 4.f * eps; 00693 00694 k = newfst; 00695 windex = wbegin + k - 1; 00696 /* Computing MAX */ 00697 i__4 = windex - 1; 00698 windmn = max(i__4,1); 00699 /* Computing MIN */ 00700 i__4 = windex + 1; 00701 windpl = min(i__4,*m); 00702 lambda = work[windex]; 00703 ++done; 00704 /* Check if eigenvector computation is to be skipped */ 00705 if (windex < *dol || windex > *dou) { 00706 eskip = TRUE_; 00707 goto L125; 00708 } else { 00709 eskip = FALSE_; 00710 } 00711 left = work[windex] - werr[windex]; 00712 right = work[windex] + werr[windex]; 00713 indeig = indexw[windex]; 00714 /* Note that since we compute the eigenpairs for a child, */ 00715 /* all eigenvalue approximations are w.r.t the same shift. */ 00716 /* In this case, the entries in WORK should be used for */ 00717 /* computing the gaps since they exhibit even very small */ 00718 /* differences in the eigenvalues, as opposed to the */ 00719 /* entries in W which might "look" the same. */ 00720 if (k == 1) { 00721 /* In the case RANGE='I' and with not much initial */ 00722 /* accuracy in LAMBDA and VL, the formula */ 00723 /* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) */ 00724 /* can lead to an overestimation of the left gap and */ 00725 /* thus to inadequately early RQI 'convergence'. */ 00726 /* Prevent this by forcing a small left gap. */ 00727 /* Computing MAX */ 00728 r__1 = dabs(left), r__2 = dabs(right); 00729 lgap = eps * dmax(r__1,r__2); 00730 } else { 00731 lgap = wgap[windmn]; 00732 } 00733 if (k == im) { 00734 /* In the case RANGE='I' and with not much initial */ 00735 /* accuracy in LAMBDA and VU, the formula */ 00736 /* can lead to an overestimation of the right gap and */ 00737 /* thus to inadequately early RQI 'convergence'. */ 00738 /* Prevent this by forcing a small right gap. */ 00739 /* Computing MAX */ 00740 r__1 = dabs(left), r__2 = dabs(right); 00741 rgap = eps * dmax(r__1,r__2); 00742 } else { 00743 rgap = wgap[windex]; 00744 } 00745 gap = dmin(lgap,rgap); 00746 if (k == 1 || k == im) { 00747 /* The eigenvector support can become wrong */ 00748 /* because significant entries could be cut off due to a */ 00749 /* large GAPTOL parameter in LAR1V. Prevent this. */ 00750 gaptol = 0.f; 00751 } else { 00752 gaptol = gap * eps; 00753 } 00754 isupmn = in; 00755 isupmx = 1; 00756 /* Update WGAP so that it holds the minimum gap */ 00757 /* to the left or the right. This is crucial in the */ 00758 /* case where bisection is used to ensure that the */ 00759 /* eigenvalue is refined up to the required precision. */ 00760 /* The correct value is restored afterwards. */ 00761 savgap = wgap[windex]; 00762 wgap[windex] = gap; 00763 /* We want to use the Rayleigh Quotient Correction */ 00764 /* as often as possible since it converges quadratically */ 00765 /* when we are close enough to the desired eigenvalue. */ 00766 /* However, the Rayleigh Quotient can have the wrong sign */ 00767 /* and lead us away from the desired eigenvalue. In this */ 00768 /* case, the best we can do is to use bisection. */ 00769 usedbs = FALSE_; 00770 usedrq = FALSE_; 00771 /* Bisection is initially turned off unless it is forced */ 00772 needbs = ! tryrqc; 00773 L120: 00774 /* Check if bisection should be used to refine eigenvalue */ 00775 if (needbs) { 00776 /* Take the bisection as new iterate */ 00777 usedbs = TRUE_; 00778 itmp1 = iwork[iindr + windex]; 00779 offset = indexw[wbegin] - 1; 00780 r__1 = eps * 2.f; 00781 slarrb_(&in, &d__[ibegin], &work[indlld + ibegin 00782 - 1], &indeig, &indeig, &c_b5, &r__1, & 00783 offset, &work[wbegin], &wgap[wbegin], & 00784 werr[wbegin], &work[indwrk], &iwork[ 00785 iindwk], pivmin, &spdiam, &itmp1, &iinfo); 00786 if (iinfo != 0) { 00787 *info = -3; 00788 return 0; 00789 } 00790 lambda = work[windex]; 00791 /* Reset twist index from inaccurate LAMBDA to */ 00792 /* force computation of true MINGMA */ 00793 iwork[iindr + windex] = 0; 00794 } 00795 /* Given LAMBDA, compute the eigenvector. */ 00796 L__1 = ! usedbs; 00797 slar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &l[ 00798 ibegin], &work[indld + ibegin - 1], &work[ 00799 indlld + ibegin - 1], pivmin, &gaptol, &z__[ 00800 ibegin + windex * z_dim1], &L__1, &negcnt, & 00801 ztz, &mingma, &iwork[iindr + windex], &isuppz[ 00802 (windex << 1) - 1], &nrminv, &resid, &rqcorr, 00803 &work[indwrk]); 00804 if (iter == 0) { 00805 bstres = resid; 00806 bstw = lambda; 00807 } else if (resid < bstres) { 00808 bstres = resid; 00809 bstw = lambda; 00810 } 00811 /* Computing MIN */ 00812 i__4 = isupmn, i__5 = isuppz[(windex << 1) - 1]; 00813 isupmn = min(i__4,i__5); 00814 /* Computing MAX */ 00815 i__4 = isupmx, i__5 = isuppz[windex * 2]; 00816 isupmx = max(i__4,i__5); 00817 ++iter; 00818 /* sin alpha <= |resid|/gap */ 00819 /* Note that both the residual and the gap are */ 00820 /* proportional to the matrix, so ||T|| doesn't play */ 00821 /* a role in the quotient */ 00822 00823 /* Convergence test for Rayleigh-Quotient iteration */ 00824 /* (omitted when Bisection has been used) */ 00825 00826 if (resid > tol * gap && dabs(rqcorr) > rqtol * dabs( 00827 lambda) && ! usedbs) { 00828 /* We need to check that the RQCORR update doesn't */ 00829 /* move the eigenvalue away from the desired one and */ 00830 /* towards a neighbor. -> protection with bisection */ 00831 if (indeig <= negcnt) { 00832 /* The wanted eigenvalue lies to the left */ 00833 sgndef = -1.f; 00834 } else { 00835 /* The wanted eigenvalue lies to the right */ 00836 sgndef = 1.f; 00837 } 00838 /* We only use the RQCORR if it improves the */ 00839 /* the iterate reasonably. */ 00840 if (rqcorr * sgndef >= 0.f && lambda + rqcorr <= 00841 right && lambda + rqcorr >= left) { 00842 usedrq = TRUE_; 00843 /* Store new midpoint of bisection interval in WORK */ 00844 if (sgndef == 1.f) { 00845 /* The current LAMBDA is on the left of the true */ 00846 /* eigenvalue */ 00847 left = lambda; 00848 /* We prefer to assume that the error estimate */ 00849 /* is correct. We could make the interval not */ 00850 /* as a bracket but to be modified if the RQCORR */ 00851 /* chooses to. In this case, the RIGHT side should */ 00852 /* be modified as follows: */ 00853 /* RIGHT = MAX(RIGHT, LAMBDA + RQCORR) */ 00854 } else { 00855 /* The current LAMBDA is on the right of the true */ 00856 /* eigenvalue */ 00857 right = lambda; 00858 /* See comment about assuming the error estimate is */ 00859 /* correct above. */ 00860 /* LEFT = MIN(LEFT, LAMBDA + RQCORR) */ 00861 } 00862 work[windex] = (right + left) * .5f; 00863 /* Take RQCORR since it has the correct sign and */ 00864 /* improves the iterate reasonably */ 00865 lambda += rqcorr; 00866 /* Update width of error interval */ 00867 werr[windex] = (right - left) * .5f; 00868 } else { 00869 needbs = TRUE_; 00870 } 00871 if (right - left < rqtol * dabs(lambda)) { 00872 /* The eigenvalue is computed to bisection accuracy */ 00873 /* compute eigenvector and stop */ 00874 usedbs = TRUE_; 00875 goto L120; 00876 } else if (iter < 10) { 00877 goto L120; 00878 } else if (iter == 10) { 00879 needbs = TRUE_; 00880 goto L120; 00881 } else { 00882 *info = 5; 00883 return 0; 00884 } 00885 } else { 00886 stp2ii = FALSE_; 00887 if (usedrq && usedbs && bstres <= resid) { 00888 lambda = bstw; 00889 stp2ii = TRUE_; 00890 } 00891 if (stp2ii) { 00892 /* improve error angle by second step */ 00893 L__1 = ! usedbs; 00894 slar1v_(&in, &c__1, &in, &lambda, &d__[ibegin] 00895 , &l[ibegin], &work[indld + ibegin - 00896 1], &work[indlld + ibegin - 1], 00897 pivmin, &gaptol, &z__[ibegin + windex 00898 * z_dim1], &L__1, &negcnt, &ztz, & 00899 mingma, &iwork[iindr + windex], & 00900 isuppz[(windex << 1) - 1], &nrminv, & 00901 resid, &rqcorr, &work[indwrk]); 00902 } 00903 work[windex] = lambda; 00904 } 00905 00906 /* Compute FP-vector support w.r.t. whole matrix */ 00907 00908 isuppz[(windex << 1) - 1] += oldien; 00909 isuppz[windex * 2] += oldien; 00910 zfrom = isuppz[(windex << 1) - 1]; 00911 zto = isuppz[windex * 2]; 00912 isupmn += oldien; 00913 isupmx += oldien; 00914 /* Ensure vector is ok if support in the RQI has changed */ 00915 if (isupmn < zfrom) { 00916 i__4 = zfrom - 1; 00917 for (ii = isupmn; ii <= i__4; ++ii) { 00918 z__[ii + windex * z_dim1] = 0.f; 00919 /* L122: */ 00920 } 00921 } 00922 if (isupmx > zto) { 00923 i__4 = isupmx; 00924 for (ii = zto + 1; ii <= i__4; ++ii) { 00925 z__[ii + windex * z_dim1] = 0.f; 00926 /* L123: */ 00927 } 00928 } 00929 i__4 = zto - zfrom + 1; 00930 sscal_(&i__4, &nrminv, &z__[zfrom + windex * z_dim1], 00931 &c__1); 00932 L125: 00933 /* Update W */ 00934 w[windex] = lambda + sigma; 00935 /* Recompute the gaps on the left and right */ 00936 /* But only allow them to become larger and not */ 00937 /* smaller (which can only happen through "bad" */ 00938 /* cancellation and doesn't reflect the theory */ 00939 /* where the initial gaps are underestimated due */ 00940 /* to WERR being too crude.) */ 00941 if (! eskip) { 00942 if (k > 1) { 00943 /* Computing MAX */ 00944 r__1 = wgap[windmn], r__2 = w[windex] - werr[ 00945 windex] - w[windmn] - werr[windmn]; 00946 wgap[windmn] = dmax(r__1,r__2); 00947 } 00948 if (windex < wend) { 00949 /* Computing MAX */ 00950 r__1 = savgap, r__2 = w[windpl] - werr[windpl] 00951 - w[windex] - werr[windex]; 00952 wgap[windex] = dmax(r__1,r__2); 00953 } 00954 } 00955 ++idone; 00956 } 00957 /* here ends the code for the current child */ 00958 00959 L139: 00960 /* Proceed to any remaining child nodes */ 00961 newfst = j + 1; 00962 L140: 00963 ; 00964 } 00965 /* L150: */ 00966 } 00967 ++ndepth; 00968 goto L40; 00969 } 00970 ibegin = iend + 1; 00971 wbegin = wend + 1; 00972 L170: 00973 ; 00974 } 00975 00976 return 0; 00977 00978 /* End of SLARRV */ 00979 00980 } /* slarrv_ */