git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlaed4.c

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