00001 /* slarrd.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 integer c__1 = 1; 00019 static integer c_n1 = -1; 00020 static integer c__3 = 3; 00021 static integer c__2 = 2; 00022 static integer c__0 = 0; 00023 00024 /* Subroutine */ int slarrd_(char *range, char *order, integer *n, real *vl, 00025 real *vu, integer *il, integer *iu, real *gers, real *reltol, real * 00026 d__, real *e, real *e2, real *pivmin, integer *nsplit, integer * 00027 isplit, integer *m, real *w, real *werr, real *wl, real *wu, integer * 00028 iblock, integer *indexw, real *work, integer *iwork, integer *info) 00029 { 00030 /* System generated locals */ 00031 integer i__1, i__2, i__3; 00032 real r__1, r__2; 00033 00034 /* Builtin functions */ 00035 double log(doublereal); 00036 00037 /* Local variables */ 00038 integer i__, j, ib, ie, je, nb; 00039 real gl; 00040 integer im, in; 00041 real gu; 00042 integer iw, jee; 00043 real eps; 00044 integer nwl; 00045 real wlu, wul; 00046 integer nwu; 00047 real tmp1, tmp2; 00048 integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc; 00049 extern logical lsame_(char *, char *); 00050 integer iinfo; 00051 real atoli; 00052 integer iwoff, itmax; 00053 real wkill, rtoli, uflow, tnorm; 00054 integer ibegin, irange, idiscl; 00055 extern doublereal slamch_(char *); 00056 integer idumma[1]; 00057 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00058 integer *, integer *); 00059 integer idiscu; 00060 extern /* Subroutine */ int slaebz_(integer *, integer *, integer *, 00061 integer *, integer *, integer *, real *, real *, real *, real *, 00062 real *, real *, integer *, real *, real *, integer *, integer *, 00063 real *, integer *, integer *); 00064 logical ncnvrg, toofew; 00065 00066 00067 /* -- LAPACK auxiliary routine (version 3.2.1) -- */ 00068 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00069 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ 00070 /* -- April 2009 -- */ 00071 00072 /* .. Scalar Arguments .. */ 00073 /* .. */ 00074 /* .. Array Arguments .. */ 00075 /* .. */ 00076 00077 /* Purpose */ 00078 /* ======= */ 00079 00080 /* SLARRD computes the eigenvalues of a symmetric tridiagonal */ 00081 /* matrix T to suitable accuracy. This is an auxiliary code to be */ 00082 /* called from SSTEMR. */ 00083 /* The user may ask for all eigenvalues, all eigenvalues */ 00084 /* in the half-open interval (VL, VU], or the IL-th through IU-th */ 00085 /* eigenvalues. */ 00086 00087 /* To avoid overflow, the matrix must be scaled so that its */ 00088 /* largest element is no greater than overflow**(1/2) * */ 00089 /* underflow**(1/4) in absolute value, and for greatest */ 00090 /* accuracy, it should not be much smaller than that. */ 00091 00092 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */ 00093 /* Matrix", Report CS41, Computer Science Dept., Stanford */ 00094 /* University, July 21, 1966. */ 00095 00096 /* Arguments */ 00097 /* ========= */ 00098 00099 /* RANGE (input) CHARACTER */ 00100 /* = 'A': ("All") all eigenvalues will be found. */ 00101 /* = 'V': ("Value") all eigenvalues in the half-open interval */ 00102 /* (VL, VU] will be found. */ 00103 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */ 00104 /* entire matrix) will be found. */ 00105 00106 /* ORDER (input) CHARACTER */ 00107 /* = 'B': ("By Block") the eigenvalues will be grouped by */ 00108 /* split-off block (see IBLOCK, ISPLIT) and */ 00109 /* ordered from smallest to largest within */ 00110 /* the block. */ 00111 /* = 'E': ("Entire matrix") */ 00112 /* the eigenvalues for the entire matrix */ 00113 /* will be ordered from smallest to */ 00114 /* largest. */ 00115 00116 /* N (input) INTEGER */ 00117 /* The order of the tridiagonal matrix T. N >= 0. */ 00118 00119 /* VL (input) REAL */ 00120 /* VU (input) REAL */ 00121 /* If RANGE='V', the lower and upper bounds of the interval to */ 00122 /* be searched for eigenvalues. Eigenvalues less than or equal */ 00123 /* to VL, or greater than VU, will not be returned. VL < VU. */ 00124 /* Not referenced if RANGE = 'A' or 'I'. */ 00125 00126 /* IL (input) INTEGER */ 00127 /* IU (input) INTEGER */ 00128 /* If RANGE='I', the indices (in ascending order) of the */ 00129 /* smallest and largest eigenvalues to be returned. */ 00130 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */ 00131 /* Not referenced if RANGE = 'A' or 'V'. */ 00132 00133 /* GERS (input) REAL array, dimension (2*N) */ 00134 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */ 00135 /* is (GERS(2*i-1), GERS(2*i)). */ 00136 00137 /* RELTOL (input) REAL */ 00138 /* The minimum relative width of an interval. When an interval */ 00139 /* is narrower than RELTOL times the larger (in */ 00140 /* magnitude) endpoint, then it is considered to be */ 00141 /* sufficiently small, i.e., converged. Note: this should */ 00142 /* always be at least radix*machine epsilon. */ 00143 00144 /* D (input) REAL array, dimension (N) */ 00145 /* The n diagonal elements of the tridiagonal matrix T. */ 00146 00147 /* E (input) REAL array, dimension (N-1) */ 00148 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */ 00149 00150 /* E2 (input) REAL array, dimension (N-1) */ 00151 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */ 00152 00153 /* PIVMIN (input) REAL */ 00154 /* The minimum pivot allowed in the Sturm sequence for T. */ 00155 00156 /* NSPLIT (input) INTEGER */ 00157 /* The number of diagonal blocks in the matrix T. */ 00158 /* 1 <= NSPLIT <= N. */ 00159 00160 /* ISPLIT (input) INTEGER array, dimension (N) */ 00161 /* The splitting points, at which T breaks up into submatrices. */ 00162 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */ 00163 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */ 00164 /* etc., and the NSPLIT-th consists of rows/columns */ 00165 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */ 00166 /* (Only the first NSPLIT elements will actually be used, but */ 00167 /* since the user cannot know a priori what value NSPLIT will */ 00168 /* have, N words must be reserved for ISPLIT.) */ 00169 00170 /* M (output) INTEGER */ 00171 /* The actual number of eigenvalues found. 0 <= M <= N. */ 00172 /* (See also the description of INFO=2,3.) */ 00173 00174 /* W (output) REAL array, dimension (N) */ 00175 /* On exit, the first M elements of W will contain the */ 00176 /* eigenvalue approximations. SLARRD computes an interval */ 00177 /* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */ 00178 /* approximation is given as the interval midpoint */ 00179 /* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */ 00180 /* WERR(j) = abs( a_j - b_j)/2 */ 00181 00182 /* WERR (output) REAL array, dimension (N) */ 00183 /* The error bound on the corresponding eigenvalue approximation */ 00184 /* in W. */ 00185 00186 /* WL (output) REAL */ 00187 /* WU (output) REAL */ 00188 /* The interval (WL, WU] contains all the wanted eigenvalues. */ 00189 /* If RANGE='V', then WL=VL and WU=VU. */ 00190 /* If RANGE='A', then WL and WU are the global Gerschgorin bounds */ 00191 /* on the spectrum. */ 00192 /* If RANGE='I', then WL and WU are computed by SLAEBZ from the */ 00193 /* index range specified. */ 00194 00195 /* IBLOCK (output) INTEGER array, dimension (N) */ 00196 /* At each row/column j where E(j) is zero or small, the */ 00197 /* matrix T is considered to split into a block diagonal */ 00198 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */ 00199 /* block (from 1 to the number of blocks) the eigenvalue W(i) */ 00200 /* belongs. (SLARRD may use the remaining N-M elements as */ 00201 /* workspace.) */ 00202 00203 /* INDEXW (output) INTEGER array, dimension (N) */ 00204 /* The indices of the eigenvalues within each block (submatrix); */ 00205 /* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */ 00206 /* i-th eigenvalue W(i) is the j-th eigenvalue in block k. */ 00207 00208 /* WORK (workspace) REAL array, dimension (4*N) */ 00209 00210 /* IWORK (workspace) INTEGER array, dimension (3*N) */ 00211 00212 /* INFO (output) INTEGER */ 00213 /* = 0: successful exit */ 00214 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00215 /* > 0: some or all of the eigenvalues failed to converge or */ 00216 /* were not computed: */ 00217 /* =1 or 3: Bisection failed to converge for some */ 00218 /* eigenvalues; these eigenvalues are flagged by a */ 00219 /* negative block number. The effect is that the */ 00220 /* eigenvalues may not be as accurate as the */ 00221 /* absolute and relative tolerances. This is */ 00222 /* generally caused by unexpectedly inaccurate */ 00223 /* arithmetic. */ 00224 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */ 00225 /* IL:IU were found. */ 00226 /* Effect: M < IU+1-IL */ 00227 /* Cause: non-monotonic arithmetic, causing the */ 00228 /* Sturm sequence to be non-monotonic. */ 00229 /* Cure: recalculate, using RANGE='A', and pick */ 00230 /* out eigenvalues IL:IU. In some cases, */ 00231 /* increasing the PARAMETER "FUDGE" may */ 00232 /* make things work. */ 00233 /* = 4: RANGE='I', and the Gershgorin interval */ 00234 /* initially used was too small. No eigenvalues */ 00235 /* were computed. */ 00236 /* Probable cause: your machine has sloppy */ 00237 /* floating-point arithmetic. */ 00238 /* Cure: Increase the PARAMETER "FUDGE", */ 00239 /* recompile, and try again. */ 00240 00241 /* Internal Parameters */ 00242 /* =================== */ 00243 00244 /* FUDGE REAL , default = 2 */ 00245 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */ 00246 /* a value of 1 should work, but on machines with sloppy */ 00247 /* arithmetic, this needs to be larger. The default for */ 00248 /* publicly released versions should be large enough to handle */ 00249 /* the worst machine around. Note that this has no effect */ 00250 /* on accuracy of the solution. */ 00251 00252 /* Based on contributions by */ 00253 /* W. Kahan, University of California, Berkeley, USA */ 00254 /* Beresford Parlett, University of California, Berkeley, USA */ 00255 /* Jim Demmel, University of California, Berkeley, USA */ 00256 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00257 /* Osni Marques, LBNL/NERSC, USA */ 00258 /* Christof Voemel, University of California, Berkeley, USA */ 00259 00260 /* ===================================================================== */ 00261 00262 /* .. Parameters .. */ 00263 /* .. */ 00264 /* .. Local Scalars .. */ 00265 /* .. */ 00266 /* .. Local Arrays .. */ 00267 /* .. */ 00268 /* .. External Functions .. */ 00269 /* .. */ 00270 /* .. External Subroutines .. */ 00271 /* .. */ 00272 /* .. Intrinsic Functions .. */ 00273 /* .. */ 00274 /* .. Executable Statements .. */ 00275 00276 /* Parameter adjustments */ 00277 --iwork; 00278 --work; 00279 --indexw; 00280 --iblock; 00281 --werr; 00282 --w; 00283 --isplit; 00284 --e2; 00285 --e; 00286 --d__; 00287 --gers; 00288 00289 /* Function Body */ 00290 *info = 0; 00291 00292 /* Decode RANGE */ 00293 00294 if (lsame_(range, "A")) { 00295 irange = 1; 00296 } else if (lsame_(range, "V")) { 00297 irange = 2; 00298 } else if (lsame_(range, "I")) { 00299 irange = 3; 00300 } else { 00301 irange = 0; 00302 } 00303 00304 /* Check for Errors */ 00305 00306 if (irange <= 0) { 00307 *info = -1; 00308 } else if (! (lsame_(order, "B") || lsame_(order, 00309 "E"))) { 00310 *info = -2; 00311 } else if (*n < 0) { 00312 *info = -3; 00313 } else if (irange == 2) { 00314 if (*vl >= *vu) { 00315 *info = -5; 00316 } 00317 } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) { 00318 *info = -6; 00319 } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) { 00320 *info = -7; 00321 } 00322 00323 if (*info != 0) { 00324 return 0; 00325 } 00326 /* Initialize error flags */ 00327 *info = 0; 00328 ncnvrg = FALSE_; 00329 toofew = FALSE_; 00330 /* Quick return if possible */ 00331 *m = 0; 00332 if (*n == 0) { 00333 return 0; 00334 } 00335 /* Simplification: */ 00336 if (irange == 3 && *il == 1 && *iu == *n) { 00337 irange = 1; 00338 } 00339 /* Get machine constants */ 00340 eps = slamch_("P"); 00341 uflow = slamch_("U"); 00342 /* Special Case when N=1 */ 00343 /* Treat case of 1x1 matrix for quick return */ 00344 if (*n == 1) { 00345 if (irange == 1 || irange == 2 && d__[1] > *vl && d__[1] <= *vu || 00346 irange == 3 && *il == 1 && *iu == 1) { 00347 *m = 1; 00348 w[1] = d__[1]; 00349 /* The computation error of the eigenvalue is zero */ 00350 werr[1] = 0.f; 00351 iblock[1] = 1; 00352 indexw[1] = 1; 00353 } 00354 return 0; 00355 } 00356 /* NB is the minimum vector length for vector bisection, or 0 */ 00357 /* if only scalar is to be done. */ 00358 nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1); 00359 if (nb <= 1) { 00360 nb = 0; 00361 } 00362 /* Find global spectral radius */ 00363 gl = d__[1]; 00364 gu = d__[1]; 00365 i__1 = *n; 00366 for (i__ = 1; i__ <= i__1; ++i__) { 00367 /* Computing MIN */ 00368 r__1 = gl, r__2 = gers[(i__ << 1) - 1]; 00369 gl = dmin(r__1,r__2); 00370 /* Computing MAX */ 00371 r__1 = gu, r__2 = gers[i__ * 2]; 00372 gu = dmax(r__1,r__2); 00373 /* L5: */ 00374 } 00375 /* Compute global Gerschgorin bounds and spectral diameter */ 00376 /* Computing MAX */ 00377 r__1 = dabs(gl), r__2 = dabs(gu); 00378 tnorm = dmax(r__1,r__2); 00379 gl = gl - tnorm * 2.f * eps * *n - *pivmin * 4.f; 00380 gu = gu + tnorm * 2.f * eps * *n + *pivmin * 4.f; 00381 /* [JAN/28/2009] remove the line below since SPDIAM variable not use */ 00382 /* SPDIAM = GU - GL */ 00383 /* Input arguments for SLAEBZ: */ 00384 /* The relative tolerance. An interval (a,b] lies within */ 00385 /* "relative tolerance" if b-a < RELTOL*max(|a|,|b|), */ 00386 rtoli = *reltol; 00387 /* Set the absolute tolerance for interval convergence to zero to force */ 00388 /* interval convergence based on relative size of the interval. */ 00389 /* This is dangerous because intervals might not converge when RELTOL is */ 00390 /* small. But at least a very small number should be selected so that for */ 00391 /* strongly graded matrices, the code can get relatively accurate */ 00392 /* eigenvalues. */ 00393 atoli = uflow * 4.f + *pivmin * 4.f; 00394 if (irange == 3) { 00395 /* RANGE='I': Compute an interval containing eigenvalues */ 00396 /* IL through IU. The initial interval [GL,GU] from the global */ 00397 /* Gerschgorin bounds GL and GU is refined by SLAEBZ. */ 00398 itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.f)) 00399 + 2; 00400 work[*n + 1] = gl; 00401 work[*n + 2] = gl; 00402 work[*n + 3] = gu; 00403 work[*n + 4] = gu; 00404 work[*n + 5] = gl; 00405 work[*n + 6] = gu; 00406 iwork[1] = -1; 00407 iwork[2] = -1; 00408 iwork[3] = *n + 1; 00409 iwork[4] = *n + 1; 00410 iwork[5] = *il - 1; 00411 iwork[6] = *iu; 00412 00413 slaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, & 00414 d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5] 00415 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo); 00416 if (iinfo != 0) { 00417 *info = iinfo; 00418 return 0; 00419 } 00420 /* On exit, output intervals may not be ordered by ascending negcount */ 00421 if (iwork[6] == *iu) { 00422 *wl = work[*n + 1]; 00423 wlu = work[*n + 3]; 00424 nwl = iwork[1]; 00425 *wu = work[*n + 4]; 00426 wul = work[*n + 2]; 00427 nwu = iwork[4]; 00428 } else { 00429 *wl = work[*n + 2]; 00430 wlu = work[*n + 4]; 00431 nwl = iwork[2]; 00432 *wu = work[*n + 3]; 00433 wul = work[*n + 1]; 00434 nwu = iwork[3]; 00435 } 00436 /* On exit, the interval [WL, WLU] contains a value with negcount NWL, */ 00437 /* and [WUL, WU] contains a value with negcount NWU. */ 00438 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) { 00439 *info = 4; 00440 return 0; 00441 } 00442 } else if (irange == 2) { 00443 *wl = *vl; 00444 *wu = *vu; 00445 } else if (irange == 1) { 00446 *wl = gl; 00447 *wu = gu; 00448 } 00449 /* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */ 00450 /* NWL accumulates the number of eigenvalues .le. WL, */ 00451 /* NWU accumulates the number of eigenvalues .le. WU */ 00452 *m = 0; 00453 iend = 0; 00454 *info = 0; 00455 nwl = 0; 00456 nwu = 0; 00457 00458 i__1 = *nsplit; 00459 for (jblk = 1; jblk <= i__1; ++jblk) { 00460 ioff = iend; 00461 ibegin = ioff + 1; 00462 iend = isplit[jblk]; 00463 in = iend - ioff; 00464 00465 if (in == 1) { 00466 /* 1x1 block */ 00467 if (*wl >= d__[ibegin] - *pivmin) { 00468 ++nwl; 00469 } 00470 if (*wu >= d__[ibegin] - *pivmin) { 00471 ++nwu; 00472 } 00473 if (irange == 1 || *wl < d__[ibegin] - *pivmin && *wu >= d__[ 00474 ibegin] - *pivmin) { 00475 ++(*m); 00476 w[*m] = d__[ibegin]; 00477 werr[*m] = 0.f; 00478 /* The gap for a single block doesn't matter for the later */ 00479 /* algorithm and is assigned an arbitrary large value */ 00480 iblock[*m] = jblk; 00481 indexw[*m] = 1; 00482 } 00483 /* Disabled 2x2 case because of a failure on the following matrix */ 00484 /* RANGE = 'I', IL = IU = 4 */ 00485 /* Original Tridiagonal, d = [ */ 00486 /* -0.150102010615740E+00 */ 00487 /* -0.849897989384260E+00 */ 00488 /* -0.128208148052635E-15 */ 00489 /* 0.128257718286320E-15 */ 00490 /* ]; */ 00491 /* e = [ */ 00492 /* -0.357171383266986E+00 */ 00493 /* -0.180411241501588E-15 */ 00494 /* -0.175152352710251E-15 */ 00495 /* ]; */ 00496 00497 /* ELSE IF( IN.EQ.2 ) THEN */ 00498 /* * 2x2 block */ 00499 /* DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */ 00500 /* TMP1 = HALF*(D(IBEGIN)+D(IEND)) */ 00501 /* L1 = TMP1 - DISC */ 00502 /* IF( WL.GE. L1-PIVMIN ) */ 00503 /* $ NWL = NWL + 1 */ 00504 /* IF( WU.GE. L1-PIVMIN ) */ 00505 /* $ NWU = NWU + 1 */ 00506 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */ 00507 /* $ L1-PIVMIN ) ) THEN */ 00508 /* M = M + 1 */ 00509 /* W( M ) = L1 */ 00510 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */ 00511 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */ 00512 /* IBLOCK( M ) = JBLK */ 00513 /* INDEXW( M ) = 1 */ 00514 /* ENDIF */ 00515 /* L2 = TMP1 + DISC */ 00516 /* IF( WL.GE. L2-PIVMIN ) */ 00517 /* $ NWL = NWL + 1 */ 00518 /* IF( WU.GE. L2-PIVMIN ) */ 00519 /* $ NWU = NWU + 1 */ 00520 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */ 00521 /* $ L2-PIVMIN ) ) THEN */ 00522 /* M = M + 1 */ 00523 /* W( M ) = L2 */ 00524 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */ 00525 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */ 00526 /* IBLOCK( M ) = JBLK */ 00527 /* INDEXW( M ) = 2 */ 00528 /* ENDIF */ 00529 } else { 00530 /* General Case - block of size IN >= 2 */ 00531 /* Compute local Gerschgorin interval and use it as the initial */ 00532 /* interval for SLAEBZ */ 00533 gu = d__[ibegin]; 00534 gl = d__[ibegin]; 00535 tmp1 = 0.f; 00536 i__2 = iend; 00537 for (j = ibegin; j <= i__2; ++j) { 00538 /* Computing MIN */ 00539 r__1 = gl, r__2 = gers[(j << 1) - 1]; 00540 gl = dmin(r__1,r__2); 00541 /* Computing MAX */ 00542 r__1 = gu, r__2 = gers[j * 2]; 00543 gu = dmax(r__1,r__2); 00544 /* L40: */ 00545 } 00546 /* [JAN/28/2009] */ 00547 /* change SPDIAM by TNORM in lines 2 and 3 thereafter */ 00548 /* line 1: remove computation of SPDIAM (not useful anymore) */ 00549 /* SPDIAM = GU - GL */ 00550 /* GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN */ 00551 /* GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN */ 00552 gl = gl - tnorm * 2.f * eps * in - *pivmin * 2.f; 00553 gu = gu + tnorm * 2.f * eps * in + *pivmin * 2.f; 00554 00555 if (irange > 1) { 00556 if (gu < *wl) { 00557 /* the local block contains none of the wanted eigenvalues */ 00558 nwl += in; 00559 nwu += in; 00560 goto L70; 00561 } 00562 /* refine search interval if possible, only range (WL,WU] matters */ 00563 gl = dmax(gl,*wl); 00564 gu = dmin(gu,*wu); 00565 if (gl >= gu) { 00566 goto L70; 00567 } 00568 } 00569 /* Find negcount of initial interval boundaries GL and GU */ 00570 work[*n + 1] = gl; 00571 work[*n + in + 1] = gu; 00572 slaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, 00573 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, & 00574 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], & 00575 w[*m + 1], &iblock[*m + 1], &iinfo); 00576 if (iinfo != 0) { 00577 *info = iinfo; 00578 return 0; 00579 } 00580 00581 nwl += iwork[1]; 00582 nwu += iwork[in + 1]; 00583 iwoff = *m - iwork[1]; 00584 /* Compute Eigenvalues */ 00585 itmax = (integer) ((log(gu - gl + *pivmin) - log(*pivmin)) / log( 00586 2.f)) + 2; 00587 slaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, 00588 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, & 00589 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1], 00590 &w[*m + 1], &iblock[*m + 1], &iinfo); 00591 if (iinfo != 0) { 00592 *info = iinfo; 00593 return 0; 00594 } 00595 00596 /* Copy eigenvalues into W and IBLOCK */ 00597 /* Use -JBLK for block number for unconverged eigenvalues. */ 00598 /* Loop over the number of output intervals from SLAEBZ */ 00599 i__2 = iout; 00600 for (j = 1; j <= i__2; ++j) { 00601 /* eigenvalue approximation is middle point of interval */ 00602 tmp1 = (work[j + *n] + work[j + in + *n]) * .5f; 00603 /* semi length of error interval */ 00604 tmp2 = (r__1 = work[j + *n] - work[j + in + *n], dabs(r__1)) * 00605 .5f; 00606 if (j > iout - iinfo) { 00607 /* Flag non-convergence. */ 00608 ncnvrg = TRUE_; 00609 ib = -jblk; 00610 } else { 00611 ib = jblk; 00612 } 00613 i__3 = iwork[j + in] + iwoff; 00614 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) { 00615 w[je] = tmp1; 00616 werr[je] = tmp2; 00617 indexw[je] = je - iwoff; 00618 iblock[je] = ib; 00619 /* L50: */ 00620 } 00621 /* L60: */ 00622 } 00623 00624 *m += im; 00625 } 00626 L70: 00627 ; 00628 } 00629 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */ 00630 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */ 00631 if (irange == 3) { 00632 idiscl = *il - 1 - nwl; 00633 idiscu = nwu - *iu; 00634 00635 if (idiscl > 0) { 00636 im = 0; 00637 i__1 = *m; 00638 for (je = 1; je <= i__1; ++je) { 00639 /* Remove some of the smallest eigenvalues from the left so that */ 00640 /* at the end IDISCL =0. Move all eigenvalues up to the left. */ 00641 if (w[je] <= wlu && idiscl > 0) { 00642 --idiscl; 00643 } else { 00644 ++im; 00645 w[im] = w[je]; 00646 werr[im] = werr[je]; 00647 indexw[im] = indexw[je]; 00648 iblock[im] = iblock[je]; 00649 } 00650 /* L80: */ 00651 } 00652 *m = im; 00653 } 00654 if (idiscu > 0) { 00655 /* Remove some of the largest eigenvalues from the right so that */ 00656 /* at the end IDISCU =0. Move all eigenvalues up to the left. */ 00657 im = *m + 1; 00658 for (je = *m; je >= 1; --je) { 00659 if (w[je] >= wul && idiscu > 0) { 00660 --idiscu; 00661 } else { 00662 --im; 00663 w[im] = w[je]; 00664 werr[im] = werr[je]; 00665 indexw[im] = indexw[je]; 00666 iblock[im] = iblock[je]; 00667 } 00668 /* L81: */ 00669 } 00670 jee = 0; 00671 i__1 = *m; 00672 for (je = im; je <= i__1; ++je) { 00673 ++jee; 00674 w[jee] = w[je]; 00675 werr[jee] = werr[je]; 00676 indexw[jee] = indexw[je]; 00677 iblock[jee] = iblock[je]; 00678 /* L82: */ 00679 } 00680 *m = *m - im + 1; 00681 } 00682 if (idiscl > 0 || idiscu > 0) { 00683 /* Code to deal with effects of bad arithmetic. (If N(w) is */ 00684 /* monotone non-decreasing, this should never happen.) */ 00685 /* Some low eigenvalues to be discarded are not in (WL,WLU], */ 00686 /* or high eigenvalues to be discarded are not in (WUL,WU] */ 00687 /* so just kill off the smallest IDISCL/largest IDISCU */ 00688 /* eigenvalues, by marking the corresponding IBLOCK = 0 */ 00689 if (idiscl > 0) { 00690 wkill = *wu; 00691 i__1 = idiscl; 00692 for (jdisc = 1; jdisc <= i__1; ++jdisc) { 00693 iw = 0; 00694 i__2 = *m; 00695 for (je = 1; je <= i__2; ++je) { 00696 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) { 00697 iw = je; 00698 wkill = w[je]; 00699 } 00700 /* L90: */ 00701 } 00702 iblock[iw] = 0; 00703 /* L100: */ 00704 } 00705 } 00706 if (idiscu > 0) { 00707 wkill = *wl; 00708 i__1 = idiscu; 00709 for (jdisc = 1; jdisc <= i__1; ++jdisc) { 00710 iw = 0; 00711 i__2 = *m; 00712 for (je = 1; je <= i__2; ++je) { 00713 if (iblock[je] != 0 && (w[je] >= wkill || iw == 0)) { 00714 iw = je; 00715 wkill = w[je]; 00716 } 00717 /* L110: */ 00718 } 00719 iblock[iw] = 0; 00720 /* L120: */ 00721 } 00722 } 00723 /* Now erase all eigenvalues with IBLOCK set to zero */ 00724 im = 0; 00725 i__1 = *m; 00726 for (je = 1; je <= i__1; ++je) { 00727 if (iblock[je] != 0) { 00728 ++im; 00729 w[im] = w[je]; 00730 werr[im] = werr[je]; 00731 indexw[im] = indexw[je]; 00732 iblock[im] = iblock[je]; 00733 } 00734 /* L130: */ 00735 } 00736 *m = im; 00737 } 00738 if (idiscl < 0 || idiscu < 0) { 00739 toofew = TRUE_; 00740 } 00741 } 00742 00743 if (irange == 1 && *m != *n || irange == 3 && *m != *iu - *il + 1) { 00744 toofew = TRUE_; 00745 } 00746 /* If ORDER='B', do nothing the eigenvalues are already sorted by */ 00747 /* block. */ 00748 /* If ORDER='E', sort the eigenvalues from smallest to largest */ 00749 if (lsame_(order, "E") && *nsplit > 1) { 00750 i__1 = *m - 1; 00751 for (je = 1; je <= i__1; ++je) { 00752 ie = 0; 00753 tmp1 = w[je]; 00754 i__2 = *m; 00755 for (j = je + 1; j <= i__2; ++j) { 00756 if (w[j] < tmp1) { 00757 ie = j; 00758 tmp1 = w[j]; 00759 } 00760 /* L140: */ 00761 } 00762 if (ie != 0) { 00763 tmp2 = werr[ie]; 00764 itmp1 = iblock[ie]; 00765 itmp2 = indexw[ie]; 00766 w[ie] = w[je]; 00767 werr[ie] = werr[je]; 00768 iblock[ie] = iblock[je]; 00769 indexw[ie] = indexw[je]; 00770 w[je] = tmp1; 00771 werr[je] = tmp2; 00772 iblock[je] = itmp1; 00773 indexw[je] = itmp2; 00774 } 00775 /* L150: */ 00776 } 00777 } 00778 00779 *info = 0; 00780 if (ncnvrg) { 00781 ++(*info); 00782 } 00783 if (toofew) { 00784 *info += 2; 00785 } 00786 return 0; 00787 00788 /* End of SLARRD */ 00789 00790 } /* slarrd_ */