#include "blaswrap.h"
/* slaed4.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"
/* Subroutine */ int slaed4_(integer *n, integer *i__, real *d__, real *z__,
real *delta, real *rho, real *dlam, integer *info)
{
/* System generated locals */
integer i__1;
real r__1;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static real a, b, c__;
static integer j;
static real w;
static integer ii;
static real dw, zz[3];
static integer ip1;
static real del, eta, phi, eps, tau, psi;
static integer iim1, iip1;
static real dphi, dpsi;
static integer iter;
static real temp, prew, temp1, dltlb, dltub, midpt;
static integer niter;
static logical swtch;
extern /* Subroutine */ int slaed5_(integer *, real *, real *, real *,
real *, real *), slaed6_(integer *, logical *, real *, real *,
real *, real *, real *, integer *);
static logical swtch3;
extern doublereal slamch_(char *);
static logical orgati;
static real erretm, rhoinv;
/* -- LAPACK routine (version 3.1) --
Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
November 2006
Purpose
=======
This subroutine computes the I-th updated eigenvalue of a symmetric
rank-one modification to a diagonal matrix whose elements are
given in the array d, and that
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 ) + 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, 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 .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th
component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5
for detail. The vector DELTA contains the information necessary
to construct the eigenvectors by SLAED3 and SLAED9.
RHO (input) REAL
The scalar in the symmetric updating formula.
DLAM (output) REAL
The computed lambda_I, the I-th updated eigenvalue.
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
=====================================================================
Since this routine is called in an inner loop, we do no argument
checking.
Quick return for N=1 and 2.
Parameter adjustments */
--delta;
--z__;
--d__;
/* Function Body */
*info = 0;
if (*n == 1) {
/* Presumably, I=1 upon entry */
*dlam = d__[1] + *rho * z__[1] * z__[1];
delta[1] = 1.f;
return 0;
}
if (*n == 2) {
slaed5_(i__, &d__[1], &z__[1], &delta[1], rho, dlam);
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 */
midpt = *rho / 2.f;
/* If ||Z||_2 is not one, then TEMP should be set to
RHO * ||Z||_2^2 / TWO */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] = d__[j] - d__[*i__] - midpt;
/* L10: */
}
psi = 0.f;
i__1 = *n - 2;
for (j = 1; j <= i__1; ++j) {
psi += z__[j] * z__[j] / delta[j];
/* L20: */
}
c__ = rhoinv + psi;
w = c__ + z__[ii] * z__[ii] / delta[ii] + z__[*n] * z__[*n] / delta[*
n];
if (w <= 0.f) {
temp = z__[*n - 1] * z__[*n - 1] / (d__[*n] - d__[*n - 1] + *rho)
+ z__[*n] * z__[*n] / *rho;
if (c__ <= temp) {
tau = *rho;
} else {
del = d__[*n] - d__[*n - 1];
a = -c__ * del + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n]
;
b = z__[*n] * z__[*n] * del;
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)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO */
dltlb = midpt;
dltub = *rho;
} else {
del = d__[*n] - d__[*n - 1];
a = -c__ * del + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
b = z__[*n] * z__[*n] * del;
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) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 */
dltlb = 0.f;
dltub = midpt;
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] = d__[j] - d__[*i__] - tau;
/* 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];
psi += z__[j] * temp;
dpsi += temp * temp;
erretm += psi;
/* L40: */
}
erretm = dabs(erretm);
/* Evaluate PHI and the derivative DPHI */
temp = z__[*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;
/* Test for convergence */
if (dabs(w) <= eps * erretm) {
*dlam = d__[*i__] + tau;
goto L250;
}
if (w <= 0.f) {
dltlb = dmax(dltlb,tau);
} else {
dltub = dmin(dltub,tau);
}
/* Calculate the new step */
++niter;
c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] * (
dpsi + dphi);
b = delta[*n - 1] * delta[*n] * w;
if (c__ < 0.f) {
c__ = dabs(c__);
}
if (c__ == 0.f) {
/* ETA = B/A
ETA = RHO - TAU */
eta = dltub - tau;
} 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 = tau + eta;
if (temp > dltub || temp < dltlb) {
if (w < 0.f) {
eta = (dltub - tau) / 2.f;
} else {
eta = (dltlb - tau) / 2.f;
}
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] -= eta;
/* L50: */
}
tau += 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] / delta[j];
psi += z__[j] * temp;
dpsi += temp * temp;
erretm += psi;
/* L60: */
}
erretm = dabs(erretm);
/* Evaluate PHI and the derivative DPHI */
temp = z__[*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 <= 30; ++niter) {
/* Test for convergence */
if (dabs(w) <= eps * erretm) {
*dlam = d__[*i__] + tau;
goto L250;
}
if (w <= 0.f) {
dltlb = dmax(dltlb,tau);
} else {
dltub = dmin(dltub,tau);
}
/* Calculate the new step */
c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] *
(dpsi + dphi);
b = delta[*n - 1] * delta[*n] * 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 = tau + eta;
if (temp > dltub || temp < dltlb) {
if (w < 0.f) {
eta = (dltub - tau) / 2.f;
} else {
eta = (dltlb - tau) / 2.f;
}
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] -= eta;
/* L70: */
}
tau += 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] / delta[j];
psi += z__[j] * temp;
dpsi += temp * temp;
erretm += psi;
/* L80: */
}
erretm = dabs(erretm);
/* Evaluate PHI and the derivative DPHI */
temp = z__[*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;
*dlam = d__[*i__] + tau;
goto L250;
/* End for the case I = N */
} else {
/* The case for I < N */
niter = 1;
ip1 = *i__ + 1;
/* Calculate initial guess */
del = d__[ip1] - d__[*i__];
midpt = del / 2.f;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] = d__[j] - d__[*i__] - midpt;
/* L100: */
}
psi = 0.f;
i__1 = *i__ - 1;
for (j = 1; j <= i__1; ++j) {
psi += z__[j] * z__[j] / delta[j];
/* L110: */
}
phi = 0.f;
i__1 = *i__ + 2;
for (j = *n; j >= i__1; --j) {
phi += z__[j] * z__[j] / delta[j];
/* L120: */
}
c__ = rhoinv + psi + phi;
w = c__ + z__[*i__] * z__[*i__] / delta[*i__] + z__[ip1] * z__[ip1] /
delta[ip1];
if (w > 0.f) {
/* d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
We choose d(i) as origin. */
orgati = TRUE_;
a = c__ * del + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
b = z__[*i__] * z__[*i__] * del;
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);
}
dltlb = 0.f;
dltub = midpt;
} else {
/* (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
We choose d(i+1) as origin. */
orgati = FALSE_;
a = c__ * del - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
b = z__[ip1] * z__[ip1] * del;
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);
}
dltlb = -midpt;
dltub = 0.f;
}
if (orgati) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] = d__[j] - d__[*i__] - tau;
/* L130: */
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] = d__[j] - d__[ip1] - tau;
/* L140: */
}
}
if (orgati) {
ii = *i__;
} else {
ii = *i__ + 1;
}
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] / 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] / 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] / 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) {
if (orgati) {
*dlam = d__[*i__] + tau;
} else {
*dlam = d__[ip1] + tau;
}
goto L250;
}
if (w <= 0.f) {
dltlb = dmax(dltlb,tau);
} else {
dltub = dmin(dltub,tau);
}
/* Calculate the new step */
++niter;
if (! swtch3) {
if (orgati) {
/* Computing 2nd power */
r__1 = z__[*i__] / delta[*i__];
c__ = w - delta[ip1] * dw - (d__[*i__] - d__[ip1]) * (r__1 *
r__1);
} else {
/* Computing 2nd power */
r__1 = z__[ip1] / delta[ip1];
c__ = w - delta[*i__] * dw - (d__[ip1] - d__[*i__]) * (r__1 *
r__1);
}
a = (delta[*i__] + delta[ip1]) * w - delta[*i__] * delta[ip1] *
dw;
b = delta[*i__] * delta[ip1] * w;
if (c__ == 0.f) {
if (a == 0.f) {
if (orgati) {
a = z__[*i__] * z__[*i__] + delta[ip1] * delta[ip1] *
(dpsi + dphi);
} else {
a = z__[ip1] * z__[ip1] + delta[*i__] * delta[*i__] *
(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 */
temp = rhoinv + psi + phi;
if (orgati) {
temp1 = z__[iim1] / delta[iim1];
temp1 *= temp1;
c__ = temp - delta[iip1] * (dpsi + dphi) - (d__[iim1] - d__[
iip1]) * temp1;
zz[0] = z__[iim1] * z__[iim1];
zz[2] = delta[iip1] * delta[iip1] * (dpsi - temp1 + dphi);
} else {
temp1 = z__[iip1] / delta[iip1];
temp1 *= temp1;
c__ = temp - delta[iim1] * (dpsi + dphi) - (d__[iip1] - d__[
iim1]) * temp1;
zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi - temp1));
zz[2] = z__[iip1] * z__[iip1];
}
zz[1] = z__[ii] * z__[ii];
slaed6_(&niter, &orgati, &c__, &delta[iim1], zz, &w, &eta, info);
if (*info != 0) {
goto L250;
}
}
/* 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;
}
temp = tau + eta;
if (temp > dltub || temp < dltlb) {
if (w < 0.f) {
eta = (dltub - tau) / 2.f;
} else {
eta = (dltlb - tau) / 2.f;
}
}
prew = w;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] -= eta;
/* L180: */
}
/* 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] / delta[j];
psi += z__[j] * temp;
dpsi += temp * temp;
erretm += psi;
/* L190: */
}
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] / delta[j];
phi += z__[j] * temp;
dphi += temp * temp;
erretm += phi;
/* L200: */
}
temp = z__[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
+ (r__1 = tau + eta, dabs(r__1)) * dw;
swtch = FALSE_;
if (orgati) {
if (-w > dabs(prew) / 10.f) {
swtch = TRUE_;
}
} else {
if (w > dabs(prew) / 10.f) {
swtch = TRUE_;
}
}
tau += eta;
/* Main loop to update the values of the array DELTA */
iter = niter + 1;
for (niter = iter; niter <= 30; ++niter) {
/* Test for convergence */
if (dabs(w) <= eps * erretm) {
if (orgati) {
*dlam = d__[*i__] + tau;
} else {
*dlam = d__[ip1] + tau;
}
goto L250;
}
if (w <= 0.f) {
dltlb = dmax(dltlb,tau);
} else {
dltub = dmin(dltub,tau);
}
/* Calculate the new step */
if (! swtch3) {
if (! swtch) {
if (orgati) {
/* Computing 2nd power */
r__1 = z__[*i__] / delta[*i__];
c__ = w - delta[ip1] * dw - (d__[*i__] - d__[ip1]) * (
r__1 * r__1);
} else {
/* Computing 2nd power */
r__1 = z__[ip1] / delta[ip1];
c__ = w - delta[*i__] * dw - (d__[ip1] - d__[*i__]) *
(r__1 * r__1);
}
} else {
temp = z__[ii] / delta[ii];
if (orgati) {
dpsi += temp * temp;
} else {
dphi += temp * temp;
}
c__ = w - delta[*i__] * dpsi - delta[ip1] * dphi;
}
a = (delta[*i__] + delta[ip1]) * w - delta[*i__] * delta[ip1]
* dw;
b = delta[*i__] * delta[ip1] * w;
if (c__ == 0.f) {
if (a == 0.f) {
if (! swtch) {
if (orgati) {
a = z__[*i__] * z__[*i__] + delta[ip1] *
delta[ip1] * (dpsi + dphi);
} else {
a = z__[ip1] * z__[ip1] + delta[*i__] * delta[
*i__] * (dpsi + dphi);
}
} else {
a = delta[*i__] * delta[*i__] * dpsi + delta[ip1]
* delta[ip1] * 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 */
temp = rhoinv + psi + phi;
if (swtch) {
c__ = temp - delta[iim1] * dpsi - delta[iip1] * dphi;
zz[0] = delta[iim1] * delta[iim1] * dpsi;
zz[2] = delta[iip1] * delta[iip1] * dphi;
} else {
if (orgati) {
temp1 = z__[iim1] / delta[iim1];
temp1 *= temp1;
c__ = temp - delta[iip1] * (dpsi + dphi) - (d__[iim1]
- d__[iip1]) * temp1;
zz[0] = z__[iim1] * z__[iim1];
zz[2] = delta[iip1] * delta[iip1] * (dpsi - temp1 +
dphi);
} else {
temp1 = z__[iip1] / delta[iip1];
temp1 *= temp1;
c__ = temp - delta[iim1] * (dpsi + dphi) - (d__[iip1]
- d__[iim1]) * temp1;
zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi -
temp1));
zz[2] = z__[iip1] * z__[iip1];
}
}
slaed6_(&niter, &orgati, &c__, &delta[iim1], zz, &w, &eta,
info);
if (*info != 0) {
goto L250;
}
}
/* 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;
}
temp = tau + eta;
if (temp > dltub || temp < dltlb) {
if (w < 0.f) {
eta = (dltub - tau) / 2.f;
} else {
eta = (dltlb - tau) / 2.f;
}
}
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
delta[j] -= eta;
/* L210: */
}
tau += eta;
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] / delta[j];
psi += z__[j] * temp;
dpsi += temp * temp;
erretm += psi;
/* L220: */
}
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] / delta[j];
phi += z__[j] * temp;
dphi += temp * temp;
erretm += phi;
/* L230: */
}
temp = z__[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;
}
/* L240: */
}
/* Return with INFO = 1, NITER = MAXIT and not converged */
*info = 1;
if (orgati) {
*dlam = d__[*i__] + tau;
} else {
*dlam = d__[ip1] + tau;
}
}
L250:
return 0;
/* End of SLAED4 */
} /* slaed4_ */