mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 15:55:37 +00:00
Remez tested
This commit is contained in:
Peter Boyle
parent
d0e4673a3f
commit
2843264bd8
@ -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);
|
||||
|
||||
@ -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
|
||||
|
||||
162
lib/algorithms/approx/bigfloat_double.h
Normal file
162
lib/algorithms/approx/bigfloat_double.h
Normal file
@ -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);
|
||||
};
|
||||
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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
|
||||
|
||||
103
tests/Grid_remez.cc
Normal file
103
tests/Grid_remez.cc
Normal file
@ -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();
|
||||
}
|
||||
@ -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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user