src/gromacs/linearalgebra/gmx_lapack/dlarfg.c

   1 #include <math.h>
   2 #include "types/simple.h"
   3
   4 #include "../gmx_blas.h"
   5 #include "../gmx_lapack.h"
   6 #include "lapack_limits.h"
   7
   8
   9 void
  10 F77_FUNC(dlarfg,DLARFG)(int   *n,
  11                         double *alpha,
  12                         double *x,
  13                         int    *incx,
  14                         double *tau)
  15 {
  16   double xnorm,t;
  17   int    ti1,knt,j;
  18   double minval,safmin,rsafmn,beta;
  19
  20   if(*n<=1) {
  21     *tau = 0;
  22     return;
  23   }
  24
  25   ti1 = *n-1;
  26
  27   xnorm = F77_FUNC(dnrm2,DNRM2)(&ti1,x,incx);
  28
  29   if(fabs(xnorm)<GMX_DOUBLE_MIN) {
  30     *tau = 0.0;
  31   } else {
  32
  33     t = F77_FUNC(dlapy2,DLAPY2)(alpha,&xnorm);
  34
  35     if(*alpha<0)
  36       beta = t;
  37     else
  38       beta = -t;
  39
  40     minval = GMX_DOUBLE_MIN;
  41
  42     safmin = minval*(1.0+GMX_DOUBLE_EPS) / GMX_DOUBLE_EPS;
  43
  44
  45     if(fabs(beta)<safmin) {
  46
  47       knt = 0;
  48       rsafmn = 1.0 / safmin;
  49
  50       while(fabs(beta)<safmin) {
  51         knt++;
  52         ti1 = *n-1;
  53         F77_FUNC(dscal,DSCAL)(&ti1,&rsafmn,x,incx);
  54         beta *= rsafmn;
  55         *alpha *= rsafmn;
  56       }
  57
  58       /* safmin <= beta <= 1 now */
  59       ti1 = *n-1;
  60       xnorm = F77_FUNC(dnrm2,DNRM2)(&ti1,x,incx);
  61       t = F77_FUNC(dlapy2,DLAPY2)(alpha,&xnorm);
  62
  63       if(*alpha<0)
  64         beta = t;
  65       else
  66         beta = -t;
  67
  68       *tau = (beta-*alpha)/beta;
  69
  70       ti1= *n-1;
  71       t = 1.0/(*alpha-beta);
  72       F77_FUNC(dscal,DSCAL)(&ti1,&t,x,incx);
  73
  74       *alpha = beta;
  75       for(j=0;j<knt;j++)
  76         *alpha *= safmin;
  77     } else {
  78       *tau = (beta-*alpha)/beta;
  79       ti1= *n-1;
  80       t = 1.0/(*alpha-beta);
  81       F77_FUNC(dscal,DSCAL)(&ti1,&t,x,incx);
  82       *alpha = beta;
  83     }
  84   }
  85
  86   return;
  87 }