src/gmxlib/gmx_lapack/slasd4.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36 #include "gmx_lapack.h"
  37 #include "lapack_limits.h"
  38
  39 #include <types/simple.h>
  40
  41 void
  42 F77_FUNC(slasd4,SLASD4)(int *n,
  43         int *i__,
  44         float *d__,
  45         float *z__,
  46         float *delta,
  47         float *rho,
  48         float *sigma,
  49         float *work,
  50         int *info)
  51 {
  52     int i__1;
  53     float d__1;
  54
  55     float a, b, c__;
  56     int j;
  57     float w, dd[3];
  58     int ii;
  59     float dw, zz[3];
  60     int ip1;
  61     float eta, phi, eps, tau, psi;
  62     int iim1, iip1;
  63     float dphi, dpsi;
  64     int iter;
  65     float temp, prew, sg2lb, sg2ub, temp1, temp2, dtiim, delsq,
  66             dtiip;
  67     int niter;
  68     float dtisq;
  69     int swtch;
  70     float dtnsq;
  71     float delsq2, dtnsq1;
  72     int swtch3;
  73     int orgati;
  74     float erretm, dtipsq, rhoinv;
  75
  76     --work;
  77     --delta;
  78     --z__;
  79     --d__;
  80
  81     *info = 0;
  82     if (*n == 1) {
  83
  84         *sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
  85         delta[1] = 1.;
  86         work[1] = 1.;
  87         return;
  88     }
  89     if (*n == 2) {
  90         F77_FUNC(slasd5,SLASD5)(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
  91         return;
  92     }
  93
  94     eps = GMX_FLOAT_EPS;
  95     rhoinv = 1. / *rho;
  96
  97     if (*i__ == *n) {
  98
  99         ii = *n - 1;
 100         niter = 1;
 101
 102         temp = *rho / 2.;
 103
 104         temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
 105         i__1 = *n;
 106         for (j = 1; j <= i__1; ++j) {
 107             work[j] = d__[j] + d__[*n] + temp1;
 108             delta[j] = d__[j] - d__[*n] - temp1;
 109         }
 110
 111         psi = 0.;
 112         i__1 = *n - 2;
 113         for (j = 1; j <= i__1; ++j) {
 114             psi += z__[j] * z__[j] / (delta[j] * work[j]);
 115         }
 116
 117         c__ = rhoinv + psi;
 118         w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[*
 119                 n] / (delta[*n] * work[*n]);
 120
 121         if (w <= 0.) {
 122             temp1 = sqrt(d__[*n] * d__[*n] + *rho);
 123             temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[*
 124                     n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] *
 125                     z__[*n] / *rho;
 126
 127             if (c__ <= temp) {
 128                 tau = *rho;
 129             } else {
 130                 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
 131                 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*
 132                         n];
 133                 b = z__[*n] * z__[*n] * delsq;
 134                 if (a < 0.) {
 135                     tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
 136                 } else {
 137                     tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
 138                 }
 139             }
 140
 141         } else {
 142             delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
 143             a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
 144             b = z__[*n] * z__[*n] * delsq;
 145
 146             if (a < 0.) {
 147                 tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
 148             } else {
 149                 tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
 150             }
 151
 152         }
 153
 154         eta = tau / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau));
 155
 156         *sigma = d__[*n] + eta;
 157         i__1 = *n;
 158         for (j = 1; j <= i__1; ++j) {
 159             delta[j] = d__[j] - d__[*i__] - eta;
 160             work[j] = d__[j] + d__[*i__] + eta;
 161         }
 162
 163         dpsi = 0.;
 164         psi = 0.;
 165         erretm = 0.;
 166         i__1 = ii;
 167         for (j = 1; j <= i__1; ++j) {
 168             temp = z__[j] / (delta[j] * work[j]);
 169             psi += z__[j] * temp;
 170             dpsi += temp * temp;
 171             erretm += psi;
 172         }
 173         erretm = fabs(erretm);
 174
 175         temp = z__[*n] / (delta[*n] * work[*n]);
 176         phi = z__[*n] * temp;
 177         dphi = temp * temp;
 178         erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
 179                 + dphi);
 180
 181         w = rhoinv + phi + psi;
 182
 183         if (fabs(w) <= eps * erretm) {
 184             goto L240;
 185         }
 186
 187         ++niter;
 188         dtnsq1 = work[*n - 1] * delta[*n - 1];
 189         dtnsq = work[*n] * delta[*n];
 190         c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
 191         a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
 192         b = dtnsq * dtnsq1 * w;
 193         if (c__ < 0.) {
 194             c__ = fabs(c__);
 195         }
 196         if ( fabs(c__)<GMX_FLOAT_MIN) {
 197             eta = *rho - *sigma * *sigma;
 198         } else if (a >= 0.) {
 199             eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__  * 2.);
 200         } else {
 201           eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
 202         }
 203
 204         if (w * eta > 0.) {
 205             eta = -w / (dpsi + dphi);
 206         }
 207         temp = eta - dtnsq;
 208         if (temp > *rho) {
 209             eta = *rho + dtnsq;
 210         }
 211
 212         tau += eta;
 213         eta /= *sigma + sqrt(eta + *sigma * *sigma);
 214         i__1 = *n;
 215         for (j = 1; j <= i__1; ++j) {
 216             delta[j] -= eta;
 217             work[j] += eta;
 218         }
 219
 220         *sigma += eta;
 221
 222         dpsi = 0.;
 223         psi = 0.;
 224         erretm = 0.;
 225         i__1 = ii;
 226         for (j = 1; j <= i__1; ++j) {
 227             temp = z__[j] / (work[j] * delta[j]);
 228             psi += z__[j] * temp;
 229             dpsi += temp * temp;
 230             erretm += psi;
 231         }
 232         erretm = fabs(erretm);
 233
 234         temp = z__[*n] / (work[*n] * delta[*n]);
 235         phi = z__[*n] * temp;
 236         dphi = temp * temp;
 237         erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
 238                 + dphi);
 239
 240         w = rhoinv + phi + psi;
 241
 242         iter = niter + 1;
 243
 244         for (niter = iter; niter <= 20; ++niter) {
 245
 246             if (fabs(w) <= eps * erretm) {
 247                 goto L240;
 248             }
 249             dtnsq1 = work[*n - 1] * delta[*n - 1];
 250             dtnsq = work[*n] * delta[*n];
 251             c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
 252             a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
 253             b = dtnsq1 * dtnsq * w;
 254             if (a >= 0.) {
 255                 eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 256             } else {
 257               eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
 258             }
 259
 260             if (w * eta > 0.) {
 261                 eta = -w / (dpsi + dphi);
 262             }
 263             temp = eta - dtnsq;
 264             if (temp <= 0.) {
 265                 eta /= 2.;
 266             }
 267
 268             tau += eta;
 269             eta /= *sigma + sqrt(eta + *sigma * *sigma);
 270             i__1 = *n;
 271             for (j = 1; j <= i__1; ++j) {
 272                 delta[j] -= eta;
 273                 work[j] += eta;
 274             }
 275
 276             *sigma += eta;
 277
 278             dpsi = 0.;
 279             psi = 0.;
 280             erretm = 0.;
 281             i__1 = ii;
 282             for (j = 1; j <= i__1; ++j) {
 283                 temp = z__[j] / (work[j] * delta[j]);
 284                 psi += z__[j] * temp;
 285                 dpsi += temp * temp;
 286                 erretm += psi;
 287             }
 288             erretm = fabs(erretm);
 289
 290             temp = z__[*n] / (work[*n] * delta[*n]);
 291             phi = z__[*n] * temp;
 292             dphi = temp * temp;
 293             erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (
 294                     dpsi + dphi);
 295
 296             w = rhoinv + phi + psi;
 297         }
 298
 299         *info = 1;
 300         goto L240;
 301
 302     } else {
 303
 304         niter = 1;
 305         ip1 = *i__ + 1;
 306
 307         delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
 308         delsq2 = delsq / 2.;
 309         temp = delsq2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + delsq2));
 310         i__1 = *n;
 311         for (j = 1; j <= i__1; ++j) {
 312             work[j] = d__[j] + d__[*i__] + temp;
 313             delta[j] = d__[j] - d__[*i__] - temp;
 314         }
 315
 316         psi = 0.;
 317         i__1 = *i__ - 1;
 318         for (j = 1; j <= i__1; ++j) {
 319             psi += z__[j] * z__[j] / (work[j] * delta[j]);
 320         }
 321
 322         phi = 0.;
 323         i__1 = *i__ + 2;
 324         for (j = *n; j >= i__1; --j) {
 325             phi += z__[j] * z__[j] / (work[j] * delta[j]);
 326         }
 327         c__ = rhoinv + psi + phi;
 328         w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[
 329                 ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
 330
 331         if (w > 0.) {
 332
 333             orgati = 1;
 334             sg2lb = 0.;
 335             sg2ub = delsq2;
 336             a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
 337             b = z__[*i__] * z__[*i__] * delsq;
 338             if (a > 0.) {
 339                 tau = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
 340             } else {
 341                 tau = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 342             }
 343             eta = tau / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau));
 344         } else {
 345
 346             orgati = 0;
 347             sg2lb = -delsq2;
 348             sg2ub = 0.;
 349             a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
 350             b = z__[ip1] * z__[ip1] * delsq;
 351             if (a < 0.) {
 352                 tau = b * 2. / (a - sqrt(fabs(a * a + b * 4. * c__)));
 353             } else {
 354                 tau = -(a + sqrt(fabs(a * a + b * 4. * c__))) / (c__ * 2.);
 355             }
 356             eta = tau / (d__[ip1] + sqrt(fabs(d__[ip1] * d__[ip1] + tau)));
 357         }
 358
 359         if (orgati) {
 360             ii = *i__;
 361             *sigma = d__[*i__] + eta;
 362             i__1 = *n;
 363             for (j = 1; j <= i__1; ++j) {
 364                 work[j] = d__[j] + d__[*i__] + eta;
 365                 delta[j] = d__[j] - d__[*i__] - eta;
 366             }
 367         } else {
 368             ii = *i__ + 1;
 369             *sigma = d__[ip1] + eta;
 370             i__1 = *n;
 371             for (j = 1; j <= i__1; ++j) {
 372                 work[j] = d__[j] + d__[ip1] + eta;
 373                 delta[j] = d__[j] - d__[ip1] - eta;
 374             }
 375         }
 376         iim1 = ii - 1;
 377         iip1 = ii + 1;
 378
 379         dpsi = 0.;
 380         psi = 0.;
 381         erretm = 0.;
 382         i__1 = iim1;
 383         for (j = 1; j <= i__1; ++j) {
 384             temp = z__[j] / (work[j] * delta[j]);
 385             psi += z__[j] * temp;
 386             dpsi += temp * temp;
 387             erretm += psi;
 388         }
 389         erretm = fabs(erretm);
 390
 391         dphi = 0.;
 392         phi = 0.;
 393         i__1 = iip1;
 394         for (j = *n; j >= i__1; --j) {
 395             temp = z__[j] / (work[j] * delta[j]);
 396             phi += z__[j] * temp;
 397             dphi += temp * temp;
 398             erretm += phi;
 399         }
 400
 401         w = rhoinv + phi + psi;
 402
 403         swtch3 = 0;
 404         if (orgati) {
 405             if (w < 0.) {
 406                 swtch3 = 1;
 407             }
 408         } else {
 409             if (w > 0.) {
 410                 swtch3 = 1;
 411             }
 412         }
 413         if (ii == 1 || ii == *n) {
 414             swtch3 = 0;
 415         }
 416
 417         temp = z__[ii] / (work[ii] * delta[ii]);
 418         dw = dpsi + dphi + temp * temp;
 419         temp = z__[ii] * temp;
 420         w += temp;
 421         erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
 422                 fabs(tau) * dw;
 423
 424         if (fabs(w) <= eps * erretm) {
 425             goto L240;
 426         }
 427
 428         if (w <= 0.) {
 429             sg2lb = (sg2lb > tau) ? sg2lb : tau;
 430         } else {
 431             sg2ub = (sg2ub < tau) ? sg2ub : tau;
 432         }
 433
 434         ++niter;
 435         if (! swtch3) {
 436             dtipsq = work[ip1] * delta[ip1];
 437             dtisq = work[*i__] * delta[*i__];
 438             if (orgati) {
 439                 d__1 = z__[*i__] / dtisq;
 440                 c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
 441             } else {
 442                 d__1 = z__[ip1] / dtipsq;
 443                 c__ = w - dtisq * dw - delsq * (d__1 * d__1);
 444             }
 445             a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
 446             b = dtipsq * dtisq * w;
 447             if ( fabs(c__)<GMX_FLOAT_MIN) {
 448                 if ( fabs(a)<GMX_FLOAT_MIN) {
 449                     if (orgati) {
 450                         a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi +
 451                                 dphi);
 452                     } else {
 453                         a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi +
 454                                 dphi);
 455                     }
 456                 }
 457                 eta = b / a;
 458             } else if (a <= 0.) {
 459                 eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 460             } else {
 461                 eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
 462             }
 463         } else {
 464
 465             dtiim = work[iim1] * delta[iim1];
 466             dtiip = work[iip1] * delta[iip1];
 467             temp = rhoinv + psi + phi;
 468             if (orgati) {
 469                 temp1 = z__[iim1] / dtiim;
 470                 temp1 *= temp1;
 471                 c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) *
 472                          (d__[iim1] + d__[iip1]) * temp1;
 473                 zz[0] = z__[iim1] * z__[iim1];
 474                 if (dpsi < temp1) {
 475                     zz[2] = dtiip * dtiip * dphi;
 476                 } else {
 477                     zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
 478                 }
 479             } else {
 480                 temp1 = z__[iip1] / dtiip;
 481                 temp1 *= temp1;
 482                 c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) *
 483                          (d__[iim1] + d__[iip1]) * temp1;
 484                 if (dphi < temp1) {
 485                     zz[0] = dtiim * dtiim * dpsi;
 486                 } else {
 487                     zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
 488                 }
 489                 zz[2] = z__[iip1] * z__[iip1];
 490             }
 491             zz[1] = z__[ii] * z__[ii];
 492             dd[0] = dtiim;
 493             dd[1] = delta[ii] * work[ii];
 494             dd[2] = dtiip;
 495             F77_FUNC(slaed6,SLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
 496             if (*info != 0) {
 497                 goto L240;
 498             }
 499         }
 500
 501         if (w * eta >= 0.) {
 502             eta = -w / dw;
 503         }
 504         if (orgati) {
 505             temp1 = work[*i__] * delta[*i__];
 506             temp = eta - temp1;
 507         } else {
 508             temp1 = work[ip1] * delta[ip1];
 509             temp = eta - temp1;
 510         }
 511         if (temp > sg2ub || temp < sg2lb) {
 512             if (w < 0.) {
 513                 eta = (sg2ub - tau) / 2.;
 514             } else {
 515                 eta = (sg2lb - tau) / 2.;
 516             }
 517         }
 518
 519         tau += eta;
 520         eta /= *sigma + sqrt(*sigma * *sigma + eta);
 521
 522         prew = w;
 523
 524         *sigma += eta;
 525         i__1 = *n;
 526         for (j = 1; j <= i__1; ++j) {
 527             work[j] += eta;
 528             delta[j] -= eta;
 529         }
 530
 531         dpsi = 0.;
 532         psi = 0.;
 533         erretm = 0.;
 534         i__1 = iim1;
 535         for (j = 1; j <= i__1; ++j) {
 536             temp = z__[j] / (work[j] * delta[j]);
 537             psi += z__[j] * temp;
 538             dpsi += temp * temp;
 539             erretm += psi;
 540         }
 541         erretm = fabs(erretm);
 542
 543         dphi = 0.;
 544         phi = 0.;
 545         i__1 = iip1;
 546         for (j = *n; j >= i__1; --j) {
 547             temp = z__[j] / (work[j] * delta[j]);
 548             phi += z__[j] * temp;
 549             dphi += temp * temp;
 550             erretm += phi;
 551         }
 552
 553         temp = z__[ii] / (work[ii] * delta[ii]);
 554         dw = dpsi + dphi + temp * temp;
 555         temp = z__[ii] * temp;
 556         w = rhoinv + phi + psi + temp;
 557         erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
 558                 fabs(tau) * dw;
 559
 560         if (w <= 0.) {
 561             sg2lb = (sg2lb > tau) ? sg2lb : tau;
 562         } else {
 563             sg2ub = (sg2ub < tau) ? sg2ub : tau;
 564         }
 565
 566         swtch = 0;
 567         if (orgati) {
 568             if (-w > fabs(prew) / 10.) {
 569                 swtch = 1;
 570             }
 571         } else {
 572             if (w > fabs(prew) / 10.) {
 573                 swtch = 1;
 574             }
 575         }
 576
 577         iter = niter + 1;
 578
 579         for (niter = iter; niter <= 20; ++niter) {
 580
 581             if (fabs(w) <= eps * erretm) {
 582                 goto L240;
 583             }
 584
 585             if (! swtch3) {
 586                 dtipsq = work[ip1] * delta[ip1];
 587                 dtisq = work[*i__] * delta[*i__];
 588                 if (! swtch) {
 589                     if (orgati) {
 590                         d__1 = z__[*i__] / dtisq;
 591                         c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
 592                     } else {
 593                         d__1 = z__[ip1] / dtipsq;
 594                         c__ = w - dtisq * dw - delsq * (d__1 * d__1);
 595                     }
 596                 } else {
 597                     temp = z__[ii] / (work[ii] * delta[ii]);
 598                     if (orgati) {
 599                         dpsi += temp * temp;
 600                     } else {
 601                         dphi += temp * temp;
 602                     }
 603                     c__ = w - dtisq * dpsi - dtipsq * dphi;
 604                 }
 605                 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
 606                 b = dtipsq * dtisq * w;
 607                 if (fabs(c__)<GMX_FLOAT_MIN) {
 608                     if (fabs(a)<GMX_FLOAT_MIN) {
 609                         if (! swtch) {
 610                             if (orgati) {
 611                                 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq *
 612                                         (dpsi + dphi);
 613                             } else {
 614                                 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (
 615                                         dpsi + dphi);
 616                             }
 617                         } else {
 618                             a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
 619                         }
 620                     }
 621                     eta = b / a;
 622                 } else if (a <= 0.) {
 623                   eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
 624                 } else {
 625                   eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
 626                 }
 627             } else {
 628
 629                 dtiim = work[iim1] * delta[iim1];
 630                 dtiip = work[iip1] * delta[iip1];
 631                 temp = rhoinv + psi + phi;
 632                 if (swtch) {
 633                     c__ = temp - dtiim * dpsi - dtiip * dphi;
 634                     zz[0] = dtiim * dtiim * dpsi;
 635                     zz[2] = dtiip * dtiip * dphi;
 636                 } else {
 637                     if (orgati) {
 638                         temp1 = z__[iim1] / dtiim;
 639                         temp1 *= temp1;
 640                         temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[
 641                                 iip1]) * temp1;
 642                         c__ = temp - dtiip * (dpsi + dphi) - temp2;
 643                         zz[0] = z__[iim1] * z__[iim1];
 644                         if (dpsi < temp1) {
 645                             zz[2] = dtiip * dtiip * dphi;
 646                         } else {
 647                             zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
 648                         }
 649                     } else {
 650                         temp1 = z__[iip1] / dtiip;
 651                         temp1 *= temp1;
 652                         temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[
 653                                 iip1]) * temp1;
 654                         c__ = temp - dtiim * (dpsi + dphi) - temp2;
 655                         if (dphi < temp1) {
 656                             zz[0] = dtiim * dtiim * dpsi;
 657                         } else {
 658                             zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
 659                         }
 660                         zz[2] = z__[iip1] * z__[iip1];
 661                     }
 662                 }
 663                 dd[0] = dtiim;
 664                 dd[1] = delta[ii] * work[ii];
 665                 dd[2] = dtiip;
 666                 F77_FUNC(slaed6,SLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
 667                 if (*info != 0) {
 668                     goto L240;
 669                 }
 670             }
 671
 672             if (w * eta >= 0.) {
 673                 eta = -w / dw;
 674             }
 675             if (orgati) {
 676                 temp1 = work[*i__] * delta[*i__];
 677                 temp = eta - temp1;
 678             } else {
 679                 temp1 = work[ip1] * delta[ip1];
 680                 temp = eta - temp1;
 681             }
 682             if (temp > sg2ub || temp < sg2lb) {
 683                 if (w < 0.) {
 684                     eta = (sg2ub - tau) / 2.;
 685                 } else {
 686                     eta = (sg2lb - tau) / 2.;
 687                 }
 688             }
 689
 690             tau += eta;
 691             eta /= *sigma + sqrt(*sigma * *sigma + eta);
 692
 693             *sigma += eta;
 694             i__1 = *n;
 695             for (j = 1; j <= i__1; ++j) {
 696                 work[j] += eta;
 697                 delta[j] -= eta;
 698             }
 699
 700             prew = w;
 701
 702             dpsi = 0.;
 703             psi = 0.;
 704             erretm = 0.;
 705             i__1 = iim1;
 706             for (j = 1; j <= i__1; ++j) {
 707                 temp = z__[j] / (work[j] * delta[j]);
 708                 psi += z__[j] * temp;
 709                 dpsi += temp * temp;
 710                 erretm += psi;
 711             }
 712             erretm = fabs(erretm);
 713
 714             dphi = 0.;
 715             phi = 0.;
 716             i__1 = iip1;
 717             for (j = *n; j >= i__1; --j) {
 718                 temp = z__[j] / (work[j] * delta[j]);
 719                 phi += z__[j] * temp;
 720                 dphi += temp * temp;
 721                 erretm += phi;
 722             }
 723
 724             temp = z__[ii] / (work[ii] * delta[ii]);
 725             dw = dpsi + dphi + temp * temp;
 726             temp = z__[ii] * temp;
 727             w = rhoinv + phi + psi + temp;
 728             erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3.
 729                     + fabs(tau) * dw;
 730             if (w * prew > 0. && fabs(w) > fabs(prew) / 10.) {
 731                 swtch = ! swtch;
 732             }
 733
 734             if (w <= 0.) {
 735                 sg2lb = (sg2lb > tau) ? sg2lb : tau;
 736             } else {
 737                 sg2ub = (sg2ub < tau) ? sg2ub : tau;
 738             }
 739         }
 740
 741         *info = 1;
 742
 743     }
 744
 745 L240:
 746     return;
 747
 748 }