[] add metering mode to CLI

1 files changed, 1010 insertions, 0 deletions

diff --git a/contrib/libs/clapack/slasd4.c b/contrib/libs/clapack/slasd4.c
new file mode 100644
index 0000000000..11ba6d23c6
--- /dev/null
+++ b/contrib/libs/clapack/slasd4.c
@@ -0,0 +1,1010 @@
+/* slasd4.f -- translated by f2c (version 20061008).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "f2c.h"
+#include "blaswrap.h"
+
+/* Subroutine */ int slasd4_(integer *n, integer *i__, real *d__, real *z__, 
+	real *delta, real *rho, real *sigma, real *work, integer *info)
+{
+    /* System generated locals */
+    integer i__1;
+    real r__1;
+
+    /* Builtin functions */
+    double sqrt(doublereal);
+
+    /* Local variables */
+    real a, b, c__;
+    integer j;
+    real w, dd[3];
+    integer ii;
+    real dw, zz[3];
+    integer ip1;
+    real eta, phi, eps, tau, psi;
+    integer iim1, iip1;
+    real dphi, dpsi;
+    integer iter;
+    real temp, prew, sg2lb, sg2ub, temp1, temp2, dtiim, delsq, dtiip;
+    integer niter;
+    real dtisq;
+    logical swtch;
+    real dtnsq;
+    extern /* Subroutine */ int slaed6_(integer *, logical *, real *, real *, 
+	    real *, real *, real *, integer *);
+    real delsq2;
+    extern /* Subroutine */ int slasd5_(integer *, real *, real *, real *, 
+	    real *, real *, real *);
+    real dtnsq1;
+    logical swtch3;
+    extern doublereal slamch_(char *);
+    logical orgati;
+    real erretm, dtipsq, rhoinv;
+
+
+/*  -- LAPACK auxiliary routine (version 3.2) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/*     November 2006 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  This subroutine computes the square root of the I-th updated */
+/*  eigenvalue of a positive symmetric rank-one modification to */
+/*  a positive diagonal matrix whose entries are given as the squares */
+/*  of the corresponding entries in the array d, and that */
+
+/*         0 <= D(i) < D(j)  for  i < j */
+
+/*  and that RHO > 0. This is arranged by the calling routine, and is */
+/*  no loss in generality.  The rank-one modified system is thus */
+
+/*         diag( D ) * diag( D ) +  RHO *  Z * Z_transpose. */
+
+/*  where we assume the Euclidean norm of Z is 1. */
+
+/*  The method consists of approximating the rational functions in the */
+/*  secular equation by simpler interpolating rational functions. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  N      (input) INTEGER */
+/*         The length of all arrays. */
+
+/*  I      (input) INTEGER */
+/*         The index of the eigenvalue to be computed.  1 <= I <= N. */
+
+/*  D      (input) REAL array, dimension ( N ) */
+/*         The original eigenvalues.  It is assumed that they are in */
+/*         order, 0 <= D(I) < D(J)  for I < J. */
+
+/*  Z      (input) REAL array, dimension (N) */
+/*         The components of the updating vector. */
+
+/*  DELTA  (output) REAL array, dimension (N) */
+/*         If N .ne. 1, DELTA contains (D(j) - sigma_I) in its  j-th */
+/*         component.  If N = 1, then DELTA(1) = 1.  The vector DELTA */
+/*         contains the information necessary to construct the */
+/*         (singular) eigenvectors. */
+
+/*  RHO    (input) REAL */
+/*         The scalar in the symmetric updating formula. */
+
+/*  SIGMA  (output) REAL */
+/*         The computed sigma_I, the I-th updated eigenvalue. */
+
+/*  WORK   (workspace) REAL array, dimension (N) */
+/*         If N .ne. 1, WORK contains (D(j) + sigma_I) in its  j-th */
+/*         component.  If N = 1, then WORK( 1 ) = 1. */
+
+/*  INFO   (output) INTEGER */
+/*         = 0:  successful exit */
+/*         > 0:  if INFO = 1, the updating process failed. */
+
+/*  Internal Parameters */
+/*  =================== */
+
+/*  Logical variable ORGATI (origin-at-i?) is used for distinguishing */
+/*  whether D(i) or D(i+1) is treated as the origin. */
+
+/*            ORGATI = .true.    origin at i */
+/*            ORGATI = .false.   origin at i+1 */
+
+/*  Logical variable SWTCH3 (switch-for-3-poles?) is for noting */
+/*  if we are working with THREE poles! */
+
+/*  MAXIT is the maximum number of iterations allowed for each */
+/*  eigenvalue. */
+
+/*  Further Details */
+/*  =============== */
+
+/*  Based on contributions by */
+/*     Ren-Cang Li, Computer Science Division, University of California */
+/*     at Berkeley, USA */
+
+/*  ===================================================================== */
+
+/*     .. Parameters .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. Local Arrays .. */
+/*     .. */
+/*     .. External Subroutines .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     Since this routine is called in an inner loop, we do no argument */
+/*     checking. */
+
+/*     Quick return for N=1 and 2. */
+
+    /* Parameter adjustments */
+    --work;
+    --delta;
+    --z__;
+    --d__;
+
+    /* Function Body */
+    *info = 0;
+    if (*n == 1) {
+
+/*        Presumably, I=1 upon entry */
+
+	*sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
+	delta[1] = 1.f;
+	work[1] = 1.f;
+	return 0;
+    }
+    if (*n == 2) {
+	slasd5_(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
+	return 0;
+    }
+
+/*     Compute machine epsilon */
+
+    eps = slamch_("Epsilon");
+    rhoinv = 1.f / *rho;
+
+/*     The case I = N */
+
+    if (*i__ == *n) {
+
+/*        Initialize some basic variables */
+
+	ii = *n - 1;
+	niter = 1;
+
+/*        Calculate initial guess */
+
+	temp = *rho / 2.f;
+
+/*        If ||Z||_2 is not one, then TEMP should be set to */
+/*        RHO * ||Z||_2^2 / TWO */
+
+	temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    work[j] = d__[j] + d__[*n] + temp1;
+	    delta[j] = d__[j] - d__[*n] - temp1;
+/* L10: */
+	}
+
+	psi = 0.f;
+	i__1 = *n - 2;
+	for (j = 1; j <= i__1; ++j) {
+	    psi += z__[j] * z__[j] / (delta[j] * work[j]);
+/* L20: */
+	}
+
+	c__ = rhoinv + psi;
+	w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[*
+		n] / (delta[*n] * work[*n]);
+
+	if (w <= 0.f) {
+	    temp1 = sqrt(d__[*n] * d__[*n] + *rho);
+	    temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[*
+		    n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] * 
+		    z__[*n] / *rho;
+
+/*           The following TAU is to approximate */
+/*           SIGMA_n^2 - D( N )*D( N ) */
+
+	    if (c__ <= temp) {
+		tau = *rho;
+	    } else {
+		delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
+		a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*
+			n];
+		b = z__[*n] * z__[*n] * delsq;
+		if (a < 0.f) {
+		    tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
+		} else {
+		    tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
+		}
+	    }
+
+/*           It can be proved that */
+/*               D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO */
+
+	} else {
+	    delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
+	    a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
+	    b = z__[*n] * z__[*n] * delsq;
+
+/*           The following TAU is to approximate */
+/*           SIGMA_n^2 - D( N )*D( N ) */
+
+	    if (a < 0.f) {
+		tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
+	    } else {
+		tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
+	    }
+
+/*           It can be proved that */
+/*           D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2 */
+
+	}
+
+/*        The following ETA is to approximate SIGMA_n - D( N ) */
+
+	eta = tau / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau));
+
+	*sigma = d__[*n] + eta;
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    delta[j] = d__[j] - d__[*i__] - eta;
+	    work[j] = d__[j] + d__[*i__] + eta;
+/* L30: */
+	}
+
+/*        Evaluate PSI and the derivative DPSI */
+
+	dpsi = 0.f;
+	psi = 0.f;
+	erretm = 0.f;
+	i__1 = ii;
+	for (j = 1; j <= i__1; ++j) {
+	    temp = z__[j] / (delta[j] * work[j]);
+	    psi += z__[j] * temp;
+	    dpsi += temp * temp;
+	    erretm += psi;
+/* L40: */
+	}
+	erretm = dabs(erretm);
+
+/*        Evaluate PHI and the derivative DPHI */
+
+	temp = z__[*n] / (delta[*n] * work[*n]);
+	phi = z__[*n] * temp;
+	dphi = temp * temp;
+	erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
+		dpsi + dphi);
+
+	w = rhoinv + phi + psi;
+
+/*        Test for convergence */
+
+	if (dabs(w) <= eps * erretm) {
+	    goto L240;
+	}
+
+/*        Calculate the new step */
+
+	++niter;
+	dtnsq1 = work[*n - 1] * delta[*n - 1];
+	dtnsq = work[*n] * delta[*n];
+	c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
+	a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
+	b = dtnsq * dtnsq1 * w;
+	if (c__ < 0.f) {
+	    c__ = dabs(c__);
+	}
+	if (c__ == 0.f) {
+	    eta = *rho - *sigma * *sigma;
+	} else if (a >= 0.f) {
+	    eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
+		    c__ * 2.f);
+	} else {
+	    eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
+		    r__1))));
+	}
+
+/*        Note, eta should be positive if w is negative, and */
+/*        eta should be negative otherwise. However, */
+/*        if for some reason caused by roundoff, eta*w > 0, */
+/*        we simply use one Newton step instead. This way */
+/*        will guarantee eta*w < 0. */
+
+	if (w * eta > 0.f) {
+	    eta = -w / (dpsi + dphi);
+	}
+	temp = eta - dtnsq;
+	if (temp > *rho) {
+	    eta = *rho + dtnsq;
+	}
+
+	tau += eta;
+	eta /= *sigma + sqrt(eta + *sigma * *sigma);
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    delta[j] -= eta;
+	    work[j] += eta;
+/* L50: */
+	}
+
+	*sigma += eta;
+
+/*        Evaluate PSI and the derivative DPSI */
+
+	dpsi = 0.f;
+	psi = 0.f;
+	erretm = 0.f;
+	i__1 = ii;
+	for (j = 1; j <= i__1; ++j) {
+	    temp = z__[j] / (work[j] * delta[j]);
+	    psi += z__[j] * temp;
+	    dpsi += temp * temp;
+	    erretm += psi;
+/* L60: */
+	}
+	erretm = dabs(erretm);
+
+/*        Evaluate PHI and the derivative DPHI */
+
+	temp = z__[*n] / (work[*n] * delta[*n]);
+	phi = z__[*n] * temp;
+	dphi = temp * temp;
+	erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
+		dpsi + dphi);
+
+	w = rhoinv + phi + psi;
+
+/*        Main loop to update the values of the array   DELTA */
+
+	iter = niter + 1;
+
+	for (niter = iter; niter <= 20; ++niter) {
+
+/*           Test for convergence */
+
+	    if (dabs(w) <= eps * erretm) {
+		goto L240;
+	    }
+
+/*           Calculate the new step */
+
+	    dtnsq1 = work[*n - 1] * delta[*n - 1];
+	    dtnsq = work[*n] * delta[*n];
+	    c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
+	    a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
+	    b = dtnsq1 * dtnsq * w;
+	    if (a >= 0.f) {
+		eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
+			 (c__ * 2.f);
+	    } else {
+		eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
+			r__1))));
+	    }
+
+/*           Note, eta should be positive if w is negative, and */
+/*           eta should be negative otherwise. However, */
+/*           if for some reason caused by roundoff, eta*w > 0, */
+/*           we simply use one Newton step instead. This way */
+/*           will guarantee eta*w < 0. */
+
+	    if (w * eta > 0.f) {
+		eta = -w / (dpsi + dphi);
+	    }
+	    temp = eta - dtnsq;
+	    if (temp <= 0.f) {
+		eta /= 2.f;
+	    }
+
+	    tau += eta;
+	    eta /= *sigma + sqrt(eta + *sigma * *sigma);
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		delta[j] -= eta;
+		work[j] += eta;
+/* L70: */
+	    }
+
+	    *sigma += eta;
+
+/*           Evaluate PSI and the derivative DPSI */
+
+	    dpsi = 0.f;
+	    psi = 0.f;
+	    erretm = 0.f;
+	    i__1 = ii;
+	    for (j = 1; j <= i__1; ++j) {
+		temp = z__[j] / (work[j] * delta[j]);
+		psi += z__[j] * temp;
+		dpsi += temp * temp;
+		erretm += psi;
+/* L80: */
+	    }
+	    erretm = dabs(erretm);
+
+/*           Evaluate PHI and the derivative DPHI */
+
+	    temp = z__[*n] / (work[*n] * delta[*n]);
+	    phi = z__[*n] * temp;
+	    dphi = temp * temp;
+	    erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * 
+		    (dpsi + dphi);
+
+	    w = rhoinv + phi + psi;
+/* L90: */
+	}
+
+/*        Return with INFO = 1, NITER = MAXIT and not converged */
+
+	*info = 1;
+	goto L240;
+
+/*        End for the case I = N */
+
+    } else {
+
+/*        The case for I < N */
+
+	niter = 1;
+	ip1 = *i__ + 1;
+
+/*        Calculate initial guess */
+
+	delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
+	delsq2 = delsq / 2.f;
+	temp = delsq2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + delsq2));
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    work[j] = d__[j] + d__[*i__] + temp;
+	    delta[j] = d__[j] - d__[*i__] - temp;
+/* L100: */
+	}
+
+	psi = 0.f;
+	i__1 = *i__ - 1;
+	for (j = 1; j <= i__1; ++j) {
+	    psi += z__[j] * z__[j] / (work[j] * delta[j]);
+/* L110: */
+	}
+
+	phi = 0.f;
+	i__1 = *i__ + 2;
+	for (j = *n; j >= i__1; --j) {
+	    phi += z__[j] * z__[j] / (work[j] * delta[j]);
+/* L120: */
+	}
+	c__ = rhoinv + psi + phi;
+	w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[
+		ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
+
+	if (w > 0.f) {
+
+/*           d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2 */
+
+/*           We choose d(i) as origin. */
+
+	    orgati = TRUE_;
+	    sg2lb = 0.f;
+	    sg2ub = delsq2;
+	    a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
+	    b = z__[*i__] * z__[*i__] * delsq;
+	    if (a > 0.f) {
+		tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
+			r__1))));
+	    } else {
+		tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
+			 (c__ * 2.f);
+	    }
+
+/*           TAU now is an estimation of SIGMA^2 - D( I )^2. The */
+/*           following, however, is the corresponding estimation of */
+/*           SIGMA - D( I ). */
+
+	    eta = tau / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau));
+	} else {
+
+/*           (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2 */
+
+/*           We choose d(i+1) as origin. */
+
+	    orgati = FALSE_;
+	    sg2lb = -delsq2;
+	    sg2ub = 0.f;
+	    a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
+	    b = z__[ip1] * z__[ip1] * delsq;
+	    if (a < 0.f) {
+		tau = b * 2.f / (a - sqrt((r__1 = a * a + b * 4.f * c__, dabs(
+			r__1))));
+	    } else {
+		tau = -(a + sqrt((r__1 = a * a + b * 4.f * c__, dabs(r__1)))) 
+			/ (c__ * 2.f);
+	    }
+
+/*           TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The */
+/*           following, however, is the corresponding estimation of */
+/*           SIGMA - D( IP1 ). */
+
+	    eta = tau / (d__[ip1] + sqrt((r__1 = d__[ip1] * d__[ip1] + tau, 
+		    dabs(r__1))));
+	}
+
+	if (orgati) {
+	    ii = *i__;
+	    *sigma = d__[*i__] + eta;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		work[j] = d__[j] + d__[*i__] + eta;
+		delta[j] = d__[j] - d__[*i__] - eta;
+/* L130: */
+	    }
+	} else {
+	    ii = *i__ + 1;
+	    *sigma = d__[ip1] + eta;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		work[j] = d__[j] + d__[ip1] + eta;
+		delta[j] = d__[j] - d__[ip1] - eta;
+/* L140: */
+	    }
+	}
+	iim1 = ii - 1;
+	iip1 = ii + 1;
+
+/*        Evaluate PSI and the derivative DPSI */
+
+	dpsi = 0.f;
+	psi = 0.f;
+	erretm = 0.f;
+	i__1 = iim1;
+	for (j = 1; j <= i__1; ++j) {
+	    temp = z__[j] / (work[j] * delta[j]);
+	    psi += z__[j] * temp;
+	    dpsi += temp * temp;
+	    erretm += psi;
+/* L150: */
+	}
+	erretm = dabs(erretm);
+
+/*        Evaluate PHI and the derivative DPHI */
+
+	dphi = 0.f;
+	phi = 0.f;
+	i__1 = iip1;
+	for (j = *n; j >= i__1; --j) {
+	    temp = z__[j] / (work[j] * delta[j]);
+	    phi += z__[j] * temp;
+	    dphi += temp * temp;
+	    erretm += phi;
+/* L160: */
+	}
+
+	w = rhoinv + phi + psi;
+
+/*        W is the value of the secular function with */
+/*        its ii-th element removed. */
+
+	swtch3 = FALSE_;
+	if (orgati) {
+	    if (w < 0.f) {
+		swtch3 = TRUE_;
+	    }
+	} else {
+	    if (w > 0.f) {
+		swtch3 = TRUE_;
+	    }
+	}
+	if (ii == 1 || ii == *n) {
+	    swtch3 = FALSE_;
+	}
+
+	temp = z__[ii] / (work[ii] * delta[ii]);
+	dw = dpsi + dphi + temp * temp;
+	temp = z__[ii] * temp;
+	w += temp;
+	erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f 
+		+ dabs(tau) * dw;
+
+/*        Test for convergence */
+
+	if (dabs(w) <= eps * erretm) {
+	    goto L240;
+	}
+
+	if (w <= 0.f) {
+	    sg2lb = dmax(sg2lb,tau);
+	} else {
+	    sg2ub = dmin(sg2ub,tau);
+	}
+
+/*        Calculate the new step */
+
+	++niter;
+	if (! swtch3) {
+	    dtipsq = work[ip1] * delta[ip1];
+	    dtisq = work[*i__] * delta[*i__];
+	    if (orgati) {
+/* Computing 2nd power */
+		r__1 = z__[*i__] / dtisq;
+		c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
+	    } else {
+/* Computing 2nd power */
+		r__1 = z__[ip1] / dtipsq;
+		c__ = w - dtisq * dw - delsq * (r__1 * r__1);
+	    }
+	    a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
+	    b = dtipsq * dtisq * w;
+	    if (c__ == 0.f) {
+		if (a == 0.f) {
+		    if (orgati) {
+			a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi + 
+				dphi);
+		    } else {
+			a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi + 
+				dphi);
+		    }
+		}
+		eta = b / a;
+	    } else if (a <= 0.f) {
+		eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
+			 (c__ * 2.f);
+	    } else {
+		eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
+			r__1))));
+	    }
+	} else {
+
+/*           Interpolation using THREE most relevant poles */
+
+	    dtiim = work[iim1] * delta[iim1];
+	    dtiip = work[iip1] * delta[iip1];
+	    temp = rhoinv + psi + phi;
+	    if (orgati) {
+		temp1 = z__[iim1] / dtiim;
+		temp1 *= temp1;
+		c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) *
+			 (d__[iim1] + d__[iip1]) * temp1;
+		zz[0] = z__[iim1] * z__[iim1];
+		if (dpsi < temp1) {
+		    zz[2] = dtiip * dtiip * dphi;
+		} else {
+		    zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
+		}
+	    } else {
+		temp1 = z__[iip1] / dtiip;
+		temp1 *= temp1;
+		c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) *
+			 (d__[iim1] + d__[iip1]) * temp1;
+		if (dphi < temp1) {
+		    zz[0] = dtiim * dtiim * dpsi;
+		} else {
+		    zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
+		}
+		zz[2] = z__[iip1] * z__[iip1];
+	    }
+	    zz[1] = z__[ii] * z__[ii];
+	    dd[0] = dtiim;
+	    dd[1] = delta[ii] * work[ii];
+	    dd[2] = dtiip;
+	    slaed6_(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
+	    if (*info != 0) {
+		goto L240;
+	    }
+	}
+
+/*        Note, eta should be positive if w is negative, and */
+/*        eta should be negative otherwise. However, */
+/*        if for some reason caused by roundoff, eta*w > 0, */
+/*        we simply use one Newton step instead. This way */
+/*        will guarantee eta*w < 0. */
+
+	if (w * eta >= 0.f) {
+	    eta = -w / dw;
+	}
+	if (orgati) {
+	    temp1 = work[*i__] * delta[*i__];
+	    temp = eta - temp1;
+	} else {
+	    temp1 = work[ip1] * delta[ip1];
+	    temp = eta - temp1;
+	}
+	if (temp > sg2ub || temp < sg2lb) {
+	    if (w < 0.f) {
+		eta = (sg2ub - tau) / 2.f;
+	    } else {
+		eta = (sg2lb - tau) / 2.f;
+	    }
+	}
+
+	tau += eta;
+	eta /= *sigma + sqrt(*sigma * *sigma + eta);
+
+	prew = w;
+
+	*sigma += eta;
+	i__1 = *n;
+	for (j = 1; j <= i__1; ++j) {
+	    work[j] += eta;
+	    delta[j] -= eta;
+/* L170: */
+	}
+
+/*        Evaluate PSI and the derivative DPSI */
+
+	dpsi = 0.f;
+	psi = 0.f;
+	erretm = 0.f;
+	i__1 = iim1;
+	for (j = 1; j <= i__1; ++j) {
+	    temp = z__[j] / (work[j] * delta[j]);
+	    psi += z__[j] * temp;
+	    dpsi += temp * temp;
+	    erretm += psi;
+/* L180: */
+	}
+	erretm = dabs(erretm);
+
+/*        Evaluate PHI and the derivative DPHI */
+
+	dphi = 0.f;
+	phi = 0.f;
+	i__1 = iip1;
+	for (j = *n; j >= i__1; --j) {
+	    temp = z__[j] / (work[j] * delta[j]);
+	    phi += z__[j] * temp;
+	    dphi += temp * temp;
+	    erretm += phi;
+/* L190: */
+	}
+
+	temp = z__[ii] / (work[ii] * delta[ii]);
+	dw = dpsi + dphi + temp * temp;
+	temp = z__[ii] * temp;
+	w = rhoinv + phi + psi + temp;
+	erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f 
+		+ dabs(tau) * dw;
+
+	if (w <= 0.f) {
+	    sg2lb = dmax(sg2lb,tau);
+	} else {
+	    sg2ub = dmin(sg2ub,tau);
+	}
+
+	swtch = FALSE_;
+	if (orgati) {
+	    if (-w > dabs(prew) / 10.f) {
+		swtch = TRUE_;
+	    }
+	} else {
+	    if (w > dabs(prew) / 10.f) {
+		swtch = TRUE_;
+	    }
+	}
+
+/*        Main loop to update the values of the array   DELTA and WORK */
+
+	iter = niter + 1;
+
+	for (niter = iter; niter <= 20; ++niter) {
+
+/*           Test for convergence */
+
+	    if (dabs(w) <= eps * erretm) {
+		goto L240;
+	    }
+
+/*           Calculate the new step */
+
+	    if (! swtch3) {
+		dtipsq = work[ip1] * delta[ip1];
+		dtisq = work[*i__] * delta[*i__];
+		if (! swtch) {
+		    if (orgati) {
+/* Computing 2nd power */
+			r__1 = z__[*i__] / dtisq;
+			c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
+		    } else {
+/* Computing 2nd power */
+			r__1 = z__[ip1] / dtipsq;
+			c__ = w - dtisq * dw - delsq * (r__1 * r__1);
+		    }
+		} else {
+		    temp = z__[ii] / (work[ii] * delta[ii]);
+		    if (orgati) {
+			dpsi += temp * temp;
+		    } else {
+			dphi += temp * temp;
+		    }
+		    c__ = w - dtisq * dpsi - dtipsq * dphi;
+		}
+		a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
+		b = dtipsq * dtisq * w;
+		if (c__ == 0.f) {
+		    if (a == 0.f) {
+			if (! swtch) {
+			    if (orgati) {
+				a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * 
+					(dpsi + dphi);
+			    } else {
+				a = z__[ip1] * z__[ip1] + dtisq * dtisq * (
+					dpsi + dphi);
+			    }
+			} else {
+			    a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
+			}
+		    }
+		    eta = b / a;
+		} else if (a <= 0.f) {
+		    eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1))
+			    )) / (c__ * 2.f);
+		} else {
+		    eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, 
+			    dabs(r__1))));
+		}
+	    } else {
+
+/*              Interpolation using THREE most relevant poles */
+
+		dtiim = work[iim1] * delta[iim1];
+		dtiip = work[iip1] * delta[iip1];
+		temp = rhoinv + psi + phi;
+		if (swtch) {
+		    c__ = temp - dtiim * dpsi - dtiip * dphi;
+		    zz[0] = dtiim * dtiim * dpsi;
+		    zz[2] = dtiip * dtiip * dphi;
+		} else {
+		    if (orgati) {
+			temp1 = z__[iim1] / dtiim;
+			temp1 *= temp1;
+			temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[
+				iip1]) * temp1;
+			c__ = temp - dtiip * (dpsi + dphi) - temp2;
+			zz[0] = z__[iim1] * z__[iim1];
+			if (dpsi < temp1) {
+			    zz[2] = dtiip * dtiip * dphi;
+			} else {
+			    zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
+			}
+		    } else {
+			temp1 = z__[iip1] / dtiip;
+			temp1 *= temp1;
+			temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[
+				iip1]) * temp1;
+			c__ = temp - dtiim * (dpsi + dphi) - temp2;
+			if (dphi < temp1) {
+			    zz[0] = dtiim * dtiim * dpsi;
+			} else {
+			    zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
+			}
+			zz[2] = z__[iip1] * z__[iip1];
+		    }
+		}
+		dd[0] = dtiim;
+		dd[1] = delta[ii] * work[ii];
+		dd[2] = dtiip;
+		slaed6_(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
+		if (*info != 0) {
+		    goto L240;
+		}
+	    }
+
+/*           Note, eta should be positive if w is negative, and */
+/*           eta should be negative otherwise. However, */
+/*           if for some reason caused by roundoff, eta*w > 0, */
+/*           we simply use one Newton step instead. This way */
+/*           will guarantee eta*w < 0. */
+
+	    if (w * eta >= 0.f) {
+		eta = -w / dw;
+	    }
+	    if (orgati) {
+		temp1 = work[*i__] * delta[*i__];
+		temp = eta - temp1;
+	    } else {
+		temp1 = work[ip1] * delta[ip1];
+		temp = eta - temp1;
+	    }
+	    if (temp > sg2ub || temp < sg2lb) {
+		if (w < 0.f) {
+		    eta = (sg2ub - tau) / 2.f;
+		} else {
+		    eta = (sg2lb - tau) / 2.f;
+		}
+	    }
+
+	    tau += eta;
+	    eta /= *sigma + sqrt(*sigma * *sigma + eta);
+
+	    *sigma += eta;
+	    i__1 = *n;
+	    for (j = 1; j <= i__1; ++j) {
+		work[j] += eta;
+		delta[j] -= eta;
+/* L200: */
+	    }
+
+	    prew = w;
+
+/*           Evaluate PSI and the derivative DPSI */
+
+	    dpsi = 0.f;
+	    psi = 0.f;
+	    erretm = 0.f;
+	    i__1 = iim1;
+	    for (j = 1; j <= i__1; ++j) {
+		temp = z__[j] / (work[j] * delta[j]);
+		psi += z__[j] * temp;
+		dpsi += temp * temp;
+		erretm += psi;
+/* L210: */
+	    }
+	    erretm = dabs(erretm);
+
+/*           Evaluate PHI and the derivative DPHI */
+
+	    dphi = 0.f;
+	    phi = 0.f;
+	    i__1 = iip1;
+	    for (j = *n; j >= i__1; --j) {
+		temp = z__[j] / (work[j] * delta[j]);
+		phi += z__[j] * temp;
+		dphi += temp * temp;
+		erretm += phi;
+/* L220: */
+	    }
+
+	    temp = z__[ii] / (work[ii] * delta[ii]);
+	    dw = dpsi + dphi + temp * temp;
+	    temp = z__[ii] * temp;
+	    w = rhoinv + phi + psi + temp;
+	    erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 
+		    3.f + dabs(tau) * dw;
+	    if (w * prew > 0.f && dabs(w) > dabs(prew) / 10.f) {
+		swtch = ! swtch;
+	    }
+
+	    if (w <= 0.f) {
+		sg2lb = dmax(sg2lb,tau);
+	    } else {
+		sg2ub = dmin(sg2ub,tau);
+	    }
+
+/* L230: */
+	}
+
+/*        Return with INFO = 1, NITER = MAXIT and not converged */
+
+	*info = 1;
+
+    }
+
+L240:
+    return 0;
+
+/*     End of SLASD4 */
+
+} /* slasd4_ */