3rdparty/lapack/slaed4.c

   1 /* slaed4.f -- translated by f2c (version 20061008).
   2    You must link the resulting object file with libf2c:
   3         on Microsoft Windows system, link with libf2c.lib;
   4         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
   5         or, if you install libf2c.a in a standard place, with -lf2c -lm
   6         -- in that order, at the end of the command line, as in
   7                 cc *.o -lf2c -lm
   8         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
   9
  10                 http://www.netlib.org/f2c/libf2c.zip
  11 */
  12
  13 #include "clapack.h"
  14
  15
  16 /* Subroutine */ int slaed4_(integer *n, integer *i__, real *d__, real *z__,
  17         real *delta, real *rho, real *dlam, integer *info)
  18 {
  19     /* System generated locals */
  20     integer i__1;
  21     real r__1;
  22
  23     /* Builtin functions */
  24     double sqrt(doublereal);
  25
  26     /* Local variables */
  27     real a, b, c__;
  28     integer j;
  29     real w;
  30     integer ii;
  31     real dw, zz[3];
  32     integer ip1;
  33     real del, eta, phi, eps, tau, psi;
  34     integer iim1, iip1;
  35     real dphi, dpsi;
  36     integer iter;
  37     real temp, prew, temp1, dltlb, dltub, midpt;
  38     integer niter;
  39     logical swtch;
  40     extern /* Subroutine */ int slaed5_(integer *, real *, real *, real *,
  41             real *, real *), slaed6_(integer *, logical *, real *, real *,
  42             real *, real *, real *, integer *);
  43     logical swtch3;
  44     extern doublereal slamch_(char *);
  45     logical orgati;
  46     real erretm, rhoinv;
  47
  48
  49 /*  -- LAPACK routine (version 3.2) -- */
  50 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  51 /*     November 2006 */
  52
  53 /*     .. Scalar Arguments .. */
  54 /*     .. */
  55 /*     .. Array Arguments .. */
  56 /*     .. */
  57
  58 /*  Purpose */
  59 /*  ======= */
  60
  61 /*  This subroutine computes the I-th updated eigenvalue of a symmetric */
  62 /*  rank-one modification to a diagonal matrix whose elements are */
  63 /*  given in the array d, and that */
  64
  65 /*             D(i) < D(j)  for  i < j */
  66
  67 /*  and that RHO > 0.  This is arranged by the calling routine, and is */
  68 /*  no loss in generality.  The rank-one modified system is thus */
  69
  70 /*             diag( D )  +  RHO *  Z * Z_transpose. */
  71
  72 /*  where we assume the Euclidean norm of Z is 1. */
  73
  74 /*  The method consists of approximating the rational functions in the */
  75 /*  secular equation by simpler interpolating rational functions. */
  76
  77 /*  Arguments */
  78 /*  ========= */
  79
  80 /*  N      (input) INTEGER */
  81 /*         The length of all arrays. */
  82
  83 /*  I      (input) INTEGER */
  84 /*         The index of the eigenvalue to be computed.  1 <= I <= N. */
  85
  86 /*  D      (input) REAL array, dimension (N) */
  87 /*         The original eigenvalues.  It is assumed that they are in */
  88 /*         order, D(I) < D(J)  for I < J. */
  89
  90 /*  Z      (input) REAL array, dimension (N) */
  91 /*         The components of the updating vector. */
  92
  93 /*  DELTA  (output) REAL array, dimension (N) */
  94 /*         If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th */
  95 /*         component.  If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5 */
  96 /*         for detail. The vector DELTA contains the information necessary */
  97 /*         to construct the eigenvectors by SLAED3 and SLAED9. */
  98
  99 /*  RHO    (input) REAL */
 100 /*         The scalar in the symmetric updating formula. */
 101
 102 /*  DLAM   (output) REAL */
 103 /*         The computed lambda_I, the I-th updated eigenvalue. */
 104
 105 /*  INFO   (output) INTEGER */
 106 /*         = 0:  successful exit */
 107 /*         > 0:  if INFO = 1, the updating process failed. */
 108
 109 /*  Internal Parameters */
 110 /*  =================== */
 111
 112 /*  Logical variable ORGATI (origin-at-i?) is used for distinguishing */
 113 /*  whether D(i) or D(i+1) is treated as the origin. */
 114
 115 /*            ORGATI = .true.    origin at i */
 116 /*            ORGATI = .false.   origin at i+1 */
 117
 118 /*   Logical variable SWTCH3 (switch-for-3-poles?) is for noting */
 119 /*   if we are working with THREE poles! */
 120
 121 /*   MAXIT is the maximum number of iterations allowed for each */
 122 /*   eigenvalue. */
 123
 124 /*  Further Details */
 125 /*  =============== */
 126
 127 /*  Based on contributions by */
 128 /*     Ren-Cang Li, Computer Science Division, University of California */
 129 /*     at Berkeley, USA */
 130
 131 /*  ===================================================================== */
 132
 133 /*     .. Parameters .. */
 134 /*     .. */
 135 /*     .. Local Scalars .. */
 136 /*     .. */
 137 /*     .. Local Arrays .. */
 138 /*     .. */
 139 /*     .. External Functions .. */
 140 /*     .. */
 141 /*     .. External Subroutines .. */
 142 /*     .. */
 143 /*     .. Intrinsic Functions .. */
 144 /*     .. */
 145 /*     .. Executable Statements .. */
 146
 147 /*     Since this routine is called in an inner loop, we do no argument */
 148 /*     checking. */
 149
 150 /*     Quick return for N=1 and 2. */
 151
 152     /* Parameter adjustments */
 153     --delta;
 154     --z__;
 155     --d__;
 156
 157     /* Function Body */
 158     *info = 0;
 159     if (*n == 1) {
 160
 161 /*         Presumably, I=1 upon entry */
 162
 163         *dlam = d__[1] + *rho * z__[1] * z__[1];
 164         delta[1] = 1.f;
 165         return 0;
 166     }
 167     if (*n == 2) {
 168         slaed5_(i__, &d__[1], &z__[1], &delta[1], rho, dlam);
 169         return 0;
 170     }
 171
 172 /*     Compute machine epsilon */
 173
 174     eps = slamch_("Epsilon");
 175     rhoinv = 1.f / *rho;
 176
 177 /*     The case I = N */
 178
 179     if (*i__ == *n) {
 180
 181 /*        Initialize some basic variables */
 182
 183         ii = *n - 1;
 184         niter = 1;
 185
 186 /*        Calculate initial guess */
 187
 188         midpt = *rho / 2.f;
 189
 190 /*        If ||Z||_2 is not one, then TEMP should be set to */
 191 /*        RHO * ||Z||_2^2 / TWO */
 192
 193         i__1 = *n;
 194         for (j = 1; j <= i__1; ++j) {
 195             delta[j] = d__[j] - d__[*i__] - midpt;
 196 /* L10: */
 197         }
 198
 199         psi = 0.f;
 200         i__1 = *n - 2;
 201         for (j = 1; j <= i__1; ++j) {
 202             psi += z__[j] * z__[j] / delta[j];
 203 /* L20: */
 204         }
 205
 206         c__ = rhoinv + psi;
 207         w = c__ + z__[ii] * z__[ii] / delta[ii] + z__[*n] * z__[*n] / delta[*
 208                 n];
 209
 210         if (w <= 0.f) {
 211             temp = z__[*n - 1] * z__[*n - 1] / (d__[*n] - d__[*n - 1] + *rho)
 212                     + z__[*n] * z__[*n] / *rho;
 213             if (c__ <= temp) {
 214                 tau = *rho;
 215             } else {
 216                 del = d__[*n] - d__[*n - 1];
 217                 a = -c__ * del + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n]
 218                         ;
 219                 b = z__[*n] * z__[*n] * del;
 220                 if (a < 0.f) {
 221                     tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
 222                 } else {
 223                     tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
 224                 }
 225             }
 226
 227 /*           It can be proved that */
 228 /*               D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO */
 229
 230             dltlb = midpt;
 231             dltub = *rho;
 232         } else {
 233             del = d__[*n] - d__[*n - 1];
 234             a = -c__ * del + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
 235             b = z__[*n] * z__[*n] * del;
 236             if (a < 0.f) {
 237                 tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
 238             } else {
 239                 tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
 240             }
 241
 242 /*           It can be proved that */
 243 /*               D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 */
 244
 245             dltlb = 0.f;
 246             dltub = midpt;
 247         }
 248
 249         i__1 = *n;
 250         for (j = 1; j <= i__1; ++j) {
 251             delta[j] = d__[j] - d__[*i__] - tau;
 252 /* L30: */
 253         }
 254
 255 /*        Evaluate PSI and the derivative DPSI */
 256
 257         dpsi = 0.f;
 258         psi = 0.f;
 259         erretm = 0.f;
 260         i__1 = ii;
 261         for (j = 1; j <= i__1; ++j) {
 262             temp = z__[j] / delta[j];
 263             psi += z__[j] * temp;
 264             dpsi += temp * temp;
 265             erretm += psi;
 266 /* L40: */
 267         }
 268         erretm = dabs(erretm);
 269
 270 /*        Evaluate PHI and the derivative DPHI */
 271
 272         temp = z__[*n] / delta[*n];
 273         phi = z__[*n] * temp;
 274         dphi = temp * temp;
 275         erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
 276                 dpsi + dphi);
 277
 278         w = rhoinv + phi + psi;
 279
 280 /*        Test for convergence */
 281
 282         if (dabs(w) <= eps * erretm) {
 283             *dlam = d__[*i__] + tau;
 284             goto L250;
 285         }
 286
 287         if (w <= 0.f) {
 288             dltlb = dmax(dltlb,tau);
 289         } else {
 290             dltub = dmin(dltub,tau);
 291         }
 292
 293 /*        Calculate the new step */
 294
 295         ++niter;
 296         c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
 297         a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] * (
 298                 dpsi + dphi);
 299         b = delta[*n - 1] * delta[*n] * w;
 300         if (c__ < 0.f) {
 301             c__ = dabs(c__);
 302         }
 303         if (c__ == 0.f) {
 304 /*          ETA = B/A */
 305 /*           ETA = RHO - TAU */
 306             eta = dltub - tau;
 307         } else if (a >= 0.f) {
 308             eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
 309                     c__ * 2.f);
 310         } else {
 311             eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
 312                     r__1))));
 313         }
 314
 315 /*        Note, eta should be positive if w is negative, and */
 316 /*        eta should be negative otherwise. However, */
 317 /*        if for some reason caused by roundoff, eta*w > 0, */
 318 /*        we simply use one Newton step instead. This way */
 319 /*        will guarantee eta*w < 0. */
 320
 321         if (w * eta > 0.f) {
 322             eta = -w / (dpsi + dphi);
 323         }
 324         temp = tau + eta;
 325         if (temp > dltub || temp < dltlb) {
 326             if (w < 0.f) {
 327                 eta = (dltub - tau) / 2.f;
 328             } else {
 329                 eta = (dltlb - tau) / 2.f;
 330             }
 331         }
 332         i__1 = *n;
 333         for (j = 1; j <= i__1; ++j) {
 334             delta[j] -= eta;
 335 /* L50: */
 336         }
 337
 338         tau += eta;
 339
 340 /*        Evaluate PSI and the derivative DPSI */
 341
 342         dpsi = 0.f;
 343         psi = 0.f;
 344         erretm = 0.f;
 345         i__1 = ii;
 346         for (j = 1; j <= i__1; ++j) {
 347             temp = z__[j] / delta[j];
 348             psi += z__[j] * temp;
 349             dpsi += temp * temp;
 350             erretm += psi;
 351 /* L60: */
 352         }
 353         erretm = dabs(erretm);
 354
 355 /*        Evaluate PHI and the derivative DPHI */
 356
 357         temp = z__[*n] / delta[*n];
 358         phi = z__[*n] * temp;
 359         dphi = temp * temp;
 360         erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
 361                 dpsi + dphi);
 362
 363         w = rhoinv + phi + psi;
 364
 365 /*        Main loop to update the values of the array   DELTA */
 366
 367         iter = niter + 1;
 368
 369         for (niter = iter; niter <= 30; ++niter) {
 370
 371 /*           Test for convergence */
 372
 373             if (dabs(w) <= eps * erretm) {
 374                 *dlam = d__[*i__] + tau;
 375                 goto L250;
 376             }
 377
 378             if (w <= 0.f) {
 379                 dltlb = dmax(dltlb,tau);
 380             } else {
 381                 dltub = dmin(dltub,tau);
 382             }
 383
 384 /*           Calculate the new step */
 385
 386             c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
 387             a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] *
 388                     (dpsi + dphi);
 389             b = delta[*n - 1] * delta[*n] * w;
 390             if (a >= 0.f) {
 391                 eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
 392                          (c__ * 2.f);
 393             } else {
 394                 eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
 395                         r__1))));
 396             }
 397
 398 /*           Note, eta should be positive if w is negative, and */
 399 /*           eta should be negative otherwise. However, */
 400 /*           if for some reason caused by roundoff, eta*w > 0, */
 401 /*           we simply use one Newton step instead. This way */
 402 /*           will guarantee eta*w < 0. */
 403
 404             if (w * eta > 0.f) {
 405                 eta = -w / (dpsi + dphi);
 406             }
 407             temp = tau + eta;
 408             if (temp > dltub || temp < dltlb) {
 409                 if (w < 0.f) {
 410                     eta = (dltub - tau) / 2.f;
 411                 } else {
 412                     eta = (dltlb - tau) / 2.f;
 413                 }
 414             }
 415             i__1 = *n;
 416             for (j = 1; j <= i__1; ++j) {
 417                 delta[j] -= eta;
 418 /* L70: */
 419             }
 420
 421             tau += eta;
 422
 423 /*           Evaluate PSI and the derivative DPSI */
 424
 425             dpsi = 0.f;
 426             psi = 0.f;
 427             erretm = 0.f;
 428             i__1 = ii;
 429             for (j = 1; j <= i__1; ++j) {
 430                 temp = z__[j] / delta[j];
 431                 psi += z__[j] * temp;
 432                 dpsi += temp * temp;
 433                 erretm += psi;
 434 /* L80: */
 435             }
 436             erretm = dabs(erretm);
 437
 438 /*           Evaluate PHI and the derivative DPHI */
 439
 440             temp = z__[*n] / delta[*n];
 441             phi = z__[*n] * temp;
 442             dphi = temp * temp;
 443             erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) *
 444                     (dpsi + dphi);
 445
 446             w = rhoinv + phi + psi;
 447 /* L90: */
 448         }
 449
 450 /*        Return with INFO = 1, NITER = MAXIT and not converged */
 451
 452         *info = 1;
 453         *dlam = d__[*i__] + tau;
 454         goto L250;
 455
 456 /*        End for the case I = N */
 457
 458     } else {
 459
 460 /*        The case for I < N */
 461
 462         niter = 1;
 463         ip1 = *i__ + 1;
 464
 465 /*        Calculate initial guess */
 466
 467         del = d__[ip1] - d__[*i__];
 468         midpt = del / 2.f;
 469         i__1 = *n;
 470         for (j = 1; j <= i__1; ++j) {
 471             delta[j] = d__[j] - d__[*i__] - midpt;
 472 /* L100: */
 473         }
 474
 475         psi = 0.f;
 476         i__1 = *i__ - 1;
 477         for (j = 1; j <= i__1; ++j) {
 478             psi += z__[j] * z__[j] / delta[j];
 479 /* L110: */
 480         }
 481
 482         phi = 0.f;
 483         i__1 = *i__ + 2;
 484         for (j = *n; j >= i__1; --j) {
 485             phi += z__[j] * z__[j] / delta[j];
 486 /* L120: */
 487         }
 488         c__ = rhoinv + psi + phi;
 489         w = c__ + z__[*i__] * z__[*i__] / delta[*i__] + z__[ip1] * z__[ip1] /
 490                 delta[ip1];
 491
 492         if (w > 0.f) {
 493
 494 /*           d(i)< the ith eigenvalue < (d(i)+d(i+1))/2 */
 495
 496 /*           We choose d(i) as origin. */
 497
 498             orgati = TRUE_;
 499             a = c__ * del + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
 500             b = z__[*i__] * z__[*i__] * del;
 501             if (a > 0.f) {
 502                 tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
 503                         r__1))));
 504             } else {
 505                 tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
 506                          (c__ * 2.f);
 507             }
 508             dltlb = 0.f;
 509             dltub = midpt;
 510         } else {
 511
 512 /*           (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1) */
 513
 514 /*           We choose d(i+1) as origin. */
 515
 516             orgati = FALSE_;
 517             a = c__ * del - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
 518             b = z__[ip1] * z__[ip1] * del;
 519             if (a < 0.f) {
 520                 tau = b * 2.f / (a - sqrt((r__1 = a * a + b * 4.f * c__, dabs(
 521                         r__1))));
 522             } else {
 523                 tau = -(a + sqrt((r__1 = a * a + b * 4.f * c__, dabs(r__1))))
 524                         / (c__ * 2.f);
 525             }
 526             dltlb = -midpt;
 527             dltub = 0.f;
 528         }
 529
 530         if (orgati) {
 531             i__1 = *n;
 532             for (j = 1; j <= i__1; ++j) {
 533                 delta[j] = d__[j] - d__[*i__] - tau;
 534 /* L130: */
 535             }
 536         } else {
 537             i__1 = *n;
 538             for (j = 1; j <= i__1; ++j) {
 539                 delta[j] = d__[j] - d__[ip1] - tau;
 540 /* L140: */
 541             }
 542         }
 543         if (orgati) {
 544             ii = *i__;
 545         } else {
 546             ii = *i__ + 1;
 547         }
 548         iim1 = ii - 1;
 549         iip1 = ii + 1;
 550
 551 /*        Evaluate PSI and the derivative DPSI */
 552
 553         dpsi = 0.f;
 554         psi = 0.f;
 555         erretm = 0.f;
 556         i__1 = iim1;
 557         for (j = 1; j <= i__1; ++j) {
 558             temp = z__[j] / delta[j];
 559             psi += z__[j] * temp;
 560             dpsi += temp * temp;
 561             erretm += psi;
 562 /* L150: */
 563         }
 564         erretm = dabs(erretm);
 565
 566 /*        Evaluate PHI and the derivative DPHI */
 567
 568         dphi = 0.f;
 569         phi = 0.f;
 570         i__1 = iip1;
 571         for (j = *n; j >= i__1; --j) {
 572             temp = z__[j] / delta[j];
 573             phi += z__[j] * temp;
 574             dphi += temp * temp;
 575             erretm += phi;
 576 /* L160: */
 577         }
 578
 579         w = rhoinv + phi + psi;
 580
 581 /*        W is the value of the secular function with */
 582 /*        its ii-th element removed. */
 583
 584         swtch3 = FALSE_;
 585         if (orgati) {
 586             if (w < 0.f) {
 587                 swtch3 = TRUE_;
 588             }
 589         } else {
 590             if (w > 0.f) {
 591                 swtch3 = TRUE_;
 592             }
 593         }
 594         if (ii == 1 || ii == *n) {
 595             swtch3 = FALSE_;
 596         }
 597
 598         temp = z__[ii] / delta[ii];
 599         dw = dpsi + dphi + temp * temp;
 600         temp = z__[ii] * temp;
 601         w += temp;
 602         erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f
 603                 + dabs(tau) * dw;
 604
 605 /*        Test for convergence */
 606
 607         if (dabs(w) <= eps * erretm) {
 608             if (orgati) {
 609                 *dlam = d__[*i__] + tau;
 610             } else {
 611                 *dlam = d__[ip1] + tau;
 612             }
 613             goto L250;
 614         }
 615
 616         if (w <= 0.f) {
 617             dltlb = dmax(dltlb,tau);
 618         } else {
 619             dltub = dmin(dltub,tau);
 620         }
 621
 622 /*        Calculate the new step */
 623
 624         ++niter;
 625         if (! swtch3) {
 626             if (orgati) {
 627 /* Computing 2nd power */
 628                 r__1 = z__[*i__] / delta[*i__];
 629                 c__ = w - delta[ip1] * dw - (d__[*i__] - d__[ip1]) * (r__1 *
 630                         r__1);
 631             } else {
 632 /* Computing 2nd power */
 633                 r__1 = z__[ip1] / delta[ip1];
 634                 c__ = w - delta[*i__] * dw - (d__[ip1] - d__[*i__]) * (r__1 *
 635                         r__1);
 636             }
 637             a = (delta[*i__] + delta[ip1]) * w - delta[*i__] * delta[ip1] *
 638                     dw;
 639             b = delta[*i__] * delta[ip1] * w;
 640             if (c__ == 0.f) {
 641                 if (a == 0.f) {
 642                     if (orgati) {
 643                         a = z__[*i__] * z__[*i__] + delta[ip1] * delta[ip1] *
 644                                 (dpsi + dphi);
 645                     } else {
 646                         a = z__[ip1] * z__[ip1] + delta[*i__] * delta[*i__] *
 647                                 (dpsi + dphi);
 648                     }
 649                 }
 650                 eta = b / a;
 651             } else if (a <= 0.f) {
 652                 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
 653                          (c__ * 2.f);
 654             } else {
 655                 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
 656                         r__1))));
 657             }
 658         } else {
 659
 660 /*           Interpolation using THREE most relevant poles */
 661
 662             temp = rhoinv + psi + phi;
 663             if (orgati) {
 664                 temp1 = z__[iim1] / delta[iim1];
 665                 temp1 *= temp1;
 666                 c__ = temp - delta[iip1] * (dpsi + dphi) - (d__[iim1] - d__[
 667                         iip1]) * temp1;
 668                 zz[0] = z__[iim1] * z__[iim1];
 669                 zz[2] = delta[iip1] * delta[iip1] * (dpsi - temp1 + dphi);
 670             } else {
 671                 temp1 = z__[iip1] / delta[iip1];
 672                 temp1 *= temp1;
 673                 c__ = temp - delta[iim1] * (dpsi + dphi) - (d__[iip1] - d__[
 674                         iim1]) * temp1;
 675                 zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi - temp1));
 676                 zz[2] = z__[iip1] * z__[iip1];
 677             }
 678             zz[1] = z__[ii] * z__[ii];
 679             slaed6_(&niter, &orgati, &c__, &delta[iim1], zz, &w, &eta, info);
 680             if (*info != 0) {
 681                 goto L250;
 682             }
 683         }
 684
 685 /*        Note, eta should be positive if w is negative, and */
 686 /*        eta should be negative otherwise. However, */
 687 /*        if for some reason caused by roundoff, eta*w > 0, */
 688 /*        we simply use one Newton step instead. This way */
 689 /*        will guarantee eta*w < 0. */
 690
 691         if (w * eta >= 0.f) {
 692             eta = -w / dw;
 693         }
 694         temp = tau + eta;
 695         if (temp > dltub || temp < dltlb) {
 696             if (w < 0.f) {
 697                 eta = (dltub - tau) / 2.f;
 698             } else {
 699                 eta = (dltlb - tau) / 2.f;
 700             }
 701         }
 702
 703         prew = w;
 704
 705         i__1 = *n;
 706         for (j = 1; j <= i__1; ++j) {
 707             delta[j] -= eta;
 708 /* L180: */
 709         }
 710
 711 /*        Evaluate PSI and the derivative DPSI */
 712
 713         dpsi = 0.f;
 714         psi = 0.f;
 715         erretm = 0.f;
 716         i__1 = iim1;
 717         for (j = 1; j <= i__1; ++j) {
 718             temp = z__[j] / delta[j];
 719             psi += z__[j] * temp;
 720             dpsi += temp * temp;
 721             erretm += psi;
 722 /* L190: */
 723         }
 724         erretm = dabs(erretm);
 725
 726 /*        Evaluate PHI and the derivative DPHI */
 727
 728         dphi = 0.f;
 729         phi = 0.f;
 730         i__1 = iip1;
 731         for (j = *n; j >= i__1; --j) {
 732             temp = z__[j] / delta[j];
 733             phi += z__[j] * temp;
 734             dphi += temp * temp;
 735             erretm += phi;
 736 /* L200: */
 737         }
 738
 739         temp = z__[ii] / delta[ii];
 740         dw = dpsi + dphi + temp * temp;
 741         temp = z__[ii] * temp;
 742         w = rhoinv + phi + psi + temp;
 743         erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f
 744                 + (r__1 = tau + eta, dabs(r__1)) * dw;
 745
 746         swtch = FALSE_;
 747         if (orgati) {
 748             if (-w > dabs(prew) / 10.f) {
 749                 swtch = TRUE_;
 750             }
 751         } else {
 752             if (w > dabs(prew) / 10.f) {
 753                 swtch = TRUE_;
 754             }
 755         }
 756
 757         tau += eta;
 758
 759 /*        Main loop to update the values of the array   DELTA */
 760
 761         iter = niter + 1;
 762
 763         for (niter = iter; niter <= 30; ++niter) {
 764
 765 /*           Test for convergence */
 766
 767             if (dabs(w) <= eps * erretm) {
 768                 if (orgati) {
 769                     *dlam = d__[*i__] + tau;
 770                 } else {
 771                     *dlam = d__[ip1] + tau;
 772                 }
 773                 goto L250;
 774             }
 775
 776             if (w <= 0.f) {
 777                 dltlb = dmax(dltlb,tau);
 778             } else {
 779                 dltub = dmin(dltub,tau);
 780             }
 781
 782 /*           Calculate the new step */
 783
 784             if (! swtch3) {
 785                 if (! swtch) {
 786                     if (orgati) {
 787 /* Computing 2nd power */
 788                         r__1 = z__[*i__] / delta[*i__];
 789                         c__ = w - delta[ip1] * dw - (d__[*i__] - d__[ip1]) * (
 790                                 r__1 * r__1);
 791                     } else {
 792 /* Computing 2nd power */
 793                         r__1 = z__[ip1] / delta[ip1];
 794                         c__ = w - delta[*i__] * dw - (d__[ip1] - d__[*i__]) *
 795                                 (r__1 * r__1);
 796                     }
 797                 } else {
 798                     temp = z__[ii] / delta[ii];
 799                     if (orgati) {
 800                         dpsi += temp * temp;
 801                     } else {
 802                         dphi += temp * temp;
 803                     }
 804                     c__ = w - delta[*i__] * dpsi - delta[ip1] * dphi;
 805                 }
 806                 a = (delta[*i__] + delta[ip1]) * w - delta[*i__] * delta[ip1]
 807                         * dw;
 808                 b = delta[*i__] * delta[ip1] * w;
 809                 if (c__ == 0.f) {
 810                     if (a == 0.f) {
 811                         if (! swtch) {
 812                             if (orgati) {
 813                                 a = z__[*i__] * z__[*i__] + delta[ip1] *
 814                                         delta[ip1] * (dpsi + dphi);
 815                             } else {
 816                                 a = z__[ip1] * z__[ip1] + delta[*i__] * delta[
 817                                         *i__] * (dpsi + dphi);
 818                             }
 819                         } else {
 820                             a = delta[*i__] * delta[*i__] * dpsi + delta[ip1]
 821                                     * delta[ip1] * dphi;
 822                         }
 823                     }
 824                     eta = b / a;
 825                 } else if (a <= 0.f) {
 826                     eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1))
 827                             )) / (c__ * 2.f);
 828                 } else {
 829                     eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__,
 830                             dabs(r__1))));
 831                 }
 832             } else {
 833
 834 /*              Interpolation using THREE most relevant poles */
 835
 836                 temp = rhoinv + psi + phi;
 837                 if (swtch) {
 838                     c__ = temp - delta[iim1] * dpsi - delta[iip1] * dphi;
 839                     zz[0] = delta[iim1] * delta[iim1] * dpsi;
 840                     zz[2] = delta[iip1] * delta[iip1] * dphi;
 841                 } else {
 842                     if (orgati) {
 843                         temp1 = z__[iim1] / delta[iim1];
 844                         temp1 *= temp1;
 845                         c__ = temp - delta[iip1] * (dpsi + dphi) - (d__[iim1]
 846                                 - d__[iip1]) * temp1;
 847                         zz[0] = z__[iim1] * z__[iim1];
 848                         zz[2] = delta[iip1] * delta[iip1] * (dpsi - temp1 +
 849                                 dphi);
 850                     } else {
 851                         temp1 = z__[iip1] / delta[iip1];
 852                         temp1 *= temp1;
 853                         c__ = temp - delta[iim1] * (dpsi + dphi) - (d__[iip1]
 854                                 - d__[iim1]) * temp1;
 855                         zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi -
 856                                 temp1));
 857                         zz[2] = z__[iip1] * z__[iip1];
 858                     }
 859                 }
 860                 slaed6_(&niter, &orgati, &c__, &delta[iim1], zz, &w, &eta,
 861                         info);
 862                 if (*info != 0) {
 863                     goto L250;
 864                 }
 865             }
 866
 867 /*           Note, eta should be positive if w is negative, and */
 868 /*           eta should be negative otherwise. However, */
 869 /*           if for some reason caused by roundoff, eta*w > 0, */
 870 /*           we simply use one Newton step instead. This way */
 871 /*           will guarantee eta*w < 0. */
 872
 873             if (w * eta >= 0.f) {
 874                 eta = -w / dw;
 875             }
 876             temp = tau + eta;
 877             if (temp > dltub || temp < dltlb) {
 878                 if (w < 0.f) {
 879                     eta = (dltub - tau) / 2.f;
 880                 } else {
 881                     eta = (dltlb - tau) / 2.f;
 882                 }
 883             }
 884
 885             i__1 = *n;
 886             for (j = 1; j <= i__1; ++j) {
 887                 delta[j] -= eta;
 888 /* L210: */
 889             }
 890
 891             tau += eta;
 892             prew = w;
 893
 894 /*           Evaluate PSI and the derivative DPSI */
 895
 896             dpsi = 0.f;
 897             psi = 0.f;
 898             erretm = 0.f;
 899             i__1 = iim1;
 900             for (j = 1; j <= i__1; ++j) {
 901                 temp = z__[j] / delta[j];
 902                 psi += z__[j] * temp;
 903                 dpsi += temp * temp;
 904                 erretm += psi;
 905 /* L220: */
 906             }
 907             erretm = dabs(erretm);
 908
 909 /*           Evaluate PHI and the derivative DPHI */
 910
 911             dphi = 0.f;
 912             phi = 0.f;
 913             i__1 = iip1;
 914             for (j = *n; j >= i__1; --j) {
 915                 temp = z__[j] / delta[j];
 916                 phi += z__[j] * temp;
 917                 dphi += temp * temp;
 918                 erretm += phi;
 919 /* L230: */
 920             }
 921
 922             temp = z__[ii] / delta[ii];
 923             dw = dpsi + dphi + temp * temp;
 924             temp = z__[ii] * temp;
 925             w = rhoinv + phi + psi + temp;
 926             erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) *
 927                     3.f + dabs(tau) * dw;
 928             if (w * prew > 0.f && dabs(w) > dabs(prew) / 10.f) {
 929                 swtch = ! swtch;
 930             }
 931
 932 /* L240: */
 933         }
 934
 935 /*        Return with INFO = 1, NITER = MAXIT and not converged */
 936
 937         *info = 1;
 938         if (orgati) {
 939             *dlam = d__[*i__] + tau;
 940         } else {
 941             *dlam = d__[ip1] + tau;
 942         }
 943
 944     }
 945
 946 L250:
 947
 948     return 0;
 949
 950 /*     End of SLAED4 */
 951
 952 } /* slaed4_ */