dlasd4.c
来自「提供矩阵类的函数库」· C语言 代码 · 共 1,087 行 · 第 1/2 页
C
1,087 行
}
/* 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. */
latime_1.ops += 20.;
orgati = FALSE_;
sg2lb = -delsq2;
sg2ub = 0.;
a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
b = z__[ip1] * z__[ip1] * delsq;
if (a < 0.) {
tau = b * 2. / (a - sqrt((d__1 = a * a + b * 4. * c__, abs(
d__1))));
} else {
tau = -(a + sqrt((d__1 = a * a + b * 4. * c__, abs(d__1)))) /
(c__ * 2.);
}
/* 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((d__1 = d__[ip1] * d__[ip1] + tau,
abs(d__1))));
}
latime_1.ops += (doublereal) ((*n << 2) + 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.;
psi = 0.;
erretm = 0.;
latime_1.ops += (doublereal) (iim1 * 7);
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 = abs(erretm);
/* Evaluate PHI and the derivative DPHI */
dphi = 0.;
phi = 0.;
latime_1.ops += (doublereal) ((*n - iip1 + 1) * 7 + 2);
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.) {
swtch3 = TRUE_;
}
} else {
if (w > 0.) {
swtch3 = TRUE_;
}
}
if (ii == 1 || ii == *n) {
swtch3 = FALSE_;
}
latime_1.ops += 17.;
temp = z__[ii] / (work[ii] * delta[ii]);
dw = dpsi + dphi + temp * temp;
temp = z__[ii] * temp;
w += temp;
erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + abs(temp) * 3. +
abs(tau) * dw;
/* Test for convergence */
if (abs(w) <= eps * erretm) {
goto L240;
}
if (w <= 0.) {
sg2lb = max(sg2lb,tau);
} else {
sg2ub = min(sg2ub,tau);
}
/* Calculate the new step */
++niter;
if (! swtch3) {
latime_1.ops += 15.;
dtipsq = work[ip1] * delta[ip1];
dtisq = work[*i__] * delta[*i__];
if (orgati) {
/* Computing 2nd power */
d__1 = z__[*i__] / dtisq;
c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
} else {
/* Computing 2nd power */
d__1 = z__[ip1] / dtipsq;
c__ = w - dtisq * dw - delsq * (d__1 * d__1);
}
a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
b = dtipsq * dtisq * w;
if (c__ == 0.) {
if (a == 0.) {
latime_1.ops += 5.;
if (orgati) {
a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi +
dphi);
} else {
a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi +
dphi);
}
}
latime_1.ops += 1.;
eta = b / a;
} else if (a <= 0.) {
latime_1.ops += 8.;
eta = (a - sqrt((d__1 = a * a - b * 4. * c__, abs(d__1)))) / (
c__ * 2.);
} else {
latime_1.ops += 8.;
eta = b * 2. / (a + sqrt((d__1 = a * a - b * 4. * c__, abs(
d__1))));
}
} else {
/* Interpolation using THREE most relevant poles */
latime_1.ops += 15.;
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) {
latime_1.ops += 2.;
zz[2] = dtiip * dtiip * dphi;
} else {
latime_1.ops += 4.;
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) {
latime_1.ops += 2.;
zz[0] = dtiim * dtiim * dpsi;
} else {
latime_1.ops += 4.;
zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
}
zz[2] = z__[iip1] * z__[iip1];
}
latime_1.ops += 2.;
zz[1] = z__[ii] * z__[ii];
dd[0] = dtiim;
dd[1] = delta[ii] * work[ii];
dd[2] = dtiip;
dlaed6_(&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. */
latime_1.ops += 1.;
if (w * eta >= 0.) {
latime_1.ops += 1.;
eta = -w / dw;
}
latime_1.ops += 8.;
if (orgati) {
temp1 = work[*i__] * delta[*i__];
temp = eta - temp1;
} else {
temp1 = work[ip1] * delta[ip1];
temp = eta - temp1;
}
if (temp > sg2ub || temp < sg2lb) {
latime_1.ops += 2.;
if (w < 0.) {
eta = (sg2ub - tau) / 2.;
} else {
eta = (sg2lb - tau) / 2.;
}
}
tau += eta;
eta /= *sigma + sqrt(*sigma * *sigma + eta);
prew = w;
latime_1.ops += (doublereal) ((*n << 1) + 1);
*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.;
psi = 0.;
erretm = 0.;
latime_1.ops += (doublereal) (iim1 * 7);
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 = abs(erretm);
/* Evaluate PHI and the derivative DPHI */
dphi = 0.;
phi = 0.;
latime_1.ops += (doublereal) ((*n - iim1 + 1) * 7);
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: */
}
latime_1.ops += 19.;
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. + erretm + rhoinv * 2. + abs(temp) * 3. +
abs(tau) * dw;
if (w <= 0.) {
sg2lb = max(sg2lb,tau);
} else {
sg2ub = min(sg2ub,tau);
}
swtch = FALSE_;
if (orgati) {
if (-w > abs(prew) / 10.) {
swtch = TRUE_;
}
} else {
if (w > abs(prew) / 10.) {
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 */
latime_1.ops += 1.;
if (abs(w) <= eps * erretm) {
goto L240;
}
/* Calculate the new step */
if (! swtch3) {
latime_1.ops += 2.;
dtipsq = work[ip1] * delta[ip1];
dtisq = work[*i__] * delta[*i__];
if (! swtch) {
latime_1.ops += 6.;
if (orgati) {
/* Computing 2nd power */
d__1 = z__[*i__] / dtisq;
c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
} else {
/* Computing 2nd power */
d__1 = z__[ip1] / dtipsq;
c__ = w - dtisq * dw - delsq * (d__1 * d__1);
}
} else {
latime_1.ops += 8.;
temp = z__[ii] / (work[ii] * delta[ii]);
if (orgati) {
dpsi += temp * temp;
} else {
dphi += temp * temp;
}
c__ = w - dtisq * dpsi - dtipsq * dphi;
}
latime_1.ops += 7.;
a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
b = dtipsq * dtisq * w;
if (c__ == 0.) {
if (a == 0.) {
latime_1.ops += 5.;
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;
}
}
latime_1.ops += 1.;
eta = b / a;
} else if (a <= 0.) {
latime_1.ops += 8.;
eta = (a - sqrt((d__1 = a * a - b * 4. * c__, abs(d__1))))
/ (c__ * 2.);
} else {
latime_1.ops += 8.;
eta = b * 2. / (a + sqrt((d__1 = a * a - b * 4. * c__,
abs(d__1))));
}
} else {
/* Interpolation using THREE most relevant poles */
latime_1.ops += 4.;
dtiim = work[iim1] * delta[iim1];
dtiip = work[iip1] * delta[iip1];
temp = rhoinv + psi + phi;
if (swtch) {
latime_1.ops += 8.;
c__ = temp - dtiim * dpsi - dtiip * dphi;
zz[0] = dtiim * dtiim * dpsi;
zz[2] = dtiip * dtiip * dphi;
} else {
if (orgati) {
latime_1.ops += 11.;
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) {
latime_1.ops += 2.;
zz[2] = dtiip * dtiip * dphi;
} else {
latime_1.ops += 4.;
zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
}
} else {
latime_1.ops += 10.;
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) {
latime_1.ops += 2.;
zz[0] = dtiim * dtiim * dpsi;
} else {
latime_1.ops += 4.;
zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
}
latime_1.ops += 1.;
zz[2] = z__[iip1] * z__[iip1];
}
}
latime_1.ops += 1.;
dd[0] = dtiim;
dd[1] = delta[ii] * work[ii];
dd[2] = dtiip;
dlaed6_(&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. */
latime_1.ops += 1.;
if (w * eta >= 0.) {
latime_1.ops += 1.;
eta = -w / dw;
}
latime_1.ops += 2.;
if (orgati) {
temp1 = work[*i__] * delta[*i__];
temp = eta - temp1;
} else {
temp1 = work[ip1] * delta[ip1];
temp = eta - temp1;
}
if (temp > sg2ub || temp < sg2lb) {
latime_1.ops += 2.;
if (w < 0.) {
eta = (sg2ub - tau) / 2.;
} else {
eta = (sg2lb - tau) / 2.;
}
}
latime_1.ops += 6.;
tau += eta;
eta /= *sigma + sqrt(*sigma * *sigma + eta);
latime_1.ops += (doublereal) ((*n << 1) + 1);
*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.;
psi = 0.;
erretm = 0.;
latime_1.ops += (doublereal) (iim1 * 7);
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 = abs(erretm);
/* Evaluate PHI and the derivative DPHI */
dphi = 0.;
phi = 0.;
latime_1.ops += (doublereal) ((iim1 - *n + 1) * 7);
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: */
}
latime_1.ops += 19.;
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. + erretm + rhoinv * 2. + abs(temp) * 3.
+ abs(tau) * dw;
if (w * prew > 0. && abs(w) > abs(prew) / 10.) {
swtch = ! swtch;
}
if (w <= 0.) {
sg2lb = max(sg2lb,tau);
} else {
sg2ub = min(sg2ub,tau);
}
/* L230: */
}
/* Return with INFO = 1, NITER = MAXIT and not converged */
*info = 1;
}
L240:
return 0;
/* End of DLASD4 */
} /* dlasd4_ */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?