src/gmxlib/gmx_lapack/dlaed6.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, 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
  37 #include "gmx_lapack.h"
  38
  39 #include <types/simple.h>
  40
  41 void
  42 F77_FUNC(dlaed6,DLAED6)(int *kniter,
  43                         int *orgati,
  44                         double *rho,
  45                         double *d__,
  46                         double *z__,
  47                         double *finit,
  48                         double *tau,
  49                         int *info)
  50 {
  51     int i__1;
  52     double r__1, r__2, r__3, r__4;
  53
  54     double a, b, c__, f;
  55     int i__;
  56     double fc, df, ddf, eta, eps, base;
  57     int iter;
  58     double temp, temp1, temp2, temp3, temp4;
  59     int scale;
  60     int niter;
  61     double small1, small2, sminv1, sminv2, dscale[3], sclfac;
  62     double zscale[3], erretm;
  63     double safemin;
  64     double sclinv = 0;
  65
  66     --z__;
  67     --d__;
  68
  69     *info = 0;
  70
  71     niter = 1;
  72     *tau = 0.f;
  73     if (*kniter == 2) {
  74         if (*orgati) {
  75             temp = (d__[3] - d__[2]) / 2.f;
  76             c__ = *rho + z__[1] / (d__[1] - d__[2] - temp);
  77             a = c__ * (d__[2] + d__[3]) + z__[2] + z__[3];
  78             b = c__ * d__[2] * d__[3] + z__[2] * d__[3] + z__[3] * d__[2];
  79         } else {
  80             temp = (d__[1] - d__[2]) / 2.f;
  81             c__ = *rho + z__[3] / (d__[3] - d__[2] - temp);
  82             a = c__ * (d__[1] + d__[2]) + z__[1] + z__[2];
  83             b = c__ * d__[1] * d__[2] + z__[1] * d__[2] + z__[2] * d__[1];
  84         }
  85         r__1 = fabs(a), r__2 = fabs(b), r__1 = ((r__1>r__2)? r__1:r__2), r__2 = fabs(c__);
  86         temp = (r__1>r__2) ? r__1 : r__2;
  87         a /= temp;
  88         b /= temp;
  89         c__ /= temp;
  90         if (c__ == 0.f) {
  91             *tau = b / a;
  92         } else if (a <= 0.f) {
  93             *tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1)))) / (
  94                     c__ * 2.f);
  95         } else {
  96             *tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1))));
  97         }
  98
  99         temp = *rho + z__[1] / (d__[1] - *tau) + z__[2] / (d__[2] - *tau) +
 100                 z__[3] / (d__[3] - *tau);
 101         if (fabs(*finit) <= fabs(temp)) {
 102             *tau = 0.f;
 103         }
 104     }
 105
 106     eps = GMX_DOUBLE_EPS;
 107     base = 2;
 108     safemin = GMX_DOUBLE_MIN*(1.0+GMX_DOUBLE_EPS);
 109     i__1 = (int) (log(safemin) / log(base) / 3.f);
 110     small1 = pow(base, i__1);
 111     sminv1 = 1.f / small1;
 112     small2 = small1 * small1;
 113     sminv2 = sminv1 * sminv1;
 114
 115     if (*orgati) {
 116         r__3 = (r__1 = d__[2] - *tau, fabs(r__1)), r__4 = (r__2 = d__[3] - *
 117                 tau, fabs(r__2));
 118         temp = (r__3<r__4) ? r__3 : r__4;
 119     } else {
 120         r__3 = (r__1 = d__[1] - *tau, fabs(r__1)), r__4 = (r__2 = d__[2] - *
 121                 tau, fabs(r__2));
 122         temp = (r__3<r__4) ? r__3 : r__4;
 123     }
 124     scale = 0;
 125     if (temp <= small1) {
 126         scale = 1;
 127         if (temp <= small2) {
 128
 129             sclfac = sminv2;
 130             sclinv = small2;
 131         } else {
 132
 133             sclfac = sminv1;
 134             sclinv = small1;
 135
 136         }
 137
 138         for (i__ = 1; i__ <= 3; ++i__) {
 139             dscale[i__ - 1] = d__[i__] * sclfac;
 140             zscale[i__ - 1] = z__[i__] * sclfac;
 141         }
 142         *tau *= sclfac;
 143     } else {
 144
 145         for (i__ = 1; i__ <= 3; ++i__) {
 146             dscale[i__ - 1] = d__[i__];
 147             zscale[i__ - 1] = z__[i__];
 148         }
 149     }
 150     fc = 0.f;
 151     df = 0.f;
 152     ddf = 0.f;
 153     for (i__ = 1; i__ <= 3; ++i__) {
 154         temp = 1.f / (dscale[i__ - 1] - *tau);
 155         temp1 = zscale[i__ - 1] * temp;
 156         temp2 = temp1 * temp;
 157         temp3 = temp2 * temp;
 158         fc += temp1 / dscale[i__ - 1];
 159         df += temp2;
 160         ddf += temp3;
 161     }
 162     f = *finit + *tau * fc;
 163
 164     if (fabs(f) <= 0.f) {
 165         goto L60;
 166     }
 167     iter = niter + 1;
 168     for (niter = iter; niter <= 20; ++niter) {
 169         if (*orgati) {
 170             temp1 = dscale[1] - *tau;
 171             temp2 = dscale[2] - *tau;
 172         } else {
 173             temp1 = dscale[0] - *tau;
 174             temp2 = dscale[1] - *tau;
 175         }
 176         a = (temp1 + temp2) * f - temp1 * temp2 * df;
 177         b = temp1 * temp2 * f;
 178         c__ = f - (temp1 + temp2) * df + temp1 * temp2 * ddf;
 179         r__1 = fabs(a), r__2 = fabs(b), r__1 = ((r__1>r__2)? r__1:r__2), r__2 = fabs(c__);
 180         temp = (r__1>r__2) ? r__1 : r__2;
 181         a /= temp;
 182         b /= temp;
 183         c__ /= temp;
 184         if (c__ == 0.f) {
 185             eta = b / a;
 186         } else if (a <= 0.f) {
 187             eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1)))) / ( c__ * 2.f);
 188         } else {
 189             eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, fabs( r__1))));
 190         }
 191         if (f * eta >= 0.f) {
 192             eta = -f / df;
 193         }
 194         temp = eta + *tau;
 195         if (*orgati) {
 196             if (eta > 0.f && temp >= dscale[2]) {
 197                 eta = (dscale[2] - *tau) / 2.f;
 198             }
 199
 200             if (eta < 0.f && temp <= dscale[1]) {
 201                 eta = (dscale[1] - *tau) / 2.f;
 202             }
 203         } else {
 204             if (eta > 0.f && temp >= dscale[1]) {
 205                 eta = (dscale[1] - *tau) / 2.f;
 206             }
 207             if (eta < 0.f && temp <= dscale[0]) {
 208                 eta = (dscale[0] - *tau) / 2.f;
 209             }
 210         }
 211         *tau += eta;
 212         fc = 0.f;
 213         erretm = 0.f;
 214         df = 0.f;
 215         ddf = 0.f;
 216         for (i__ = 1; i__ <= 3; ++i__) {
 217             temp = 1.f / (dscale[i__ - 1] - *tau);
 218             temp1 = zscale[i__ - 1] * temp;
 219             temp2 = temp1 * temp;
 220             temp3 = temp2 * temp;
 221             temp4 = temp1 / dscale[i__ - 1];
 222             fc += temp4;
 223             erretm += fabs(temp4);
 224             df += temp2;
 225             ddf += temp3;
 226         }
 227         f = *finit + *tau * fc;
 228         erretm = (fabs(*finit) + fabs(*tau) * erretm) * 8.f + fabs(*tau) * df;
 229         if (fabs(f) <= eps * erretm) {
 230             goto L60;
 231         }
 232     }
 233     *info = 1;
 234 L60:
 235     if (scale) {
 236         *tau *= sclinv;
 237     }
 238     return;
 239 }
 240
 241