#include "blaswrap.h" #include "f2c.h" /* Subroutine */ int slarrj_(integer *n, real *d__, real *e2, integer *ifirst, integer *ilast, real *rtol, integer *offset, real *w, real *werr, real *work, integer *iwork, real *pivmin, real *spdiam, integer *info) { /* -- LAPACK auxiliary routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= Given the initial eigenvalue approximations of T, SLARRJ does bisection to refine the eigenvalues of T, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial guesses for these eigenvalues are input in W, the corresponding estimate of the error in these guesses in WERR. During bisection, intervals [left, right] are maintained by storing their mid-points and semi-widths in the arrays W and WERR respectively. Arguments ========= N (input) INTEGER The order of the matrix. D (input) REAL array, dimension (N) The N diagonal elements of T. E2 (input) REAL array, dimension (N-1) The Squares of the (N-1) subdiagonal elements of T. IFIRST (input) INTEGER The index of the first eigenvalue to be computed. ILAST (input) INTEGER The index of the last eigenvalue to be computed. RTOL (input) REAL Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|). OFFSET (input) INTEGER Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are to be used. W (input/output) REAL array, dimension (N) On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. On output, these estimates are refined. WERR (input/output) REAL array, dimension (N) On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined. WORK (workspace) REAL array, dimension (2*N) Workspace. IWORK (workspace) INTEGER array, dimension (2*N) Workspace. PIVMIN (input) DOUBLE PRECISION The minimum pivot in the Sturm sequence for T. SPDIAM (input) DOUBLE PRECISION The spectral diameter of T. INFO (output) INTEGER Error flag. Further Details =============== Based on contributions by Beresford Parlett, University of California, Berkeley, USA Jim Demmel, University of California, Berkeley, USA Inderjit Dhillon, University of Texas, Austin, USA Osni Marques, LBNL/NERSC, USA Christof Voemel, University of California, Berkeley, USA ===================================================================== Parameter adjustments */ /* System generated locals */ integer i__1, i__2; real r__1, r__2; /* Builtin functions */ double log(doublereal); /* Local variables */ static integer i__, j, k, p; static real s; static integer i1, i2, ii; static real fac, mid; static integer cnt; static real tmp, left; static integer iter, nint, prev, next, savi1; static real right, width, dplus; static integer olnint, maxitr; --iwork; --work; --werr; --w; --e2; --d__; /* Function Body */ *info = 0; maxitr = (integer) ((log(*spdiam + *pivmin) - log(*pivmin)) / log(2.f)) + 2; /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) for an unconverged interval is set to the index of the next unconverged interval, and is -1 or 0 for a converged interval. Thus a linked list of unconverged intervals is set up. */ i1 = *ifirst; i2 = *ilast; /* The number of unconverged intervals */ nint = 0; /* The last unconverged interval found */ prev = 0; i__1 = i2; for (i__ = i1; i__ <= i__1; ++i__) { k = i__ << 1; ii = i__ - *offset; left = w[ii] - werr[ii]; mid = w[ii]; right = w[ii] + werr[ii]; width = right - mid; /* Computing MAX */ r__1 = dabs(left), r__2 = dabs(right); tmp = dmax(r__1,r__2); /* The following test prevents the test of converged intervals */ if (width < *rtol * tmp) { /* This interval has already converged and does not need refinement. (Note that the gaps might change through refining the eigenvalues, however, they can only get bigger.) Remove it from the list. */ iwork[k - 1] = -1; /* Make sure that I1 always points to the first unconverged interval */ if (i__ == i1 && i__ < i2) { i1 = i__ + 1; } if (prev >= i1 && i__ <= i2) { iwork[(prev << 1) - 1] = i__ + 1; } } else { /* unconverged interval found */ prev = i__; /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue Do while( CNT(LEFT).GT.I-1 ) */ fac = 1.f; L20: cnt = 0; s = left; dplus = d__[1] - s; if (dplus < 0.f) { ++cnt; } i__2 = *n; for (j = 2; j <= i__2; ++j) { dplus = d__[j] - s - e2[j - 1] / dplus; if (dplus < 0.f) { ++cnt; } /* L30: */ } if (cnt > i__ - 1) { left -= werr[ii] * fac; fac *= 2.f; goto L20; } /* Do while( CNT(RIGHT).LT.I ) */ fac = 1.f; L50: cnt = 0; s = right; dplus = d__[1] - s; if (dplus < 0.f) { ++cnt; } i__2 = *n; for (j = 2; j <= i__2; ++j) { dplus = d__[j] - s - e2[j - 1] / dplus; if (dplus < 0.f) { ++cnt; } /* L60: */ } if (cnt < i__) { right += werr[ii] * fac; fac *= 2.f; goto L50; } ++nint; iwork[k - 1] = i__ + 1; iwork[k] = cnt; } work[k - 1] = left; work[k] = right; /* L75: */ } savi1 = i1; /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals and while (ITER.LT.MAXITR) */ iter = 0; L80: prev = i1 - 1; i__ = i1; olnint = nint; i__1 = olnint; for (p = 1; p <= i__1; ++p) { k = i__ << 1; ii = i__ - *offset; next = iwork[k - 1]; left = work[k - 1]; right = work[k]; mid = (left + right) * .5f; /* semiwidth of interval */ width = right - mid; /* Computing MAX */ r__1 = dabs(left), r__2 = dabs(right); tmp = dmax(r__1,r__2); if (width < *rtol * tmp || iter == maxitr) { /* reduce number of unconverged intervals */ --nint; /* Mark interval as converged. */ iwork[k - 1] = 0; if (i1 == i__) { i1 = next; } else { /* Prev holds the last unconverged interval previously examined */ if (prev >= i1) { iwork[(prev << 1) - 1] = next; } } i__ = next; goto L100; } prev = i__; /* Perform one bisection step */ cnt = 0; s = mid; dplus = d__[1] - s; if (dplus < 0.f) { ++cnt; } i__2 = *n; for (j = 2; j <= i__2; ++j) { dplus = d__[j] - s - e2[j - 1] / dplus; if (dplus < 0.f) { ++cnt; } /* L90: */ } if (cnt <= i__ - 1) { work[k - 1] = mid; } else { work[k] = mid; } i__ = next; L100: ; } ++iter; /* do another loop if there are still unconverged intervals However, in the last iteration, all intervals are accepted since this is the best we can do. */ if (nint > 0 && iter <= maxitr) { goto L80; } /* At this point, all the intervals have converged */ i__1 = *ilast; for (i__ = savi1; i__ <= i__1; ++i__) { k = i__ << 1; ii = i__ - *offset; /* All intervals marked by '0' have been refined. */ if (iwork[k - 1] == 0) { w[ii] = (work[k - 1] + work[k]) * .5f; werr[ii] = work[k] - w[ii]; } /* L110: */ } return 0; /* End of SLARRJ */ } /* slarrj_ */