--- /dev/null
+#include "clapack.h"
+
+/* Table of constant values */
+
+static integer c__1 = 1;
+
+/* Subroutine */ int dlarrf_(integer *n, doublereal *d__, doublereal *l,
+ doublereal *ld, integer *clstrt, integer *clend, doublereal *w,
+ doublereal *wgap, doublereal *werr, doublereal *spdiam, doublereal *
+ clgapl, doublereal *clgapr, doublereal *pivmin, doublereal *sigma,
+ doublereal *dplus, doublereal *lplus, doublereal *work, integer *info)
+{
+ /* System generated locals */
+ integer i__1;
+ doublereal d__1, d__2, d__3;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ integer i__;
+ doublereal s, bestshift, smlgrowth, eps, tmp, max1, max2, rrr1, rrr2,
+ znm2, growthbound, fail, fact, oldp;
+ integer indx;
+ doublereal prod;
+ integer ktry;
+ doublereal fail2, avgap, ldmax, rdmax;
+ integer shift;
+ extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
+ doublereal *, integer *);
+ logical dorrr1;
+ extern doublereal dlamch_(char *);
+ doublereal ldelta;
+ logical nofail;
+ doublereal mingap, lsigma, rdelta;
+ extern logical disnan_(doublereal *);
+ logical forcer;
+ doublereal rsigma, clwdth;
+ logical sawnan1, sawnan2, tryrrr1;
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+/* * */
+/* .. Scalar Arguments .. */
+/* .. */
+/* .. Array Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* Given the initial representation L D L^T and its cluster of close */
+/* eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... */
+/* W( CLEND ), DLARRF finds a new relatively robust representation */
+/* L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the */
+/* eigenvalues of L(+) D(+) L(+)^T is relatively isolated. */
+
+/* Arguments */
+/* ========= */
+
+/* N (input) INTEGER */
+/* The order of the matrix (subblock, if the matrix splitted). */
+
+/* D (input) DOUBLE PRECISION array, dimension (N) */
+/* The N diagonal elements of the diagonal matrix D. */
+
+/* L (input) DOUBLE PRECISION array, dimension (N-1) */
+/* The (N-1) subdiagonal elements of the unit bidiagonal */
+/* matrix L. */
+
+/* LD (input) DOUBLE PRECISION array, dimension (N-1) */
+/* The (N-1) elements L(i)*D(i). */
+
+/* CLSTRT (input) INTEGER */
+/* The index of the first eigenvalue in the cluster. */
+
+/* CLEND (input) INTEGER */
+/* The index of the last eigenvalue in the cluster. */
+
+/* W (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1) */
+/* The eigenvalue APPROXIMATIONS of L D L^T in ascending order. */
+/* W( CLSTRT ) through W( CLEND ) form the cluster of relatively */
+/* close eigenalues. */
+
+/* WGAP (input/output) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1) */
+/* The separation from the right neighbor eigenvalue in W. */
+
+/* WERR (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1) */
+/* WERR contain the semiwidth of the uncertainty */
+/* interval of the corresponding eigenvalue APPROXIMATION in W */
+
+/* SPDIAM (input) estimate of the spectral diameter obtained from the */
+/* Gerschgorin intervals */
+
+/* CLGAPL, CLGAPR (input) absolute gap on each end of the cluster. */
+/* Set by the calling routine to protect against shifts too close */
+/* to eigenvalues outside the cluster. */
+
+/* PIVMIN (input) DOUBLE PRECISION */
+/* The minimum pivot allowed in the Sturm sequence. */
+
+/* SIGMA (output) DOUBLE PRECISION */
+/* The shift used to form L(+) D(+) L(+)^T. */
+
+/* DPLUS (output) DOUBLE PRECISION array, dimension (N) */
+/* The N diagonal elements of the diagonal matrix D(+). */
+
+/* LPLUS (output) DOUBLE PRECISION array, dimension (N-1) */
+/* The first (N-1) elements of LPLUS contain the subdiagonal */
+/* elements of the unit bidiagonal matrix L(+). */
+
+/* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */
+/* Workspace. */
+
+/* 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 */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. External Subroutines .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+ /* Parameter adjustments */
+ --work;
+ --lplus;
+ --dplus;
+ --werr;
+ --wgap;
+ --w;
+ --ld;
+ --l;
+ --d__;
+
+ /* Function Body */
+ *info = 0;
+ fact = 2.;
+ eps = dlamch_("Precision");
+ shift = 0;
+ forcer = FALSE_;
+/* Note that we cannot guarantee that for any of the shifts tried, */
+/* the factorization has a small or even moderate element growth. */
+/* There could be Ritz values at both ends of the cluster and despite */
+/* backing off, there are examples where all factorizations tried */
+/* (in IEEE mode, allowing zero pivots & infinities) have INFINITE */
+/* element growth. */
+/* For this reason, we should use PIVMIN in this subroutine so that at */
+/* least the L D L^T factorization exists. It can be checked afterwards */
+/* whether the element growth caused bad residuals/orthogonality. */
+/* Decide whether the code should accept the best among all */
+/* representations despite large element growth or signal INFO=1 */
+ nofail = TRUE_;
+
+/* Compute the average gap length of the cluster */
+ clwdth = (d__1 = w[*clend] - w[*clstrt], abs(d__1)) + werr[*clend] + werr[
+ *clstrt];
+ avgap = clwdth / (doublereal) (*clend - *clstrt);
+ mingap = min(*clgapl,*clgapr);
+/* Initial values for shifts to both ends of cluster */
+/* Computing MIN */
+ d__1 = w[*clstrt], d__2 = w[*clend];
+ lsigma = min(d__1,d__2) - werr[*clstrt];
+/* Computing MAX */
+ d__1 = w[*clstrt], d__2 = w[*clend];
+ rsigma = max(d__1,d__2) + werr[*clend];
+/* Use a small fudge to make sure that we really shift to the outside */
+ lsigma -= abs(lsigma) * 4. * eps;
+ rsigma += abs(rsigma) * 4. * eps;
+/* Compute upper bounds for how much to back off the initial shifts */
+ ldmax = mingap * .25 + *pivmin * 2.;
+ rdmax = mingap * .25 + *pivmin * 2.;
+/* Computing MAX */
+ d__1 = avgap, d__2 = wgap[*clstrt];
+ ldelta = max(d__1,d__2) / fact;
+/* Computing MAX */
+ d__1 = avgap, d__2 = wgap[*clend - 1];
+ rdelta = max(d__1,d__2) / fact;
+
+/* Initialize the record of the best representation found */
+
+ s = dlamch_("S");
+ smlgrowth = 1. / s;
+ fail = (doublereal) (*n - 1) * mingap / (*spdiam * eps);
+ fail2 = (doublereal) (*n - 1) * mingap / (*spdiam * sqrt(eps));
+ bestshift = lsigma;
+
+/* while (KTRY <= KTRYMAX) */
+ ktry = 0;
+ growthbound = *spdiam * 8.;
+L5:
+ sawnan1 = FALSE_;
+ sawnan2 = FALSE_;
+/* Ensure that we do not back off too much of the initial shifts */
+ ldelta = min(ldmax,ldelta);
+ rdelta = min(rdmax,rdelta);
+/* Compute the element growth when shifting to both ends of the cluster */
+/* accept the shift if there is no element growth at one of the two ends */
+/* Left end */
+ s = -lsigma;
+ dplus[1] = d__[1] + s;
+ if (abs(dplus[1]) < *pivmin) {
+ dplus[1] = -(*pivmin);
+/* Need to set SAWNAN1 because refined RRR test should not be used */
+/* in this case */
+ sawnan1 = TRUE_;
+ }
+ max1 = abs(dplus[1]);
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ lplus[i__] = ld[i__] / dplus[i__];
+ s = s * lplus[i__] * l[i__] - lsigma;
+ dplus[i__ + 1] = d__[i__ + 1] + s;
+ if ((d__1 = dplus[i__ + 1], abs(d__1)) < *pivmin) {
+ dplus[i__ + 1] = -(*pivmin);
+/* Need to set SAWNAN1 because refined RRR test should not be used */
+/* in this case */
+ sawnan1 = TRUE_;
+ }
+/* Computing MAX */
+ d__2 = max1, d__3 = (d__1 = dplus[i__ + 1], abs(d__1));
+ max1 = max(d__2,d__3);
+/* L6: */
+ }
+ sawnan1 = sawnan1 || disnan_(&max1);
+ if (forcer || max1 <= growthbound && ! sawnan1) {
+ *sigma = lsigma;
+ shift = 1;
+ goto L100;
+ }
+/* Right end */
+ s = -rsigma;
+ work[1] = d__[1] + s;
+ if (abs(work[1]) < *pivmin) {
+ work[1] = -(*pivmin);
+/* Need to set SAWNAN2 because refined RRR test should not be used */
+/* in this case */
+ sawnan2 = TRUE_;
+ }
+ max2 = abs(work[1]);
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ work[*n + i__] = ld[i__] / work[i__];
+ s = s * work[*n + i__] * l[i__] - rsigma;
+ work[i__ + 1] = d__[i__ + 1] + s;
+ if ((d__1 = work[i__ + 1], abs(d__1)) < *pivmin) {
+ work[i__ + 1] = -(*pivmin);
+/* Need to set SAWNAN2 because refined RRR test should not be used */
+/* in this case */
+ sawnan2 = TRUE_;
+ }
+/* Computing MAX */
+ d__2 = max2, d__3 = (d__1 = work[i__ + 1], abs(d__1));
+ max2 = max(d__2,d__3);
+/* L7: */
+ }
+ sawnan2 = sawnan2 || disnan_(&max2);
+ if (forcer || max2 <= growthbound && ! sawnan2) {
+ *sigma = rsigma;
+ shift = 2;
+ goto L100;
+ }
+/* If we are at this point, both shifts led to too much element growth */
+/* Record the better of the two shifts (provided it didn't lead to NaN) */
+ if (sawnan1 && sawnan2) {
+/* both MAX1 and MAX2 are NaN */
+ goto L50;
+ } else {
+ if (! sawnan1) {
+ indx = 1;
+ if (max1 <= smlgrowth) {
+ smlgrowth = max1;
+ bestshift = lsigma;
+ }
+ }
+ if (! sawnan2) {
+ if (sawnan1 || max2 <= max1) {
+ indx = 2;
+ }
+ if (max2 <= smlgrowth) {
+ smlgrowth = max2;
+ bestshift = rsigma;
+ }
+ }
+ }
+/* If we are here, both the left and the right shift led to */
+/* element growth. If the element growth is moderate, then */
+/* we may still accept the representation, if it passes a */
+/* refined test for RRR. This test supposes that no NaN occurred. */
+/* Moreover, we use the refined RRR test only for isolated clusters. */
+ if (clwdth < mingap / 128. && min(max1,max2) < fail2 && ! sawnan1 && !
+ sawnan2) {
+ dorrr1 = TRUE_;
+ } else {
+ dorrr1 = FALSE_;
+ }
+ tryrrr1 = TRUE_;
+ if (tryrrr1 && dorrr1) {
+ if (indx == 1) {
+ tmp = (d__1 = dplus[*n], abs(d__1));
+ znm2 = 1.;
+ prod = 1.;
+ oldp = 1.;
+ for (i__ = *n - 1; i__ >= 1; --i__) {
+ if (prod <= eps) {
+ prod = dplus[i__ + 1] * work[*n + i__ + 1] / (dplus[i__] *
+ work[*n + i__]) * oldp;
+ } else {
+ prod *= (d__1 = work[*n + i__], abs(d__1));
+ }
+ oldp = prod;
+/* Computing 2nd power */
+ d__1 = prod;
+ znm2 += d__1 * d__1;
+/* Computing MAX */
+ d__2 = tmp, d__3 = (d__1 = dplus[i__] * prod, abs(d__1));
+ tmp = max(d__2,d__3);
+/* L15: */
+ }
+ rrr1 = tmp / (*spdiam * sqrt(znm2));
+ if (rrr1 <= 8.) {
+ *sigma = lsigma;
+ shift = 1;
+ goto L100;
+ }
+ } else if (indx == 2) {
+ tmp = (d__1 = work[*n], abs(d__1));
+ znm2 = 1.;
+ prod = 1.;
+ oldp = 1.;
+ for (i__ = *n - 1; i__ >= 1; --i__) {
+ if (prod <= eps) {
+ prod = work[i__ + 1] * lplus[i__ + 1] / (work[i__] *
+ lplus[i__]) * oldp;
+ } else {
+ prod *= (d__1 = lplus[i__], abs(d__1));
+ }
+ oldp = prod;
+/* Computing 2nd power */
+ d__1 = prod;
+ znm2 += d__1 * d__1;
+/* Computing MAX */
+ d__2 = tmp, d__3 = (d__1 = work[i__] * prod, abs(d__1));
+ tmp = max(d__2,d__3);
+/* L16: */
+ }
+ rrr2 = tmp / (*spdiam * sqrt(znm2));
+ if (rrr2 <= 8.) {
+ *sigma = rsigma;
+ shift = 2;
+ goto L100;
+ }
+ }
+ }
+L50:
+ if (ktry < 1) {
+/* If we are here, both shifts failed also the RRR test. */
+/* Back off to the outside */
+/* Computing MAX */
+ d__1 = lsigma - ldelta, d__2 = lsigma - ldmax;
+ lsigma = max(d__1,d__2);
+/* Computing MIN */
+ d__1 = rsigma + rdelta, d__2 = rsigma + rdmax;
+ rsigma = min(d__1,d__2);
+ ldelta *= 2.;
+ rdelta *= 2.;
+ ++ktry;
+ goto L5;
+ } else {
+/* None of the representations investigated satisfied our */
+/* criteria. Take the best one we found. */
+ if (smlgrowth < fail || nofail) {
+ lsigma = bestshift;
+ rsigma = bestshift;
+ forcer = TRUE_;
+ goto L5;
+ } else {
+ *info = 1;
+ return 0;
+ }
+ }
+L100:
+ if (shift == 1) {
+ } else if (shift == 2) {
+/* store new L and D back into DPLUS, LPLUS */
+ dcopy_(n, &work[1], &c__1, &dplus[1], &c__1);
+ i__1 = *n - 1;
+ dcopy_(&i__1, &work[*n + 1], &c__1, &lplus[1], &c__1);
+ }
+ return 0;
+
+/* End of DLARRF */
+
+} /* dlarrf_ */