1 /* slarrj.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
10 http://www.netlib.org/f2c/libf2c.zip
16 /* Subroutine */ int slarrj_(integer *n, real *d__, real *e2, integer *ifirst,
17 integer *ilast, real *rtol, integer *offset, real *w, real *werr,
18 real *work, integer *iwork, real *pivmin, real *spdiam, integer *info)
20 /* System generated locals */
24 /* Builtin functions */
25 double log(doublereal);
34 integer iter, nint, prev, next, savi1;
35 real right, width, dplus;
36 integer olnint, maxitr;
39 /* -- LAPACK auxiliary routine (version 3.2) -- */
40 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
43 /* .. Scalar Arguments .. */
45 /* .. Array Arguments .. */
51 /* Given the initial eigenvalue approximations of T, SLARRJ */
52 /* does bisection to refine the eigenvalues of T, */
53 /* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial */
54 /* guesses for these eigenvalues are input in W, the corresponding estimate */
55 /* of the error in these guesses in WERR. During bisection, intervals */
56 /* [left, right] are maintained by storing their mid-points and */
57 /* semi-widths in the arrays W and WERR respectively. */
62 /* N (input) INTEGER */
63 /* The order of the matrix. */
65 /* D (input) REAL array, dimension (N) */
66 /* The N diagonal elements of T. */
68 /* E2 (input) REAL array, dimension (N-1) */
69 /* The Squares of the (N-1) subdiagonal elements of T. */
71 /* IFIRST (input) INTEGER */
72 /* The index of the first eigenvalue to be computed. */
74 /* ILAST (input) INTEGER */
75 /* The index of the last eigenvalue to be computed. */
77 /* RTOL (input) REAL */
78 /* Tolerance for the convergence of the bisection intervals. */
79 /* An interval [LEFT,RIGHT] has converged if */
80 /* RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|). */
82 /* OFFSET (input) INTEGER */
83 /* Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET */
84 /* through ILAST-OFFSET elements of these arrays are to be used. */
86 /* W (input/output) REAL array, dimension (N) */
87 /* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are */
88 /* estimates of the eigenvalues of L D L^T indexed IFIRST through */
90 /* On output, these estimates are refined. */
92 /* WERR (input/output) REAL array, dimension (N) */
93 /* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are */
94 /* the errors in the estimates of the corresponding elements in W. */
95 /* On output, these errors are refined. */
97 /* WORK (workspace) REAL array, dimension (2*N) */
100 /* IWORK (workspace) INTEGER array, dimension (2*N) */
103 /* PIVMIN (input) DOUBLE PRECISION */
104 /* The minimum pivot in the Sturm sequence for T. */
106 /* SPDIAM (input) DOUBLE PRECISION */
107 /* The spectral diameter of T. */
109 /* INFO (output) INTEGER */
112 /* Further Details */
113 /* =============== */
115 /* Based on contributions by */
116 /* Beresford Parlett, University of California, Berkeley, USA */
117 /* Jim Demmel, University of California, Berkeley, USA */
118 /* Inderjit Dhillon, University of Texas, Austin, USA */
119 /* Osni Marques, LBNL/NERSC, USA */
120 /* Christof Voemel, University of California, Berkeley, USA */
122 /* ===================================================================== */
124 /* .. Parameters .. */
126 /* .. Local Scalars .. */
129 /* .. Intrinsic Functions .. */
131 /* .. Executable Statements .. */
133 /* Parameter adjustments */
144 maxitr = (integer) ((log(*spdiam + *pivmin) - log(*pivmin)) / log(2.f)) +
147 /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. */
148 /* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while */
149 /* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) */
150 /* for an unconverged interval is set to the index of the next unconverged */
151 /* interval, and is -1 or 0 for a converged interval. Thus a linked */
152 /* list of unconverged intervals is set up. */
156 /* The number of unconverged intervals */
158 /* The last unconverged interval found */
161 for (i__ = i1; i__ <= i__1; ++i__) {
164 left = w[ii] - werr[ii];
166 right = w[ii] + werr[ii];
169 r__1 = dabs(left), r__2 = dabs(right);
170 tmp = dmax(r__1,r__2);
171 /* The following test prevents the test of converged intervals */
172 if (width < *rtol * tmp) {
173 /* This interval has already converged and does not need refinement. */
174 /* (Note that the gaps might change through refining the */
175 /* eigenvalues, however, they can only get bigger.) */
176 /* Remove it from the list. */
178 /* Make sure that I1 always points to the first unconverged interval */
179 if (i__ == i1 && i__ < i2) {
182 if (prev >= i1 && i__ <= i2) {
183 iwork[(prev << 1) - 1] = i__ + 1;
186 /* unconverged interval found */
188 /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue */
190 /* Do while( CNT(LEFT).GT.I-1 ) */
201 for (j = 2; j <= i__2; ++j) {
202 dplus = d__[j] - s - e2[j - 1] / dplus;
209 left -= werr[ii] * fac;
214 /* Do while( CNT(RIGHT).LT.I ) */
225 for (j = 2; j <= i__2; ++j) {
226 dplus = d__[j] - s - e2[j - 1] / dplus;
233 right += werr[ii] * fac;
238 iwork[k - 1] = i__ + 1;
247 /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals */
248 /* and while (ITER.LT.MAXITR) */
256 for (p = 1; p <= i__1; ++p) {
262 mid = (left + right) * .5f;
263 /* semiwidth of interval */
266 r__1 = dabs(left), r__2 = dabs(right);
267 tmp = dmax(r__1,r__2);
268 if (width < *rtol * tmp || iter == maxitr) {
269 /* reduce number of unconverged intervals */
271 /* Mark interval as converged. */
276 /* Prev holds the last unconverged interval previously examined */
278 iwork[(prev << 1) - 1] = next;
286 /* Perform one bisection step */
295 for (j = 2; j <= i__2; ++j) {
296 dplus = d__[j] - s - e2[j - 1] / dplus;
302 if (cnt <= i__ - 1) {
312 /* do another loop if there are still unconverged intervals */
313 /* However, in the last iteration, all intervals are accepted */
314 /* since this is the best we can do. */
315 if (nint > 0 && iter <= maxitr) {
320 /* At this point, all the intervals have converged */
322 for (i__ = savi1; i__ <= i__1; ++i__) {
325 /* All intervals marked by '0' have been refined. */
326 if (iwork[k - 1] == 0) {
327 w[ii] = (work[k - 1] + work[k]) * .5f;
328 werr[ii] = work[k] - w[ii];