src/gromacs/linearalgebra/gmx_lapack/dlasd4.c

   1 #include <math.h>
   2 #include "../gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 #include "types/simple.h"
   6
   7 void
   8 F77_FUNC(dlasd4,DLASD4)(int *n,
   9         int *i__,
  10         double *d__,
  11         double *z__,
  12         double *delta,
  13         double *rho,
  14         double *sigma,
  15         double *work,
  16         int *info)
  17 {
  18     int i__1;
  19     double d__1;
  20
  21     double a, b, c__;
  22     int j;
  23     double w, dd[3];
  24     int ii;
  25     double dw, zz[3];
  26     int ip1;
  27     double eta, phi, eps, tau, psi;
  28     int iim1, iip1;
  29     double dphi, dpsi;
  30     int iter;
  31     double temp, prew, sg2lb, sg2ub, temp1, temp2, dtiim, delsq,
  32             dtiip;
  33     int niter;
  34     double dtisq;
  35     int swtch;
  36     double dtnsq;
  37     double delsq2, dtnsq1;
  38     int swtch3;
  39     int orgati;
  40     double erretm, dtipsq, rhoinv;
  41
  42     --work;
  43     --delta;
  44     --z__;
  45     --d__;
  46
  47     *info = 0;
  48     if (*n == 1) {
  49
  50         *sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
  51         delta[1] = 1.;
  52         work[1] = 1.;
  53         return;
  54     }
  55     if (*n == 2) {
  56         F77_FUNC(dlasd5,DLASD5)(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
  57         return;
  58     }
  59
  60     eps = GMX_DOUBLE_EPS;
  61     rhoinv = 1. / *rho;
  62
  63     if (*i__ == *n) {
  64
  65         ii = *n - 1;
  66         niter = 1;
  67
  68         temp = *rho / 2.;
  69
  70         temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
  71         i__1 = *n;
  72         for (j = 1; j <= i__1; ++j) {
  73             work[j] = d__[j] + d__[*n] + temp1;
  74             delta[j] = d__[j] - d__[*n] - temp1;
  75         }
  76
  77         psi = 0.;
  78         i__1 = *n - 2;
  79         for (j = 1; j <= i__1; ++j) {
  80             psi += z__[j] * z__[j] / (delta[j] * work[j]);
  81         }
  82
  83         c__ = rhoinv + psi;
  84         w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[*
  85                 n] / (delta[*n] * work[*n]);
  86
  87         if (w <= 0.) {
  88             temp1 = sqrt(d__[*n] * d__[*n] + *rho);
  89             temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[*
  90                     n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] *
  91                     z__[*n] / *rho;
  92
  93             if (c__ <= temp) {
  94                 tau = *rho;
  95             } else {
  96                 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
  97                 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*
  98                         n];
  99                 b = z__[*n] * z__[*n] * delsq;
 100                 if (a < 0.) {
 101                     tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
 102                 } else {
 103                     tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
 104                 }
 105             }
 106
 107         } else {
 108             delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
 109             a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
 110             b = z__[*n] * z__[*n] * delsq;
 111
 112             if (a < 0.) {
 113                 tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
 114             } else {
 115                 tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
 116             }
 117
 118         }
 119
 120         eta = tau / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau));
 121
 122         *sigma = d__[*n] + eta;
 123         i__1 = *n;
 124         for (j = 1; j <= i__1; ++j) {
 125             delta[j] = d__[j] - d__[*i__] - eta;
 126             work[j] = d__[j] + d__[*i__] + eta;
 127         }
 128
 129         dpsi = 0.;
 130         psi = 0.;
 131         erretm = 0.;
 132         i__1 = ii;
 133         for (j = 1; j <= i__1; ++j) {
 134             temp = z__[j] / (delta[j] * work[j]);
 135             psi += z__[j] * temp;
 136             dpsi += temp * temp;
 137             erretm += psi;
 138         }
 139         erretm = fabs(erretm);
 140
 141         temp = z__[*n] / (delta[*n] * work[*n]);
 142         phi = z__[*n] * temp;
 143         dphi = temp * temp;
 144         erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
 145                 + dphi);
 146
 147         w = rhoinv + phi + psi;
 148
 149         if (fabs(w) <= eps * erretm) {
 150             goto L240;
 151         }
 152
 153         ++niter;
 154         dtnsq1 = work[*n - 1] * delta[*n - 1];
 155         dtnsq = work[*n] * delta[*n];
 156         c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
 157         a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
 158         b = dtnsq * dtnsq1 * w;
 159         if (c__ < 0.) {
 160             c__ = fabs(c__);
 161         }
 162         if ( fabs(c__)<GMX_DOUBLE_MIN) {
 163             eta = *rho - *sigma * *sigma;
 164         } else if (a >= 0.) {
 165             eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__  * 2.);
 166         } else {
 167           eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
 168         }
 169
 170         if (w * eta > 0.) {
 171             eta = -w / (dpsi + dphi);
 172         }
 173         temp = eta - dtnsq;
 174         if (temp > *rho) {
 175             eta = *rho + dtnsq;
 176         }
 177
 178         tau += eta;
 179         eta /= *sigma + sqrt(eta + *sigma * *sigma);
 180         i__1 = *n;
 181         for (j = 1; j <= i__1; ++j) {
 182             delta[j] -= eta;
 183             work[j] += eta;
 184         }
 185
 186         *sigma += eta;
 187
 188         dpsi = 0.;
 189         psi = 0.;
 190         erretm = 0.;
 191         i__1 = ii;
 192         for (j = 1; j <= i__1; ++j) {
 193             temp = z__[j] / (work[j] * delta[j]);
 194             psi += z__[j] * temp;
 195             dpsi += temp * temp;
 196             erretm += psi;
 197         }
 198         erretm = fabs(erretm);
 199
 200         temp = z__[*n] / (work[*n] * delta[*n]);
 201         phi = z__[*n] * temp;
 202         dphi = temp * temp;
 203         erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
 204                 + dphi);
 205
 206         w = rhoinv + phi + psi;
 207
 208         iter = niter + 1;
 209
 210         for (niter = iter; niter <= 20; ++niter) {
 211
 212             if (fabs(w) <= eps * erretm) {
 213                 goto L240;
 214             }
 215             dtnsq1 = work[*n - 1] * delta[*n - 1];
 216             dtnsq = work[*n] * delta[*n];
 217             c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
 218             a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
 219             b = dtnsq1 * dtnsq * w;
 220             if (a >= 0.) {
 221                 eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 222             } else {
 223               eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
 224             }
 225
 226             if (w * eta > 0.) {
 227                 eta = -w / (dpsi + dphi);
 228             }
 229             temp = eta - dtnsq;
 230             if (temp <= 0.) {
 231                 eta /= 2.;
 232             }
 233
 234             tau += eta;
 235             eta /= *sigma + sqrt(eta + *sigma * *sigma);
 236             i__1 = *n;
 237             for (j = 1; j <= i__1; ++j) {
 238                 delta[j] -= eta;
 239                 work[j] += eta;
 240             }
 241
 242             *sigma += eta;
 243
 244             dpsi = 0.;
 245             psi = 0.;
 246             erretm = 0.;
 247             i__1 = ii;
 248             for (j = 1; j <= i__1; ++j) {
 249                 temp = z__[j] / (work[j] * delta[j]);
 250                 psi += z__[j] * temp;
 251                 dpsi += temp * temp;
 252                 erretm += psi;
 253             }
 254             erretm = fabs(erretm);
 255
 256             temp = z__[*n] / (work[*n] * delta[*n]);
 257             phi = z__[*n] * temp;
 258             dphi = temp * temp;
 259             erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (
 260                     dpsi + dphi);
 261
 262             w = rhoinv + phi + psi;
 263         }
 264
 265         *info = 1;
 266         goto L240;
 267
 268     } else {
 269
 270         niter = 1;
 271         ip1 = *i__ + 1;
 272
 273         delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
 274         delsq2 = delsq / 2.;
 275         temp = delsq2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + delsq2));
 276         i__1 = *n;
 277         for (j = 1; j <= i__1; ++j) {
 278             work[j] = d__[j] + d__[*i__] + temp;
 279             delta[j] = d__[j] - d__[*i__] - temp;
 280         }
 281
 282         psi = 0.;
 283         i__1 = *i__ - 1;
 284         for (j = 1; j <= i__1; ++j) {
 285             psi += z__[j] * z__[j] / (work[j] * delta[j]);
 286         }
 287
 288         phi = 0.;
 289         i__1 = *i__ + 2;
 290         for (j = *n; j >= i__1; --j) {
 291             phi += z__[j] * z__[j] / (work[j] * delta[j]);
 292         }
 293         c__ = rhoinv + psi + phi;
 294         w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[
 295                 ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
 296
 297         if (w > 0.) {
 298
 299             orgati = 1;
 300             sg2lb = 0.;
 301             sg2ub = delsq2;
 302             a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
 303             b = z__[*i__] * z__[*i__] * delsq;
 304             if (a > 0.) {
 305                 tau = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
 306             } else {
 307                 tau = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 308             }
 309             eta = tau / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau));
 310         } else {
 311
 312             orgati = 0;
 313             sg2lb = -delsq2;
 314             sg2ub = 0.;
 315             a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
 316             b = z__[ip1] * z__[ip1] * delsq;
 317             if (a < 0.) {
 318                 tau = b * 2. / (a - sqrt(fabs(a * a + b * 4. * c__)));
 319             } else {
 320                 tau = -(a + sqrt(fabs(a * a + b * 4. * c__))) / (c__ * 2.);
 321             }
 322             eta = tau / (d__[ip1] + sqrt(fabs(d__[ip1] * d__[ip1] + tau)));
 323         }
 324
 325         if (orgati) {
 326             ii = *i__;
 327             *sigma = d__[*i__] + eta;
 328             i__1 = *n;
 329             for (j = 1; j <= i__1; ++j) {
 330                 work[j] = d__[j] + d__[*i__] + eta;
 331                 delta[j] = d__[j] - d__[*i__] - eta;
 332             }
 333         } else {
 334             ii = *i__ + 1;
 335             *sigma = d__[ip1] + eta;
 336             i__1 = *n;
 337             for (j = 1; j <= i__1; ++j) {
 338                 work[j] = d__[j] + d__[ip1] + eta;
 339                 delta[j] = d__[j] - d__[ip1] - eta;
 340             }
 341         }
 342         iim1 = ii - 1;
 343         iip1 = ii + 1;
 344
 345         dpsi = 0.;
 346         psi = 0.;
 347         erretm = 0.;
 348         i__1 = iim1;
 349         for (j = 1; j <= i__1; ++j) {
 350             temp = z__[j] / (work[j] * delta[j]);
 351             psi += z__[j] * temp;
 352             dpsi += temp * temp;
 353             erretm += psi;
 354         }
 355         erretm = fabs(erretm);
 356
 357         dphi = 0.;
 358         phi = 0.;
 359         i__1 = iip1;
 360         for (j = *n; j >= i__1; --j) {
 361             temp = z__[j] / (work[j] * delta[j]);
 362             phi += z__[j] * temp;
 363             dphi += temp * temp;
 364             erretm += phi;
 365         }
 366
 367         w = rhoinv + phi + psi;
 368
 369         swtch3 = 0;
 370         if (orgati) {
 371             if (w < 0.) {
 372                 swtch3 = 1;
 373             }
 374         } else {
 375             if (w > 0.) {
 376                 swtch3 = 1;
 377             }
 378         }
 379         if (ii == 1 || ii == *n) {
 380             swtch3 = 0;
 381         }
 382
 383         temp = z__[ii] / (work[ii] * delta[ii]);
 384         dw = dpsi + dphi + temp * temp;
 385         temp = z__[ii] * temp;
 386         w += temp;
 387         erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
 388                 fabs(tau) * dw;
 389
 390         if (fabs(w) <= eps * erretm) {
 391             goto L240;
 392         }
 393
 394         if (w <= 0.) {
 395             sg2lb = (sg2lb > tau) ? sg2lb : tau;
 396         } else {
 397             sg2ub = (sg2ub < tau) ? sg2ub : tau;
 398         }
 399
 400         ++niter;
 401         if (! swtch3) {
 402             dtipsq = work[ip1] * delta[ip1];
 403             dtisq = work[*i__] * delta[*i__];
 404             if (orgati) {
 405                 d__1 = z__[*i__] / dtisq;
 406                 c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
 407             } else {
 408                 d__1 = z__[ip1] / dtipsq;
 409                 c__ = w - dtisq * dw - delsq * (d__1 * d__1);
 410             }
 411             a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
 412             b = dtipsq * dtisq * w;
 413             if ( fabs(c__)<GMX_DOUBLE_MIN) {
 414                 if ( fabs(a)<GMX_DOUBLE_MIN) {
 415                     if (orgati) {
 416                         a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi +
 417                                 dphi);
 418                     } else {
 419                         a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi +
 420                                 dphi);
 421                     }
 422                 }
 423                 eta = b / a;
 424             } else if (a <= 0.) {
 425                 eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 426             } else {
 427                 eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
 428             }
 429         } else {
 430
 431             dtiim = work[iim1] * delta[iim1];
 432             dtiip = work[iip1] * delta[iip1];
 433             temp = rhoinv + psi + phi;
 434             if (orgati) {
 435                 temp1 = z__[iim1] / dtiim;
 436                 temp1 *= temp1;
 437                 c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) *
 438                          (d__[iim1] + d__[iip1]) * temp1;
 439                 zz[0] = z__[iim1] * z__[iim1];
 440                 if (dpsi < temp1) {
 441                     zz[2] = dtiip * dtiip * dphi;
 442                 } else {
 443                     zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
 444                 }
 445             } else {
 446                 temp1 = z__[iip1] / dtiip;
 447                 temp1 *= temp1;
 448                 c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) *
 449                          (d__[iim1] + d__[iip1]) * temp1;
 450                 if (dphi < temp1) {
 451                     zz[0] = dtiim * dtiim * dpsi;
 452                 } else {
 453                     zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
 454                 }
 455                 zz[2] = z__[iip1] * z__[iip1];
 456             }
 457             zz[1] = z__[ii] * z__[ii];
 458             dd[0] = dtiim;
 459             dd[1] = delta[ii] * work[ii];
 460             dd[2] = dtiip;
 461             F77_FUNC(dlaed6,DLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
 462             if (*info != 0) {
 463                 goto L240;
 464             }
 465         }
 466
 467         if (w * eta >= 0.) {
 468             eta = -w / dw;
 469         }
 470         if (orgati) {
 471             temp1 = work[*i__] * delta[*i__];
 472             temp = eta - temp1;
 473         } else {
 474             temp1 = work[ip1] * delta[ip1];
 475             temp = eta - temp1;
 476         }
 477         if (temp > sg2ub || temp < sg2lb) {
 478             if (w < 0.) {
 479                 eta = (sg2ub - tau) / 2.;
 480             } else {
 481                 eta = (sg2lb - tau) / 2.;
 482             }
 483         }
 484
 485         tau += eta;
 486         eta /= *sigma + sqrt(*sigma * *sigma + eta);
 487
 488         prew = w;
 489
 490         *sigma += eta;
 491         i__1 = *n;
 492         for (j = 1; j <= i__1; ++j) {
 493             work[j] += eta;
 494             delta[j] -= eta;
 495         }
 496
 497         dpsi = 0.;
 498         psi = 0.;
 499         erretm = 0.;
 500         i__1 = iim1;
 501         for (j = 1; j <= i__1; ++j) {
 502             temp = z__[j] / (work[j] * delta[j]);
 503             psi += z__[j] * temp;
 504             dpsi += temp * temp;
 505             erretm += psi;
 506         }
 507         erretm = fabs(erretm);
 508
 509         dphi = 0.;
 510         phi = 0.;
 511         i__1 = iip1;
 512         for (j = *n; j >= i__1; --j) {
 513             temp = z__[j] / (work[j] * delta[j]);
 514             phi += z__[j] * temp;
 515             dphi += temp * temp;
 516             erretm += phi;
 517         }
 518
 519         temp = z__[ii] / (work[ii] * delta[ii]);
 520         dw = dpsi + dphi + temp * temp;
 521         temp = z__[ii] * temp;
 522         w = rhoinv + phi + psi + temp;
 523         erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
 524                 fabs(tau) * dw;
 525
 526         if (w <= 0.) {
 527             sg2lb = (sg2lb > tau) ? sg2lb : tau;
 528         } else {
 529             sg2ub = (sg2ub < tau) ? sg2ub : tau;
 530         }
 531
 532         swtch = 0;
 533         if (orgati) {
 534             if (-w > fabs(prew) / 10.) {
 535                 swtch = 1;
 536             }
 537         } else {
 538             if (w > fabs(prew) / 10.) {
 539                 swtch = 1;
 540             }
 541         }
 542
 543         iter = niter + 1;
 544
 545         for (niter = iter; niter <= 20; ++niter) {
 546
 547             if (fabs(w) <= eps * erretm) {
 548                 goto L240;
 549             }
 550
 551             if (! swtch3) {
 552                 dtipsq = work[ip1] * delta[ip1];
 553                 dtisq = work[*i__] * delta[*i__];
 554                 if (! swtch) {
 555                     if (orgati) {
 556                         d__1 = z__[*i__] / dtisq;
 557                         c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
 558                     } else {
 559                         d__1 = z__[ip1] / dtipsq;
 560                         c__ = w - dtisq * dw - delsq * (d__1 * d__1);
 561                     }
 562                 } else {
 563                     temp = z__[ii] / (work[ii] * delta[ii]);
 564                     if (orgati) {
 565                         dpsi += temp * temp;
 566                     } else {
 567                         dphi += temp * temp;
 568                     }
 569                     c__ = w - dtisq * dpsi - dtipsq * dphi;
 570                 }
 571                 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
 572                 b = dtipsq * dtisq * w;
 573                 if (fabs(c__)<GMX_DOUBLE_MIN) {
 574                     if (fabs(a)<GMX_DOUBLE_MIN) {
 575                         if (! swtch) {
 576                             if (orgati) {
 577                                 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq *
 578                                         (dpsi + dphi);
 579                             } else {
 580                                 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (
 581                                         dpsi + dphi);
 582                             }
 583                         } else {
 584                             a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
 585                         }
 586                     }
 587                     eta = b / a;
 588                 } else if (a <= 0.) {
 589                   eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 590                 } else {
 591                   eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
 592                 }
 593             } else {
 594
 595                 dtiim = work[iim1] * delta[iim1];
 596                 dtiip = work[iip1] * delta[iip1];
 597                 temp = rhoinv + psi + phi;
 598                 if (swtch) {
 599                     c__ = temp - dtiim * dpsi - dtiip * dphi;
 600                     zz[0] = dtiim * dtiim * dpsi;
 601                     zz[2] = dtiip * dtiip * dphi;
 602                 } else {
 603                     if (orgati) {
 604                         temp1 = z__[iim1] / dtiim;
 605                         temp1 *= temp1;
 606                         temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[
 607                                 iip1]) * temp1;
 608                         c__ = temp - dtiip * (dpsi + dphi) - temp2;
 609                         zz[0] = z__[iim1] * z__[iim1];
 610                         if (dpsi < temp1) {
 611                             zz[2] = dtiip * dtiip * dphi;
 612                         } else {
 613                             zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
 614                         }
 615                     } else {
 616                         temp1 = z__[iip1] / dtiip;
 617                         temp1 *= temp1;
 618                         temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[
 619                                 iip1]) * temp1;
 620                         c__ = temp - dtiim * (dpsi + dphi) - temp2;
 621                         if (dphi < temp1) {
 622                             zz[0] = dtiim * dtiim * dpsi;
 623                         } else {
 624                             zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
 625                         }
 626                         zz[2] = z__[iip1] * z__[iip1];
 627                     }
 628                 }
 629                 dd[0] = dtiim;
 630                 dd[1] = delta[ii] * work[ii];
 631                 dd[2] = dtiip;
 632                 F77_FUNC(dlaed6,DLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
 633                 if (*info != 0) {
 634                     goto L240;
 635                 }
 636             }
 637
 638             if (w * eta >= 0.) {
 639                 eta = -w / dw;
 640             }
 641             if (orgati) {
 642                 temp1 = work[*i__] * delta[*i__];
 643                 temp = eta - temp1;
 644             } else {
 645                 temp1 = work[ip1] * delta[ip1];
 646                 temp = eta - temp1;
 647             }
 648             if (temp > sg2ub || temp < sg2lb) {
 649                 if (w < 0.) {
 650                     eta = (sg2ub - tau) / 2.;
 651                 } else {
 652                     eta = (sg2lb - tau) / 2.;
 653                 }
 654             }
 655
 656             tau += eta;
 657             eta /= *sigma + sqrt(*sigma * *sigma + eta);
 658
 659             *sigma += eta;
 660             i__1 = *n;
 661             for (j = 1; j <= i__1; ++j) {
 662                 work[j] += eta;
 663                 delta[j] -= eta;
 664             }
 665
 666             prew = w;
 667
 668             dpsi = 0.;
 669             psi = 0.;
 670             erretm = 0.;
 671             i__1 = iim1;
 672             for (j = 1; j <= i__1; ++j) {
 673                 temp = z__[j] / (work[j] * delta[j]);
 674                 psi += z__[j] * temp;
 675                 dpsi += temp * temp;
 676                 erretm += psi;
 677             }
 678             erretm = fabs(erretm);
 679
 680             dphi = 0.;
 681             phi = 0.;
 682             i__1 = iip1;
 683             for (j = *n; j >= i__1; --j) {
 684                 temp = z__[j] / (work[j] * delta[j]);
 685                 phi += z__[j] * temp;
 686                 dphi += temp * temp;
 687                 erretm += phi;
 688             }
 689
 690             temp = z__[ii] / (work[ii] * delta[ii]);
 691             dw = dpsi + dphi + temp * temp;
 692             temp = z__[ii] * temp;
 693             w = rhoinv + phi + psi + temp;
 694             erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3.
 695                     + fabs(tau) * dw;
 696             if (w * prew > 0. && fabs(w) > fabs(prew) / 10.) {
 697                 swtch = ! swtch;
 698             }
 699
 700             if (w <= 0.) {
 701                 sg2lb = (sg2lb > tau) ? sg2lb : tau;
 702             } else {
 703                 sg2ub = (sg2ub < tau) ? sg2ub : tau;
 704             }
 705         }
 706
 707         *info = 1;
 708
 709     }
 710
 711 L240:
 712     return;
 713
 714 }