3 /* Subroutine */ int slaed4_(integer *n, integer *i__, real *d__, real *z__,
4 real *delta, real *rho, real *dlam, integer *info)
6 /* System generated locals */
10 /* Builtin functions */
11 double sqrt(doublereal);
20 real del, eta, phi, eps, tau, psi;
24 real temp, prew, temp1, dltlb, dltub, midpt;
27 extern /* Subroutine */ int slaed5_(integer *, real *, real *, real *,
28 real *, real *), slaed6_(integer *, logical *, real *, real *,
29 real *, real *, real *, integer *);
31 extern doublereal slamch_(char *);
36 /* -- LAPACK routine (version 3.1) -- */
37 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
40 /* .. Scalar Arguments .. */
42 /* .. Array Arguments .. */
48 /* This subroutine computes the I-th updated eigenvalue of a symmetric */
49 /* rank-one modification to a diagonal matrix whose elements are */
50 /* given in the array d, and that */
52 /* D(i) < D(j) for i < j */
54 /* and that RHO > 0. This is arranged by the calling routine, and is */
55 /* no loss in generality. The rank-one modified system is thus */
57 /* diag( D ) + RHO * Z * Z_transpose. */
59 /* where we assume the Euclidean norm of Z is 1. */
61 /* The method consists of approximating the rational functions in the */
62 /* secular equation by simpler interpolating rational functions. */
67 /* N (input) INTEGER */
68 /* The length of all arrays. */
70 /* I (input) INTEGER */
71 /* The index of the eigenvalue to be computed. 1 <= I <= N. */
73 /* D (input) REAL array, dimension (N) */
74 /* The original eigenvalues. It is assumed that they are in */
75 /* order, D(I) < D(J) for I < J. */
77 /* Z (input) REAL array, dimension (N) */
78 /* The components of the updating vector. */
80 /* DELTA (output) REAL array, dimension (N) */
81 /* If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th */
82 /* component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5 */
83 /* for detail. The vector DELTA contains the information necessary */
84 /* to construct the eigenvectors by SLAED3 and SLAED9. */
86 /* RHO (input) REAL */
87 /* The scalar in the symmetric updating formula. */
89 /* DLAM (output) REAL */
90 /* The computed lambda_I, the I-th updated eigenvalue. */
92 /* INFO (output) INTEGER */
93 /* = 0: successful exit */
94 /* > 0: if INFO = 1, the updating process failed. */
96 /* Internal Parameters */
97 /* =================== */
99 /* Logical variable ORGATI (origin-at-i?) is used for distinguishing */
100 /* whether D(i) or D(i+1) is treated as the origin. */
102 /* ORGATI = .true. origin at i */
103 /* ORGATI = .false. origin at i+1 */
105 /* Logical variable SWTCH3 (switch-for-3-poles?) is for noting */
106 /* if we are working with THREE poles! */
108 /* MAXIT is the maximum number of iterations allowed for each */
111 /* Further Details */
112 /* =============== */
114 /* Based on contributions by */
115 /* Ren-Cang Li, Computer Science Division, University of California */
116 /* at Berkeley, USA */
118 /* ===================================================================== */
120 /* .. Parameters .. */
122 /* .. Local Scalars .. */
124 /* .. Local Arrays .. */
126 /* .. External Functions .. */
128 /* .. External Subroutines .. */
130 /* .. Intrinsic Functions .. */
132 /* .. Executable Statements .. */
134 /* Since this routine is called in an inner loop, we do no argument */
137 /* Quick return for N=1 and 2. */
139 /* Parameter adjustments */
148 /* Presumably, I=1 upon entry */
150 *dlam = d__[1] + *rho * z__[1] * z__[1];
155 slaed5_(i__, &d__[1], &z__[1], &delta[1], rho, dlam);
159 /* Compute machine epsilon */
161 eps = slamch_("Epsilon");
168 /* Initialize some basic variables */
173 /* Calculate initial guess */
177 /* If ||Z||_2 is not one, then TEMP should be set to */
178 /* RHO * ||Z||_2^2 / TWO */
181 for (j = 1; j <= i__1; ++j) {
182 delta[j] = d__[j] - d__[*i__] - midpt;
188 for (j = 1; j <= i__1; ++j) {
189 psi += z__[j] * z__[j] / delta[j];
194 w = c__ + z__[ii] * z__[ii] / delta[ii] + z__[*n] * z__[*n] / delta[*
198 temp = z__[*n - 1] * z__[*n - 1] / (d__[*n] - d__[*n - 1] + *rho)
199 + z__[*n] * z__[*n] / *rho;
203 del = d__[*n] - d__[*n - 1];
204 a = -c__ * del + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n]
206 b = z__[*n] * z__[*n] * del;
208 tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
210 tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
214 /* It can be proved that */
215 /* D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO */
220 del = d__[*n] - d__[*n - 1];
221 a = -c__ * del + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
222 b = z__[*n] * z__[*n] * del;
224 tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
226 tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
229 /* It can be proved that */
230 /* D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 */
237 for (j = 1; j <= i__1; ++j) {
238 delta[j] = d__[j] - d__[*i__] - tau;
242 /* Evaluate PSI and the derivative DPSI */
248 for (j = 1; j <= i__1; ++j) {
249 temp = z__[j] / delta[j];
250 psi += z__[j] * temp;
255 erretm = dabs(erretm);
257 /* Evaluate PHI and the derivative DPHI */
259 temp = z__[*n] / delta[*n];
260 phi = z__[*n] * temp;
262 erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
265 w = rhoinv + phi + psi;
267 /* Test for convergence */
269 if (dabs(w) <= eps * erretm) {
270 *dlam = d__[*i__] + tau;
275 dltlb = dmax(dltlb,tau);
277 dltub = dmin(dltub,tau);
280 /* Calculate the new step */
283 c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
284 a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] * (
286 b = delta[*n - 1] * delta[*n] * w;
292 /* ETA = RHO - TAU */
294 } else if (a >= 0.f) {
295 eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
298 eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
302 /* Note, eta should be positive if w is negative, and */
303 /* eta should be negative otherwise. However, */
304 /* if for some reason caused by roundoff, eta*w > 0, */
305 /* we simply use one Newton step instead. This way */
306 /* will guarantee eta*w < 0. */
309 eta = -w / (dpsi + dphi);
312 if (temp > dltub || temp < dltlb) {
314 eta = (dltub - tau) / 2.f;
316 eta = (dltlb - tau) / 2.f;
320 for (j = 1; j <= i__1; ++j) {
327 /* Evaluate PSI and the derivative DPSI */
333 for (j = 1; j <= i__1; ++j) {
334 temp = z__[j] / delta[j];
335 psi += z__[j] * temp;
340 erretm = dabs(erretm);
342 /* Evaluate PHI and the derivative DPHI */
344 temp = z__[*n] / delta[*n];
345 phi = z__[*n] * temp;
347 erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
350 w = rhoinv + phi + psi;
352 /* Main loop to update the values of the array DELTA */
356 for (niter = iter; niter <= 30; ++niter) {
358 /* Test for convergence */
360 if (dabs(w) <= eps * erretm) {
361 *dlam = d__[*i__] + tau;
366 dltlb = dmax(dltlb,tau);
368 dltub = dmin(dltub,tau);
371 /* Calculate the new step */
373 c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
374 a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] *
376 b = delta[*n - 1] * delta[*n] * w;
378 eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
381 eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
385 /* Note, eta should be positive if w is negative, and */
386 /* eta should be negative otherwise. However, */
387 /* if for some reason caused by roundoff, eta*w > 0, */
388 /* we simply use one Newton step instead. This way */
389 /* will guarantee eta*w < 0. */
392 eta = -w / (dpsi + dphi);
395 if (temp > dltub || temp < dltlb) {
397 eta = (dltub - tau) / 2.f;
399 eta = (dltlb - tau) / 2.f;
403 for (j = 1; j <= i__1; ++j) {
410 /* Evaluate PSI and the derivative DPSI */
416 for (j = 1; j <= i__1; ++j) {
417 temp = z__[j] / delta[j];
418 psi += z__[j] * temp;
423 erretm = dabs(erretm);
425 /* Evaluate PHI and the derivative DPHI */
427 temp = z__[*n] / delta[*n];
428 phi = z__[*n] * temp;
430 erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) *
433 w = rhoinv + phi + psi;
437 /* Return with INFO = 1, NITER = MAXIT and not converged */
440 *dlam = d__[*i__] + tau;
443 /* End for the case I = N */
447 /* The case for I < N */
452 /* Calculate initial guess */
454 del = d__[ip1] - d__[*i__];
457 for (j = 1; j <= i__1; ++j) {
458 delta[j] = d__[j] - d__[*i__] - midpt;
464 for (j = 1; j <= i__1; ++j) {
465 psi += z__[j] * z__[j] / delta[j];
471 for (j = *n; j >= i__1; --j) {
472 phi += z__[j] * z__[j] / delta[j];
475 c__ = rhoinv + psi + phi;
476 w = c__ + z__[*i__] * z__[*i__] / delta[*i__] + z__[ip1] * z__[ip1] /
481 /* d(i)< the ith eigenvalue < (d(i)+d(i+1))/2 */
483 /* We choose d(i) as origin. */
486 a = c__ * del + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
487 b = z__[*i__] * z__[*i__] * del;
489 tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
492 tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
499 /* (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1) */
501 /* We choose d(i+1) as origin. */
504 a = c__ * del - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
505 b = z__[ip1] * z__[ip1] * del;
507 tau = b * 2.f / (a - sqrt((r__1 = a * a + b * 4.f * c__, dabs(
510 tau = -(a + sqrt((r__1 = a * a + b * 4.f * c__, dabs(r__1))))
519 for (j = 1; j <= i__1; ++j) {
520 delta[j] = d__[j] - d__[*i__] - tau;
525 for (j = 1; j <= i__1; ++j) {
526 delta[j] = d__[j] - d__[ip1] - tau;
538 /* Evaluate PSI and the derivative DPSI */
544 for (j = 1; j <= i__1; ++j) {
545 temp = z__[j] / delta[j];
546 psi += z__[j] * temp;
551 erretm = dabs(erretm);
553 /* Evaluate PHI and the derivative DPHI */
558 for (j = *n; j >= i__1; --j) {
559 temp = z__[j] / delta[j];
560 phi += z__[j] * temp;
566 w = rhoinv + phi + psi;
568 /* W is the value of the secular function with */
569 /* its ii-th element removed. */
581 if (ii == 1 || ii == *n) {
585 temp = z__[ii] / delta[ii];
586 dw = dpsi + dphi + temp * temp;
587 temp = z__[ii] * temp;
589 erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f
592 /* Test for convergence */
594 if (dabs(w) <= eps * erretm) {
596 *dlam = d__[*i__] + tau;
598 *dlam = d__[ip1] + tau;
604 dltlb = dmax(dltlb,tau);
606 dltub = dmin(dltub,tau);
609 /* Calculate the new step */
614 /* Computing 2nd power */
615 r__1 = z__[*i__] / delta[*i__];
616 c__ = w - delta[ip1] * dw - (d__[*i__] - d__[ip1]) * (r__1 *
619 /* Computing 2nd power */
620 r__1 = z__[ip1] / delta[ip1];
621 c__ = w - delta[*i__] * dw - (d__[ip1] - d__[*i__]) * (r__1 *
624 a = (delta[*i__] + delta[ip1]) * w - delta[*i__] * delta[ip1] *
626 b = delta[*i__] * delta[ip1] * w;
630 a = z__[*i__] * z__[*i__] + delta[ip1] * delta[ip1] *
633 a = z__[ip1] * z__[ip1] + delta[*i__] * delta[*i__] *
638 } else if (a <= 0.f) {
639 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
642 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
647 /* Interpolation using THREE most relevant poles */
649 temp = rhoinv + psi + phi;
651 temp1 = z__[iim1] / delta[iim1];
653 c__ = temp - delta[iip1] * (dpsi + dphi) - (d__[iim1] - d__[
655 zz[0] = z__[iim1] * z__[iim1];
656 zz[2] = delta[iip1] * delta[iip1] * (dpsi - temp1 + dphi);
658 temp1 = z__[iip1] / delta[iip1];
660 c__ = temp - delta[iim1] * (dpsi + dphi) - (d__[iip1] - d__[
662 zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi - temp1));
663 zz[2] = z__[iip1] * z__[iip1];
665 zz[1] = z__[ii] * z__[ii];
666 slaed6_(&niter, &orgati, &c__, &delta[iim1], zz, &w, &eta, info);
672 /* Note, eta should be positive if w is negative, and */
673 /* eta should be negative otherwise. However, */
674 /* if for some reason caused by roundoff, eta*w > 0, */
675 /* we simply use one Newton step instead. This way */
676 /* will guarantee eta*w < 0. */
678 if (w * eta >= 0.f) {
682 if (temp > dltub || temp < dltlb) {
684 eta = (dltub - tau) / 2.f;
686 eta = (dltlb - tau) / 2.f;
693 for (j = 1; j <= i__1; ++j) {
698 /* Evaluate PSI and the derivative DPSI */
704 for (j = 1; j <= i__1; ++j) {
705 temp = z__[j] / delta[j];
706 psi += z__[j] * temp;
711 erretm = dabs(erretm);
713 /* Evaluate PHI and the derivative DPHI */
718 for (j = *n; j >= i__1; --j) {
719 temp = z__[j] / delta[j];
720 phi += z__[j] * temp;
726 temp = z__[ii] / delta[ii];
727 dw = dpsi + dphi + temp * temp;
728 temp = z__[ii] * temp;
729 w = rhoinv + phi + psi + temp;
730 erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f
731 + (r__1 = tau + eta, dabs(r__1)) * dw;
735 if (-w > dabs(prew) / 10.f) {
739 if (w > dabs(prew) / 10.f) {
746 /* Main loop to update the values of the array DELTA */
750 for (niter = iter; niter <= 30; ++niter) {
752 /* Test for convergence */
754 if (dabs(w) <= eps * erretm) {
756 *dlam = d__[*i__] + tau;
758 *dlam = d__[ip1] + tau;
764 dltlb = dmax(dltlb,tau);
766 dltub = dmin(dltub,tau);
769 /* Calculate the new step */
774 /* Computing 2nd power */
775 r__1 = z__[*i__] / delta[*i__];
776 c__ = w - delta[ip1] * dw - (d__[*i__] - d__[ip1]) * (
779 /* Computing 2nd power */
780 r__1 = z__[ip1] / delta[ip1];
781 c__ = w - delta[*i__] * dw - (d__[ip1] - d__[*i__]) *
785 temp = z__[ii] / delta[ii];
791 c__ = w - delta[*i__] * dpsi - delta[ip1] * dphi;
793 a = (delta[*i__] + delta[ip1]) * w - delta[*i__] * delta[ip1]
795 b = delta[*i__] * delta[ip1] * w;
800 a = z__[*i__] * z__[*i__] + delta[ip1] *
801 delta[ip1] * (dpsi + dphi);
803 a = z__[ip1] * z__[ip1] + delta[*i__] * delta[
804 *i__] * (dpsi + dphi);
807 a = delta[*i__] * delta[*i__] * dpsi + delta[ip1]
812 } else if (a <= 0.f) {
813 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1))
816 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__,
821 /* Interpolation using THREE most relevant poles */
823 temp = rhoinv + psi + phi;
825 c__ = temp - delta[iim1] * dpsi - delta[iip1] * dphi;
826 zz[0] = delta[iim1] * delta[iim1] * dpsi;
827 zz[2] = delta[iip1] * delta[iip1] * dphi;
830 temp1 = z__[iim1] / delta[iim1];
832 c__ = temp - delta[iip1] * (dpsi + dphi) - (d__[iim1]
833 - d__[iip1]) * temp1;
834 zz[0] = z__[iim1] * z__[iim1];
835 zz[2] = delta[iip1] * delta[iip1] * (dpsi - temp1 +
838 temp1 = z__[iip1] / delta[iip1];
840 c__ = temp - delta[iim1] * (dpsi + dphi) - (d__[iip1]
841 - d__[iim1]) * temp1;
842 zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi -
844 zz[2] = z__[iip1] * z__[iip1];
847 slaed6_(&niter, &orgati, &c__, &delta[iim1], zz, &w, &eta,
854 /* Note, eta should be positive if w is negative, and */
855 /* eta should be negative otherwise. However, */
856 /* if for some reason caused by roundoff, eta*w > 0, */
857 /* we simply use one Newton step instead. This way */
858 /* will guarantee eta*w < 0. */
860 if (w * eta >= 0.f) {
864 if (temp > dltub || temp < dltlb) {
866 eta = (dltub - tau) / 2.f;
868 eta = (dltlb - tau) / 2.f;
873 for (j = 1; j <= i__1; ++j) {
881 /* Evaluate PSI and the derivative DPSI */
887 for (j = 1; j <= i__1; ++j) {
888 temp = z__[j] / delta[j];
889 psi += z__[j] * temp;
894 erretm = dabs(erretm);
896 /* Evaluate PHI and the derivative DPHI */
901 for (j = *n; j >= i__1; --j) {
902 temp = z__[j] / delta[j];
903 phi += z__[j] * temp;
909 temp = z__[ii] / delta[ii];
910 dw = dpsi + dphi + temp * temp;
911 temp = z__[ii] * temp;
912 w = rhoinv + phi + psi + temp;
913 erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) *
914 3.f + dabs(tau) * dw;
915 if (w * prew > 0.f && dabs(w) > dabs(prew) / 10.f) {
922 /* Return with INFO = 1, NITER = MAXIT and not converged */
926 *dlam = d__[*i__] + tau;
928 *dlam = d__[ip1] + tau;