src/gromacs/linearalgebra/gmx_lapack/dlasq1.cpp

   1 #include <cmath>
   2 #include "../gmx_blas.h"
   3 #include "../gmx_lapack.h"
   4 #include "lapack_limits.h"
   5
   6 #include "gromacs/utility/real.h"
   7
   8 void
   9 F77_FUNC(dlasq1,DLASQ1)(int *n,
  10         double *d,
  11         double *e,
  12         double *work,
  13         int *info)
  14 {
  15   double sigmx = 0.0;
  16   int i,j,k,iinfo;
  17   double minval,safemin;
  18   double dtemp,scale;
  19   double eps;
  20
  21   eps = GMX_DOUBLE_EPS;
  22   minval = GMX_DOUBLE_MIN;
  23   safemin = minval*(1.0+GMX_DOUBLE_EPS);
  24   *info = 0;
  25
  26   if(*n<0) {
  27     *info = -2;
  28     return;
  29   }
  30
  31   for(i=0;i<*n-1;i++) {
  32     d[i] = std::abs(d[i]);
  33     dtemp = std::abs(e[i]);
  34     if(dtemp>sigmx)
  35       sigmx=dtemp;
  36   }
  37   d[*n-1] = std::abs(d[*n-1]);
  38
  39   if(std::abs(sigmx)<GMX_DOUBLE_MIN) {
  40     F77_FUNC(dlasrt,DLASRT)("D",n,d,&iinfo);
  41     return;
  42   }
  43
  44   for(i=0;i<*n;i++) {
  45     if(d[i]>sigmx)
  46       sigmx=d[i];
  47   }
  48
  49   /* Copy d and e into work (z format) and scale.
  50    * Squaring input data makes scaling by a power of the
  51    * radix pointless.
  52    */
  53   scale =  std::sqrt(eps/safemin);
  54   i = 1;
  55   j = 2;
  56   F77_FUNC(dcopy,DCOPY)(n,d,&i,work,&j);
  57   k = *n-1;
  58   F77_FUNC(dcopy,DCOPY)(&k,e,&i,work+1,&j);
  59   i = 0;
  60   j = 2*(*n)-1;
  61   k = 1;
  62   F77_FUNC(dlascl,DLASCL)("G",&i,&i,&sigmx,&scale,&j,&k,work,&j,&iinfo);
  63
  64
  65   /* Compute q and e elements */
  66   for(i=0;i<2*(*n)-1;i++)
  67     work[i] = work[i]*work[i];
  68
  69   work[2*(*n)-1] = 0.0;
  70
  71   F77_FUNC(dlasq2,DLASQ2)(n,work,info);
  72
  73   j = 0;
  74   k = 1;
  75   if(*info==0) {
  76     for(i=0;i<*n;i++)
  77       d[i]= std::sqrt(work[i]);
  78     F77_FUNC(dlascl,DLASCL)("G",&j,&j,&scale,&sigmx,n,&k,d,n,&iinfo);
  79   }
  80   return;
  81 }