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

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