Remez tested

lib/algorithms/approx/Remez.cc

@@ -116,7 +116,7 @@ double AlgRemez::generateApprox(int num_degree, int den_degree,
 				int a_len, double *a_param, int *a_pow)
 {
   std::cout<<"Degree of the approximation is ("<<num_degree<<","<<den_degree<<")\n";
-  std::cout<<"Approximating the function x^("<<pnum<<"/"<<pden<<"\n";
+  std::cout<<"Approximating the function x^("<<pnum<<"/"<<pden<<")\n";
 
   // Reallocate arrays, since degree has changed
   if (num_degree != n || den_degree != d) allocate(num_degree,den_degree);
@@ -149,6 +149,10 @@ double AlgRemez::generateApprox(int num_degree, int den_degree,
       std::cout<<"Iteration " <<iter-1<<" spread "<<(double)spread<<" delta "<<(double)delta<<std::endl; 
 
     equations();
+    if (delta < tolerance) {
+      std::cout<<"Delta too small, try increasing precision\n";
+      assert(0);
+    };    
     assert( delta>= tolerance);
 
     search(step);

lib/algorithms/approx/bigfloat.h

@@ -10,7 +10,7 @@
 #ifndef INCLUDED_BIGFLOAT_H
 #define INCLUDED_BIGFLOAT_H
 
-#ifdef HAVE_GMP
+
 #include <gmp.h>
 #include <mpf2mpfr.h>
 #include <mpfr.h>
@@ -202,167 +202,5 @@ public:
   friend int sgn(const bigfloat&);
 
 };
-#else
-
-typedef double mfloat; 
-class bigfloat {
-private:
-
-  mfloat x;
-
-public:
-
-  bigfloat() { }
-  bigfloat(const bigfloat& y) { x=y.x; }
-  bigfloat(const unsigned long u) { x=u; }
-  bigfloat(const long i) { x=i; }
-  bigfloat(const int i) { x=i;}
-  bigfloat(const float d) { x=d;}
-  bigfloat(const double d) {  x=d;}
-  bigfloat(const char *str) { x=std::stod(std::string(str));}
-  ~bigfloat(void) { }
-  operator double (void) const { return (double)x; }
-  static void setDefaultPrecision(unsigned long dprec) {
-  }
-
-  void setPrecision(unsigned long dprec) {
-  }
-  
-  unsigned long getPrecision(void) const { return 64; }
-  unsigned long getDefaultPrecision(void) const { return 64; }
-
-  bigfloat& operator=(const bigfloat& y)     { x=y.x;    return *this;  }
-  bigfloat& operator=(const unsigned long y) { x=y; return *this; }
-  bigfloat& operator=(const signed long y)   { x=y; return *this; }
-  bigfloat& operator=(const float y)    { x=y; return *this; }
-  bigfloat& operator=(const double y)   { x=y; return *this; }
-
-  size_t write(void);
-  size_t read(void);
-
-  /* Arithmetic Functions */
-
-  bigfloat& operator+=(const bigfloat& y) { return *this = *this + y; }
-  bigfloat& operator-=(const bigfloat& y) { return *this = *this - y; }
-  bigfloat& operator*=(const bigfloat& y) { return *this = *this * y; }
-  bigfloat& operator/=(const bigfloat& y) { return *this = *this / y; }
-
-  friend bigfloat operator+(const bigfloat& x, const bigfloat& y) { 
-    bigfloat a;
-    a.x=x.x+y.x;
-    return a;
-  }
-
-  friend bigfloat operator+(const bigfloat& x, const unsigned long y) {
-    bigfloat a;
-    a.x=x.x+y;
-    return a;
-  }
-
-  friend bigfloat operator-(const bigfloat& x, const bigfloat& y) {
-    bigfloat a;
-    a.x=x.x-y.x;
-    return a;
-  }
-  
-  friend bigfloat operator-(const unsigned long x, const bigfloat& y) {
-    bigfloat bx(x);
-    return bx-y;
-  }
-  
-  friend bigfloat operator-(const bigfloat& x, const unsigned long y) {
-    bigfloat by(y);
-    return x-by;
-  }
-
-  friend bigfloat operator-(const bigfloat& x) {
-    bigfloat a;
-    a.x=-x.x;
-    return a;
-  }
-
-  friend bigfloat operator*(const bigfloat& x, const bigfloat& y) {
-    bigfloat a;
-    a.x=x.x*y.x;
-    return a;
-  }
-
-  friend bigfloat operator*(const bigfloat& x, const unsigned long y) {
-    bigfloat a;
-    a.x=x.x*y;
-    return a;
-  }
-
-  friend bigfloat operator/(const bigfloat& x, const bigfloat& y){
-    bigfloat a;
-    a.x=x.x/y.x;
-    return a;
-  }
-
-  friend bigfloat operator/(const unsigned long x, const bigfloat& y){
-    bigfloat bx(x);
-    return bx/y;
-  }
-
-  friend bigfloat operator/(const bigfloat& x, const unsigned long y){
-    bigfloat by(y);
-    return x/by;
-  }
-
-  friend bigfloat sqrt_bf(const bigfloat& x){
-    bigfloat a;
-    a.x= sqrt(x.x);
-    return a;
-  }
-
-  friend bigfloat sqrt_bf(const unsigned long x){
-    bigfloat a(x);
-    return sqrt_bf(a);
-  }
-
-  friend bigfloat abs_bf(const bigfloat& x){
-    bigfloat a;
-    a.x=abs(x.x);
-    return a;
-  }
-
-  friend bigfloat pow_bf(const bigfloat& a, long power) {
-    bigfloat b;
-    b.x=pow(a.x,power);
-    return b;
-  }
-
-  friend bigfloat pow_bf(const bigfloat& a, bigfloat &power) {
-    bigfloat b;
-    b.x=pow(a.x,power.x);
-    return b;
-  }
-
-  friend bigfloat exp_bf(const bigfloat& a) {
-    bigfloat b;
-    b.x=exp(a.x);
-    return b;
-  }
-
-  /* Comparison Functions */
-  friend int operator>(const bigfloat& x, const bigfloat& y) {
-    return x.x>y.x;
-  }
-
-  friend int operator<(const bigfloat& x, const bigfloat& y) {
-    return x.x<y.x;
-  }
-
-  friend int sgn(const bigfloat& x) {
-    if ( x.x>=0 )  return 1;   
-    else return 0;
-  }
-
-  /* Miscellaneous Functions */
-
-  //  friend bigfloat& random(void);
-};
-
-#endif
 
 #endif

lib/algorithms/approx/bigfloat_double.h

@@ -0,0 +1,162 @@
+#include <math.h>
+
+typedef double mfloat; 
+class bigfloat {
+private:
+
+  mfloat x;
+
+public:
+
+  bigfloat() { }
+  bigfloat(const bigfloat& y) { x=y.x; }
+  bigfloat(const unsigned long u) { x=u; }
+  bigfloat(const long i) { x=i; }
+  bigfloat(const int i) { x=i;}
+  bigfloat(const float d) { x=d;}
+  bigfloat(const double d) {  x=d;}
+  bigfloat(const char *str) { x=std::stod(std::string(str));}
+  ~bigfloat(void) { }
+  operator double (void) const { return (double)x; }
+  static void setDefaultPrecision(unsigned long dprec) {
+  }
+
+  void setPrecision(unsigned long dprec) {
+  }
+  
+  unsigned long getPrecision(void) const { return 64; }
+  unsigned long getDefaultPrecision(void) const { return 64; }
+
+  bigfloat& operator=(const bigfloat& y)     { x=y.x;    return *this;  }
+  bigfloat& operator=(const unsigned long y) { x=y; return *this; }
+  bigfloat& operator=(const signed long y)   { x=y; return *this; }
+  bigfloat& operator=(const float y)    { x=y; return *this; }
+  bigfloat& operator=(const double y)   { x=y; return *this; }
+
+  size_t write(void);
+  size_t read(void);
+
+  /* Arithmetic Functions */
+
+  bigfloat& operator+=(const bigfloat& y) { return *this = *this + y; }
+  bigfloat& operator-=(const bigfloat& y) { return *this = *this - y; }
+  bigfloat& operator*=(const bigfloat& y) { return *this = *this * y; }
+  bigfloat& operator/=(const bigfloat& y) { return *this = *this / y; }
+
+  friend bigfloat operator+(const bigfloat& x, const bigfloat& y) { 
+    bigfloat a;
+    a.x=x.x+y.x;
+    return a;
+  }
+
+  friend bigfloat operator+(const bigfloat& x, const unsigned long y) {
+    bigfloat a;
+    a.x=x.x+y;
+    return a;
+  }
+
+  friend bigfloat operator-(const bigfloat& x, const bigfloat& y) {
+    bigfloat a;
+    a.x=x.x-y.x;
+    return a;
+  }
+  
+  friend bigfloat operator-(const unsigned long x, const bigfloat& y) {
+    bigfloat bx(x);
+    return bx-y;
+  }
+  
+  friend bigfloat operator-(const bigfloat& x, const unsigned long y) {
+    bigfloat by(y);
+    return x-by;
+  }
+
+  friend bigfloat operator-(const bigfloat& x) {
+    bigfloat a;
+    a.x=-x.x;
+    return a;
+  }
+
+  friend bigfloat operator*(const bigfloat& x, const bigfloat& y) {
+    bigfloat a;
+    a.x=x.x*y.x;
+    return a;
+  }
+
+  friend bigfloat operator*(const bigfloat& x, const unsigned long y) {
+    bigfloat a;
+    a.x=x.x*y;
+    return a;
+  }
+
+  friend bigfloat operator/(const bigfloat& x, const bigfloat& y){
+    bigfloat a;
+    a.x=x.x/y.x;
+    return a;
+  }
+
+  friend bigfloat operator/(const unsigned long x, const bigfloat& y){
+    bigfloat bx(x);
+    return bx/y;
+  }
+
+  friend bigfloat operator/(const bigfloat& x, const unsigned long y){
+    bigfloat by(y);
+    return x/by;
+  }
+
+  friend bigfloat sqrt_bf(const bigfloat& x){
+    bigfloat a;
+    a.x= sqrt(x.x);
+    return a;
+  }
+
+  friend bigfloat sqrt_bf(const unsigned long x){
+    bigfloat a(x);
+    return sqrt_bf(a);
+  }
+
+  friend bigfloat abs_bf(const bigfloat& x){
+    bigfloat a;
+    a.x=fabs(x.x);
+    return a;
+  }
+
+  friend bigfloat pow_bf(const bigfloat& a, long power) {
+    bigfloat b;
+    b.x=pow(a.x,power);
+    return b;
+  }
+
+  friend bigfloat pow_bf(const bigfloat& a, bigfloat &power) {
+    bigfloat b;
+    b.x=pow(a.x,power.x);
+    return b;
+  }
+
+  friend bigfloat exp_bf(const bigfloat& a) {
+    bigfloat b;
+    b.x=exp(a.x);
+    return b;
+  }
+
+  /* Comparison Functions */
+  friend int operator>(const bigfloat& x, const bigfloat& y) {
+    return x.x>y.x;
+  }
+
+  friend int operator<(const bigfloat& x, const bigfloat& y) {
+    return x.x<y.x;
+  }
+
+  friend int sgn(const bigfloat& x) {
+    if ( x.x>=0 )  return 1;   
+    else return 0;
+  }
+
+  /* Miscellaneous Functions */
+
+  //  friend bigfloat& random(void);
+};
+
+

lib/parallelIO/GridNerscIO.h

@@ -481,7 +481,7 @@ inline void writeNerscObject(Lattice<vobj> &Umu,std::string file,munger munge,in
   // Now write the body
   //////////////////////////////////////////////////
   {
-    std::ofstream fout(file,std::ios::binary|std::ios::in);
+    std::ofstream fout(file,std::ios::binary|std::ios::out);
     fout.seekp(offset);
 
     Umu = zero;

scripts/configure-commands

@@ -17,46 +17,46 @@ echo -e $YELLOW
 
 case $WD in
 g++-sse4)
-  CXX=g++-5 ../../configure --enable-simd=SSE4 CXXFLAGS="-msse4 -O3 -std=c++11" --enable-comms=none
+  CXX=g++-5 ../../configure --enable-simd=SSE4 CXXFLAGS="-msse4 -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 g++-avx)
-  CXX=g++-5 ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" --enable-comms=none
+  CXX=g++-5 ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 icpc-avx)
-  CXX=icpc ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" --enable-comms=none
+  CXX=icpc ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 icpc-avx2)
-  CXX=icpc ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -std=c++11" --enable-comms=none
+  CXX=icpc ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -std=c++11" LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 icpc-avx512)
-  CXX=icpc ../../configure --enable-simd=AVX512 CXXFLAGS="-xCOMMON-AVX512 -O3  -std=c++11" --host=none  --enable-comms=none
+  CXX=icpc ../../configure --enable-simd=AVX512 CXXFLAGS="-xCOMMON-AVX512 -O3  -std=c++11" --host=none  LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 icpc-mic)
-  CXX=icpc ../../configure --host=none --enable-simd=AVX512 CXXFLAGS="-mmic -O3  -std=c++11" LDFLAGS=-mmic --enable-comms=none
+  CXX=icpc ../../configure --host=none --enable-simd=AVX512 CXXFLAGS="-mmic -O3  -std=c++11" LDFLAGS=-mmic LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 clang-avx)
-CXX=clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11"  --enable-comms=none
+CXX=clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -std=c++11"  LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 clang-avx2)
-CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -std=c++11"  --enable-comms=none
+CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -std=c++11"  LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 clang-avx-openmp)
-CXX=clang-omp++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -std=c++11" LDFLAGS="-fopenmp"  --enable-comms=none
+CXX=clang-omp++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -std=c++11" LDFLAGS="-fopenmp"  LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 clang-avx2-openmp)
-CXX=clang-omp++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -fopenmp -O3 -std=c++11" LDFLAGS="-fopenmp"  --enable-comms=none
+CXX=clang-omp++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -fopenmp -O3 -std=c++11" LDFLAGS="-fopenmp"  LIBS="-lgmp -lmpfr" --enable-comms=none
   ;;
 clang-avx-openmp-mpi)
