src/gromacs/linearalgebra/gmx_lapack/slas2.c

   1 #include <math.h>
   2 #include "types/simple.h"
   3
   4 #include "../gmx_lapack.h"
   5
   6 void
   7 F77_FUNC(slas2,SLAS2)(float *f,
   8        float *g,
   9        float *h,
  10        float *ssmin,
  11        float *ssmax)
  12 {
  13   float fa = fabs(*f);
  14   float ga = fabs(*g);
  15   float ha = fabs(*h);
  16   float fhmin,fhmax,tmax,tmin,tmp1,tmp2;
  17   float as,at,au,c;
  18
  19   fhmin = (fa<ha) ? fa : ha;
  20   fhmax = (fa>ha) ? fa : ha;
  21
  22   if(fabs(fhmin)<GMX_FLOAT_MIN) {
  23     *ssmin = 0.0;
  24     if(fabs(fhmax)<GMX_FLOAT_MIN)
  25       *ssmax = ga;
  26     else {
  27       tmax = (fhmax>ga) ? fhmax : ga;
  28       tmin = (fhmax<ga) ? fhmax : ga;
  29       tmp1 = tmin / tmax;
  30       tmp1 = tmp1 * tmp1;
  31       *ssmax = tmax*sqrt(1.0 + tmp1);
  32     }
  33   } else {
  34     if(ga<fhmax) {
  35       as = 1.0 + fhmin / fhmax;
  36       at = (fhmax-fhmin) / fhmax;
  37       au = (ga/fhmax);
  38       au = au * au;
  39       c = 2.0 / ( sqrt(as*as+au) + sqrt(at*at+au) );
  40       *ssmin = fhmin * c;
  41       *ssmax = fhmax / c;
  42     } else {
  43       au = fhmax / ga;
  44       if(fabs(au)<GMX_FLOAT_MIN) {
  45         *ssmin = (fhmin*fhmax)/ga;
  46         *ssmax = ga;
  47       } else {
  48         as = 1.0 + fhmin / fhmax;
  49         at = (fhmax-fhmin)/fhmax;
  50         tmp1 = as*au;
  51         tmp2 = at*au;
  52         c = 1.0 / ( sqrt(1.0+tmp1*tmp1) + sqrt(1.0+tmp2*tmp2));
  53         *ssmin = (fhmin*c)*au;
  54         *ssmin = *ssmin + *ssmin;
  55         *ssmax = ga / (c+c);
  56       }
  57     }
  58   }
  59   return;
  60 }