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++.h418
1 files changed, 209 insertions, 209 deletions
diff --git a/graphics/asymptote/fftw++.h b/graphics/asymptote/fftw++.h
index 01e1f55152..467387d34f 100644
--- a/graphics/asymptote/fftw++.h
+++ b/graphics/asymptote/fftw++.h
@@ -2,7 +2,7 @@
Copyright (C) 2004-16
John C. Bowman, University of Alberta
Malcolm Roberts, University of Strasbourg
-
+
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
the Free Software Foundation; either version 3 of the License, or
@@ -39,22 +39,22 @@
#include <omp.h>
#endif
-inline int get_thread_num()
+inline int get_thread_num()
{
#ifdef FFTWPP_SINGLE_THREAD
return 0;
#else
return omp_get_thread_num();
-#endif
+#endif
}
-inline int get_max_threads()
+inline int get_max_threads()
{
#ifdef FFTWPP_SINGLE_THREAD
return 1;
#else
return omp_get_max_threads();
-#endif
+#endif
}
#ifndef FFTWPP_SINGLE_THREAD
@@ -69,7 +69,7 @@ inline int get_max_threads()
#define PARALLEL(code) \
{ \
code \
- }
+ }
#endif
#ifndef __Complex_h__
@@ -117,12 +117,12 @@ class ThreadBase
protected:
unsigned int threads;
unsigned int innerthreads;
-public:
+public:
ThreadBase();
ThreadBase(unsigned int threads) : threads(threads) {}
void Threads(unsigned int nthreads) {threads=nthreads;}
unsigned int Threads() {return threads;}
-
+
void multithread(unsigned int nx) {
if(nx >= threads) {
innerthreads=1;
@@ -137,17 +137,17 @@ inline unsigned int realsize(unsigned int n, Complex *in, Complex *out=NULL)
{
return (!out || in == out) ? 2*(n/2+1) : n;
}
-
+
inline unsigned int realsize(unsigned int n, Complex *in, double *out)
{
- return realsize(n,in,(Complex *) out);
+ return realsize(n,in,(Complex *) out);
}
-
+
inline unsigned int realsize(unsigned int n, double *in, Complex *out)
{
- return realsize(n,(Complex *) in,out);
+ return realsize(n,(Complex *) in,out);
}
-
+
// Base clase for fft routines
//
class fftw : public ThreadBase {
@@ -159,24 +159,24 @@ protected:
fftw_plan plan;
bool inplace;
-
+
unsigned int Dist(unsigned int 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 fftw_plan (*planner)(fftw *f, Complex *in, Complex *out);
-
+
virtual unsigned int Threads() {return threads;}
-
+
static const char *oddshift;
-
+
// Inplace 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) {
@@ -260,41 +260,41 @@ public:
exit(1);
}
}
-
+
fftw() : plan(NULL) {}
fftw(unsigned int doubles, int sign, unsigned int threads,
unsigned int n=0) :
- doubles(doubles), sign(sign), threads(threads),
+ doubles(doubles), sign(sign), threads(threads),
norm(1.0/(n ? n : doubles/2)), plan(NULL) {
#ifndef FFTWPP_SINGLE_THREAD
fftw_init_threads();
-#endif
+#endif
}
-
+
virtual ~fftw() {
if(plan) fftw_destroy_plan(plan);
}
-
+
virtual fftw_plan Plan(Complex *in, Complex *out) {return NULL;};
-
+
inline void CheckAlign(Complex *p, const char *s) {
if((size_t) p % sizeof(Complex) == 0) return;
- std::cerr << "WARNING: " << s << " array is not " << sizeof(Complex)
+ std::cerr << "WARNING: " << s << " array is not " << sizeof(Complex)
<< "-byte aligned: address " << p << std::endl;
}
-
+
void noplan() {
std::cerr << "Unable to construct FFTW plan" << std::endl;
exit(1);
}
-
+
static void planThreads(unsigned int threads) {
#ifndef FFTWPP_SINGLE_THREAD
omp_set_num_threads(threads);
fftw_plan_with_nthreads(threads);
-#endif
+#endif
}
-
+
threaddata time(fftw_plan plan1, fftw_plan planT, Complex *in, Complex *out,
unsigned int Threads) {
utils::statistics S,ST;
@@ -339,30 +339,30 @@ public:
}
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
+#ifndef NO_CHECK_ALIGN
CheckAlign(in,constructor ? "constructor input" : "input");
if(out) CheckAlign(out,constructor ? "constructor output" : "output");
else out=in;
#else
if(!out) out=in;
-#endif
+#endif
return out;
}
-
+
threaddata 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);
@@ -370,13 +370,13 @@ public:
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);
@@ -384,11 +384,11 @@ public:
store(inplace,threaddata(threads,data.mean,data.stdev));
}
}
-
+
if(alloc) Array::deleteAlign(in,(doubles+1)/2);
return data;
}
-
+
threaddata Setup(Complex *in, double *out) {
return Setup(in,(Complex *) out);
}
@@ -396,11 +396,11 @@ public:
threaddata Setup(double *in, Complex *out=NULL) {
return Setup((Complex *) in,out);
}
-
+
virtual void Execute(Complex *in, Complex *out, bool=false) {
fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
}
-
+
Complex *Setout(Complex *in, Complex *out) {
out=CheckAlign(in,out,false);
if(inplace ^ (out == in)) {
@@ -409,33 +409,33 @@ public:
}
return out;
}
-
+
void fft(Complex *in, Complex *out=NULL) {
out=Setout(in,out);
Execute(in,out);
}
-
+
void fft(double *in, Complex *out=NULL) {
fft((Complex *) in,out);
}
-
+
void fft(Complex *in, double *out) {
fft(in,(Complex *) out);
}
-
+
void fft0(Complex *in, Complex *out=NULL) {
out=Setout(in,out);
Execute(in,out,true);
}
-
+
void fft0(double *in, Complex *out=NULL) {
fft0((Complex *) in,out);
}
-
+
void fft0(Complex *in, double *out) {
fft0(in,(Complex *) out);
}
-
+
void Normalize(Complex *out) {
unsigned int stop=doubles/2;
#ifndef FFTWPP_SINGLE_THREAD
@@ -450,29 +450,29 @@ public:
#endif
for(unsigned int i=0; i < doubles; i++) out[i] *= norm;
}
-
- virtual void fftNormalized(Complex *in, Complex *out=NULL, bool shift=false)
+
+ virtual void fftNormalized(Complex *in, Complex *out=NULL, bool shift=false)
{
out=Setout(in,out);
Execute(in,out,shift);
Normalize(out);
}
-
+
virtual void fftNormalized(Complex *in, double *out, bool shift=false) {
out=(double *) Setout(in,(Complex *) out);
Execute(in,(Complex *) out,shift);
Normalize(out);
}
-
+
virtual void fftNormalized(double *in, Complex *out, bool shift=false) {
fftNormalized((Complex *) in,out,shift);
}
-
+
template<class I, class O>
void fft0Normalized(I in, O out) {
fftNormalized(in,out,true);
}
-
+
template<class O>
void Normalize(unsigned int nx, unsigned int M, size_t ostride,
size_t odist, O *out) {
@@ -488,7 +488,7 @@ 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) {
@@ -508,7 +508,7 @@ public:
T *in, T *out=NULL, unsigned int threads=fftw::maxthreads) {
unsigned int size=sizeof(T);
if(size % sizeof(double) != 0) {
- std::cerr << "ERROR: Transpose is not implemented for type of size "
+ std::cerr << "ERROR: Transpose is not implemented for type of size "
<< size;
exit(1);
}
@@ -520,14 +520,14 @@ public:
if(!out) out=in;
inplace=(out==in);
fftw::planThreads(threads);
-
+
fftw_iodim dims[3];
- dims[0].n=rows;
+ dims[0].n=rows;
dims[0].is=cols*length;
dims[0].os=length;
-
- dims[1].n=cols;
+
+ dims[1].n=cols;
dims[1].is=length;
dims[1].os=rows*length;
@@ -543,7 +543,7 @@ public:
~Transpose() {
if(plan) fftw_destroy_plan(plan);
}
-
+
template<class T>
void transpose(T *in, T *out=NULL) {
if(!plan) return;
@@ -565,7 +565,7 @@ public:
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;
}
@@ -575,10 +575,10 @@ struct keytype1 {
unsigned int nx;
unsigned int threads;
bool inplace;
- 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 &&
@@ -593,10 +593,10 @@ struct keytype2 {
unsigned int threads;
bool inplace;
keytype2(unsigned int nx, unsigned int ny, unsigned int threads,
- bool inplace) :
+ 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 &&
@@ -614,13 +614,13 @@ struct keytype3 {
unsigned int threads;
bool inplace;
keytype3(unsigned int nx, unsigned int ny, unsigned int nz,
- unsigned int threads, bool inplace) :
+ 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 &&
+ 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 &&
@@ -634,7 +634,7 @@ struct keyless3 {
// Before calling fft(), the arrays in and out (which may coincide) must be
// allocated as Complex[n].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// fft1d Forward(n,-1,in,out);
// Forward.fft(in,out);
@@ -656,31 +656,31 @@ struct keyless3 {
class fft1d : public fftw, public Threadtable<keytype1,keyless1> {
unsigned int nx;
static Table threadtable;
-public:
+public:
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);}
-
+ : 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)
- : fftw(2*in.Nx(),sign,threads), nx(in.Nx()) {Setup(in,out);}
-#endif
-
+ 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);
}
};
-
+
template<class I, class O>
class fftwblock : public virtual fftw {
public:
@@ -699,19 +699,19 @@ public:
T=1;
Q=M;
R=0;
-
+
threaddata S1=Setup(in,out);
fftw_plan planT1=plan;
-
+
if(fftw::maxthreads > 1) {
if(Threads > 1) {
T=std::min(M,Threads);
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;
@@ -735,23 +735,23 @@ public:
} else
Setup(in,out); // Synchronize wisdom
}
- }
-
+ }
+
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);
}
-
+
fftw_plan Plan(int Q, double *in, fftw_complex *out) {
return fftw_plan_many_dft_r2c(1,&nx,Q,in,NULL,istride,idist,
out,NULL,ostride,odist,effort);
}
-
+
fftw_plan Plan(int Q, fftw_complex *in, double *out) {
return fftw_plan_many_dft_c2r(1,&nx,Q,in,NULL,istride,idist,
out,NULL,ostride,odist,effort);
}
-
+
fftw_plan Plan(Complex *in, Complex *out) {
if(R > 0) {
plan2=Plan(Q+1,(I *) in,(O *) out);
@@ -760,11 +760,11 @@ public:
}
return Plan(Q,(I *) in,(O *) out);
}
-
+
void Execute(fftw_plan plan, fftw_complex *in, fftw_complex *out) {
fftw_execute_dft(plan,in,out);
}
-
+
void Execute(fftw_plan plan, double *in, fftw_complex *out) {
fftw_execute_dft_r2c(plan,in,out);
}
@@ -792,20 +792,20 @@ public:
}
}
}
-
+
unsigned int Threads() {return std::max(T,threads);}
-
+
~fftwblock() {
if(plan2) fftw_destroy_plan(plan2);
}
};
-
+
// Compute the complex Fourier transform of M complex vectors, each of
// length n.
// Before calling fft(), the arrays in and out (which may coincide) must be
// allocated as Complex[M*n].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// mfft1d Forward(n,-1,M,stride,dist,in,out);
// Forward.fft(in,out);
@@ -823,14 +823,14 @@ public:
class mfft1d : public fftwblock<fftw_complex,fftw_complex>,
public Threadtable<keytype3,keyless3> {
static Table threadtable;
-public:
+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,
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) {}
-
+ (nx,M,stride,stride,dist,dist,in,out,threads) {}
+
mfft1d(unsigned int nx, int sign, unsigned int M,
size_t istride, size_t ostride, size_t idist, size_t odist,
Complex *in=NULL, Complex *out=NULL, unsigned int threads=maxthreads):
@@ -838,8 +838,8 @@ public:
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) {}
-
+ out,threads) {}
+
threaddata lookup(bool inplace, unsigned int threads) {
return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace));
}
@@ -847,13 +847,13 @@ public:
Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data);
}
};
-
+
// Compute the complex Fourier transform of n real values, using phase sign -1.
// Before calling fft(), the array in must be allocated as double[n] and
// the array out must be allocated as Complex[n/2+1]. The arrays in and out
// may coincide, allocated as Complex[n/2+1].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// rcfft1d Forward(n,in,out);
// Forward.fft(in,out);
@@ -862,7 +862,7 @@ public:
//
// rcfft1d Forward(n);
// Forward.fft(out);
-//
+//
// Notes:
// in contains the n real values stored as a Complex array;
// out contains the first n/2+1 Complex Fourier values.
@@ -870,36 +870,36 @@ public:
class rcfft1d : public fftw, public Threadtable<keytype1,keyless1> {
unsigned int nx;
static Table threadtable;
-public:
- rcfft1d(unsigned int nx, Complex *out=NULL, unsigned int threads=maxthreads)
+public:
+ 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,
- unsigned int threads=maxthreads)
+ unsigned int 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);
}
-
+
void Execute(Complex *in, Complex *out, bool=false) {
fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
}
};
-
+
// Compute the real inverse Fourier transform of the n/2+1 Complex values
// corresponding to the non-negative part of the frequency spectrum, using
// phase sign +1.
// Before calling fft(), the array in must be allocated as Complex[n/2+1]
// and the array out must be allocated as double[n]. The arrays in and out
-// may coincide, allocated as Complex[n/2+1].
+// may coincide, allocated as Complex[n/2+1].
//
// Out-of-place usage (input destroyed):
//
@@ -910,7 +910,7 @@ public:
//
// crfft1d Backward(n);
// Backward.fft(in);
-//
+//
// Notes:
// in contains the first n/2+1 Complex Fourier values.
// out contains the n real values stored as a Complex array;
@@ -918,25 +918,25 @@ public:
class crfft1d : public fftw, public Threadtable<keytype1,keyless1> {
unsigned int nx;
static Table threadtable;
-public:
- 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,
+public:
+ 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,
unsigned int threads=maxthreads)
- : fftw(realsize(nx,in,out),1,threads,nx), nx(nx) {Setup(in,out);}
-
+ : 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);
}
-
+
fftw_plan Plan(Complex *in, Complex *out) {
return fftw_plan_dft_c2r_1d(nx,(fftw_complex *) in,(double *) out,effort);
}
-
+
void Execute(Complex *in, Complex *out, bool=false) {
fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
}
@@ -948,7 +948,7 @@ public:
// Complex[M*(n/2+1)]. The arrays in and out may coincide,
// allocated as Complex[M*(n/2+1)].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// mrcfft1d Forward(n,M,istride,ostride,idist,odist,in,out);
// Forward.fft(in,out);
@@ -957,7 +957,7 @@ public:
//
// mrcfft1d Forward(n,M,istride,ostride,idist,odist);
// Forward.fft(out);
-//
+//
// Notes:
// istride is the spacing between the elements of each real vector;
// ostride is the spacing between the elements of each Complex vector;
@@ -974,28 +974,28 @@ public:
size_t istride, size_t ostride,
size_t idist, size_t odist,
double *in=NULL, Complex *out=NULL,
- unsigned int threads=maxthreads)
+ 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) {}
-
+
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);
}
-
+
void Normalize(Complex *out) {
fftw::Normalize<Complex>(nx/2+1,M,ostride,odist,out);
}
-
+
void fftNormalized(double *in, Complex *out=NULL, bool shift=false) {
fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out,false);
}
-
+
void fft0Normalized(double *in, Complex *out=NULL) {
fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out,true);
}
@@ -1006,7 +1006,7 @@ public:
// spectra, using phase sign +1. Before calling fft(), the array in must be
// allocated as Complex[M*(n/2+1)] and the array out must be allocated as
// double[M*n]. The arrays in and out may coincide,
-// allocated as Complex[M*(n/2+1)].
+// allocated as Complex[M*(n/2+1)].
//
// Out-of-place usage (input destroyed):
//
@@ -1017,7 +1017,7 @@ public:
//
// mcrfft1d Backward(n,M,istride,ostride,idist,odist);
// Backward.fft(out);
-//
+//
// Notes:
// stride is the spacing between the elements of each Complex vector;
// dist is the spacing between the first elements of the vectors;
@@ -1030,38 +1030,38 @@ class mcrfft1d : public fftwblock<fftw_complex,double>,
public:
mcrfft1d(unsigned int nx, unsigned int M, size_t istride, size_t ostride,
size_t idist, size_t odist, Complex *in=NULL, double *out=NULL,
- unsigned int threads=maxthreads)
+ 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) {}
-
+
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);
}
-
+
void Normalize(double *out) {
fftw::Normalize<double>(nx,M,ostride,odist,out);
}
-
+
void fftNormalized(Complex *in, double *out=NULL, bool shift=false) {
fftw::fftNormalized<Complex,double>(nx,M,ostride,odist,in,out,false);
}
-
+
void fft0Normalized(Complex *in, double *out=NULL) {
fftw::fftNormalized<Complex,double>(nx,M,ostride,odist,in,out,true);
}
};
-
+
// Compute the complex two-dimensional Fourier transform of nx times ny
// complex values. Before calling fft(), the arrays in and out (which may
// coincide) must be allocated as Complex[nx*ny].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// fft2d Forward(nx,ny,-1,in,out);
// Forward.fft(in,out);
@@ -1087,32 +1087,32 @@ class fft2d : public fftw, public Threadtable<keytype2,keyless2> {
unsigned int nx;
unsigned int ny;
static Table threadtable;
-public:
+public:
fft2d(unsigned int nx, unsigned int ny, int sign, Complex *in=NULL,
- Complex *out=NULL, unsigned int threads=maxthreads)
- : fftw(2*nx*ny,sign,threads), 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,
- unsigned int threads=maxthreads)
+ 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
-
+#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);
}
-
+
void Execute(Complex *in, Complex *out, bool=false) {
fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
}
@@ -1122,9 +1122,9 @@ public:
// values, using phase sign -1.
// Before calling fft(), the array in must be allocated as double[nx*ny] and
// the array out must be allocated as Complex[nx*(ny/2+1)]. The arrays in
-// and out may coincide, allocated as Complex[nx*(ny/2+1)].
+// and out may coincide, allocated as Complex[nx*(ny/2+1)].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// rcfft2d Forward(nx,ny,in,out);
// Forward.fft(in,out); // Origin of Fourier domain at (0,0)
@@ -1136,7 +1136,7 @@ public:
// rcfft2d Forward(nx,ny);
// Forward.fft(in); // Origin of Fourier domain at (0,0)
// Forward.fft0(in); // Origin of Fourier domain at (nx/2,0)
-//
+//
// Notes:
// in contains the nx*ny real values stored as a Complex array;
// out contains the upper-half portion (ky >= 0) of the Complex transform.
@@ -1144,22 +1144,22 @@ public:
class rcfft2d : public fftw {
unsigned int nx;
unsigned int ny;
-public:
+public:
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);}
-
+ 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,
- unsigned int threads=maxthreads)
+ unsigned int threads=maxthreads)
: fftw(2*nx*(ny/2+1),-1,threads,nx*ny), nx(nx), ny(ny) {
Setup(in,out);
- }
-
+ }
+
fftw_plan Plan(Complex *in, Complex *out) {
return fftw_plan_dft_r2c_2d(nx,ny,(double *) in,(fftw_complex *) out,
effort);
}
-
+
void Execute(Complex *in, Complex *out, bool shift=false) {
if(shift) {
if(inplace) Shift(in,nx,ny,threads);
@@ -1167,7 +1167,7 @@ public:
}
fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
}
-
+
// Set Nyquist modes of even shifted transforms to zero.
void deNyquist(Complex *f) {
unsigned int nyp=ny/2+1;
@@ -1185,14 +1185,14 @@ public:
f[(i+1)*nyp-1]=0.0;
}
};
-
+
// Compute the real two-dimensional inverse Fourier transform of the
// 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/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)].
+// allocated as Complex[nx*(ny/2+1)].
//
// Out-of-place usage (input destroyed):
//
@@ -1205,7 +1205,7 @@ public:
// crfft2d Backward(nx,ny);
// Backward.fft(in); // Origin of Fourier domain at (0,0)
// Backward.fft0(in); // Origin of Fourier domain at (nx/2,0)
-//
+//
// Notes:
// in contains the upper-half portion (ky >= 0) of the Complex transform;
// out contains the nx*ny real values stored as a Complex array.
@@ -1213,22 +1213,22 @@ public:
class crfft2d : public fftw {
unsigned int nx;
unsigned int ny;
-public:
+public:
crfft2d(unsigned int nx, unsigned int ny, double *out=NULL,
unsigned int threads=maxthreads) :
- fftw(2*nx*(ny/2+1),1,threads,nx*ny), nx(nx), ny(ny) {Setup(out);}
-
+ 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) {
Setup(in,out);
- }
-
+ }
+
fftw_plan Plan(Complex *in, Complex *out) {
return fftw_plan_dft_c2r_2d(nx,ny,(fftw_complex *) in,(double *) out,
effort);
}
-
+
void Execute(Complex *in, Complex *out, bool shift=false) {
fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
if(shift) {
@@ -1236,7 +1236,7 @@ public:
else Shift((double *) out,nx,ny,threads);
}
}
-
+
// Set Nyquist modes of even shifted transforms to zero.
void deNyquist(Complex *f) {
unsigned int nyp=ny/2+1;
@@ -1255,11 +1255,11 @@ public:
}
};
-// Compute the complex three-dimensional Fourier transform of
+// Compute the complex three-dimensional Fourier transform of
// nx times ny times nz complex values. Before calling fft(), the arrays in
// and out (which may coincide) must be allocated as Complex[nx*ny*nz].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// fft3d Forward(nx,ny,nz,-1,in,out);
// Forward.fft(in,out);
@@ -1286,20 +1286,20 @@ class fft3d : public fftw {
unsigned int nx;
unsigned int ny;
unsigned int nz;
-public:
+public:
fft3d(unsigned int nx, unsigned int ny, unsigned int nz,
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);}
-
+ 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,
- unsigned int threads=maxthreads)
- : fftw(2*in.Size(),sign,threads), nx(in.Nx()), ny(in.Ny()), nz(in.Nz())
+ unsigned int threads=maxthreads)
+ : fftw(2*in.Size(),sign,threads), nx(in.Nx()), ny(in.Ny()), nz(in.Nz())
{Setup(in,out);}
-#endif
-
+#endif
+
fftw_plan Plan(Complex *in, Complex *out) {
return fftw_plan_dft_3d(nx,ny,nz,(fftw_complex *) in,
(fftw_complex *) out, sign, effort);
@@ -1310,9 +1310,9 @@ public:
// nx times ny times nz real values, using phase sign -1.
// Before calling fft(), the array in must be allocated as double[nx*ny*nz]
// and the array out must be allocated as Complex[nx*ny*(nz/2+1)]. The
-// arrays in and out may coincide, allocated as Complex[nx*ny*(nz/2+1)].
+// arrays in and out may coincide, allocated as Complex[nx*ny*(nz/2+1)].
//
-// Out-of-place usage:
+// Out-of-place usage:
//
// rcfft3d Forward(nx,ny,nz,in,out);
// Forward.fft(in,out); // Origin of Fourier domain at (0,0)
@@ -1323,7 +1323,7 @@ public:
// rcfft3d Forward(nx,ny,nz);
// Forward.fft(in); // Origin of Fourier domain at (0,0)
// Forward.fft0(in); // Origin of Fourier domain at (nx/2,ny/2,0)
-//
+//
// Notes:
// in contains the nx*ny*nz real values stored as a Complex array;
// out contains the upper-half portion (kz >= 0) of the Complex transform.
@@ -1332,23 +1332,23 @@ class rcfft3d : public fftw {
unsigned int nx;
unsigned int ny;
unsigned int nz;
-public:
+public:
rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *out=NULL,
unsigned int 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)
+ Complex *out=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,out);}
-
+ nx(nx), ny(ny), nz(nz) {Setup(in,out);}
+
fftw_plan Plan(Complex *in, Complex *out) {
return fftw_plan_dft_r2c_3d(nx,ny,nz,(double *) in,(fftw_complex *) out,
effort);
}
-
+
void Execute(Complex *in, Complex *out, bool shift=false) {
if(shift) {
if(inplace) Shift(in,nx,ny,nz,threads);
@@ -1356,7 +1356,7 @@ public:
}
fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
}
-
+
// Set Nyquist modes of even shifted transforms to zero.
void deNyquist(Complex *f) {
unsigned int nzp=nz/2+1;
@@ -1368,7 +1368,7 @@ public:
for(unsigned int k=0; k < yz; ++k)
f[k]=0.0;
}
-
+
if(ny % 2 == 0) {
#ifndef FFTWPP_SINGLE_THREAD
#pragma omp parallel for num_threads(threads)
@@ -1379,7 +1379,7 @@ public:
f[iyz+k]=0.0;
}
}
-
+
if(nz % 2 == 0)
#ifndef FFTWPP_SINGLE_THREAD
#pragma omp parallel for num_threads(threads)
@@ -1389,14 +1389,14 @@ public:
f[i*yz+(j+1)*nzp-1]=0.0;
}
};
-
+
// Compute the real two-dimensional inverse Fourier transform of the
// nx*ny*(nz/2+1) Complex values corresponding to the spectral values in the
// half-plane kz >= 0, using phase sign +1.
// Before calling fft(), the array in must be allocated as
// Complex[nx*ny*(nz+1)/2] and the array out must be allocated as
// double[nx*ny*nz]. The arrays in and out may coincide,
-// allocated as Complex[nx*ny*(nz/2+1)].
+// allocated as Complex[nx*ny*(nz/2+1)].
//
// Out-of-place usage (input destroyed):
//
@@ -1409,7 +1409,7 @@ public:
// crfft3d Backward(nx,ny,nz);
// Backward.fft(in); // Origin of Fourier domain at (0,0)
// Backward.fft0(in); // Origin of Fourier domain at (nx/2,ny/2,0)
-//
+//
// Notes:
// in contains the upper-half portion (kz >= 0) of the Complex transform;
// out contains the nx*ny*nz real values stored as a Complex array.
@@ -1418,22 +1418,22 @@ class crfft3d : public fftw {
unsigned int nx;
unsigned int ny;
unsigned int nz;
-public:
+public:
crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, double *out=NULL,
- unsigned int threads=maxthreads)
+ unsigned int threads=maxthreads)
: fftw(2*nx*ny*(nz/2+1),1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz)
- {Setup(out);}
-
+ {Setup(out);}
+
crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in,
- double *out=NULL, unsigned int threads=maxthreads)
+ 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);}
-
+ 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,
effort);
}
-
+
void Execute(Complex *in, Complex *out, bool shift=false) {
fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
if(shift) {
@@ -1441,7 +1441,7 @@ public:
else Shift((double *) out,nx,ny,nz,threads);
}
}
-
+
// Set Nyquist modes of even shifted transforms to zero.
void deNyquist(Complex *f) {
unsigned int nzp=nz/2+1;
@@ -1453,7 +1453,7 @@ public:
for(unsigned int k=0; k < yz; ++k)
f[k]=0.0;
}
-
+
if(ny % 2 == 0) {
#ifndef FFTWPP_SINGLE_THREAD
#pragma omp parallel for num_threads(threads)
@@ -1464,7 +1464,7 @@ public:
f[iyz+k]=0.0;
}
}
-
+
if(nz % 2 == 0)
#ifndef FFTWPP_SINGLE_THREAD
#pragma omp parallel for num_threads(threads)