src/gromacs/linearalgebra/gmx_lapack/dlartg.c

   1 #include <math.h>
   2 #include "../gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 #include "types/simple.h"
   6
   7 void
   8 F77_FUNC(dlartg,DLARTG)(double *f,
   9         double *g,
  10         double *cs,
  11         double *sn,
  12         double *r)
  13 {
  14   double minval,safemin, safemin2, safemx2, eps;
  15   double f1,g1,f1a,g1a,scale;
  16   int i,n,count;
  17
  18   eps = GMX_DOUBLE_EPS;
  19   minval = GMX_DOUBLE_MIN;
  20   safemin = minval*(1.0+eps);
  21   n = 0.5*log( safemin/eps ) / log(2);
  22   safemin2 = pow(2,n);
  23
  24   safemx2 = 1.0 / safemin2;
  25
  26   if(fabs(*g)<GMX_DOUBLE_MIN) {
  27     *cs = 1.0;
  28     *sn = 0.0;
  29     *r = *f;
  30   } else if (fabs(*f)<GMX_DOUBLE_MIN) {
  31     *cs = 0.0;
  32     *sn = 1.0;
  33     *r = *g;
  34   } else {
  35     f1 = *f;
  36     g1 = *g;
  37     f1a = fabs(f1);
  38     g1a = fabs(g1);
  39     scale = (f1a > g1a) ? f1a : g1a;
  40     if(scale >= safemx2) {
  41       count = 0;
  42       while(scale >= safemx2) {
  43         count++;
  44         f1 *= safemin2;
  45         g1 *= safemin2;
  46         f1a = fabs(f1);
  47         g1a = fabs(g1);
  48         scale = (f1a > g1a) ? f1a : g1a;
  49       }
  50       *r = sqrt(f1*f1 + g1*g1);
  51       *cs = f1 / *r;
  52       *sn = g1 / *r;
  53       for(i=0;i<count;i++)
  54         *r *= safemx2;
  55     } else if (scale<=safemin2) {
  56       count = 0;
  57       while(scale <= safemin2) {
  58         count++;
  59         f1 *= safemx2;
  60         g1 *= safemx2;
  61         f1a = fabs(f1);
  62         g1a = fabs(g1);
  63         scale = (f1a > g1a) ? f1a : g1a;
  64       }
  65       *r = sqrt(f1*f1 + g1*g1);
  66       *cs = f1 / *r;
  67       *sn = g1 / *r;
  68       for(i=0;i<count;i++)
  69         *r *= safemin2;
  70     } else {
  71       *r = sqrt(f1*f1 + g1*g1);
  72       *cs = f1 / *r;
  73       *sn = g1 / *r;
  74     }
  75     if(fabs(*f)>fabs(*g) && *cs<0.0) {
  76       *cs *= -1.0;
  77       *sn *= -1.0;
  78       *r  *= -1.0;
  79     }
  80   }
  81   return;
  82 }
  83