diff options
Diffstat (limited to 'graphics/asymptote/fftw++.h')
-rw-r--r-- | graphics/asymptote/fftw++.h | 1046 |
1 files changed, 557 insertions, 489 deletions
diff --git a/graphics/asymptote/fftw++.h b/graphics/asymptote/fftw++.h index bf4102e06d..162954623e 100644 --- a/graphics/asymptote/fftw++.h +++ b/graphics/asymptote/fftw++.h @@ -1,5 +1,5 @@ /* Fast Fourier transform C++ header class for the FFTW3 Library - Copyright (C) 2004-16 + Copyright (C) 2004-2022 John C. Bowman, University of Alberta Malcolm Roberts, University of Strasbourg @@ -20,7 +20,7 @@ #ifndef __fftwpp_h__ #define __fftwpp_h__ 1 -#define __FFTWPP_H_VERSION__ 2.10 +#define __FFTWPP_H_VERSION__ 2.11 #include <cstdlib> #include <fstream> @@ -29,56 +29,16 @@ #include <cerrno> #include <map> #include <typeinfo> +#include <climits> -#ifndef _OPENMP -#ifndef FFTWPP_SINGLE_THREAD -#define FFTWPP_SINGLE_THREAD -#endif -#endif - -#ifndef FFTWPP_SINGLE_THREAD -#include <omp.h> -#endif - -inline int get_thread_num() -{ -#ifdef FFTWPP_SINGLE_THREAD - return 0; -#else - return omp_get_thread_num(); -#endif -} - -inline int get_max_threads() -{ -#ifdef FFTWPP_SINGLE_THREAD - return 1; -#else - return omp_get_max_threads(); -#endif -} - -#ifndef FFTWPP_SINGLE_THREAD -#define PARALLEL(code) \ - if(threads > 1) { \ - _Pragma("omp parallel for num_threads(threads)") \ - code \ - } else { \ - code \ - } -#else -#define PARALLEL(code) \ - { \ - code \ - } -#endif +#include "seconds.h" +#include "parallel.h" #ifndef __Complex_h__ #include <complex> typedef std::complex<double> Complex; #endif -#include "seconds.h" #include "statistics.h" #include "align.h" @@ -97,35 +57,24 @@ extern "C" size_t fftw_alignment(); class fftw; extern "C" fftw_plan Planner(fftw *F, Complex *in, Complex *out); -void LoadWisdom(); -void SaveWisdom(); +void loadWisdom(); +void saveWisdom(); +extern std::string wisdomName; extern const char *inout; -struct threaddata { - unsigned int threads; - double mean; - double stdev; - threaddata() : threads(0), mean(0.0), stdev(0.0) {} - threaddata(unsigned int threads, double mean, double stdev) : - threads(threads), mean(mean), stdev(stdev) {} -}; - -class fftw; - class ThreadBase { -protected: - unsigned int threads; - unsigned int innerthreads; public: + size_t threads; + size_t innerthreads; ThreadBase(); - ThreadBase(unsigned int threads) : threads(threads) {} - void Threads(unsigned int nthreads) {threads=nthreads;} - unsigned int Threads() {return threads;} - unsigned int Innerthreads() {return innerthreads;} + ThreadBase(size_t threads) : threads(threads) {} + void Threads(size_t nthreads) {threads=nthreads;} + size_t Threads() {return threads;} + size_t Innerthreads() {return innerthreads;} - void multithread(unsigned int n) { + void multithread(size_t n) { if(n >= threads) { innerthreads=1; } else { @@ -133,65 +82,87 @@ public: threads=1; } } + + int get_thread_num0() { + return threads > 1 ? parallel::get_thread_num() : 0; + } }; -inline unsigned int realsize(unsigned int n, Complex *in, Complex *out=NULL) +inline size_t realsize(size_t n, bool inplace) +{ + return inplace ? 2*(n/2+1) : n; +} + +inline size_t Inplace(Complex *in, Complex *out=NULL) { - return (!out || in == out) ? 2*(n/2+1) : n; + return !out || in == out; } -inline unsigned int realsize(unsigned int n, Complex *in, double *out) +inline size_t Inplace(Complex *in, double *out) { - return realsize(n,in,(Complex *) out); + return Inplace(in,(Complex *) out); } -inline unsigned int realsize(unsigned int n, double *in, Complex *out) +inline size_t Inplace(double *in, Complex *out) { - return realsize(n,(Complex *) in,out); + return Inplace((Complex *) in,out); } +class Doubles { +public: + size_t rsize,csize; + + Doubles(size_t nx, size_t M, + size_t istride, size_t ostride, + size_t idist, size_t odist, bool inplace) { + rsize=(realsize(nx,inplace)-2)*istride+(M-1)*idist+2; + csize=2*(nx/2*ostride+(M-1)*odist+1); + if(inplace) + rsize=csize=std::max(rsize,csize); + } +}; + // Base clase for fft routines // class fftw : public ThreadBase { protected: - unsigned int doubles; // number of double precision values in dataset + size_t doubles; // number of double precision values in output int sign; - unsigned int threads; + size_t threads; double norm; fftw_plan plan; bool inplace; - unsigned int Dist(unsigned int n, size_t stride, size_t dist) { + size_t Dist(size_t n, size_t stride, size_t dist) { return dist ? dist : ((stride == 1) ? n : 1); } static const double twopi; public: - static unsigned int effort; - static unsigned int maxthreads; - static double testseconds; - static const char *WisdomName; + static size_t effort; + static size_t maxthreads; static fftw_plan (*planner)(fftw *f, Complex *in, Complex *out); + static bool wiser; - virtual unsigned int Threads() {return threads;} + virtual size_t Threads() {return threads;} static const char *oddshift; // 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; - unsigned int stop=nx*nyp; + static void Shift(Complex *data, size_t nx, size_t ny, + size_t threads) { + size_t nyp=ny/2+1; + size_t stop=nx*nyp; if(nx % 2 == 0) { - unsigned int inc=2*nyp; + size_t inc=2*nyp; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=nyp; i < stop; i += inc) { + for(size_t i=nyp; i < stop; i += inc) { Complex *p=data+i; - for(unsigned int j=0; j < nyp; j++) p[j]=-p[j]; + for(size_t j=0; j < nyp; j++) p[j]=-p[j]; } } else { std::cerr << oddshift << std::endl; @@ -200,17 +171,17 @@ public: } // Out-of-place shift of Fourier origin to (nx/2,0) for even nx. - static void Shift(double *data, unsigned int nx, unsigned int ny, - unsigned int threads) { + static void Shift(double *data, size_t nx, size_t ny, + size_t threads) { if(nx % 2 == 0) { - unsigned int stop=nx*ny; - unsigned int inc=2*ny; + size_t stop=nx*ny; + size_t inc=2*ny; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=ny; i < stop; i += inc) { + for(size_t i=ny; i < stop; i += inc) { double *p=data+i; - for(unsigned int j=0; j < ny; j++) p[j]=-p[j]; + for(size_t j=0; j < ny; j++) p[j]=-p[j]; } } else { std::cerr << oddshift << std::endl; @@ -219,20 +190,20 @@ public: } // 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; - unsigned int nyzp=ny*nzp; + static void Shift(Complex *data, size_t nx, size_t ny, + size_t nz, size_t threads) { + size_t nzp=nz/2+1; + size_t nyzp=ny*nzp; if(nx % 2 == 0 && ny % 2 == 0) { - unsigned int pinc=2*nzp; + size_t pinc=2*nzp; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; i++) { + for(size_t i=0; i < nx; i++) { Complex *pstart=data+i*nyzp; Complex *pstop=pstart+nyzp; for(Complex *p=pstart+(1-(i % 2))*nzp; p < pstop; p += pinc) { - for(unsigned int k=0; k < nzp; k++) p[k]=-p[k]; + for(size_t k=0; k < nzp; k++) p[k]=-p[k]; } } } else { @@ -242,19 +213,19 @@ public: } // Out-of-place shift of Fourier origin to (nx/2,ny/2,0) for even nx and ny. - static void Shift(double *data, unsigned int nx, unsigned int ny, - unsigned int nz, unsigned int threads) { - unsigned int nyz=ny*nz; + static void Shift(double *data, size_t nx, size_t ny, + size_t nz, size_t threads) { + size_t nyz=ny*nz; if(nx % 2 == 0 && ny % 2 == 0) { - unsigned int pinc=2*nz; + size_t pinc=2*nz; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; i++) { + for(size_t i=0; i < nx; i++) { double *pstart=data+i*nyz; double *pstop=pstart+nyz; for(double *p=pstart+(1-(i % 2))*nz; p < pstop; p += pinc) { - for(unsigned int k=0; k < nz; k++) p[k]=-p[k]; + for(size_t k=0; k < nz; k++) p[k]=-p[k]; } } } else { @@ -264,8 +235,8 @@ public: } fftw() : plan(NULL) {} - fftw(unsigned int doubles, int sign, unsigned int threads, - unsigned int n=0) : + fftw(size_t doubles, int sign, size_t threads, + size_t n=0) : doubles(doubles), sign(sign), threads(threads), norm(1.0/(n ? n : doubles/2)), plan(NULL) { #ifndef FFTWPP_SINGLE_THREAD @@ -274,7 +245,8 @@ public: } virtual ~fftw() { - if(plan) fftw_destroy_plan(plan); + if(plan) + fftw_destroy_plan(plan); } virtual fftw_plan Plan(Complex *in, Complex *out) {return NULL;}; @@ -290,72 +262,13 @@ public: exit(1); } - static void planThreads(unsigned int threads) { + static void planThreads(size_t threads) { #ifndef FFTWPP_SINGLE_THREAD omp_set_num_threads(threads); fftw_plan_with_nthreads(threads); #endif } - threaddata time(fftw_plan plan1, fftw_plan planT, Complex *in, Complex *out, - unsigned int Threads) { - utils::statistics S,ST; - double stop=utils::totalseconds()+testseconds; - threads=1; - plan=plan1; - fft(in,out); - threads=Threads; - 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 < 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 < 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; - S.clear(); - ST.clear(); - } - if(S.count() >= 10) { - double error=S.stdev(); - double diff=ST.mean()-S.mean(); - if(diff >= 0.0 || t > stop) { - threads=1; - plan=plan1; - fftw_destroy_plan(planT); - break; - } - if(diff < -error) { - threads=Threads; - fftw_destroy_plan(plan1); - break; - } - } - } - return threaddata(threads,S.mean(),S.stdev()); - } - - virtual threaddata lookup(bool inplace, unsigned int threads) { - return threaddata(); - } - virtual void store(bool inplace, const threaddata& data) {} - inline Complex *CheckAlign(Complex *in, Complex *out, bool constructor=true) { #ifndef NO_CHECK_ALIGN @@ -368,51 +281,41 @@ public: return out; } - threaddata Setup(Complex *in, Complex *out=NULL) { + void Setup(Complex *in, Complex *out=NULL) { bool alloc=!in; if(alloc) in=utils::ComplexAlign((doubles+1)/2); out=CheckAlign(in,out); inplace=(out==in); - threaddata data; - unsigned int Threads=threads; - - if(threads > 1) data=lookup(inplace,threads); - else data=threaddata(1,0.0,0.0); + parallel::Threshold(threads); + if(doubles < 2*threshold) + threads=1; - threads=data.threads > 0 ? data.threads : 1; planThreads(threads); plan=(*planner)(this,in,out); if(!plan) noplan(); - fftw_plan planT; - if(fftw::maxthreads > 1) { - threads=Threads; - planThreads(threads); - planT=(*planner)(this,in,out); - - if(data.threads == 0) { - if(planT) - data=time(plan,planT,in,out,threads); - else noplan(); - store(inplace,threaddata(threads,data.mean,data.stdev)); - } - } - if(alloc) Array::deleteAlign(in,(doubles+1)/2); #ifdef FFTWPP_VERBOSE if(threads > 1) std::cout << "Using " << threads << " threads." << std::endl; #endif - return data; } - threaddata Setup(Complex *in, double *out) { - return Setup(in,(Complex *) out); + void Setup(Complex *in, double *out) { + parallel::Threshold(threads); + if(doubles < 4*threshold) + threads=1; + + Setup(in,(Complex *) out); } - threaddata Setup(double *in, Complex *out=NULL) { - return Setup((Complex *) in,out); + void Setup(double *in, Complex *out=NULL) { + parallel::Threshold(threads); + if(doubles < 4*threshold) + threads=1; + + Setup((Complex *) in,out); } virtual void Execute(Complex *in, Complex *out, bool=false) { @@ -455,51 +358,67 @@ public: } void Normalize(Complex *out) { - unsigned int stop=doubles/2; + size_t stop=doubles/2; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < stop; i++) out[i] *= norm; + for(size_t i=0; i < stop; i++) out[i] *= norm; } void Normalize(double *out) { #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < doubles; i++) out[i] *= norm; + for(size_t i=0; i < doubles; i++) out[i] *= norm; } - virtual void fftNormalized(Complex *in, Complex *out=NULL, bool shift=false) + void fftNormalized(Complex *in, Complex *out=NULL) { out=Setout(in,out); - Execute(in,out,shift); + Execute(in,out); Normalize(out); } - virtual void fftNormalized(Complex *in, double *out, bool shift=false) { + void fftNormalized(Complex *in, double *out) { out=(double *) Setout(in,(Complex *) out); - Execute(in,(Complex *) out,shift); + Execute(in,(Complex *) out); Normalize(out); } - virtual void fftNormalized(double *in, Complex *out, bool shift=false) { - fftNormalized((Complex *) in,out,shift); + void fftNormalized(double *in, Complex *out) { + out=Setout((Complex *) in,out); + Execute((Complex *) in,out); + Normalize(out); } - template<class I, class O> - void fft0Normalized(I in, O out) { - fftNormalized(in,out,true); + void fft0Normalized(Complex *in, Complex *out=NULL) + { + out=Setout(in,out); + Execute(in,out,true); + Normalize(out); + } + + void fft0Normalized(Complex *in, double *out) { + out=(double *) Setout(in,(Complex *) out); + Execute(in,(Complex *) out,true); + Normalize(out); + } + + void fft0Normalized(double *in, Complex *out) { + out=Setout((Complex *) in,out); + Execute((Complex *) in,out,true); + Normalize(out); } template<class O> - void Normalize(unsigned int nx, unsigned int M, size_t ostride, + void Normalize(size_t nx, size_t M, size_t ostride, size_t odist, O *out) { - unsigned int stop=nx*ostride; + size_t stop=nx*ostride; O *outMdist=out+M*odist; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < stop; i += ostride) { + for(size_t i=0; i < stop; i += ostride) { O *pstop=outMdist+i; for(O *p=out+i; p < pstop; p += odist) { *p *= norm; @@ -508,10 +427,10 @@ public: } template<class I, class O> - void fftNormalized(unsigned int nx, unsigned int M, size_t ostride, - size_t odist, I *in, O *out=NULL, bool shift=false) { + void fftNormalized(size_t nx, size_t M, size_t ostride, + size_t odist, I *in, O *out) { out=(O *) Setout((Complex *) in,(Complex *) out); - Execute((Complex *) in,(Complex *) out,shift); + Execute((Complex *) in,(Complex *) out); Normalize(nx,M,ostride,odist,out); } @@ -522,9 +441,9 @@ class Transpose { bool inplace; public: template<class T> - Transpose(unsigned int rows, unsigned int cols, unsigned int length, - T *in, T *out=NULL, unsigned int threads=fftw::maxthreads) { - unsigned int size=sizeof(T); + Transpose(size_t rows, size_t cols, size_t length, + T *in, T *out=NULL, size_t threads=fftw::maxthreads) { + size_t size=sizeof(T); if(size % sizeof(double) != 0) { std::cerr << "ERROR: Transpose is not implemented for type of size " << size; @@ -537,6 +456,10 @@ public: size /= sizeof(double); length *= size; + parallel::Threshold(threads); + if(length*rows*cols/2 < threshold) + threads=1; + fftw::planThreads(threads); fftw_iodim dims[3]; @@ -559,7 +482,8 @@ public: } ~Transpose() { - if(plan) fftw_destroy_plan(plan); + if(plan) + fftw_destroy_plan(plan); } template<class T> @@ -577,74 +501,35 @@ public: template<class T, class L> class Threadtable { public: - typedef std::map<T,threaddata,L> Table; + typedef std::map<T,size_t,L> Table; - threaddata Lookup(Table& table, T key) { + size_t Lookup(Table& table, T key) { typename Table::iterator p=table.find(key); - return p == table.end() ? threaddata() : p->second; + return p == table.end() ? 0 : p->second; } - void Store(Table& threadtable, T key, const threaddata& data) { - threadtable[key]=data; + void Store(Table& threadtable, T key, size_t t) { + threadtable[key]=t; } }; -struct keytype1 { - unsigned int nx; - unsigned int threads; +struct keytype { + size_t nx; + size_t M; + size_t threads; bool inplace; - keytype1(unsigned int nx, unsigned int threads, bool inplace) : - nx(nx), threads(threads), inplace(inplace) {} + keytype(size_t nx, size_t M, size_t threads, + bool inplace) : + nx(nx), M(M), threads(threads), inplace(inplace) {} }; -struct keyless1 { - bool operator()(const keytype1& a, const keytype1& b) const { +struct keyless { + bool operator()(const keytype& a, const keytype& b) const { return a.nx < b.nx || (a.nx == b.nx && - (a.threads < b.threads || (a.threads == b.threads && - a.inplace < b.inplace))); - } -}; - -struct keytype2 { - unsigned int nx; - unsigned int ny; - unsigned int threads; - bool inplace; - keytype2(unsigned int nx, unsigned int ny, unsigned int threads, - bool inplace) : - nx(nx), ny(ny), threads(threads), inplace(inplace) {} -}; - -struct keyless2 { - bool operator()(const keytype2& a, const keytype2& b) const { - return a.nx < b.nx || (a.nx == b.nx && - (a.ny < b.ny || (a.ny == b.ny && - (a.threads < b.threads || - (a.threads == b.threads && - a.inplace < b.inplace))))); - } -}; - -struct keytype3 { - unsigned int nx; - unsigned int ny; - unsigned int nz; - unsigned int threads; - bool inplace; - keytype3(unsigned int nx, unsigned int ny, unsigned int nz, - unsigned int threads, bool inplace) : - nx(nx), ny(ny), nz(nz), threads(threads), inplace(inplace) {} -}; - -struct keyless3 { - bool operator()(const keytype3& a, const keytype3& b) const { - return a.nx < b.nx || (a.nx == b.nx && - (a.ny < b.ny || (a.ny == b.ny && - (a.nz < b.nz || - (a.nz == b.nz && + (a.M < b.M || (a.M == b.M && (a.threads < b.threads || (a.threads == b.threads && - a.inplace < b.inplace))))))); + a.inplace < b.inplace))))); } }; @@ -671,28 +556,20 @@ struct keyless3 { // fft1d Backward(n,1); // Backward.fft(in); // -class fft1d : public fftw, public Threadtable<keytype1,keyless1> { - unsigned int nx; - static Table threadtable; +class fft1d : public fftw { + size_t nx; public: - fft1d(unsigned int nx, int sign, Complex *in=NULL, Complex *out=NULL, - unsigned int threads=maxthreads) + fft1d(size_t nx, int sign, Complex *in=NULL, Complex *out=NULL, + size_t threads=maxthreads) : fftw(2*nx,sign,threads), nx(nx) {Setup(in,out);} #ifdef __Array_h__ fft1d(int sign, const Array::array1<Complex>& in, const Array::array1<Complex>& out=Array::NULL1, - unsigned int threads=maxthreads) + size_t threads=maxthreads) : fftw(2*in.Nx(),sign,threads), nx(in.Nx()) {Setup(in,out);} #endif - threaddata lookup(bool inplace, unsigned int threads) { - return this->Lookup(threadtable,keytype1(nx,threads,inplace)); - } - void store(bool inplace, const threaddata& data) { - this->Store(threadtable,keytype1(nx,data.threads,inplace),data); - } - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_1d(nx,(fftw_complex *) in,(fftw_complex *) out, sign,effort); @@ -700,62 +577,126 @@ public: }; template<class I, class O> -class fftwblock : public virtual fftw { +class fftwblock : public virtual fftw, + public virtual Threadtable<keytype,keyless> { public: int nx; - unsigned int M; + size_t M; size_t istride,ostride; size_t idist,odist; fftw_plan plan1,plan2; - unsigned int T,Q,R; - fftwblock(unsigned int nx, unsigned int M, - size_t istride, size_t ostride, size_t idist, size_t odist, - Complex *in, Complex *out, unsigned int Threads) + size_t T,Q,R; + fftwblock() : plan1(NULL), plan2(NULL) {} + + fftwblock(size_t nx, size_t M, + size_t istride, size_t ostride, size_t idist, size_t odist) : fftw(), nx(nx), M(M), istride(istride), ostride(ostride), idist(Dist(nx,istride,idist)), odist(Dist(nx,ostride,odist)), - plan1(NULL), plan2(NULL) { + plan1(NULL), plan2(NULL) {} + + void init(Complex *in, Complex *out, size_t Threads, + Table& threadtable) { T=1; Q=M; R=0; - threaddata S1=Setup(in,out); - fftw_plan planT1=plan; - threads=S1.threads; - bool hermitian=typeid(I) == typeid(double) || typeid(O) == typeid(double); + if(Threads > M && M > 1) Threads=M; - if(fftw::maxthreads > 1 && (!hermitian || ostride*(nx/2+1) < idist)) { - if(Threads > 1) { - T=std::min(M,Threads); + threads=Threads; + Setup(in,out); + Threads=threads; + + size_t T0=Threads; + if(T0 > 1) { + size_t nxp=nx/2+1; + size_t olength=0; + size_t ilength=0; + if(typeid(I) == typeid(double)) { + ilength=nx; + olength=nxp; + } + if(typeid(O) == typeid(double)) { + ilength=nxp; + olength=nx; + } + if(!inplace || + (ostride*olength*sizeof(O) <= idist*sizeof(I) && + odist*sizeof(O) >= istride*ilength*sizeof(I))) { + T=T0; Q=T > 0 ? M/T : 0; R=M-Q*T; - threads=Threads; - threaddata ST=Setup(in,out); - - if(R > 0 && threads == 1 && plan1 != plan2) { - fftw_destroy_plan(plan2); - plan2=plan1; + size_t data=Lookup(threadtable,keytype(nx,M,Threads,inplace)); + if(data == 1) + T0=1; + else { + fftw_plan planFFTW=plan; + threads=1; + Setup(in,out); + plan1=plan; + if(data == T) { + plan=NULL; + return; + } + plan=planFFTW; } + } else T0=1; + } - if(ST.mean > S1.mean-S1.stdev) { // Use FFTW's multi-threading - fftw_destroy_plan(plan); - if(R > 0) { - fftw_destroy_plan(plan2); - plan2=NULL; - } - T=1; - Q=M; - R=0; - plan=planT1; - } else { // Do the multi-threading ourselves - fftw_destroy_plan(planT1); - threads=ST.threads; + if(T0 == 1 || time(in,out)) { // Use FFTW's multithreading + T=1; + if(plan1) { + fftw_destroy_plan(plan1); + plan1=NULL; + if(plan2) { + fftw_destroy_plan(plan2); + plan2=NULL; } - } else - Setup(in,out); // Synchronize wisdom + threads=Threads; + Store(threadtable,keytype(nx,M,Threads,inplace),T); + } + } else { // Do the multithreading ourselves + T=T0; + threads=T; + Store(threadtable,keytype(nx,M,Threads,inplace),T); } } + bool time(Complex *in, Complex *out) { + utils::statistics S(true),ST(true); + utils::statistics medianS(true),medianST(true); + + double eps=0.02; + + size_t T0=T; + + do { + T=1; // FFTW + utils::cpuTimer C; + inplace ? fftNormalized(in,out) : fft(in,out); + S.add(C.nanoseconds()); + + T=T0; // BLOCK + utils::cpuTimer CT; + inplace ? fftNormalized(in,out) : fft(in,out); + ST.add(CT.nanoseconds()); + + if(S.count() >= 4 && ST.min() >= S.max()) + return true; + + if(S.count() >= 4 && S.min() >= ST.max()) + return false; + + medianS.add(S.median()); + medianST.add(ST.median()); + + } while(S.count() < 5 || medianS.stderror() > eps*medianS.mean() || + medianST.stderror() > eps*medianST.mean()); + + return S.median() <= ST.median(); + } + + fftw_plan Plan(int Q, fftw_complex *in, fftw_complex *out) { return fftw_plan_many_dft(1,&nx,Q,in,NULL,istride,idist, out,NULL,ostride,odist,sign,effort); @@ -775,7 +716,6 @@ public: if(R > 0) { plan2=Plan(Q+1,(I *) in,(O *) out); if(!plan2) return NULL; - if(threads == 1) plan1=plan2; } return Plan(Q,(I *) in,(O *) out); } @@ -796,26 +736,62 @@ public: if(T == 1) Execute(plan,(I *) in,(O *) out); else { - unsigned int extra=T-R; + size_t extra=T-R; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(T) #endif - for(unsigned int i=0; i < T; ++i) { - unsigned int iQ=i*Q; + for(size_t i=0; i < T; ++i) { + size_t iQ=i*Q; if(i < extra) - Execute(plan,(I *) in+iQ*idist,(O *) out+iQ*odist); + Execute(plan1,(I *) in+iQ*idist,(O *) out+iQ*odist); else { - unsigned int offset=iQ+i-extra; + size_t offset=iQ+i-extra; Execute(plan2,(I *) in+offset*idist,(O *) out+offset*odist); } } } } - unsigned int Threads() {return std::max(T,threads);} + size_t Threads() {return std::max(T,threads);} ~fftwblock() { - if(plan2) fftw_destroy_plan(plan2); + if(plan1) + fftw_destroy_plan(plan1); + if(plan2) + fftw_destroy_plan(plan2); + } +}; + +class Mfft1d : public fftwblock<fftw_complex,fftw_complex>, + public virtual Threadtable<keytype,keyless> { + static Table threadtable; +public: + Mfft1d(size_t nx, int sign, size_t M=1, + Complex *in=NULL, Complex *out=NULL, + size_t threads=maxthreads) : + fftw(2*((nx-1)+(M-1)*nx+1),sign,threads,nx), + fftwblock<fftw_complex,fftw_complex>(nx,M,1,1,nx,nx) { + init(in,out,threads,threadtable); + } + + Mfft1d(size_t nx, int sign, size_t M, size_t stride=1, + size_t dist=0, Complex *in=NULL, Complex *out=NULL, + size_t threads=maxthreads) : + fftw(2*((nx-1)*stride+(M-1)*Dist(nx,stride,dist)+1),sign,threads,nx), + fftwblock<fftw_complex,fftw_complex> + (nx,M,stride,stride,dist,dist) { + init(in,out,threads,threadtable); + } + + Mfft1d(size_t nx, int sign, size_t M, + size_t istride, size_t ostride, size_t idist, size_t odist, + Complex *in, Complex *out, size_t 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), + fftwblock<fftw_complex,fftw_complex>(nx,M,istride,ostride,idist,odist) + { + init(in,out,threads,threadtable); } }; @@ -844,38 +820,74 @@ public: // dist is the spacing between the first elements of the vectors. // // -class mfft1d : public fftwblock<fftw_complex,fftw_complex>, - public Threadtable<keytype3,keyless3> { - static Table threadtable; +class mfft1d { + bool single; + fft1d *fft1; + Mfft1d *fftm; public: - mfft1d(unsigned int nx, int sign, unsigned int M=1, + mfft1d(size_t nx, int sign, size_t 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) {} + size_t threads=fftw::maxthreads) : single(M == 1) { + if(single) + fft1=new fft1d(nx,sign,in,out,threads); + else + fftm=new Mfft1d(nx,sign,M,in,out,threads); + } - mfft1d(unsigned int nx, int sign, unsigned int M, size_t stride=1, + mfft1d(size_t nx, int sign, size_t 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), - fftwblock<fftw_complex,fftw_complex> - (nx,M,stride,stride,dist,dist,in,out,threads) {} + size_t threads=fftw::maxthreads) : + single(M == 1 && stride == 1) { + if(single) + fft1=new fft1d(nx,sign,in,out,threads); + else + fftm=new Mfft1d(nx,sign,M,stride,dist,in,out,threads); + } - mfft1d(unsigned int nx, int sign, unsigned int M, + mfft1d(size_t nx, int sign, size_t M, size_t istride, size_t ostride, size_t idist, size_t odist, - 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), - fftwblock<fftw_complex,fftw_complex>(nx,M,istride,ostride,idist,odist,in, - out,threads) {} + Complex *in, Complex *out, size_t threads=fftw::maxthreads) : + single(M == 1 && istride == 1 && ostride == 1) { + if(single) + fft1=new fft1d(nx,sign,in,out,threads); + else + fftm=new Mfft1d(nx,sign,M,istride,ostride,idist,odist,in,out,threads); + } + + size_t Threads() { + return single ? fft1->Threads() : fftm->Threads(); + } - threaddata lookup(bool inplace, unsigned int threads) { - return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace)); + template<class I> + void fft(I in) { + single ? fft1->fft(in) : fftm->fft(in); } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data); + + template<class I, class O> + void fft(I in, O out) { + single ? fft1->fft(in,out) : fftm->fft(in,out); + } + + template<class I> + void fftNormalized(I in) { + single ? fft1->fftNormalized(in) : fftm->fftNormalized(in); + } + + template<class I, class O> + void fftNormalized(I in, O out) { + single ? fft1->fftNormalized(in,out) : fftm->fftNormalized(in,out); + } + + template<class O> + void Normalize(O out) { + single ? fft1->Normalize(out) : fftm->Normalize(out); + } + + ~mfft1d() { + if(single) + delete fft1; + else + delete fftm; } }; @@ -898,24 +910,16 @@ public: // in contains the n real values stored as a Complex array; // out contains the first n/2+1 Complex Fourier values. // -class rcfft1d : public fftw, public Threadtable<keytype1,keyless1> { - unsigned int nx; - static Table threadtable; +class rcfft1d : public fftw { + size_t nx; public: - rcfft1d(unsigned int nx, Complex *out=NULL, unsigned int threads=maxthreads) + rcfft1d(size_t nx, Complex *out=NULL, size_t threads=maxthreads) : fftw(2*(nx/2+1),-1,threads,nx), nx(nx) {Setup(out,(double*) NULL);} - rcfft1d(unsigned int nx, double *in, Complex *out=NULL, - unsigned int threads=maxthreads) + rcfft1d(size_t nx, double *in, Complex *out=NULL, + size_t threads=maxthreads) : fftw(2*(nx/2+1),-1,threads,nx), nx(nx) {Setup(in,out);} - threaddata lookup(bool inplace, unsigned int threads) { - return Lookup(threadtable,keytype1(nx,threads,inplace)); - } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype1(nx,data.threads,inplace),data); - } - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_r2c_1d(nx,(double *) in,(fftw_complex *) out, effort); } @@ -946,23 +950,15 @@ public: // in contains the first n/2+1 Complex Fourier values. // out contains the n real values stored as a Complex array; // -class crfft1d : public fftw, public Threadtable<keytype1,keyless1> { - unsigned int nx; - static Table threadtable; +class crfft1d : public fftw { + size_t nx; public: - crfft1d(unsigned int nx, double *out=NULL, unsigned int threads=maxthreads) + crfft1d(size_t nx, double *out=NULL, size_t threads=maxthreads) : fftw(2*(nx/2+1),1,threads,nx), nx(nx) {Setup(out);} - crfft1d(unsigned int nx, Complex *in, double *out=NULL, - unsigned int threads=maxthreads) - : fftw(realsize(nx,in,out),1,threads,nx), nx(nx) {Setup(in,out);} - - threaddata lookup(bool inplace, unsigned int threads) { - return Lookup(threadtable,keytype1(nx,threads,inplace)); - } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype1(nx,data.threads,inplace),data); - } + crfft1d(size_t nx, Complex *in, double *out=NULL, + size_t threads=maxthreads) + : fftw(realsize(nx,Inplace(in,out)),1,threads,nx), nx(nx) {Setup(in,out);} fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_c2r_1d(nx,(fftw_complex *) in,(double *) out,effort); @@ -973,6 +969,29 @@ public: } }; +class Mrcfft1d : public fftwblock<double,fftw_complex>, + public virtual Threadtable<keytype,keyless> { + static Table threadtable; +public: + Mrcfft1d(size_t nx, size_t M, size_t istride, size_t ostride, + size_t idist, size_t odist, double *in=NULL, Complex *out=NULL, + size_t threads=maxthreads) + : fftw(Doubles(nx,M,istride,ostride,idist,odist,Inplace(in,out)).csize, + -1,threads,nx), + fftwblock<double,fftw_complex> + (nx,M,istride,ostride,idist,odist) { + init((Complex *) in,out,threads,threadtable); + } + + void Normalize(Complex *out) { + fftw::Normalize<Complex>(nx/2+1,M,ostride,odist,out); + } + + void fftNormalized(double *in, Complex *out=NULL) { + fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out); + } +}; + // Compute the real Fourier transform of M real vectors, each of length n, // using phase sign -1. Before calling fft(), the array in must be // allocated as double[M*n] and the array out must be allocated as @@ -997,38 +1016,77 @@ public: // in contains the n real values stored as a Complex array; // out contains the first n/2+1 Complex Fourier values. // -class mrcfft1d : public fftwblock<double,fftw_complex>, - public Threadtable<keytype3,keyless3> { - static Table threadtable; +class mrcfft1d { + bool single; + rcfft1d *fft1; + Mrcfft1d *fftm; public: - mrcfft1d(unsigned int nx, unsigned int M, - size_t istride, size_t ostride, - size_t idist, size_t odist, - double *in=NULL, Complex *out=NULL, - unsigned int threads=maxthreads) - : fftw(std::max((realsize(nx,in,out)-2)*istride+(M-1)*idist+2, - 2*(nx/2*ostride+(M-1)*odist+1)),-1,threads,nx), - fftwblock<double,fftw_complex> - (nx,M,istride,ostride,idist,odist,(Complex *) in,out,threads) {} + mrcfft1d(size_t nx, size_t M, size_t istride, size_t ostride, + size_t idist, size_t odist, double *in=NULL, Complex *out=NULL, + size_t threads=fftw::maxthreads) : + single(M == 1 && istride == 1 && ostride == 1) { + if(single) + fft1=new rcfft1d(nx,in,out,threads); + else + fftm=new Mrcfft1d(nx,M,istride,ostride,idist,odist,in,out,threads); + } - threaddata lookup(bool inplace, unsigned int threads) { - return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace)); + size_t Threads() { + return single ? fft1->Threads() : fftm->Threads(); } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data); + template<class I> + void fft(I in) { + single ? fft1->fft(in) : fftm->fft(in); + } + + template<class I, class O> + void fft(I in, O out) { + single ? fft1->fft(in,out) : fftm->fft(in,out); } void Normalize(Complex *out) { - fftw::Normalize<Complex>(nx/2+1,M,ostride,odist,out); + single ? fft1->Normalize(out) : fftm->Normalize(out); + } + + template<class I> + void fftNormalized(I in) { + single ? fft1->fftNormalized(in) : fftm->fftNormalized(in); } - void fftNormalized(double *in, Complex *out=NULL, bool shift=false) { - fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out,false); + template<class I, class O> + void fftNormalized(I in, O out=NULL) { + single ? fft1->fftNormalized(in,out) : fftm->fftNormalized(in,out); } - void fft0Normalized(double *in, Complex *out=NULL) { - fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out,true); + ~mrcfft1d() { + if(single) + delete fft1; + else + delete fftm; + } +}; + +class Mcrfft1d : public fftwblock<fftw_complex,double>, + public virtual Threadtable<keytype,keyless> { + static Table threadtable; +public: + Mcrfft1d(size_t nx, size_t M, size_t istride, size_t ostride, + size_t idist, size_t odist, Complex *in=NULL, double *out=NULL, + size_t threads=maxthreads) + : fftw(Doubles(nx,M,ostride,istride,odist,idist,Inplace(in,out)).rsize, + 1,threads,nx), + fftwblock<fftw_complex,double> + (nx,M,istride,ostride,idist,odist) { + init(in,(Complex *) out,threads,threadtable); + } + + void Normalize(double *out) { + fftw::Normalize<double>(nx,M,ostride,odist,out); + } + + void fftNormalized(Complex *in, double *out=NULL) { + fftw::fftNormalized<Complex,double>(nx,M,ostride,odist,in,out); } }; @@ -1055,36 +1113,54 @@ public: // in contains the first n/2+1 Complex Fourier values; // out contains the n real values stored as a Complex array. // -class mcrfft1d : public fftwblock<fftw_complex,double>, - public Threadtable<keytype3,keyless3> { - static Table threadtable; +class mcrfft1d { + bool single; + crfft1d *fft1; + Mcrfft1d *fftm; public: - mcrfft1d(unsigned int nx, unsigned int M, size_t istride, size_t ostride, + mcrfft1d(size_t nx, size_t M, size_t istride, size_t ostride, size_t idist, size_t odist, Complex *in=NULL, double *out=NULL, - unsigned int threads=maxthreads) - : fftw(std::max(2*(nx/2*istride+(M-1)*idist+1), - (realsize(nx,in,out)-2)*ostride+(M-1)*odist+2),1,threads,nx), - fftwblock<fftw_complex,double> - (nx,M,istride,ostride,idist,odist,in,(Complex *) out,threads) {} + size_t threads=fftw::maxthreads) : + single(M == 1 && istride == 1 && ostride == 1) { + if(single) + fft1=new crfft1d(nx,in,out,threads); + else + fftm=new Mcrfft1d(nx,M,istride,ostride,idist,odist,in,out,threads); + } + + size_t Threads() { + return single ? fft1->Threads() : fftm->Threads(); + } - threaddata lookup(bool inplace, unsigned int threads) { - return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace)); + template<class I> + void fft(I in) { + single ? fft1->fft(in) : fftm->fft(in); } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data); + template<class I, class O> + void fft(I in, O out) { + single ? fft1->fft(in,out) : fftm->fft(in,out); } void Normalize(double *out) { - fftw::Normalize<double>(nx,M,ostride,odist,out); + single ? fft1->Normalize(out) : fftm->Normalize(out); } - void fftNormalized(Complex *in, double *out=NULL, bool shift=false) { - fftw::fftNormalized<Complex,double>(nx,M,ostride,odist,in,out,false); + template<class I> + void fftNormalized(I in) { + single ? fft1->fftNormalized(in) : fftm->fftNormalized(in); } - void fft0Normalized(Complex *in, double *out=NULL) { - fftw::fftNormalized<Complex,double>(nx,M,ostride,odist,in,out,true); + template<class I, class O> + void fftNormalized(I in, O out=NULL) { + single ? fft1->fftNormalized(in,out) : fftm->fftNormalized(in,out); + } + + ~mcrfft1d() { + if(single) + delete fft1; + else + delete fftm; } }; @@ -1114,31 +1190,23 @@ public: // Note: // in[ny*i+j] contains the ny Complex values for each i=0,...,nx-1. // -class fft2d : public fftw, public Threadtable<keytype2,keyless2> { - unsigned int nx; - unsigned int ny; - static Table threadtable; +class fft2d : public fftw { + size_t nx; + size_t ny; public: - fft2d(unsigned int nx, unsigned int ny, int sign, Complex *in=NULL, - Complex *out=NULL, unsigned int threads=maxthreads) + fft2d(size_t nx, size_t ny, int sign, Complex *in=NULL, + Complex *out=NULL, size_t threads=maxthreads) : fftw(2*nx*ny,sign,threads), nx(nx), ny(ny) {Setup(in,out);} #ifdef __Array_h__ fft2d(int sign, const Array::array2<Complex>& in, const Array::array2<Complex>& out=Array::NULL2, - unsigned int threads=maxthreads) + size_t threads=maxthreads) : fftw(2*in.Size(),sign,threads), nx(in.Nx()), ny(in.Ny()) { Setup(in,out); } #endif - threaddata lookup(bool inplace, unsigned int threads) { - return this->Lookup(threadtable,keytype2(nx,ny,threads,inplace)); - } - void store(bool inplace, const threaddata& data) { - this->Store(threadtable,keytype2(nx,ny,data.threads,inplace),data); - } - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_2d(nx,ny,(fftw_complex *) in,(fftw_complex *) out, sign,effort); @@ -1173,15 +1241,15 @@ public: // out contains the upper-half portion (ky >= 0) of the Complex transform. // class rcfft2d : public fftw { - unsigned int nx; - unsigned int ny; + size_t nx; + size_t ny; public: - rcfft2d(unsigned int nx, unsigned int ny, Complex *out=NULL, - unsigned int threads=maxthreads) + rcfft2d(size_t nx, size_t ny, Complex *out=NULL, + size_t threads=maxthreads) : fftw(2*nx*(ny/2+1),-1,threads,nx*ny), nx(nx), ny(ny) {Setup(out);} - rcfft2d(unsigned int nx, unsigned int ny, double *in, Complex *out=NULL, - unsigned int threads=maxthreads) + rcfft2d(size_t nx, size_t ny, double *in, Complex *out=NULL, + size_t threads=maxthreads) : fftw(2*nx*(ny/2+1),-1,threads,nx*ny), nx(nx), ny(ny) { Setup(in,out); } @@ -1201,18 +1269,18 @@ public: // Set Nyquist modes of even shifted transforms to zero. void deNyquist(Complex *f) { - unsigned int nyp=ny/2+1; + size_t nyp=ny/2+1; if(nx % 2 == 0) #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int j=0; j < nyp; ++j) + for(size_t j=0; j < nyp; ++j) f[j]=0.0; if(ny % 2 == 0) #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; ++i) + for(size_t i=0; i < nx; ++i) f[(i+1)*nyp-1]=0.0; } }; @@ -1242,16 +1310,16 @@ public: // out contains the nx*ny real values stored as a Complex array. // class crfft2d : public fftw { - unsigned int nx; - unsigned int ny; + size_t nx; + size_t ny; public: - crfft2d(unsigned int nx, unsigned int ny, double *out=NULL, - unsigned int threads=maxthreads) : + crfft2d(size_t nx, size_t ny, double *out=NULL, + size_t threads=maxthreads) : fftw(2*nx*(ny/2+1),1,threads,nx*ny), nx(nx), ny(ny) {Setup(out);} - crfft2d(unsigned int nx, unsigned int ny, Complex *in, double *out=NULL, - unsigned int threads=maxthreads) - : fftw(nx*realsize(ny,in,out),1,threads,nx*ny), nx(nx), ny(ny) { + crfft2d(size_t nx, size_t ny, Complex *in, double *out=NULL, + size_t threads=maxthreads) + : fftw(nx*realsize(ny,Inplace(in,out)),1,threads,nx*ny), nx(nx), ny(ny) { Setup(in,out); } @@ -1270,18 +1338,18 @@ public: // Set Nyquist modes of even shifted transforms to zero. void deNyquist(Complex *f) { - unsigned int nyp=ny/2+1; + size_t nyp=ny/2+1; if(nx % 2 == 0) #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int j=0; j < nyp; ++j) + for(size_t j=0; j < nyp; ++j) f[j]=0.0; if(ny % 2 == 0) #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; ++i) + for(size_t i=0; i < nx; ++i) f[(i+1)*nyp-1]=0.0; } }; @@ -1314,19 +1382,19 @@ public: // indexed by i=0,...,nx-1, j=0,...,ny-1, and k=0,...,nz-1. // class fft3d : public fftw { - unsigned int nx; - unsigned int ny; - unsigned int nz; + size_t nx; + size_t ny; + size_t nz; public: - fft3d(unsigned int nx, unsigned int ny, unsigned int nz, + fft3d(size_t nx, size_t ny, size_t nz, int sign, Complex *in=NULL, Complex *out=NULL, - unsigned int threads=maxthreads) + size_t threads=maxthreads) : fftw(2*nx*ny*nz,sign,threads), nx(nx), ny(ny), nz(nz) {Setup(in,out);} #ifdef __Array_h__ fft3d(int sign, const Array::array3<Complex>& in, const Array::array3<Complex>& out=Array::NULL3, - unsigned int threads=maxthreads) + size_t threads=maxthreads) : fftw(2*in.Size(),sign,threads), nx(in.Nx()), ny(in.Ny()), nz(in.Nz()) {Setup(in,out);} #endif @@ -1360,18 +1428,18 @@ public: // out contains the upper-half portion (kz >= 0) of the Complex transform. // class rcfft3d : public fftw { - unsigned int nx; - unsigned int ny; - unsigned int nz; + size_t nx; + size_t ny; + size_t nz; public: - rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *out=NULL, - unsigned int threads=maxthreads) + rcfft3d(size_t nx, size_t ny, size_t nz, Complex *out=NULL, + size_t threads=maxthreads) : fftw(2*nx*ny*(nz/2+1),-1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz) { Setup(out); } - rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, double *in, - Complex *out=NULL, unsigned int threads=maxthreads) + rcfft3d(size_t nx, size_t ny, size_t nz, double *in, + Complex *out=NULL, size_t threads=maxthreads) : fftw(2*nx*ny*(nz/2+1),-1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz) {Setup(in,out);} @@ -1390,13 +1458,13 @@ public: // Set Nyquist modes of even shifted transforms to zero. void deNyquist(Complex *f) { - unsigned int nzp=nz/2+1; - unsigned int yz=ny*nzp; + size_t nzp=nz/2+1; + size_t yz=ny*nzp; if(nx % 2 == 0) { #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int k=0; k < yz; ++k) + for(size_t k=0; k < yz; ++k) f[k]=0.0; } @@ -1404,9 +1472,9 @@ public: #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; ++i) { - unsigned int iyz=i*yz; - for(unsigned int k=0; k < nzp; ++k) + for(size_t i=0; i < nx; ++i) { + size_t iyz=i*yz; + for(size_t k=0; k < nzp; ++k) f[iyz+k]=0.0; } } @@ -1415,8 +1483,8 @@ public: #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; ++i) - for(unsigned int j=0; j < ny; ++j) + for(size_t i=0; i < nx; ++i) + for(size_t j=0; j < ny; ++j) f[i*yz+(j+1)*nzp-1]=0.0; } }; @@ -1446,19 +1514,19 @@ public: // out contains the nx*ny*nz real values stored as a Complex array. // class crfft3d : public fftw { - unsigned int nx; - unsigned int ny; - unsigned int nz; + size_t nx; + size_t ny; + size_t nz; public: - crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, double *out=NULL, - unsigned int threads=maxthreads) + crfft3d(size_t nx, size_t ny, size_t nz, double *out=NULL, + size_t threads=maxthreads) : fftw(2*nx*ny*(nz/2+1),1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz) {Setup(out);} - crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in, - double *out=NULL, unsigned int threads=maxthreads) - : fftw(nx*ny*(realsize(nz,in,out)),1,threads,nx*ny*nz), nx(nx), ny(ny), - nz(nz) {Setup(in,out);} + crfft3d(size_t nx, size_t ny, size_t nz, Complex *in, + double *out=NULL, size_t threads=maxthreads) + : fftw(nx*ny*(realsize(nz,Inplace(in,out))),1,threads,nx*ny*nz), + nx(nx), ny(ny), nz(nz) {Setup(in,out);} fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_c2r_3d(nx,ny,nz,(fftw_complex *) in,(double *) out, @@ -1475,13 +1543,13 @@ public: // Set Nyquist modes of even shifted transforms to zero. void deNyquist(Complex *f) { - unsigned int nzp=nz/2+1; - unsigned int yz=ny*nzp; + size_t nzp=nz/2+1; + size_t yz=ny*nzp; if(nx % 2 == 0) { #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int k=0; k < yz; ++k) + for(size_t k=0; k < yz; ++k) f[k]=0.0; } @@ -1489,9 +1557,9 @@ public: #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; ++i) { - unsigned int iyz=i*yz; - for(unsigned int k=0; k < nzp; ++k) + for(size_t i=0; i < nx; ++i) { + size_t iyz=i*yz; + for(size_t k=0; k < nzp; ++k) f[iyz+k]=0.0; } } @@ -1500,8 +1568,8 @@ public: #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif - for(unsigned int i=0; i < nx; ++i) - for(unsigned int j=0; j < ny; ++j) + for(size_t i=0; i < nx; ++i) + for(size_t j=0; j < ny; ++j) f[i*yz+(j+1)*nzp-1]=0.0; } }; |