3 /* Subroutine */ int slarrk_(integer *n, integer *iw, real *gl, real *gu,
4 real *d__, real *e2, real *pivmin, real *reltol, real *w, real *werr,
7 /* System generated locals */
11 /* Builtin functions */
12 double log(doublereal);
16 real mid, eps, tmp1, tmp2, left, atoli, right;
19 extern doublereal slamch_(char *);
23 /* -- LAPACK auxiliary routine (version 3.1) -- */
24 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
27 /* .. Scalar Arguments .. */
29 /* .. Array Arguments .. */
35 /* SLARRK computes one eigenvalue of a symmetric tridiagonal */
36 /* matrix T to suitable accuracy. This is an auxiliary code to be */
37 /* called from SSTEMR. */
39 /* To avoid overflow, the matrix must be scaled so that its */
40 /* largest element is no greater than overflow**(1/2) * */
41 /* underflow**(1/4) in absolute value, and for greatest */
42 /* accuracy, it should not be much smaller than that. */
44 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
45 /* Matrix", Report CS41, Computer Science Dept., Stanford */
46 /* University, July 21, 1966. */
51 /* N (input) INTEGER */
52 /* The order of the tridiagonal matrix T. N >= 0. */
54 /* IW (input) INTEGER */
55 /* The index of the eigenvalues to be returned. */
59 /* An upper and a lower bound on the eigenvalue. */
61 /* D (input) REAL array, dimension (N) */
62 /* The n diagonal elements of the tridiagonal matrix T. */
64 /* E2 (input) REAL array, dimension (N-1) */
65 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
67 /* PIVMIN (input) REAL */
68 /* The minimum pivot allowed in the Sturm sequence for T. */
70 /* RELTOL (input) REAL */
71 /* The minimum relative width of an interval. When an interval */
72 /* is narrower than RELTOL times the larger (in */
73 /* magnitude) endpoint, then it is considered to be */
74 /* sufficiently small, i.e., converged. Note: this should */
75 /* always be at least radix*machine epsilon. */
79 /* WERR (output) REAL */
80 /* The error bound on the corresponding eigenvalue approximation */
83 /* INFO (output) INTEGER */
84 /* = 0: Eigenvalue converged */
85 /* = -1: Eigenvalue did NOT converge */
87 /* Internal Parameters */
88 /* =================== */
90 /* FUDGE REAL , default = 2 */
91 /* A "fudge factor" to widen the Gershgorin intervals. */
93 /* ===================================================================== */
95 /* .. Parameters .. */
97 /* .. Local Scalars .. */
99 /* .. External Functions .. */
101 /* .. Intrinsic Functions .. */
103 /* .. Executable Statements .. */
105 /* Get machine constants */
106 /* Parameter adjustments */
113 r__1 = dabs(*gl), r__2 = dabs(*gu);
114 tnorm = dmax(r__1,r__2);
116 atoli = *pivmin * 4.f;
117 itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.f)) + 2;
119 left = *gl - tnorm * 2.f * eps * *n - *pivmin * 4.f;
120 right = *gu + tnorm * 2.f * eps * *n + *pivmin * 4.f;
124 /* Check if interval converged or maximum number of iterations reached */
126 tmp1 = (r__1 = right - left, dabs(r__1));
128 r__1 = dabs(right), r__2 = dabs(left);
129 tmp2 = dmax(r__1,r__2);
131 r__1 = max(atoli,*pivmin), r__2 = rtoli * tmp2;
132 if (tmp1 < dmax(r__1,r__2)) {
140 /* Count number of negative pivots for mid-point */
143 mid = (left + right) * .5f;
146 if (dabs(tmp1) < *pivmin) {
154 for (i__ = 2; i__ <= i__1; ++i__) {
155 tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid;
156 if (dabs(tmp1) < *pivmin) {
172 /* Converged or maximum number of iterations reached */
174 *w = (left + right) * .5f;
175 *werr = (r__1 = right - left, dabs(r__1)) * .5f;