diff options
Diffstat (limited to 'Build/source/utils/asymptote/fftw++.h')
-rw-r--r-- | Build/source/utils/asymptote/fftw++.h | 130 |
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);} |