src/gmxlib/gmx_lapack/slasd5.c

   1 #include <math.h>
   2 #include "gmx_lapack.h"
   3
   4 void
   5 F77_FUNC(slasd5,SLASD5)(int *i__,
   6         float *d__,
   7         float *z__,
   8         float *delta,
   9         float *rho,
  10         float *dsigma,
  11         float *work)
  12 {
  13     float b, c__, w, del, tau, delsq;
  14
  15     --work;
  16     --delta;
  17     --z__;
  18     --d__;
  19
  20     del = d__[2] - d__[1];
  21     delsq = del * (d__[2] + d__[1]);
  22     if (*i__ == 1) {
  23         w = *rho * 4. * (z__[2] * z__[2] / (d__[1] + d__[2] * 3.) - z__[1] *
  24                 z__[1] / (d__[1] * 3. + d__[2])) / del + 1.;
  25         if (w > 0.) {
  26             b = delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  27             c__ = *rho * z__[1] * z__[1] * delsq;
  28
  29             tau = c__ * 2. / (b + sqrt(fabs(b * b - c__ * 4.)));
  30
  31             tau /= d__[1] + sqrt(d__[1] * d__[1] + tau);
  32             *dsigma = d__[1] + tau;
  33             delta[1] = -tau;
  34             delta[2] = del - tau;
  35             work[1] = d__[1] * 2. + tau;
  36             work[2] = d__[1] + tau + d__[2];
  37         } else {
  38             b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  39             c__ = *rho * z__[2] * z__[2] * delsq;
  40
  41             if (b > 0.) {
  42                 tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
  43             } else {
  44                 tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
  45             }
  46
  47             tau /= d__[2] + sqrt(fabs(d__[2] * d__[2] + tau));
  48             *dsigma = d__[2] + tau;
  49             delta[1] = -(del + tau);
  50             delta[2] = -tau;
  51             work[1] = d__[1] + tau + d__[2];
  52             work[2] = d__[2] * 2. + tau;
  53         }
  54     } else {
  55
  56         b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  57         c__ = *rho * z__[2] * z__[2] * delsq;
  58
  59         if (b > 0.) {
  60             tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
  61         } else {
  62             tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
  63         }
  64         tau /= d__[2] + sqrt(d__[2] * d__[2] + tau);
  65         *dsigma = d__[2] + tau;
  66         delta[1] = -(del + tau);
  67         delta[2] = -tau;
  68         work[1] = d__[1] + tau + d__[2];
  69         work[2] = d__[2] * 2. + tau;
  70     }
  71     return;
  72
  73 }