summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/fftw++.h
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/fftw++.h')
-rw-r--r--Build/source/utils/asymptote/fftw++.h130
1 files changed, 64 insertions, 66 deletions
diff --git a/Build/source/utils/asymptote/fftw++.h b/Build/source/utils/asymptote/fftw++.h
index be68c899290..addb91cff05 100644
--- a/Build/source/utils/asymptote/fftw++.h
+++ b/Build/source/utils/asymptote/fftw++.h
@@ -18,7 +18,7 @@
#ifndef __fftwpp_h__
#define __fftwpp_h__ 1
-#define __FFTWPP_H_VERSION__ 1.12
+#define __FFTWPP_H_VERSION__ 1.14
#include <cstdlib>
#include <fstream>
@@ -30,21 +30,21 @@
#include <omp.h>
#endif
-inline int get_max_threads()
+inline int get_thread_num()
{
#ifdef FFTWPP_SINGLE_THREAD
- return 1;
+ return 0;
#else
- return omp_get_max_threads();
+ return omp_get_thread_num();
#endif
}
-inline int get_thread_num()
+inline int get_max_threads()
{
#ifdef FFTWPP_SINGLE_THREAD
- return 0;
+ return 1;
#else
- return omp_get_thread_num();
+ return omp_get_max_threads();
#endif
}
@@ -198,7 +198,7 @@ protected:
fftw_plan plan;
bool inplace;
- unsigned int Dist(unsigned int n, unsigned int stride, unsigned int dist) {
+ unsigned int Dist(unsigned int n, size_t stride, size_t dist) {
return dist ? dist : ((stride == 1) ? n : 1);
}
@@ -214,13 +214,12 @@ protected:
return realsize(n,(Complex *) in,out);
}
- static std::ifstream ifWisdom;
- static std::ofstream ofWisdom;
static bool Wise;
static const double twopi;
unsigned int threads;
public:
+ static unsigned int effort;
static unsigned int maxthreads;
static double testseconds;
@@ -296,18 +295,15 @@ public:
}
}
- static unsigned int effort;
- static const char *WisdomName;
-
+ static bool autothreads;
fftw(unsigned int doubles, int sign, unsigned int n=0) :
doubles(doubles), sign(sign), norm(1.0/(n ? n : (doubles+1)/2)),
plan(NULL) {
- static bool initialize=true;
- if(initialize) {
+ if(autothreads) {
#ifndef FFTWPP_SINGLE_THREAD
fftw_init_threads();
#endif
- initialize=false;
+ autothreads=false;
}
}
@@ -330,6 +326,13 @@ public:
exit(1);
}
+ void planThreads(unsigned int threads) {
+#ifndef FFTWPP_SINGLE_THREAD
+ omp_set_num_threads(threads);
+ fftw_plan_with_nthreads(threads);
+#endif
+ }
+
void Setup(Complex *in, Complex *out=NULL) {
if(!Wise) LoadWisdom();
bool alloc=!in;
@@ -343,48 +346,53 @@ public:
#endif
inplace=(out==in);
-#ifndef FFTWPP_SINGLE_THREAD
- fftw_plan_with_nthreads(1);
-#endif
+ planThreads(1);
fftw_plan plan1=Plan(in,out);
if(!plan1) noplan();
plan=plan1;
if(maxthreads > 1) {
- double sum=0.0;
double sum2=0.0;
unsigned int N=1;
- const unsigned int microseconds=1000000;
- unsigned int limit=(int) (testseconds*microseconds);
- for(; N < limit; ++N) {
- seconds();
+ double begin=totalseconds();
+ double lastseconds=begin;
+ double stop=begin+testseconds;
+ for(;;++N) {
fft(in,out);
- double t=seconds();
- sum += t;
- sum2 += t*t;
- if(sum > testseconds)
+ double t=totalseconds();
+ double seconds=t-lastseconds;
+ sum2 += seconds*seconds;
+ lastseconds=t;
+ if(t > stop)
break;
}
+
+ double end=totalseconds();
+ double sum=end-begin;
double mean1=sum/N;
double stdev1=stdev(N,sum,sum2);
- threads=get_max_threads();
- if(maxthreads > threads) maxthreads=threads;
-#ifndef FFTWPP_SINGLE_THREAD
- fftw_plan_with_nthreads(maxthreads);
-#endif
+ threads=maxthreads;
+ planThreads(threads);
plan=Plan(in,out);
if(!plan) noplan();
if(plan) {
- sum=0.0;
- seconds();
- for(unsigned int i=1; i <= N; ++i) {
- seconds();
+ double begin=totalseconds();
+ double lastseconds=begin;
+ double stop=begin+testseconds;
+ N=1;
+ for(;;++N) {
fft(in,out);
- double t=seconds();
- sum += t;
+ double t=totalseconds();
+ double seconds=t-lastseconds;
+ sum2 += seconds*seconds;
+ lastseconds=t;
+ if(t > stop)
+ break;
}
+ double end=totalseconds();
+ double sum=end-begin;
double mean2=sum/N;
if(mean2 > mean1-stdev1) {
@@ -406,18 +414,8 @@ public:
void Setup(Complex *in, double *out) {Setup(in,(Complex *) out);}
void Setup(double *in, Complex *out=NULL) {Setup((Complex *) in,out);}
- void LoadWisdom() {
- ifWisdom.open(WisdomName);
- fftwpp_import_wisdom(GetWisdom,ifWisdom);
- ifWisdom.close();
- Wise=true;
- }
-
- void SaveWisdom() {
- ofWisdom.open(WisdomName);
- fftwpp_export_wisdom(PutWisdom,ofWisdom);
- ofWisdom.close();
- }
+ void LoadWisdom();
+ void SaveWisdom();
virtual void Execute(Complex *in, Complex *out, bool=false) {
fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
@@ -507,7 +505,7 @@ public:
void fftNormalized(Complex *in, Complex *out,
unsigned int nx, unsigned int M,
- unsigned int stride, unsigned int dist) {
+ size_t stride, size_t dist) {
if(stride == 1 && dist == nx) fftw::fftNormalized(in,out);
else if(stride == nx && dist == 1) fftw::fftNormalized(in,out);
else {
@@ -590,11 +588,11 @@ public:
class mfft1d : public fftw {
unsigned int nx;
unsigned int M;
- unsigned int stride;
- unsigned int dist;
+ size_t stride;
+ size_t dist;
public:
- mfft1d(unsigned int nx, int sign, unsigned int M=1, unsigned int stride=1,
- unsigned int dist=0, Complex *in=NULL, Complex *out=NULL)
+ mfft1d(unsigned int nx, int sign, unsigned int M=1, size_t stride=1,
+ size_t dist=0, Complex *in=NULL, Complex *out=NULL)
: fftw(2*((nx-1)*stride+(M-1)*Dist(nx,stride,dist)+1),sign,nx),
nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist))
{Setup(in,out);}
@@ -734,16 +732,16 @@ public:
class mrcfft1d : public fftw {
unsigned int nx;
unsigned int M;
- unsigned int stride;
- unsigned int dist;
+ size_t stride;
+ size_t dist;
public:
- mrcfft1d(unsigned int nx, unsigned int M=1, unsigned int stride=1,
- unsigned int dist=0, Complex *out=NULL)
+ mrcfft1d(unsigned int nx, unsigned int M=1, size_t stride=1,
+ size_t dist=0, Complex *out=NULL)
: fftw(2*(nx/2*stride+(M-1)*Dist(nx,stride,dist)+1),-1,nx), nx(nx), M(M),
stride(stride), dist(Dist(nx,stride,dist)) {Setup(out);}
- mrcfft1d(unsigned int nx, unsigned int M=1, unsigned int stride=1,
- unsigned int dist=0, double *in=NULL, Complex *out=NULL)
+ mrcfft1d(unsigned int nx, unsigned int M=1, size_t stride=1,
+ size_t dist=0, double *in=NULL, Complex *out=NULL)
: fftw(2*(nx/2*stride+(M-1)*Dist(nx,stride,dist)+1),-1,nx), nx(nx), M(M),
stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
@@ -790,11 +788,11 @@ public:
class mcrfft1d : public fftw {
unsigned int nx;
unsigned int M;
- unsigned int stride;
- unsigned int dist;
+ size_t stride;
+ size_t dist;
public:
- mcrfft1d(unsigned int nx, unsigned int M=1, unsigned int stride=1,
- unsigned int dist=0, Complex *in=NULL, double *out=NULL)
+ mcrfft1d(unsigned int nx, unsigned int M=1, size_t stride=1,
+ size_t dist=0, Complex *in=NULL, double *out=NULL)
: fftw((realsize(nx,in,out)-2)*stride+2*(M-1)*Dist(nx,stride,dist)+2,1,nx),
nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}