src/gromacs/linearalgebra/gmx_lapack/slaed6.c

   1 #include <math.h>
   2
   3 #include "../gmx_lapack.h"
   4
   5 #include "gromacs/utility/real.h"
   6
   7 void
   8 F77_FUNC(slaed6,SLAED6)(int *kniter,
   9                         int *orgati,
  10                         float *rho,
  11                         float *d__,
  12                         float *z__,
  13                         float *finit,
  14                         float *tau,
  15                         int *info)
  16 {
  17     int i__1;
  18     float r__1, r__2, r__3, r__4;
  19
  20     float a, b, c__, f;
  21     int i__;
  22     float fc, df, ddf, eta, eps, base;
  23     int iter;
  24     float temp, temp1, temp2, temp3, temp4;
  25     int scale;
  26     int niter;
  27     float small1, small2, sminv1, sminv2, dscale[3], sclfac;
  28     float zscale[3], erretm;
  29     float safemin;
  30     float sclinv = 0;
  31
  32     --z__;
  33     --d__;
  34
  35     *info = 0;
  36
  37     niter = 1;
  38     *tau = 0.f;
  39     if (*kniter == 2) {
  40         if (*orgati) {
  41             temp = (d__[3] - d__[2]) / 2.f;
  42             c__ = *rho + z__[1] / (d__[1] - d__[2] - temp);
  43             a = c__ * (d__[2] + d__[3]) + z__[2] + z__[3];
  44             b = c__ * d__[2] * d__[3] + z__[2] * d__[3] + z__[3] * d__[2];
  45         } else {
  46             temp = (d__[1] - d__[2]) / 2.f;
  47             c__ = *rho + z__[3] / (d__[3] - d__[2] - temp);
  48             a = c__ * (d__[1] + d__[2]) + z__[1] + z__[2];
  49             b = c__ * d__[1] * d__[2] + z__[1] * d__[2] + z__[2] * d__[1];
  50         }
  51         r__1 = fabs(a), r__2 = fabs(b), r__1 = ((r__1>r__2)? r__1:r__2), r__2 = fabs(c__);
  52         temp = (r__1>r__2) ? r__1 : r__2;
  53         a /= temp;
  54         b /= temp;
  55         c__ /= temp;
  56         if (c__ == 0.f) {
  57             *tau = b / a;
  58         } else if (a <= 0.f) {
  59             *tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1)))) / (
  60                     c__ * 2.f);
  61         } else {
  62             *tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1))));
  63         }
  64
  65         temp = *rho + z__[1] / (d__[1] - *tau) + z__[2] / (d__[2] - *tau) +
  66                 z__[3] / (d__[3] - *tau);
  67         if (fabs(*finit) <= fabs(temp)) {
  68             *tau = 0.f;
  69         }
  70     }
  71
  72     eps = GMX_FLOAT_EPS;
  73     base = 2;
  74     safemin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
  75     i__1 = (int) (log(safemin) / log(base) / 3.f);
  76     small1 = pow(base, i__1);
  77     sminv1 = 1.f / small1;
  78     small2 = small1 * small1;
  79     sminv2 = sminv1 * sminv1;
  80
  81     if (*orgati) {
  82         r__3 = (r__1 = d__[2] - *tau, fabs(r__1)), r__4 = (r__2 = d__[3] - *
  83                 tau, fabs(r__2));
  84         temp = (r__3<r__4) ? r__3 : r__4;
  85     } else {
  86         r__3 = (r__1 = d__[1] - *tau, fabs(r__1)), r__4 = (r__2 = d__[2] - *
  87                 tau, fabs(r__2));
  88         temp = (r__3<r__4) ? r__3 : r__4;
  89     }
  90     scale = 0;
  91     if (temp <= small1) {
  92         scale = 1;
  93         if (temp <= small2) {
  94
  95             sclfac = sminv2;
  96             sclinv = small2;
  97         } else {
  98
  99             sclfac = sminv1;
 100             sclinv = small1;
 101
 102         }
 103
 104         for (i__ = 1; i__ <= 3; ++i__) {
 105             dscale[i__ - 1] = d__[i__] * sclfac;
 106             zscale[i__ - 1] = z__[i__] * sclfac;
 107         }
 108         *tau *= sclfac;
 109     } else {
 110
 111         for (i__ = 1; i__ <= 3; ++i__) {
 112             dscale[i__ - 1] = d__[i__];
 113             zscale[i__ - 1] = z__[i__];
 114         }
 115     }
 116     fc = 0.f;
 117     df = 0.f;
 118     ddf = 0.f;
 119     for (i__ = 1; i__ <= 3; ++i__) {
 120         temp = 1.f / (dscale[i__ - 1] - *tau);
 121         temp1 = zscale[i__ - 1] * temp;
 122         temp2 = temp1 * temp;
 123         temp3 = temp2 * temp;
 124         fc += temp1 / dscale[i__ - 1];
 125         df += temp2;
 126         ddf += temp3;
 127     }
 128     f = *finit + *tau * fc;
 129
 130     if (fabs(f) <= 0.f) {
 131         goto L60;
 132     }
 133     iter = niter + 1;
 134     for (niter = iter; niter <= 20; ++niter) {
 135         if (*orgati) {
 136             temp1 = dscale[1] - *tau;
 137             temp2 = dscale[2] - *tau;
 138         } else {
 139             temp1 = dscale[0] - *tau;
 140             temp2 = dscale[1] - *tau;
 141         }
 142         a = (temp1 + temp2) * f - temp1 * temp2 * df;
 143         b = temp1 * temp2 * f;
 144         c__ = f - (temp1 + temp2) * df + temp1 * temp2 * ddf;
 145         r__1 = fabs(a), r__2 = fabs(b), r__1 = ((r__1>r__2)? r__1:r__2), r__2 = fabs(c__);
 146         temp = (r__1>r__2) ? r__1 : r__2;
 147         a /= temp;
 148         b /= temp;
 149         c__ /= temp;
 150         if (c__ == 0.f) {
 151             eta = b / a;
 152         } else if (a <= 0.f) {
 153             eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1)))) / ( c__ * 2.f);
 154         } else {
 155             eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, fabs( r__1))));
 156         }
 157         if (f * eta >= 0.f) {
 158             eta = -f / df;
 159         }
 160         temp = eta + *tau;
 161         if (*orgati) {
 162             if (eta > 0.f && temp >= dscale[2]) {
 163                 eta = (dscale[2] - *tau) / 2.f;
 164             }
 165
 166             if (eta < 0.f && temp <= dscale[1]) {
 167                 eta = (dscale[1] - *tau) / 2.f;
 168             }
 169         } else {
 170             if (eta > 0.f && temp >= dscale[1]) {
 171                 eta = (dscale[1] - *tau) / 2.f;
 172             }
 173             if (eta < 0.f && temp <= dscale[0]) {
 174                 eta = (dscale[0] - *tau) / 2.f;
 175             }
 176         }
 177         *tau += eta;
 178         fc = 0.f;
 179         erretm = 0.f;
 180         df = 0.f;
 181         ddf = 0.f;
 182         for (i__ = 1; i__ <= 3; ++i__) {
 183             temp = 1.f / (dscale[i__ - 1] - *tau);
 184             temp1 = zscale[i__ - 1] * temp;
 185             temp2 = temp1 * temp;
 186             temp3 = temp2 * temp;
 187             temp4 = temp1 / dscale[i__ - 1];
 188             fc += temp4;
 189             erretm += fabs(temp4);
 190             df += temp2;
 191             ddf += temp3;
 192         }
 193         f = *finit + *tau * fc;
 194         erretm = (fabs(*finit) + fabs(*tau) * erretm) * 8.f + fabs(*tau) * df;
 195         if (fabs(f) <= eps * erretm) {
 196             goto L60;
 197         }
 198     }
 199     *info = 1;
 200 L60:
 201     if (scale) {
 202         *tau *= sclinv;
 203     }
 204     return;
 205 }
 206
 207