summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/fftw++.h
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2010-06-08 13:46:05 +0000
committerKarl Berry <karl@freefriends.org>2010-06-08 13:46:05 +0000
commita960e44eb527236f39aec81babc0474911a86078 (patch)
tree9950eca71791d90820a80a521a7cc252c0955db5 /Build/source/utils/asymptote/fftw++.h
parent6443467452320c296faa1f43f0606a9457bd4463 (diff)
asy 1.96
git-svn-id: svn://tug.org/texlive/trunk@18817 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/fftw++.h')
-rw-r--r--Build/source/utils/asymptote/fftw++.h207
1 files changed, 114 insertions, 93 deletions
diff --git a/Build/source/utils/asymptote/fftw++.h b/Build/source/utils/asymptote/fftw++.h
index 658768551be..1b4373a08b6 100644
--- a/Build/source/utils/asymptote/fftw++.h
+++ b/Build/source/utils/asymptote/fftw++.h
@@ -18,7 +18,7 @@
#ifndef __fftwpp_h__
#define __fftwpp_h__ 1
-#define __FFTWPP_H_VERSION__ 1.05
+#define __FFTWPP_H_VERSION__ 1.08svn
#include <cstdlib>
#include <fstream>
@@ -31,10 +31,6 @@
typedef std::complex<double> Complex;
#endif
-#ifndef M_PI
-#define M_PI acos(-1.0)
-#endif
-
#ifndef HAVE_POSIX_MEMALIGN
#ifdef __GLIBC_PREREQ
@@ -50,17 +46,12 @@ typedef std::complex<double> Complex;
#endif
#ifdef __Array_h__
-using Array::array1;
-using Array::array2;
-using Array::array3;
-
-static array1<Complex> NULL1;
-static array2<Complex> NULL2;
-static array3<Complex> NULL3;
-static array1<double> NULL1D;
-static array2<double> NULL2D;
-static array3<double> NULL3D;
+namespace Array {
+static const array1<Complex> NULL1;
+static const array2<Complex> NULL2;
+static const array3<Complex> NULL3;
+}
#else
@@ -69,6 +60,8 @@ static array3<double> NULL3D;
extern "C" int posix_memalign(void **memptr, size_t alignment, size_t size);
#endif
#else
+namespace Array {
+
// Adapted from FFTW aligned malloc/free. Assumes that malloc is at least
// sizeof(void*)-aligned. Allocated memory must be freed with free0.
inline int posix_memalign0(void **memptr, size_t alignment, size_t size)
@@ -87,8 +80,12 @@ inline void free0(void *p)
{
if(p) free(*((void **) p-1));
}
+
+}
#endif
+namespace Array {
+
template<class T>
inline void newAlign(T *&v, size_t len, size_t align)
{
@@ -117,20 +114,23 @@ inline void deleteAlign(T *v, size_t len)
free0(v);
#endif
}
+}
#endif
+namespace fftwpp {
+
inline Complex *ComplexAlign(size_t size)
{
Complex *v;
- newAlign(v,size,sizeof(Complex));
+ Array::newAlign(v,size,sizeof(Complex));
return v;
}
inline double *doubleAlign(size_t size)
{
double *v;
- newAlign(v,size,sizeof(Complex));
+ Array::newAlign(v,size,sizeof(Complex));
return v;
}
@@ -140,7 +140,7 @@ inline void deleteAlign(T *p)
#ifdef HAVE_POSIX_MEMALIGN
free(p);
#else
- free0(p);
+ Array::free0(p);
#endif
}
@@ -193,6 +193,7 @@ protected:
static std::ifstream ifWisdom;
static std::ofstream ofWisdom;
static bool Wise;
+ static const double twopi;
public:
// Shift the Fourier origin to (nx/2,0).
@@ -208,7 +209,7 @@ public:
} else {
if(sign) {
unsigned int c=nx/2;
- double arg=2.0*M_PI*c/nx;
+ double arg=twopi*c/nx;
for(unsigned int i=0; i < nx; i++) {
double iarg=i*arg;
Complex zeta(cos(iarg),sign*sin(iarg));
@@ -243,9 +244,8 @@ public:
} else {
if(sign) {
unsigned int cx=nx/2;
- unsigned int cy=ny/2;
- double twopi=2.0*M_PI;
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;
@@ -302,7 +302,7 @@ public:
exit(1);
}
- if(alloc) deleteAlign(in,(doubles+1)/2);
+ if(alloc) Array::deleteAlign(in,(doubles+1)/2);
SaveWisdom();
}
@@ -326,7 +326,7 @@ public:
fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
}
- void Setout(Complex *in, Complex *&out) {
+ Complex *Setout(Complex *in, Complex *out) {
#ifndef NO_CHECK_ALIGN
CheckAlign(in,"input");
if(out) CheckAlign(out,"output");
@@ -338,10 +338,11 @@ public:
std::cerr << "ERROR: fft constructor and call must be both in place or both out of place" << std::endl;
exit(1);
}
+ return out;
}
void fft(Complex *in, Complex *out=NULL) {
- Setout(in,out);
+ out=Setout(in,out);
Execute(in,out);
}
@@ -354,7 +355,7 @@ public:
}
void fft0(Complex *in, Complex *out=NULL) {
- Setout(in,out);
+ out=Setout(in,out);
Execute(in,out,true);
}
@@ -376,13 +377,15 @@ public:
}
virtual void fftNormalized(Complex *in, Complex *out=NULL) {
- Setout(in,out);
+ out=Setout(in,out);
Execute(in,out);
Normalize(out);
}
void fftNormalized(Complex *in, double *out) {
- fftNormalized(in,(Complex *) out);
+ out=(double *) Setout(in,(Complex *) out);
+ Execute(in,(Complex *) out);
+ Normalize(out);
}
void fftNormalized(double *in, Complex *out) {
@@ -390,13 +393,15 @@ public:
}
void fft0Normalized(Complex *in, Complex *out=NULL) {
- Setout(in,out);
+ out=Setout(in,out);
Execute(in,out,true);
Normalize(out);
}
void fft0Normalized(Complex *in, double *out) {
- fft0Normalized(in,(Complex *) out);
+ out=(double *) Setout(in,(Complex *) out);
+ Execute(in,(Complex *) out,true);
+ Normalize(out);
}
void fft0Normalized(double *in, Complex *out) {
@@ -404,14 +409,14 @@ public:
}
void fftNormalized(Complex *in, Complex *out,
- unsigned int nx, unsigned int m,
+ unsigned int nx, unsigned int M,
unsigned int stride, unsigned int dist) {
if(stride == 1 && dist == nx) fftw::fftNormalized(in,out);
else if(stride == nx && dist == 1) fftw::fftNormalized(in,out);
else {
- Setout(in,out);
+ out=Setout(in,out);
Execute(in,out);
- for(unsigned int k=0; k < m; k++) {
+ for(unsigned int k=0; k < M; k++) {
for(unsigned int j=0; j < nx; j++) {
out[j*stride+k*dist] *= norm;
}
@@ -451,7 +456,8 @@ public:
: fftw(2*nx,sign), nx(nx) {Setup(in,out);}
#ifdef __Array_h__
- fft1d(int sign, const array1<Complex>& in, const array1<Complex>& out=NULL1)
+ 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);}
#endif
@@ -461,19 +467,19 @@ public:
}
};
-// Compute the complex Fourier transform of m complex vectors, each of
+// 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].
+// allocated as Complex[M*n].
//
// Out-of-place usage:
//
-// mfft1d Forward(n,-1,m,stride,dist,in,out);
+// mfft1d Forward(n,-1,M,stride,dist,in,out);
// Forward.fft(in,out);
//
// In-place usage:
//
-// mfft1d Forward(n,-1,m,stride,dist);
+// mfft1d Forward(n,-1,M,stride,dist);
// Forward.fft(in);
//
// Notes:
@@ -483,26 +489,26 @@ public:
//
class mfft1d : public fftw {
unsigned int nx;
- unsigned int m;
+ unsigned int M;
unsigned int stride;
unsigned int dist;
public:
- mfft1d(unsigned int nx, int sign, unsigned int m=1, unsigned int stride=1,
+ mfft1d(unsigned int nx, int sign, unsigned int M=1, unsigned int stride=1,
unsigned int 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))
+ : 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);}
fftw_plan Plan(Complex *in, Complex *out) {
int n[1]={nx};
- return fftw_plan_many_dft(1,n,m,
+ return fftw_plan_many_dft(1,n,M,
(fftw_complex *) in,NULL,stride,dist,
(fftw_complex *) out,NULL,stride,dist,
sign,effort);
}
void fftNormalized(Complex *in, Complex *out=NULL) {
- fftw::fftNormalized(in,out,nx,m,stride,dist);
+ fftw::fftNormalized(in,out,nx,M,stride,dist);
}
};
@@ -535,11 +541,11 @@ public:
: fftw(realsize(nx,in,out),-1,nx), nx(nx) {Setup(in,out);}
#ifdef __Array_h__
- rcfft1d(unsigned int nx, const array1<Complex>& out)
+ rcfft1d(unsigned int nx, const Array::array1<Complex>& out)
: fftw(out.Size(),-1,nx), nx(nx) {Setup(out);}
- rcfft1d(unsigned int nx, const array1<double>& in,
- const array1<Complex>& out=NULL1)
+ 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);}
#endif
@@ -583,12 +589,15 @@ public:
: fftw(realsize(nx,in,out),1,nx), nx(nx) {Setup(in,out);}
#ifdef __Array_h__
- crfft1d(unsigned int nx, const array1<double>& out)
+ crfft1d(unsigned int nx, const Array::array1<double>& out)
: fftw(out.Size(),1,nx), nx(nx) {Setup(out);}
- crfft1d(unsigned int nx, const array1<Complex>& in,
- const array1<double>& out=NULL1D)
- : fftw(2*in.Size(),1,nx), nx(nx) {Setup(in,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,
+ const Array::array1<double>& out)
+ : fftw(out.Size(),1,nx), nx(nx) {Setup(in,out);}
#endif
fftw_plan Plan(Complex *in, Complex *out) {
@@ -600,20 +609,20 @@ public:
}
};
-// Compute the real Fourier transform of m real vectors, each of length n,
+// 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
-// Complex[m*(n/2+1)]. The arrays in and out may coincide,
-// allocated as Complex[m*(n/2+1)].
+// allocated as double[M*n] and the array out must be allocated as
+// Complex[M*(n/2+1)]. The arrays in and out may coincide,
+// allocated as Complex[M*(n/2+1)].
//
// Out-of-place usage:
//
-// mrcfft1d Forward(n,m,stride,dist,in,out);
+// mrcfft1d Forward(n,M,stride,dist,in,out);
// Forward.fft(in,out);
//
// In-place usage:
//
-// mrcfft1d Forward(n,m,stride,dist);
+// mrcfft1d Forward(n,M,stride,dist);
// Forward.fft(out);
//
// Notes:
@@ -624,23 +633,23 @@ public:
//
class mrcfft1d : public fftw {
unsigned int nx;
- unsigned int m;
+ unsigned int M;
unsigned int stride;
unsigned int dist;
public:
- mrcfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
+ mrcfft1d(unsigned int nx, unsigned int M=1, unsigned int stride=1,
unsigned int dist=0, Complex *out=NULL)
- : fftw(2*(nx/2*stride+(m-1)*Dist(nx,stride,dist)+1),-1,nx), nx(nx), m(m),
+ : 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);}
- mrcfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
+ mrcfft1d(unsigned int nx, unsigned int M=1, unsigned int stride=1,
unsigned int 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),
+ : 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);}
fftw_plan Plan(Complex *in, Complex *out) {
const int n[1]={nx};
- return fftw_plan_many_dft_r2c(1,n,m,
+ return fftw_plan_many_dft_r2c(1,n,M,
(double *) in,NULL,stride,2*dist,
(fftw_complex *) out,NULL,stride,dist,
effort);
@@ -651,25 +660,25 @@ public:
}
void fftNormalized(Complex *in, Complex *out=NULL) {
- fftw::fftNormalized(in,out,nx/2+1,m,stride,dist);
+ fftw::fftNormalized(in,out,nx/2+1,M,stride,dist);
}
};
-// Compute the real inverse Fourier transform of m complex vectors, each of
+// Compute the real inverse Fourier transform of M complex vectors, each of
// length n/2+1, corresponding to the non-negative parts of the frequency
// 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)] 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)].
//
// Out-of-place usage (input destroyed):
//
-// mcrfft1d Backward(n,m,stride,dist,in,out);
+// mcrfft1d Backward(n,M,stride,dist,in,out);
// Backward.fft(in,out);
//
// In-place usage:
//
-// mcrfft1d Backward(n,m,stride,dist);
+// mcrfft1d Backward(n,M,stride,dist);
// Backward.fft(out);
//
// Notes:
@@ -680,18 +689,18 @@ public:
//
class mcrfft1d : public fftw {
unsigned int nx;
- unsigned int m;
+ unsigned int M;
unsigned int stride;
unsigned int dist;
public:
- mcrfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
+ mcrfft1d(unsigned int nx, unsigned int M=1, unsigned int stride=1,
unsigned int 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),
- nx(nx), m(m), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
+ : fftw((realsize(nx,in,out)-2)*stride+2*(M-1)*Dist(nx,stride,dist)+2,1,nx),
+ nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
fftw_plan Plan(Complex *in, Complex *out) {
const int n[1]={nx};
- return fftw_plan_many_dft_c2r(1,n,m,
+ return fftw_plan_many_dft_c2r(1,n,M,
(fftw_complex *) in,NULL,stride,dist,
(double *) out,NULL,stride,2*dist,
effort);
@@ -702,7 +711,7 @@ public:
}
void fftNormalized(Complex *in, Complex *out=NULL) {
- fftw::fftNormalized(in,out,(nx/2+1),m,stride,dist);
+ fftw::fftNormalized(in,out,(nx/2+1),M,stride,dist);
}
};
@@ -741,7 +750,8 @@ public:
: fftw(2*nx*ny,sign), nx(nx), ny(ny) {Setup(in,out);}
#ifdef __Array_h__
- fft2d(int sign, const array2<Complex>& in, const array2<Complex>& out=NULL2)
+ 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);}
#endif
@@ -787,11 +797,11 @@ public:
: fftw(nx*realsize(ny,in,out),-1,nx*ny), nx(nx), ny(ny) {Setup(in,out);}
#ifdef __Array_h__
- rcfft2d(unsigned int ny, const array2<Complex>& out)
+ 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 array2<double>& in,
- const array2<Complex>& out=NULL2)
+ 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);}
#endif
@@ -842,12 +852,15 @@ public:
fftw(nx*realsize(ny,in,out),1,nx*ny), nx(nx), ny(ny) {Setup(in,out);}
#ifdef __Array_h__
- crfft2d(unsigned int ny, const array2<double>& out)
+ 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 array2<Complex>& in,
- const array2<double>& out=NULL2D)
- : fftw(2*in.Size(),1,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in,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<Complex>& in,
+ const Array::array2<double>& out)
+ : fftw(out.Size(),1,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in,out);}
#endif
fftw_plan Plan(Complex *in, Complex *out) {
@@ -898,7 +911,8 @@ public:
: fftw(2*nx*ny*nz,sign), nx(nx), ny(ny), nz(nz) {Setup(in,out);}
#ifdef __Array_h__
- fft3d(int sign, const array3<Complex>& in, const array3<Complex>& out=NULL3)
+ 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())
{Setup(in,out);}
#endif
@@ -943,12 +957,12 @@ public:
nx(nx), ny(ny), nz(nz) {Setup(in,out);}
#ifdef __Array_h__
- rcfft3d(unsigned int nz, const array3<Complex>& out)
+ rcfft3d(unsigned int nz, const Array::array3<Complex>& out)
: fftw(out.Size(),-1,out.Nx()*out.Ny()*nz),
nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);}
- rcfft3d(unsigned int nz, const array3<double>& in,
- const array3<Complex>& out=NULL3)
+ 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()),
nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);}
#endif
@@ -1001,15 +1015,20 @@ public:
double *out=NULL) :
fftw(nx*ny*(realsize(nz,in,out)),1,nx*ny*nz), nx(nx), ny(ny), nz(nz)
{Setup(in,out);}
+
#ifdef __Array_h__
- crfft3d(unsigned int nz, const array3<double>& out) :
- fftw(out.Size(),1,out.Nx()*out.Ny()*nz),
+ crfft3d(unsigned int nz, const Array::array3<double>& out)
+ : fftw(out.Size(),1,out.Nx()*out.Ny()*nz),
nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);}
- crfft3d(unsigned int nz, const array3<Complex>& in,
- const array3<double>& out=NULL3D) :
- fftw(2*in.Size(),1,in.Nx()*in.Ny()*nz),
- nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);}
+ crfft3d(unsigned int nz, const Array::array3<Complex>& in)
+ : fftw(2*in.Size(),1,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),
+ nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);}
#endif
fftw_plan Plan(Complex *in, Complex *out) {
@@ -1023,4 +1042,6 @@ public:
}
};
+}
+
#endif