summaryrefslogtreecommitdiff
path: root/graphics/asymptote/fftw++.h
diff options
context:
space:
mode:
Diffstat (limited to 'graphics/asymptote/fftw++.h')
-rw-r--r--graphics/asymptote/fftw++.h1046
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;
}
};