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

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