3 /* Subroutine */ int dlarrk_(integer *n, integer *iw, doublereal *gl,
4 doublereal *gu, doublereal *d__, doublereal *e2, doublereal *pivmin,
5 doublereal *reltol, doublereal *w, doublereal *werr, integer *info)
7 /* System generated locals */
11 /* Builtin functions */
12 double log(doublereal);
16 doublereal mid, eps, tmp1, tmp2, left, atoli, right;
18 doublereal rtoli, tnorm;
19 extern doublereal dlamch_(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 /* DLARRK computes one eigenvalue of a symmetric tridiagonal */
36 /* matrix T to suitable accuracy. This is an auxiliary code to be */
37 /* called from DSTEMR. */
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. */
57 /* GL (input) DOUBLE PRECISION */
58 /* GU (input) DOUBLE PRECISION */
59 /* An upper and a lower bound on the eigenvalue. */
61 /* D (input) DOUBLE PRECISION array, dimension (N) */
62 /* The n diagonal elements of the tridiagonal matrix T. */
64 /* E2 (input) DOUBLE PRECISION array, dimension (N-1) */
65 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
67 /* PIVMIN (input) DOUBLE PRECISION */
68 /* The minimum pivot allowed in the Sturm sequence for T. */
70 /* RELTOL (input) DOUBLE PRECISION */
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. */
77 /* W (output) DOUBLE PRECISION */
79 /* WERR (output) DOUBLE PRECISION */
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 DOUBLE PRECISION, 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 d__1 = abs(*gl), d__2 = abs(*gu);
114 tnorm = max(d__1,d__2);
116 atoli = *pivmin * 4.;
117 itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.)) + 2;
119 left = *gl - tnorm * 2. * eps * *n - *pivmin * 4.;
120 right = *gu + tnorm * 2. * eps * *n + *pivmin * 4.;
124 /* Check if interval converged or maximum number of iterations reached */
126 tmp1 = (d__1 = right - left, abs(d__1));
128 d__1 = abs(right), d__2 = abs(left);
129 tmp2 = max(d__1,d__2);
131 d__1 = max(atoli,*pivmin), d__2 = rtoli * tmp2;
132 if (tmp1 < max(d__1,d__2)) {
140 /* Count number of negative pivots for mid-point */
143 mid = (left + right) * .5;
146 if (abs(tmp1) < *pivmin) {
154 for (i__ = 2; i__ <= i__1; ++i__) {
155 tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid;
156 if (abs(tmp1) < *pivmin) {
172 /* Converged or maximum number of iterations reached */
174 *w = (left + right) * .5;
175 *werr = (d__1 = right - left, abs(d__1)) * .5;