diff options
Diffstat (limited to 'graphics/asymptote/fftw++.h')
-rw-r--r-- | graphics/asymptote/fftw++.h | 59 |
1 files changed, 45 insertions, 14 deletions
diff --git a/graphics/asymptote/fftw++.h b/graphics/asymptote/fftw++.h index 467387d34f..bf4102e06d 100644 --- a/graphics/asymptote/fftw++.h +++ b/graphics/asymptote/fftw++.h @@ -20,7 +20,7 @@ #ifndef __fftwpp_h__ #define __fftwpp_h__ 1 -#define __FFTWPP_H_VERSION__ 2.08 +#define __FFTWPP_H_VERSION__ 2.10 #include <cstdlib> #include <fstream> @@ -28,6 +28,7 @@ #include <fftw3.h> #include <cerrno> #include <map> +#include <typeinfo> #ifndef _OPENMP #ifndef FFTWPP_SINGLE_THREAD @@ -122,9 +123,10 @@ public: ThreadBase(unsigned int threads) : threads(threads) {} void Threads(unsigned int nthreads) {threads=nthreads;} unsigned int Threads() {return threads;} + unsigned int Innerthreads() {return innerthreads;} - void multithread(unsigned int nx) { - if(nx >= threads) { + void multithread(unsigned int n) { + if(n >= threads) { innerthreads=1; } else { innerthreads=threads; @@ -177,7 +179,7 @@ public: static const char *oddshift; - // Inplace shift of Fourier origin to (nx/2,0) for even nx. + // In-place shift of Fourier origin to (nx/2,0) for even nx. static void Shift(Complex *data, unsigned int nx, unsigned int ny, unsigned int threads) { unsigned int nyp=ny/2+1; @@ -216,7 +218,7 @@ public: } } - // Inplace shift of Fourier origin to (nx/2,ny/2,0) for even nx and ny. + // In-place shift of Fourier origin to (nx/2,ny/2,0) for even nx and ny. static void Shift(Complex *data, unsigned int nx, unsigned int ny, unsigned int nz, unsigned int threads) { unsigned int nzp=nz/2+1; @@ -306,21 +308,30 @@ public: plan=planT; fft(in,out); unsigned int N=1; + unsigned int ndoubles=doubles/2; for(;;) { double t0=utils::totalseconds(); threads=1; plan=plan1; - for(unsigned int i=0; i < N; ++i) + for(unsigned int i=0; i < N; ++i) { + for(unsigned int i=0; i < ndoubles; ++i) out[i]=i; fft(in,out); + } double t1=utils::totalseconds(); threads=Threads; plan=planT; - for(unsigned int i=0; i < N; ++i) + for(unsigned int i=0; i < N; ++i) { + for(unsigned int i=0; i < ndoubles; ++i) out[i]=i; fft(in,out); + } double t=utils::totalseconds(); S.add(t1-t0); ST.add(t-t1); - if(S.mean() < 100.0/CLOCKS_PER_SEC) N *= 2; + if(S.mean() < 100.0/CLOCKS_PER_SEC) { + N *= 2; + S.clear(); + ST.clear(); + } if(S.count() >= 10) { double error=S.stdev(); double diff=ST.mean()-S.mean(); @@ -365,7 +376,10 @@ public: threaddata data; unsigned int Threads=threads; + if(threads > 1) data=lookup(inplace,threads); + else data=threaddata(1,0.0,0.0); + threads=data.threads > 0 ? data.threads : 1; planThreads(threads); plan=(*planner)(this,in,out); @@ -386,6 +400,10 @@ public: } if(alloc) Array::deleteAlign(in,(doubles+1)/2); +#ifdef FFTWPP_VERBOSE + if(threads > 1) + std::cout << "Using " << threads << " threads." << std::endl; +#endif return data; } @@ -513,12 +531,12 @@ public: exit(1); } plan=NULL; + if(!out) out=in; + inplace=(out==in); if(rows == 0 || cols == 0) return; size /= sizeof(double); length *= size; - if(!out) out=in; - inplace=(out==in); fftw::planThreads(threads); fftw_iodim dims[3]; @@ -702,8 +720,10 @@ public: threaddata S1=Setup(in,out); fftw_plan planT1=plan; + threads=S1.threads; + bool hermitian=typeid(I) == typeid(double) || typeid(O) == typeid(double); - if(fftw::maxthreads > 1) { + if(fftw::maxthreads > 1 && (!hermitian || ostride*(nx/2+1) < idist)) { if(Threads > 1) { T=std::min(M,Threads); Q=T > 0 ? M/T : 0; @@ -727,7 +747,6 @@ public: Q=M; R=0; plan=planT1; - threads=S1.threads; } else { // Do the multi-threading ourselves fftw_destroy_plan(planT1); threads=ST.threads; @@ -810,11 +829,16 @@ public: // mfft1d Forward(n,-1,M,stride,dist,in,out); // Forward.fft(in,out); // +// mfft1d Forward(n,-1,M,istride,ostride,idist,odist,in,out); +// Forward.fft(in,out); +// // In-place usage: // // mfft1d Forward(n,-1,M,stride,dist); // Forward.fft(in); // +// +// // Notes: // stride is the spacing between the elements of each Complex vector; // dist is the spacing between the first elements of the vectors. @@ -824,7 +848,14 @@ class mfft1d : public fftwblock<fftw_complex,fftw_complex>, public Threadtable<keytype3,keyless3> { static Table threadtable; public: - mfft1d(unsigned int nx, int sign, unsigned int M=1, size_t stride=1, + mfft1d(unsigned int nx, int sign, unsigned int M=1, + Complex *in=NULL, Complex *out=NULL, + unsigned int threads=maxthreads) : + fftw(2*((nx-1)+(M-1)*nx+1),sign,threads,nx), + fftwblock<fftw_complex,fftw_complex> + (nx,M,1,1,nx,nx,in,out,threads) {} + + mfft1d(unsigned int nx, int sign, unsigned int M, size_t stride=1, size_t dist=0, Complex *in=NULL, Complex *out=NULL, unsigned int threads=maxthreads) : fftw(2*((nx-1)*stride+(M-1)*Dist(nx,stride,dist)+1),sign,threads,nx), @@ -833,7 +864,7 @@ public: mfft1d(unsigned int nx, int sign, unsigned int M, size_t istride, size_t ostride, size_t idist, size_t odist, - Complex *in=NULL, Complex *out=NULL, unsigned int threads=maxthreads): + Complex *in, Complex *out, unsigned int threads=maxthreads): fftw(std::max(2*((nx-1)*istride+(M-1)*Dist(nx,istride,idist)+1), 2*((nx-1)*ostride+(M-1)*Dist(nx,ostride,odist)+1)),sign, threads, nx), |