-CXX=clang-omp++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp"  --enable-comms=mpi
+CXX=clang-omp++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp"  LIBS="-lgmp -lmpfr" --enable-comms=mpi
 ;;
 clang-avx2-openmp-mpi)
-CXX=clang-omp++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp"  --enable-comms=mpi
+CXX=clang-omp++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -fopenmp -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx -fopenmp"  LIBS="-lgmp -lmpfr" --enable-comms=mpi
 ;;
 clang-avx-mpi)
-CXX=clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx "  --enable-comms=mpi
+CXX=clang++ ../../configure --enable-simd=AVX CXXFLAGS="-mavx -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx "  LIBS="-lgmp -lmpfr" --enable-comms=mpi
 ;;
 clang-avx2-mpi)
-CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx "  --enable-comms=mpi
+CXX=clang++ ../../configure --enable-simd=AVX2 CXXFLAGS="-mavx2 -mfma -O3 -I/opt/local/include/openmpi-mp/ -std=c++11" LDFLAGS=-L/opt/local/lib/openmpi-mp/ LIBS="-lmpi -lmpi_cxx "  LIBS="-lgmp -lmpfr" --enable-comms=mpi
 ;;
 esac
 echo -e $NORMAL

tests/Grid_remez.cc

@@ -0,0 +1,103 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+class MultiShiftFunction {
+public:
+  std::vector<double> poles;
+  std::vector<double> residues;
+  double norm;
+  double lo,hi;
+  MultiShiftFunction(int n,double _lo,double _hi): poles(n), residues(n), lo(_lo), hi(_hi) {;};
+  double approx(double x);
+  void csv(std::ostream &out);
+  void gnuplot(std::ostream &out);
+};
+double MultiShiftFunction::approx(double x)
+{
+  double a = norm;
+  for(int n=0;n<poles.size();n++){
+    a = a + residues[n]/(x+poles[n]);
+  }
+  return a;
+}
+void MultiShiftFunction::gnuplot(std::ostream &out)
+{
+  out<<"f(x) = "<<norm<<"";
+  for(int n=0;n<poles.size();n++){
+    out<<"+("<<residues[n]<<"/(x+"<<poles[n]<<"))";
+  }
+  out<<";"<<std::endl;
+}
+void MultiShiftFunction::csv(std::ostream &out)
+{
+  for (double x=lo; x<hi; x*=1.05) {
+    double f = approx(x);
+    double r = sqrt(x);
+    out<< x<<","<<r<<","<<f<<","<<r-f<<std::endl;
+  }
+  return;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::cout << "Testing Remez"<<std::endl;
+
+  double     lo=0.01;
+  double     hi=1.0;
+  int precision=64;
+  int    degree=10;
+  AlgRemez remez(0.001,1.0,precision);
+
+  ////////////////////////////////////////
+  // sqrt and inverse sqrt
+  ////////////////////////////////////////
+  MultiShiftFunction Sqrt(degree,lo,hi);
+  MultiShiftFunction InvSqrt(degree,lo,hi);
+
+  MultiShiftFunction SqrtSqrt(degree,lo,hi);
+  MultiShiftFunction InvSqrtSqrt(degree,lo,hi);
+
+
+  std::cout << "Generating degree "<<degree<<" for x^(1/2)"<<std::endl;
+  remez.generateApprox(degree,1,2);
+  remez.getPFE (&   Sqrt.residues[0],&   Sqrt.poles[0],&   Sqrt.norm);
+  remez.getIPFE(&InvSqrt.residues[0],&InvSqrt.poles[0],&InvSqrt.norm);
+
+  std::cout << "Generating degree "<<degree<<" for x^(1/4)"<<std::endl;
+  remez.generateApprox(degree,1,4);
+  remez.getPFE (&SqrtSqrt.residues[0],&SqrtSqrt.poles[0],&SqrtSqrt.norm);
+  remez.getIPFE(&InvSqrtSqrt.residues[0],&InvSqrtSqrt.poles[0],&InvSqrtSqrt.norm);
+  
+  ofstream gnuplot(std::string("Sqrt.gnu"),std::ios::out|std::ios::trunc);
+  Sqrt.gnuplot(gnuplot);
+
+  ofstream gnuplot_inv(std::string("InvSqrt.gnu"),std::ios::out|std::ios::trunc);
+  InvSqrt.gnuplot(gnuplot);
+
+  double x=0.6789;
+  double sx=sqrt(x);
+  double ssx=sqrt(sx);
+  double isx=1.0/sx;
+  double issx=1.0/ssx;
+
+  double asx  =Sqrt.approx(x);
+  double assx =SqrtSqrt.approx(x);
+  double aisx =InvSqrt.approx(x);
+  double aissx=InvSqrtSqrt.approx(x);
+
+  std::cout << "x^(1/2) : "<<sx<<" "<<asx<<std::endl;
+  std::cout << "x^(1/4) : "<<ssx<<" "<<assx<<std::endl;
+  std::cout << "x^(-1/2): "<<isx<<" "<<aisx<<std::endl;
+  std::cout << "x^(-1/4): "<<issx<<" "<<aissx<<std::endl;
+  assert(fabs(sx-asx)<1.0e-6);
+  assert(fabs(ssx-assx)<1.0e-6);
+  assert(fabs(isx-aisx)<1.0e-6);
+  assert(fabs(issx-aissx)<1.0e-6);
+
+  Grid_finalize();
+}

tests/Makefile.am

@@ -5,7 +5,7 @@ AM_LDFLAGS = -L$(top_builddir)/lib
 #
 # Test code
 #
-bin_PROGRAMS = Grid_main Grid_stencil Grid_nersc_io Grid_cshift Grid_gamma  Grid_simd Grid_rng
+bin_PROGRAMS = Grid_main Grid_stencil Grid_nersc_io Grid_cshift Grid_gamma  Grid_simd Grid_rng Grid_remez
 
 Grid_main_SOURCES = Grid_main.cc
 Grid_main_LDADD = -lGrid
@@ -13,6 +13,9 @@ Grid_main_LDADD = -lGrid
 Grid_rng_SOURCES = Grid_rng.cc
 Grid_rng_LDADD = -lGrid
 
+Grid_remez_SOURCES = Grid_remez.cc
+Grid_remez_LDADD = -lGrid
+
 Grid_nersc_io_SOURCES = Grid_nersc_io.cc
 Grid_nersc_io_LDADD = -lGrid