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++.h923
1 files changed, 692 insertions, 231 deletions
diff --git a/Build/source/utils/asymptote/fftw++.h b/Build/source/utils/asymptote/fftw++.h
index addb91cff05..35dfc6f249b 100644
--- a/Build/source/utils/asymptote/fftw++.h
+++ b/Build/source/utils/asymptote/fftw++.h
@@ -1,5 +1,5 @@
/* Fast Fourier transform C++ header class for the FFTW3 Library
- Copyright (C) 2004-12 John C. Bowman, University of Alberta
+ Copyright (C) 2004-13 John C. Bowman, University of Alberta
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
@@ -18,13 +18,14 @@
#ifndef __fftwpp_h__
#define __fftwpp_h__ 1
-#define __FFTWPP_H_VERSION__ 1.14
+#define __FFTWPP_H_VERSION__ 1.14pre
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <fftw3.h>
#include <cerrno>
+#include <map>
#ifndef FFTWPP_SINGLE_THREAD
#include <omp.h>
@@ -130,8 +131,7 @@ inline void newAlign(T *&v, size_t len, size_t align)
template<class T>
inline void deleteAlign(T *v, size_t len)
{
- for(size_t i=len-1; i > 0; i--) v[i].~T();
- v[0].~T();
+ for(size_t i=len; i-- > 0;) v[i].~T();
#ifdef HAVE_POSIX_MEMALIGN
free(v);
#else
@@ -187,12 +187,24 @@ inline int fftwpp_import_wisdom(int (*g)(std::ifstream& s), std::ifstream &s)
inline void PutWisdom(char c, std::ofstream& s) {s.put(c);}
inline int GetWisdom(std::ifstream& s) {return s.get();}
+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) {}
+};
+
// Base clase for fft routines
//
class fftw {
protected:
unsigned int doubles; // number of double precision values in dataset
int sign;
+ unsigned int threads;
double norm;
fftw_plan plan;
@@ -214,54 +226,71 @@ protected:
return realsize(n,(Complex *) in,out);
}
- static bool Wise;
+ static unsigned int Wise;
static const double twopi;
- unsigned int threads;
public:
static unsigned int effort;
static unsigned int maxthreads;
static double testseconds;
+ static const char *WisdomName;
+ static bool mpi; // TODO: Remove when FFTW transpose routines are replaced
+
+ virtual unsigned int Threads() {return threads;}
- unsigned int Threads() {return threads;}
+ static const char *oddshift;
- // Shift the Fourier origin to (nx/2,0).
+ // Inplace shift of Fourier origin to (nx/2,0) for even nx.
static void Shift(Complex *data, unsigned int nx, unsigned int ny,
- int sign=0) {
+ unsigned int threads) {
const unsigned int nyp=ny/2+1;
- Complex *pstop=data+nx*nyp;
+ unsigned int stop=nx*nyp;
if(nx % 2 == 0) {
- int pinc=2*nyp;
- for(Complex *p=data+nyp; p < pstop; p += pinc) {
+ unsigned int inc=2*nyp;
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(threads)
+#endif
+ for(unsigned int i=nyp; i < stop; i += inc) {
+ Complex *p=data+i;
for(unsigned int j=0; j < nyp; j++) p[j]=-p[j];
}
} else {
- if(sign) {
- unsigned int c=nx/2;
- double arg=twopi*c/nx;
- for(unsigned int i=0; i < nx; i++) {
- double iarg=i*arg;
- Complex zeta(cos(iarg),sign*sin(iarg));
- Complex *datai=data+i*nyp;
- for(unsigned int j=0; j < nyp; j++) datai[j] *= zeta;
- }
- } else {
- std::cerr << "Shift for odd nx must be signed and interleaved"
- << std::endl;
- exit(1);
+ std::cerr << oddshift << std::endl;
+ exit(1);
+ }
+ }
+
+ // 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) {
+ if(nx % 2 == 0) {
+ unsigned int stop=nx*ny;
+ unsigned int inc=2*ny;
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(threads)
+#endif
+ for(unsigned int i=ny; i < stop; i += inc) {
+ double *p=data+i;
+ for(unsigned int j=0; j < ny; j++) p[j]=-p[j];
}
+ } else {
+ std::cerr << oddshift << std::endl;
+ exit(1);
}
}
- // Shift the Fourier origin to (nx/2,ny/2,0).
+ // Inplace 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, int sign=0) {
+ unsigned int nz, unsigned int threads) {
const unsigned int nzp=nz/2+1;
const unsigned int nyzp=ny*nzp;
if(nx % 2 == 0 && ny % 2 == 0) {
const unsigned int pinc=2*nzp;
Complex *pstop=data;
Complex *p=data;
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(threads)
+#endif
for(unsigned i=0; i < nx; i++) {
if(i % 2) p -= nzp;
else p += nzp;
@@ -271,43 +300,52 @@ public:
}
}
} else {
- if(sign) {
- unsigned int cx=nx/2;
- double argx=twopi*cx/nx;
- unsigned int cy=ny/2;
- double argy=twopi*cy/ny;
- for(unsigned i=0; i < nx; i++) {
- double iarg=i*argx;
- Complex zetax(cos(iarg),sign*sin(iarg));
- Complex *datai=data+nyzp*i;
- for(unsigned j=0; j < ny; j++) {
- double jarg=j*argy;
- Complex zeta=zetax*Complex(cos(jarg),sign*sin(jarg));
- Complex *dataij=datai+nzp*j;
- for(unsigned int k=0; k < nzp; k++) dataij[k] *= zeta;
- }
+ std::cerr << oddshift << "or odd ny" << std::endl;
+ exit(1);
+ }
+ }
+
+ // 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) {
+ const unsigned int nyz=ny*nz;
+ if(nx % 2 == 0 && ny % 2 == 0) {
+ const unsigned int pinc=2*nz;
+ double *pstop=data;
+ double *p=data;
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(threads)
+#endif
+ for(unsigned i=0; i < nx; i++) {
+ if(i % 2) p -= nz;
+ else p += nz;
+ pstop += nyz;
+ for(; p < pstop; p += pinc) {
+ for(unsigned int k=0; k < nz; k++) p[k]=-p[k];
}
- } else {
- std::cerr << "Shift for odd nx or ny must be signed and interleaved"
- << std::endl;
- exit(1);
}
+ } else {
+ std::cerr << oddshift << "or odd ny" << std::endl;
+ exit(1);
}
}
-
- 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) {
- if(autothreads) {
+
+ fftw(unsigned int doubles, int sign, unsigned int threads,
+ unsigned int n=0) :
+ doubles(doubles), sign(sign), threads(threads),
+ norm(1.0/(n ? n : doubles/2)), plan(NULL) {
#ifndef FFTWPP_SINGLE_THREAD
- fftw_init_threads();
+ if(!mpi) fftw_init_threads();
#endif
- autothreads=false;
- }
}
- virtual ~fftw() {if(plan) fftw_destroy_plan(plan);}
+ virtual ~fftw() {
+ if(Wise > 0)
+ --Wise;
+ if(Wise == 0)
+ SaveWisdom();
+ if(plan) fftw_destroy_plan(plan);
+ }
virtual fftw_plan Plan(Complex *in, Complex *out)=0;
@@ -317,24 +355,112 @@ public:
<< "-byte aligned: address " << p << std::endl;
}
- double stdev(unsigned int N, double sum, double sum2) {
- return N > 1 ? sqrt((sum2-sum*sum/N)/(N-1)) : 0.0;
- }
-
void noplan() {
std::cerr << "Unable to construct FFTW plan" << std::endl;
exit(1);
}
- void planThreads(unsigned int threads) {
+ static 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();
+ class statistics {
+ unsigned int N;
+ double A;
+ double varL;
+ double varH;
+ public:
+ statistics() : N(0), A(0.0), varL(0.0), varH(0.0) {}
+ double count() {return N;}
+ double mean() {return A;}
+ void add(double t) {
+ ++N;
+ double diff=t-A;
+ A += diff/N;
+ double v=diff*(t-A);
+ if(diff < 0.0)
+ varL += v;
+ else
+ varH += v;
+ }
+ double stdev(double var, double f) {
+ double factor=N > f ? f/(N-f) : 0.0;
+ return sqrt(var*factor);
+ }
+ double stdev() {
+ return stdev(varL+varH,1.0);
+ }
+ double stdevL() {
+ return stdev(varL,2.0);
+ }
+ double stdevH() {
+ return stdev(varH,2.0);
+ }
+ void output(const char *text, unsigned int m) {
+ std::cout << text << ":\n"
+ << m << "\t"
+ << A << "\t"
+ << stdevL() << "\t"
+ << stdevH() << std::endl;
+ }
+ };
+
+ threaddata time(fftw_plan plan1, fftw_plan planT, Complex *in, Complex *out,
+ unsigned int Threads) {
+ statistics S,ST;
+ double stop=totalseconds()+testseconds;
+ threads=1;
+ plan=plan1;
+ fft(in,out);
+ threads=Threads;
+ plan=planT;
+ fft(in,out);
+ unsigned int N=1;
+ for(;;) {
+ double t0=totalseconds();
+ threads=1;
+ plan=plan1;
+ for(unsigned int i=0; i < N; ++i)
+ fft(in,out);
+ double t1=totalseconds();
+ threads=Threads;
+ plan=planT;
+ for(unsigned int i=0; i < N; ++i)
+ fft(in,out);
+ double t=totalseconds();
+ S.add(t1-t0);
+ ST.add(t-t1);
+ if(S.mean() < 100.0/CLOCKS_PER_SEC) N *= 2;
+ 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) {}
+
+ threaddata Setup(Complex *in, Complex *out=NULL) {
+ if(!Wise) {LoadWisdom(); ++Wise;}
+
bool alloc=!in;
if(alloc) in=ComplexAlign((doubles+1)/2);
#ifndef NO_CHECK_ALIGN
@@ -346,76 +472,36 @@ public:
#endif
inplace=(out==in);
- planThreads(1);
- fftw_plan plan1=Plan(in,out);
- if(!plan1) noplan();
- plan=plan1;
+ threaddata data;
+ unsigned int Threads=threads;
+ if(threads > 1) data=lookup(inplace,threads);
+ threads=data.threads > 0 ? data.threads : 1;
+ planThreads(threads);
+ plan=Plan(in,out);
+ if(!plan) noplan();
- if(maxthreads > 1) {
- double sum2=0.0;
- unsigned int N=1;
- double begin=totalseconds();
- double lastseconds=begin;
- double stop=begin+testseconds;
- for(;;++N) {
- fft(in,out);
- 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=maxthreads;
+ if(Threads > 1 && data.threads == 0) {
+ threads=Threads;
planThreads(threads);
- plan=Plan(in,out);
- if(!plan) noplan();
-
- if(plan) {
- double begin=totalseconds();
- double lastseconds=begin;
- double stop=begin+testseconds;
- N=1;
- for(;;++N) {
- fft(in,out);
- 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) {
- threads=1;
- fftw_destroy_plan(plan);
- plan=plan1;
- } else {
- fftw_destroy_plan(plan1);
- threads=maxthreads;
- }
- }
- } else
- threads=1;
+ fftw_plan planT=Plan(in,out);
+ 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);
- SaveWisdom();
+ return data;
}
- void Setup(Complex *in, double *out) {Setup(in,(Complex *) out);}
- void Setup(double *in, Complex *out=NULL) {Setup((Complex *) in,out);}
+ void Setup(Complex *in, double *out) {
+ Setup(in,(Complex *) out);}
+ void Setup(double *in, Complex *out=NULL) {
+ Setup((Complex *) in,out);
+ }
- void LoadWisdom();
- void SaveWisdom();
+ static void LoadWisdom();
+ static void SaveWisdom();
virtual void Execute(Complex *in, Complex *out, bool=false) {
fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
@@ -430,7 +516,7 @@ public:
if(!out) out=in;
#endif
if(inplace ^ (out == in)) {
- std::cerr << "ERROR: fft constructor and call must be both in place or both out of place" << std::endl;
+ std::cerr << "ERROR: fft " << inout << std::endl;
exit(1);
}
return out;
@@ -464,10 +550,16 @@ public:
void Normalize(Complex *out) {
unsigned int stop=(doubles+1)/2;
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(threads)
+#endif
for(unsigned int 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;
}
@@ -507,21 +599,184 @@ public:
unsigned int nx, unsigned int M,
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 {
out=Setout(in,out);
Execute(in,out);
+
+ unsigned int stop=nx*stride;
+ Complex *outMdist=out+M*dist;
#ifndef FFTWPP_SINGLE_THREAD
#pragma omp parallel for num_threads(threads)
#endif
- for(unsigned int k=0; k < M; k++) {
- for(unsigned int j=0; j < nx; j++) {
- out[j*stride+k*dist] *= norm;
+ for(unsigned int i=0; i < stop; i += stride) {
+ Complex *pstop=outMdist+i;
+ for(Complex *p=out+i; p < pstop; p += dist) {
+ *p *= norm;
}
}
}
}
+
+}; // class fftw
+
+class Transpose {
+ fftw_plan plan;
+ unsigned int a,b;
+ unsigned int nlength,mlength;
+ unsigned int instride,outstride;
+ unsigned int threads;
+ bool inplace;
+ unsigned int size;
+public:
+ template<class T>
+ Transpose(unsigned int rows, unsigned int cols, unsigned int length,
+ T *in, T *out=NULL, unsigned int threads=fftw::maxthreads) :
+ threads(threads) {
+ size=sizeof(T);
+ if(size % sizeof(double) != 0) {
+ std::cerr << "ERROR: Transpose is not implemented for type of size "
+ << size;
+ exit(1);
+ }
+ size /= sizeof(double);
+ length *= size;
+ if(!out) out=in;
+ inplace=(out==in);
+ if(inplace) threads=1;
+
+ fftw_iodim dims[3];
+
+ a=std::min(rows,threads);
+ b=std::min(cols,threads/a);
+ unsigned int n=rows/a;
+ unsigned int m=cols/b;
+ nlength=n*length;
+ mlength=m*length;
+
+ instride=cols;
+ outstride=rows;
+
+ dims[0].n=n;
+ dims[0].is=instride*length;
+ dims[0].os=length;
+
+ dims[1].n=m;
+ dims[1].is=length;
+ dims[1].os=outstride*length;
+
+ dims[2].n=length;
+ dims[2].is=1;
+ dims[2].os=1;
+
+ fftw::planThreads(inplace ? threads : 1);
+ plan=fftw_plan_guru_r2r(0,NULL,3,dims,(double *) in,
+ (double *) out,NULL,fftw::effort);
+ }
+
+ ~Transpose() {if(plan) fftw_destroy_plan(plan);}
+
+ template<class T>
+ void transpose(T *in, T *out=NULL) {
+ if(!out) out=in;
+ if(inplace ^ (out == in)) {
+ std::cerr << "ERROR: Transpose " << inout << std::endl;
+ exit(1);
+ }
+#ifndef FFTWPP_SINGLE_THREAD
+ if(threads == 1)
+#endif
+ fftw_execute_r2r(plan,(double *) in,(double*) out);
+#ifndef FFTWPP_SINGLE_THREAD
+ else {
+ int A=a, B=b;
+#pragma omp parallel for num_threads(A)
+ for(unsigned int i=0; i < a; ++i) {
+ unsigned int I=i*nlength;
+#pragma omp parallel for num_threads(B)
+ for(unsigned int j=0; j < b; ++j) {
+ unsigned int J=j*mlength;
+ fftw_execute_r2r(plan,(double *) in+instride*I+J,
+ (double *) out+outstride*J+I);
+ }
+ }
+ }
+#endif
+ }
+};
+
+template<class T, class L>
+class Threadtable {
+public:
+ typedef std::map<T,threaddata,L> Table;
+
+ threaddata Lookup(Table& table, T key) {
+ typename Table::iterator p=table.find(key);
+ return p == table.end() ? threaddata() : p->second;
+ }
+
+ void Store(Table& threadtable, T key, const threaddata& data) {
+ threadtable[key]=data;
+ }
+};
+
+struct keytype1 {
+ unsigned int nx;
+ unsigned int threads;
+ bool inplace;
+ keytype1(unsigned int nx, unsigned int threads, bool inplace) :
+ nx(nx), threads(threads), inplace(inplace) {}
+};
+
+struct keyless1 {
+ bool operator()(const keytype1& a, const keytype1& 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.threads < b.threads ||
+ (a.threads == b.threads &&
+ a.inplace < b.inplace)))))));
+ }
};
// Compute the complex Fourier transform of n complex values.
@@ -547,18 +802,28 @@ public:
// fft1d Backward(n,1);
// Backward.fft(in);
//
-class fft1d : public fftw {
+class fft1d : public fftw, public Threadtable<keytype1,keyless1> {
unsigned int nx;
+ static Table threadtable;
public:
- fft1d(unsigned int nx, int sign, Complex *in=NULL, Complex *out=NULL)
- : fftw(2*nx,sign), nx(nx) {Setup(in,out);}
+ fft1d(unsigned int nx, int sign, Complex *in=NULL, Complex *out=NULL,
+ unsigned int 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)
- : fftw(2*in.Nx(),sign), nx(in.Nx()) {Setup(in,out);}
+ const Array::array1<Complex>& out=Array::NULL1,
+ unsigned int 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);
@@ -582,29 +847,110 @@ public:
//
// Notes:
// stride is the spacing between the elements of each Complex vector;
-// dist is the spacing between the first elements of the vectors;
+// dist is the spacing between the first elements of the vectors.
//
//
-class mfft1d : public fftw {
+class mfft1d : public fftw, public Threadtable<keytype3,keyless3> {
unsigned int nx;
unsigned int M;
+ unsigned int T,Q,R;
size_t stride;
size_t dist;
+ fftw_plan plan1;
+ fftw_plan plan2;
+ static Table threadtable;
public:
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);}
+ 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),
+ nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)),
+ plan1(NULL), plan2(NULL)
+ {
+ T=1;
+ Q=M;
+ R=0;
+ threaddata S1=Setup(in,out);
+ fftw_plan planT1=plan;
+
+ T=std::min(M,Threads);
+ if(T > 1) {
+ Q=M/T;
+ R=M-Q*T;
+ threads=Threads;
+ threaddata ST=Setup(in,out);
+
+ if(R > 0 && threads == 1 && plan1 != plan2) {
+ fftw_destroy_plan(plan2);
+ plan2=plan1;
+ }
+
+ if(ST.mean > S1.mean-S1.stdev) {
+ fftw_destroy_plan(plan);
+ if(R > 0) {
+ fftw_destroy_plan(plan2);
+ plan2=NULL;
+ }
+ T=1;
+ Q=M;
+ R=0;
+ plan=planT1;
+ threads=S1.threads;
+ } else {
+ fftw_destroy_plan(planT1);
+ threads=ST.threads;
+ }
+ }
+ }
+
+ unsigned int Threads() {return std::max(T,threads);}
+
+ threaddata lookup(bool inplace, unsigned int threads) {
+ return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace));
+ }
+ void store(bool inplace, const threaddata& data) {
+ Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data);
+ }
fftw_plan Plan(Complex *in, Complex *out) {
int n=(int) nx;
- return fftw_plan_many_dft(1,&n,M,
+ if(R > 0) {
+ plan2=fftw_plan_many_dft(1,&n,Q+1,
+ (fftw_complex *) in,NULL,stride,dist,
+ (fftw_complex *) out,NULL,stride,dist,
+ sign,effort);
+ if(!plan2) noplan();
+ if(threads == 1) plan1=plan2;
+ }
+ return fftw_plan_many_dft(1,&n,Q,
(fftw_complex *) in,NULL,stride,dist,
(fftw_complex *) out,NULL,stride,dist,
sign,effort);
}
+ ~mfft1d() {
+ if(plan2) fftw_destroy_plan(plan2);
+ }
+
+ void Execute(Complex *in, Complex *out, bool=false) {
+ if(T == 1) {
+ fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
+ } else {
+ unsigned int Tdist=T*dist;
+ unsigned int extra=(T-R)*dist;
+
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(T)
+#endif
+ for(unsigned int i=0; i < Tdist; i += dist) {
+ bool normal=i < extra;
+ unsigned int offset=normal ? Q*i : Q*i+i-extra;
+ fftw_execute_dft(normal ? plan : plan2,(fftw_complex *) in+offset,
+ (fftw_complex *) out+offset);
+ }
+ }
+ }
+
void fftNormalized(Complex *in, Complex *out=NULL) {
fftw::fftNormalized(in,out,nx,M,stride,dist);
}
@@ -629,24 +975,35 @@ 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 {
+class rcfft1d : public fftw, public Threadtable<keytype1,keyless1> {
unsigned int nx;
+ static Table threadtable;
public:
- rcfft1d(unsigned int nx, Complex *out=NULL)
- : fftw(2*(nx/2+1),-1,nx), nx(nx) {Setup(out);}
+ rcfft1d(unsigned int nx, Complex *out=NULL, unsigned int 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)
- : fftw(realsize(nx,in,out),-1,nx), nx(nx) {Setup(in,out);}
+ rcfft1d(unsigned int nx, double *in, Complex *out=NULL,
+ unsigned int threads=maxthreads)
+ : fftw(realsize(nx,in,out),-1,threads,nx), nx(nx) {Setup(in,out);}
#ifdef __Array_h__
- rcfft1d(unsigned int nx, const Array::array1<Complex>& out)
- : fftw(out.Size(),-1,nx), nx(nx) {Setup(out);}
+ rcfft1d(unsigned int nx, const Array::array1<Complex>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),-1,threads,nx), nx(nx) {Setup(out,(double*) NULL);}
rcfft1d(unsigned int nx, const Array::array1<double>& in,
- const Array::array1<Complex>& out=Array::NULL1)
- : fftw(realsize(nx,in(),out()),-1,nx), nx(nx) {Setup(in,out);}
+ const Array::array1<Complex>& out=Array::NULL1,
+ unsigned int threads=maxthreads)
+ : fftw(realsize(nx,in(),out()),-1,threads,nx), nx(nx) {Setup(in,out);}
#endif
+ 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);
}
@@ -677,27 +1034,38 @@ 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 {
+class crfft1d : public fftw, public Threadtable<keytype1,keyless1> {
unsigned int nx;
+ static Table threadtable;
public:
- crfft1d(unsigned int nx, double *out=NULL)
- : fftw(2*(nx/2+1),1,nx), nx(nx) {Setup(out);}
+ crfft1d(unsigned int nx, double *out=NULL, unsigned int threads=maxthreads)
+ : fftw(2*(nx/2+1),1,threads,nx), nx(nx) {Setup(out);}
- crfft1d(unsigned int nx, Complex *in, double *out=NULL)
- : fftw(realsize(nx,in,out),1,nx), nx(nx) {Setup(in,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);}
#ifdef __Array_h__
- crfft1d(unsigned int nx, const Array::array1<double>& out)
- : fftw(out.Size(),1,nx), nx(nx) {Setup(out);}
+ crfft1d(unsigned int nx, const Array::array1<double>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),1,threads,nx), nx(nx) {Setup(out);}
- crfft1d(unsigned int nx, const Array::array1<Complex>& in)
- : fftw(2*in.Size(),1,nx), nx(nx) {Setup(in);}
+ crfft1d(unsigned int nx, const Array::array1<Complex>& in,
+ unsigned int threads=maxthreads)
+ : fftw(2*in.Size(),1,threads,nx), nx(nx) {Setup(in);}
crfft1d(unsigned int nx, const Array::array1<Complex>& in,
const Array::array1<double>& out)
- : fftw(out.Size(),1,nx), nx(nx) {Setup(in,out);}
+ : fftw(out.Size(),1,threads,nx), nx(nx) {Setup(in,out);}
#endif
+ 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_c2r_1d(nx,(fftw_complex *) in,(double *) out,effort);
}
@@ -729,27 +1097,42 @@ 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 fftw {
+class mrcfft1d : public fftw, public Threadtable<keytype2,keyless2> {
unsigned int nx;
unsigned int M;
size_t stride;
size_t dist;
+ static Table threadtable;
public:
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);}
+ size_t dist=0, Complex *out=NULL,
+ unsigned int threads=maxthreads)
+ : fftw(2*(nx/2*stride+(M-1)*Dist(nx,stride,dist)+1),-1,threads,nx), nx(nx),
+ M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(out);}
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);}
+ size_t dist=0, double *in=NULL, Complex *out=NULL,
+ unsigned int threads=maxthreads)
+ : fftw(2*(nx/2*stride+(M-1)*Dist(nx,stride,dist)+1),-1,threads,nx), nx(nx),
+ M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
+
+ threaddata lookup(bool inplace, unsigned int threads) {
+ return Lookup(threadtable,keytype2(nx,M,threads,inplace));
+ }
+ void store(bool inplace, const threaddata& data) {
+ Store(threadtable,keytype2(nx,M,data.threads,inplace),data);
+ }
fftw_plan Plan(Complex *in, Complex *out) {
int n=(int) nx;
+ if(in != out)
+ return fftw_plan_many_dft_r2c(1,&n,M,
+ (double *) in,NULL,stride,dist,
+ (fftw_complex *) out,NULL,stride,dist/2+1,
+ effort);
return fftw_plan_many_dft_r2c(1,&n,M,
- (double *) in,NULL,stride,2*dist,
- (fftw_complex *) out,NULL,stride,dist,
+ (double *) in,NULL,stride,dist,
+ (fftw_complex *) out,NULL,stride,dist/2,
effort);
}
@@ -782,22 +1165,37 @@ public:
// Notes:
// stride is the spacing between the elements of each Complex vector;
// dist is the spacing between the first elements of the vectors;
-// in contains the first n/2+1 Complex Fourier values.
-// out contains the n real values stored as a Complex array;
+// in contains the first n/2+1 Complex Fourier values;
+// out contains the n real values stored as a Complex array.
//
-class mcrfft1d : public fftw {
+class mcrfft1d : public fftw, public Threadtable<keytype2,keyless2> {
unsigned int nx;
unsigned int M;
size_t stride;
size_t dist;
+ static Table threadtable;
public:
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),
+ size_t dist=0, Complex *in=NULL, double *out=NULL,
+ unsigned int threads=maxthreads)
+ : fftw((realsize(nx,in,out)-2)*stride+2*(M-1)*Dist(nx,stride,dist)+2,1,
+ threads,nx),
nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
+ threaddata lookup(bool inplace, unsigned int threads) {
+ return Lookup(threadtable,keytype2(nx,M,threads,inplace));
+ }
+ void store(bool inplace, const threaddata& data) {
+ Store(threadtable,keytype2(nx,M,data.threads,inplace),data);
+ }
+
fftw_plan Plan(Complex *in, Complex *out) {
int n=(int) nx;
+ if(in != out)
+ return fftw_plan_many_dft_c2r(1,&n,M,
+ (fftw_complex *) in,NULL,stride,dist,
+ (double *) out,NULL,stride,2*(dist-1),
+ effort);
return fftw_plan_many_dft_c2r(1,&n,M,
(fftw_complex *) in,NULL,stride,dist,
(double *) out,NULL,stride,2*dist,
@@ -808,8 +1206,23 @@ public:
fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
}
- void fftNormalized(Complex *in, Complex *out=NULL) {
- fftw::fftNormalized(in,out,(nx/2+1),M,stride,dist);
+ void fftNormalized(Complex *in, double *out=NULL) {
+ // TODO: improve threading (cf: mfft1d).
+ fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
+
+ unsigned int stop=nx*stride;
+ unsigned int odist=in == (Complex *) out ? 2*dist : 2*(dist-1);
+ double *outMdist=out+M*odist;
+
+#ifndef FFTWPP_SINGLE_THREAD
+#pragma omp parallel for num_threads(threads)
+#endif
+ for(unsigned int i=0; i < stop; i += stride) {
+ double *pstop=outMdist+i;
+ for(double *p=out+i; p < pstop; p += odist) {
+ *p *= norm;
+ }
+ }
}
};
@@ -839,20 +1252,29 @@ public:
// Note:
// in[ny*i+j] contains the ny Complex values for each i=0,...,nx-1.
//
-class fft2d : public fftw {
+class fft2d : public fftw, public Threadtable<keytype2,keyless2> {
unsigned int nx;
unsigned int ny;
+ static Table threadtable;
public:
fft2d(unsigned int nx, unsigned int ny, int sign, Complex *in=NULL,
- Complex *out=NULL)
- : fftw(2*nx*ny,sign), nx(nx), ny(ny) {Setup(in,out);}
+ Complex *out=NULL, unsigned int 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)
- : fftw(2*in.Size(),sign), nx(in.Nx()), ny(in.Ny()) {Setup(in,out);}
+ const Array::array2<Complex>& out=Array::NULL2,
+ unsigned int 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);
@@ -873,6 +1295,8 @@ public:
//
// rcfft2d Forward(nx,ny,in,out);
// Forward.fft(in,out); // Origin of Fourier domain at (0,0)
+// Forward.fft0(in,out); // Origin of Fourier domain at (nx/2,0);
+// input destroyed.
//
// In-place usage:
//
@@ -888,20 +1312,28 @@ class rcfft2d : public fftw {
unsigned int nx;
unsigned int ny;
public:
- rcfft2d(unsigned int nx, unsigned int ny, Complex *out=NULL)
- : fftw(2*nx*(ny/2+1),-1,nx*ny), nx(nx), ny(ny) {Setup(out);}
+ rcfft2d(unsigned int nx, unsigned int ny, Complex *out=NULL,
+ unsigned int 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)
- : fftw(nx*realsize(ny,in,out),-1,nx*ny), nx(nx), ny(ny) {Setup(in,out);}
+ rcfft2d(unsigned int nx, unsigned int ny, double *in, Complex *out=NULL,
+ unsigned int threads=maxthreads)
+ : fftw(nx*realsize(ny,in,out),-1,threads,nx*ny), nx(nx), ny(ny) {
+ Setup(in,out);
+ }
#ifdef __Array_h__
- rcfft2d(unsigned int ny, const Array::array2<Complex>& out)
- : fftw(out.Size(),-1,out.Nx()*ny), nx(out.Nx()), ny(ny) {Setup(out);}
+ rcfft2d(unsigned int ny, const Array::array2<Complex>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),-1,threads,out.Nx()*ny), nx(out.Nx()), ny(ny) {
+ Setup(out);
+ }
rcfft2d(unsigned int ny, const Array::array2<double>& in,
- const Array::array2<Complex>& out=Array::NULL2)
- : fftw(in.Nx()*realsize(ny,in(),out()),-1,in.Nx()*ny), nx(in.Nx()), ny(ny)
- {Setup(in,out);}
+ const Array::array2<Complex>& out=Array::NULL2,
+ unsigned int threads=maxthreads)
+ : fftw(in.Nx()*realsize(ny,in(),out()),-1,threads,in.Nx()*ny),
+ nx(in.Nx()), ny(ny) {Setup(in,out);}
#endif
fftw_plan Plan(Complex *in, Complex *out) {
@@ -910,7 +1342,10 @@ public:
}
void Execute(Complex *in, Complex *out, bool shift=false) {
- if(shift && inplace) Shift(in,nx,ny);
+ if(shift) {
+ if(inplace) Shift(in,nx,ny,threads);
+ else Shift((double *) in,nx,ny,threads);
+ }
fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
}
};
@@ -919,7 +1354,7 @@ public:
// nx*(ny/2+1) Complex values corresponding to the spectral values in the
// half-plane ky >= 0, using phase sign +1.
// Before calling fft(), the array in must be allocated as
-// Complex[nx*(ny+1)/2] and the array out must be allocated as
+// Complex[nx*(ny/2+1)] and the array out must be allocated as
// double[nx*ny]. The arrays in and out may coincide,
// allocated as Complex[nx*(ny/2+1)].
//
@@ -943,22 +1378,29 @@ class crfft2d : public fftw {
unsigned int nx;
unsigned int ny;
public:
- crfft2d(unsigned int nx, unsigned int ny, Complex *in=NULL) :
- fftw(2*nx*(ny/2+1),1,nx*ny), nx(nx), ny(ny) {Setup(in);}
+ crfft2d(unsigned int nx, unsigned int ny, Complex *in=NULL,
+ unsigned int threads=maxthreads) :
+ fftw(2*nx*(ny/2+1),1,threads,nx*ny), nx(nx), ny(ny) {Setup(in);}
- crfft2d(unsigned int nx, unsigned int ny, Complex *in, double *out) :
- fftw(nx*realsize(ny,in,out),1,nx*ny), nx(nx), ny(ny) {Setup(in,out);}
+ crfft2d(unsigned int nx, unsigned int ny, Complex *in, double *out,
+ unsigned int threads=maxthreads)
+ : fftw(nx*realsize(ny,in,out),1,threads,nx*ny), nx(nx), ny(ny) {
+ Setup(in,out);
+ }
#ifdef __Array_h__
- crfft2d(unsigned int ny, const Array::array2<double>& out)
- : fftw(out.Size(),1,out.Nx()*ny), nx(out.Nx()), ny(ny) {Setup(out);}
-
- crfft2d(unsigned int ny, const Array::array2<Complex>& in)
- : fftw(2*in.Size(),1,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in);}
+ crfft2d(unsigned int ny, const Array::array2<double>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),1,threads,out.Nx()*ny), nx(out.Nx()), ny(ny) {Setup(out);}
+
+ crfft2d(unsigned int ny, const Array::array2<Complex>& in,
+ unsigned int threads=maxthreads)
+ : fftw(2*in.Size(),1,threads,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in);}
crfft2d(unsigned int ny, const Array::array2<Complex>& in,
- const Array::array2<double>& out)
- : fftw(out.Size(),1,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in,out);}
+ const Array::array2<double>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),1,threads,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in,out);}
#endif
fftw_plan Plan(Complex *in, Complex *out) {
@@ -968,7 +1410,10 @@ public:
void Execute(Complex *in, Complex *out, bool shift=false) {
fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
- if(shift) Shift(out,nx,ny);
+ if(shift) {
+ if(inplace) Shift(out,nx,ny,threads);
+ else Shift((double *) out,nx,ny,threads);
+ }
}
};
@@ -1005,13 +1450,15 @@ class fft3d : public fftw {
unsigned int nz;
public:
fft3d(unsigned int nx, unsigned int ny, unsigned int nz,
- int sign, Complex *in=NULL, Complex *out=NULL)
- : fftw(2*nx*ny*nz,sign), nx(nx), ny(ny), nz(nz) {Setup(in,out);}
+ int sign, Complex *in=NULL, Complex *out=NULL,
+ unsigned int 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)
- : fftw(2*in.Size(),sign), nx(in.Nx()), ny(in.Ny()), nz(in.Nz())
+ const Array::array3<Complex>& out=Array::NULL3,
+ unsigned int threads=maxthreads)
+ : fftw(2*in.Size(),sign,threads), nx(in.Nx()), ny(in.Ny()), nz(in.Nz())
{Setup(in,out);}
#endif
@@ -1031,7 +1478,8 @@ public:
//
// rcfft3d Forward(nx,ny,nz,in,out);
// Forward.fft(in,out); // Origin of Fourier domain at (0,0)
-//
+// Forward.fft0(in,out); // Origin of Fourier domain at (nx/2,ny/2,0);
+// input destroyed
// In-place usage:
//
// rcfft3d Forward(nx,ny,nz);
@@ -1048,20 +1496,23 @@ class rcfft3d : public fftw {
unsigned int nz;
public:
rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *out=NULL)
- : fftw(2*nx*ny*(nz/2+1),-1,nx*ny*nz), nx(nx), ny(ny), nz(nz) {Setup(out);}
+ : 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) : fftw(nx*ny*realsize(nz,in,out),-1,nx*ny*nz),
- nx(nx), ny(ny), nz(nz) {Setup(in,out);}
+ Complex *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);}
#ifdef __Array_h__
- rcfft3d(unsigned int nz, const Array::array3<Complex>& out)
- : fftw(out.Size(),-1,out.Nx()*out.Ny()*nz),
+ rcfft3d(unsigned int nz, const Array::array3<Complex>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),-1,threads,out.Nx()*out.Ny()*nz),
nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);}
rcfft3d(unsigned int nz, const Array::array3<double>& in,
- const Array::array3<Complex>& out=Array::NULL3)
- : fftw(in.Nx()*in.Ny()*realsize(nz,in(),out()),-1,in.Size()),
+ const Array::array3<Complex>& out=Array::NULL3,
+ unsigned int threads=maxthreads)
+ : fftw(in.Nx()*in.Ny()*realsize(nz,in(),out()),-1,threads,in.Size()),
nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);}
#endif
@@ -1071,7 +1522,10 @@ public:
}
void Execute(Complex *in, Complex *out, bool shift=false) {
- if(shift && inplace) Shift(in,nx,ny,nz);
+ if(shift) {
+ if(inplace) Shift(in,nx,ny,nz,threads);
+ else Shift((double *) in,nx,ny,nz,threads);
+ }
fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
}
};
@@ -1105,27 +1559,31 @@ class crfft3d : public fftw {
unsigned int ny;
unsigned int nz;
public:
- crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in=NULL)
- : fftw(2*nx*ny*(nz/2+1),1,nx*ny*nz), nx(nx), ny(ny), nz(nz)
+ crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in=NULL,
+ unsigned int threads=maxthreads)
+ : fftw(2*nx*ny*(nz/2+1),1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz)
{Setup(in);}
crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in,
- double *out=NULL) :
- fftw(nx*ny*(realsize(nz,in,out)),1,nx*ny*nz), nx(nx), ny(ny), nz(nz)
- {Setup(in,out);}
+ 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);}
#ifdef __Array_h__
- crfft3d(unsigned int nz, const Array::array3<double>& out)
- : fftw(out.Size(),1,out.Nx()*out.Ny()*nz),
+ crfft3d(unsigned int nz, const Array::array3<double>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),1,threads,out.Nx()*out.Ny()*nz),
nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);}
- crfft3d(unsigned int nz, const Array::array3<Complex>& in)
- : fftw(2*in.Size(),1,in.Nx()*in.Ny()*nz),
+ crfft3d(unsigned int nz, const Array::array3<Complex>& in,
+ unsigned int threads=maxthreads)
+ : fftw(2*in.Size(),1,threads,in.Nx()*in.Ny()*nz),
nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in);}
crfft3d(unsigned int nz, const Array::array3<Complex>& in,
- const Array::array3<double>& out)
- : fftw(out.Size(),1,in.Nx()*in.Ny()*nz),
+ const Array::array3<double>& out,
+ unsigned int threads=maxthreads)
+ : fftw(out.Size(),1,threads,in.Nx()*in.Ny()*nz),
nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);}
#endif
@@ -1136,7 +1594,10 @@ public:
void Execute(Complex *in, Complex *out, bool shift=false) {
fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
- if(shift) Shift(out,nx,ny,nz);
+ if(shift) {
+ if(inplace) Shift(out,nx,ny,nz,threads);
+ else Shift((double *) out,nx,ny,nz,threads);
+ }
}
};