diff options
Diffstat (limited to 'Build/source/utils/asymptote/fftw++.h')
-rw-r--r-- | Build/source/utils/asymptote/fftw++.h | 1005 |
1 files changed, 465 insertions, 540 deletions
diff --git a/Build/source/utils/asymptote/fftw++.h b/Build/source/utils/asymptote/fftw++.h index 35dfc6f249b..9ce6d4430e8 100644 --- a/Build/source/utils/asymptote/fftw++.h +++ b/Build/source/utils/asymptote/fftw++.h @@ -1,6 +1,7 @@ /* Fast Fourier transform C++ header class for the FFTW3 Library - Copyright (C) 2004-13 John C. Bowman, University of Alberta - + Copyright (C) 2004-15 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 @@ -18,7 +19,7 @@ #ifndef __fftwpp_h__ #define __fftwpp_h__ 1 -#define __FFTWPP_H_VERSION__ 1.14pre +#define __FFTWPP_H_VERSION__ 2.00 #include <cstdlib> #include <fstream> @@ -27,6 +28,12 @@ #include <cerrno> #include <map> +#ifndef _OPENMP +#ifndef FFTWPP_SINGLE_THREAD +#define FFTWPP_SINGLE_THREAD +#endif +#endif + #ifndef FFTWPP_SINGLE_THREAD #include <omp.h> #endif @@ -49,143 +56,42 @@ inline int 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 + #ifndef __Complex_h__ #include <complex> typedef std::complex<double> Complex; #endif #include "seconds.h" - -#ifndef HAVE_POSIX_MEMALIGN - -#ifdef __GLIBC_PREREQ -#if __GLIBC_PREREQ(2,3) -#define HAVE_POSIX_MEMALIGN -#endif -#else -#ifdef _POSIX_SOURCE -#define HAVE_POSIX_MEMALIGN -#endif -#endif - -#endif - -#ifdef __Array_h__ - -namespace Array { -static const array1<Complex> NULL1; -static const array2<Complex> NULL2; -static const array3<Complex> NULL3; -} - -#else - -#ifdef HAVE_POSIX_MEMALIGN -#ifdef _AIX -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) -{ - if(alignment % sizeof (void *) != 0 || (alignment & (alignment - 1)) != 0) - return EINVAL; - void *p0=malloc(size+alignment); - if(!p0) return ENOMEM; - void *p=(void *)(((size_t) p0+alignment)&~(alignment-1)); - *((void **) p-1)=p0; - *memptr=p; - return 0; -} - -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) -{ - void *mem=NULL; - const char *invalid="Invalid alignment requested"; - const char *nomem="Memory limits exceeded"; -#ifdef HAVE_POSIX_MEMALIGN - int rc=posix_memalign(&mem,align,len*sizeof(T)); -#else - int rc=posix_memalign0(&mem,align,len*sizeof(T)); -#endif - if(rc == EINVAL) std::cerr << invalid << std::endl; - if(rc == ENOMEM) std::cerr << nomem << std::endl; - v=(T *) mem; - for(size_t i=0; i < len; i++) new(v+i) T; -} - -template<class T> -inline void deleteAlign(T *v, size_t len) -{ - for(size_t i=len; i-- > 0;) v[i].~T(); -#ifdef HAVE_POSIX_MEMALIGN - free(v); -#else - free0(v); -#endif -} -} - -#endif +#include "statistics.h" +#include "align.h" namespace fftwpp { -inline Complex *ComplexAlign(size_t size) -{ - Complex *v; - Array::newAlign(v,size,sizeof(Complex)); - return v; -} - -inline double *doubleAlign(size_t size) -{ - double *v; - Array::newAlign(v,size,sizeof(Complex)); - return v; -} - -template<class T> -inline void deleteAlign(T *p) -{ -#ifdef HAVE_POSIX_MEMALIGN - free(p); -#else - Array::free0(p); -#endif -} - // Obsolete names: #define FFTWComplex ComplexAlign #define FFTWdouble doubleAlign #define FFTWdelete deleteAlign -inline void fftwpp_export_wisdom(void (*emitter)(char c, std::ofstream& s), - std::ofstream& s) -{ - fftw_export_wisdom((void (*) (char, void *)) emitter,(void *) &s); -} +class fftw; -inline int fftwpp_import_wisdom(int (*g)(std::ifstream& s), std::ifstream &s) -{ - return fftw_import_wisdom((int (*) (void *)) g,(void *) &s); -} - -inline void PutWisdom(char c, std::ofstream& s) {s.put(c);} -inline int GetWisdom(std::ifstream& s) {return s.get();} +extern "C" fftw_plan Planner(fftw *F, Complex *in, Complex *out); +void LoadWisdom(); +void SaveWisdom(); extern const char *inout; @@ -198,9 +104,32 @@ struct threaddata { threads(threads), mean(mean), stdev(stdev) {} }; +class fftw; + +class ThreadBase +{ +protected: + unsigned int threads; + unsigned int innerthreads; +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; + } else { + innerthreads=threads; + threads=1; + } + } +}; + // Base clase for fft routines // -class fftw { +class fftw : public ThreadBase { protected: unsigned int doubles; // number of double precision values in dataset int sign; @@ -226,7 +155,6 @@ protected: return realsize(n,(Complex *) in,out); } - static unsigned int Wise; static const double twopi; public: @@ -234,7 +162,7 @@ public: static unsigned int maxthreads; static double testseconds; static const char *WisdomName; - static bool mpi; // TODO: Remove when FFTW transpose routines are replaced + static fftw_plan (*planner)(fftw *f, Complex *in, Complex *out); virtual unsigned int Threads() {return threads;} @@ -243,7 +171,7 @@ public: // 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) { - const unsigned int nyp=ny/2+1; + unsigned int nyp=ny/2+1; unsigned int stop=nx*nyp; if(nx % 2 == 0) { unsigned int inc=2*nyp; @@ -282,10 +210,10 @@ public: // Inplace shift of Fourier origin to (nx/2,ny/2,0) for even nx and ny. static void Shift(Complex *data, unsigned int nx, unsigned int ny, unsigned int nz, unsigned int threads) { - const unsigned int nzp=nz/2+1; - const unsigned int nyzp=ny*nzp; + unsigned int nzp=nz/2+1; + unsigned int nyzp=ny*nzp; if(nx % 2 == 0 && ny % 2 == 0) { - const unsigned int pinc=2*nzp; + unsigned int pinc=2*nzp; Complex *pstop=data; Complex *p=data; #ifndef FFTWPP_SINGLE_THREAD @@ -300,7 +228,7 @@ public: } } } else { - std::cerr << oddshift << "or odd ny" << std::endl; + std::cerr << oddshift << " or odd ny" << std::endl; exit(1); } } @@ -308,9 +236,9 @@ 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) { - const unsigned int nyz=ny*nz; + unsigned int nyz=ny*nz; if(nx % 2 == 0 && ny % 2 == 0) { - const unsigned int pinc=2*nz; + unsigned int pinc=2*nz; double *pstop=data; double *p=data; #ifndef FFTWPP_SINGLE_THREAD @@ -325,29 +253,26 @@ public: } } } else { - std::cerr << oddshift << "or odd ny" << std::endl; + std::cerr << oddshift << " or odd ny" << std::endl; exit(1); } } + fftw() : plan(NULL) {} fftw(unsigned int doubles, int sign, unsigned int threads, unsigned int n=0) : doubles(doubles), sign(sign), threads(threads), norm(1.0/(n ? n : doubles/2)), plan(NULL) { #ifndef FFTWPP_SINGLE_THREAD - if(!mpi) fftw_init_threads(); + fftw_init_threads(); #endif } virtual ~fftw() { - if(Wise > 0) - --Wise; - if(Wise == 0) - SaveWisdom(); if(plan) fftw_destroy_plan(plan); } - virtual fftw_plan Plan(Complex *in, Complex *out)=0; + 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; @@ -367,51 +292,10 @@ public: #endif } - class statistics { - unsigned int N; - double A; - double varL; - double varH; - public: - statistics() : N(0), A(0.0), varL(0.0), varH(0.0) {} - double count() {return N;} - double mean() {return A;} - void add(double t) { - ++N; - double diff=t-A; - A += diff/N; - double v=diff*(t-A); - if(diff < 0.0) - varL += v; - else - varH += v; - } - double stdev(double var, double f) { - double factor=N > f ? f/(N-f) : 0.0; - return sqrt(var*factor); - } - double stdev() { - return stdev(varL+varH,1.0); - } - double stdevL() { - return stdev(varL,2.0); - } - double stdevH() { - return stdev(varH,2.0); - } - void output(const char *text, unsigned int m) { - std::cout << text << ":\n" - << m << "\t" - << A << "\t" - << stdevL() << "\t" - << stdevH() << std::endl; - } - }; - threaddata time(fftw_plan plan1, fftw_plan planT, Complex *in, Complex *out, unsigned int Threads) { - statistics S,ST; - double stop=totalseconds()+testseconds; + utils::statistics S,ST; + double stop=utils::totalseconds()+testseconds; threads=1; plan=plan1; fft(in,out); @@ -420,17 +304,17 @@ public: fft(in,out); unsigned int N=1; for(;;) { - double t0=totalseconds(); + double t0=utils::totalseconds(); threads=1; plan=plan1; for(unsigned int i=0; i < N; ++i) fft(in,out); - double t1=totalseconds(); + double t1=utils::totalseconds(); threads=Threads; plan=planT; for(unsigned int i=0; i < N; ++i) fft(in,out); - double t=totalseconds(); + double t=utils::totalseconds(); S.add(t1-t0); ST.add(t-t1); if(S.mean() < 100.0/CLOCKS_PER_SEC) N *= 2; @@ -446,7 +330,7 @@ public: if(diff < -error) { threads=Threads; fftw_destroy_plan(plan1); - break; + break; } } } @@ -458,18 +342,22 @@ public: } virtual void store(bool inplace, const threaddata& data) {} - threaddata Setup(Complex *in, Complex *out=NULL) { - if(!Wise) {LoadWisdom(); ++Wise;} - - bool alloc=!in; - if(alloc) in=ComplexAlign((doubles+1)/2); + inline Complex *CheckAlign(Complex *in, Complex *out, bool constructor=true) + { #ifndef NO_CHECK_ALIGN - CheckAlign(in,"constructor input"); - if(out) CheckAlign(out,"constructor output"); + CheckAlign(in,constructor ? "constructor input" : "input"); + if(out) CheckAlign(out,constructor ? "constructor output" : "output"); else out=in; #else if(!out) out=in; #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; @@ -477,44 +365,41 @@ public: if(threads > 1) data=lookup(inplace,threads); threads=data.threads > 0 ? data.threads : 1; planThreads(threads); - plan=Plan(in,out); + plan=(*planner)(this,in,out); if(!plan) noplan(); - if(Threads > 1 && data.threads == 0) { + fftw_plan planT; + if(Threads > 1) { threads=Threads; planThreads(threads); - fftw_plan planT=Plan(in,out); - if(planT) - data=time(plan,planT,in,out,threads); - else noplan(); - store(inplace,threaddata(threads,data.mean,data.stdev)); + 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); return data; } - void Setup(Complex *in, double *out) { - Setup(in,(Complex *) out);} - void Setup(double *in, Complex *out=NULL) { - Setup((Complex *) in,out); + threaddata Setup(Complex *in, double *out) { + return Setup(in,(Complex *) out); + } + + threaddata Setup(double *in, Complex *out=NULL) { + return Setup((Complex *) in,out); } - - static void LoadWisdom(); - static void SaveWisdom(); virtual void Execute(Complex *in, Complex *out, bool=false) { fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out); } Complex *Setout(Complex *in, Complex *out) { -#ifndef NO_CHECK_ALIGN - CheckAlign(in,"input"); - if(out) CheckAlign(out,"output"); - else out=in; -#else - if(!out) out=in; -#endif + out=CheckAlign(in,out,false); if(inplace ^ (out == in)) { std::cerr << "ERROR: fft " << inout << std::endl; exit(1); @@ -549,7 +434,7 @@ public: } void Normalize(Complex *out) { - unsigned int stop=(doubles+1)/2; + unsigned int stop=doubles/2; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(threads) #endif @@ -563,67 +448,56 @@ public: for(unsigned int i=0; i < doubles; i++) out[i] *= norm; } - virtual void fftNormalized(Complex *in, Complex *out=NULL) { + virtual void fftNormalized(Complex *in, Complex *out=NULL, bool shift=false) + { out=Setout(in,out); - Execute(in,out); + Execute(in,out,shift); Normalize(out); } - void fftNormalized(Complex *in, double *out) { + void fftNormalized(Complex *in, double *out, bool shift=false) { out=(double *) Setout(in,(Complex *) out); - Execute(in,(Complex *) out); + Execute(in,(Complex *) out,shift); Normalize(out); } - void fftNormalized(double *in, Complex *out) { - fftNormalized((Complex *) in,out); - } - - 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 fftNormalized(double *in, Complex *out, bool shift=false) { + fftNormalized((Complex *) in,out,shift); } - void fft0Normalized(double *in, Complex *out) { - fft0Normalized((Complex *) in,out); + template<class I, class O> + void fft0Normalized(I in, O out) { + fftNormalized(in,out,true); } - void fftNormalized(Complex *in, Complex *out, - unsigned int nx, unsigned int M, - size_t stride, size_t dist) { - if(stride == 1 && dist == nx) fftw::fftNormalized(in,out); - else { - out=Setout(in,out); - Execute(in,out); + 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) { + out=(O *) Setout((Complex *) in,(Complex *) out); + Execute((Complex *) in,(Complex *) out,shift); - unsigned int stop=nx*stride; - Complex *outMdist=out+M*dist; + unsigned int 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 += stride) { - Complex *pstop=outMdist+i; - for(Complex *p=out+i; p < pstop; p += dist) { - *p *= norm; - } + for(unsigned int i=0; i < stop; i += ostride) { + O *pstop=outMdist+i; + for(O *p=out+i; p < pstop; p += odist) { + *p *= norm; } } } - + }; // class fftw class Transpose { fftw_plan plan; + fftw_plan plan2; unsigned int a,b; unsigned int nlength,mlength; - unsigned int instride,outstride; + unsigned int ilast,jlast; + unsigned int rows,cols; unsigned int threads; bool inplace; unsigned int size; @@ -631,50 +505,77 @@ public: template<class T> Transpose(unsigned int rows, unsigned int cols, unsigned int length, T *in, T *out=NULL, unsigned int threads=fftw::maxthreads) : - threads(threads) { + rows(rows), cols(cols), threads(threads) { size=sizeof(T); if(size % sizeof(double) != 0) { std::cerr << "ERROR: Transpose is not implemented for type of size " << size; exit(1); } + if(rows == 0 || cols == 0) return; size /= sizeof(double); length *= size; if(!out) out=in; inplace=(out==in); - if(inplace) threads=1; + if(inplace) threads=1; // TODO: Generalize to inplace fftw_iodim dims[3]; a=std::min(rows,threads); b=std::min(cols,threads/a); - unsigned int n=rows/a; - unsigned int m=cols/b; + + unsigned int n=utils::ceilquotient(rows,a); + unsigned int m=utils::ceilquotient(cols,b); + + // If rows <= threads then a=rows and n=1. + // If rows >= threads then b=1 and m=cols. + nlength=n*length; mlength=m*length; - instride=cols; - outstride=rows; - - dims[0].n=n; - dims[0].is=instride*length; + dims[0].n=n; + dims[0].is=cols*length; dims[0].os=length; - - dims[1].n=m; + + dims[1].n=m; dims[1].is=length; - dims[1].os=outstride*length; + dims[1].os=rows*length; dims[2].n=length; dims[2].is=1; dims[2].os=1; - fftw::planThreads(inplace ? threads : 1); - plan=fftw_plan_guru_r2r(0,NULL,3,dims,(double *) in, - (double *) out,NULL,fftw::effort); + fftw::planThreads(1); + + // A plan with rank=0 is a transpose. + plan=fftw_plan_guru_r2r(0,NULL,3,dims,(double *) in,(double *) out, + NULL,fftw::effort); + plan2=NULL; + ilast=a; + jlast=b; + + if(n*a > rows) { // Only happens when rows > threads. + a=utils::ceilquotient(rows,n); + ilast=a-1; + dims[0].n=rows-n*ilast; + plan2=fftw_plan_guru_r2r(0,NULL,3,dims,(double *) in,(double *) out, + NULL,fftw::effort); + } else { // Only happens when rows < threads. + if(m*b > cols) { + b=utils::ceilquotient(cols,m); + jlast=b-1; + dims[1].n=cols-m*jlast; + plan2=fftw_plan_guru_r2r(0,NULL,3,dims,(double *) in,(double *) out, + NULL,fftw::effort); + } + } } - ~Transpose() {if(plan) fftw_destroy_plan(plan);} + ~Transpose() { + if(plan) fftw_destroy_plan(plan); + if(plan2) fftw_destroy_plan(plan2); + } template<class T> void transpose(T *in, T *out=NULL) { @@ -684,24 +585,40 @@ public: exit(1); } #ifndef FFTWPP_SINGLE_THREAD - if(threads == 1) -#endif - fftw_execute_r2r(plan,(double *) in,(double*) out); -#ifndef FFTWPP_SINGLE_THREAD - else { - int A=a, B=b; + if(a > 1) { + if(b > 1) { + int A=a, B=b; #pragma omp parallel for num_threads(A) - for(unsigned int i=0; i < a; ++i) { - unsigned int I=i*nlength; + for(unsigned int i=0; i < a; ++i) { + unsigned int I=i*nlength; +#pragma omp parallel for num_threads(B) + for(unsigned int j=0; j < b; ++j) { + unsigned int J=j*mlength; + fftw_execute_r2r((i < ilast && j < jlast) ? plan : plan2, + (double *) in+cols*I+J, + (double *) out+rows*J+I); + } + } + } else { + int A=a; +#pragma omp parallel for num_threads(A) + for(unsigned int i=0; i < a; ++i) { + unsigned int I=i*nlength; + fftw_execute_r2r(i < ilast ? plan : plan2, + (double *) in+cols*I,(double *) out+I); + } + } + } else if(b > 1) { + int B=b; #pragma omp parallel for num_threads(B) for(unsigned int j=0; j < b; ++j) { unsigned int J=j*mlength; - fftw_execute_r2r(plan,(double *) in+instride*I+J, - (double *) out+outstride*J+I); + fftw_execute_r2r(j < jlast ? plan : plan2, + (double *) in+J,(double *) out+rows*J); } - } - } + } else #endif + fftw_execute_r2r(plan,(double *) in,(double*) out); } }; @@ -830,52 +747,31 @@ public: } }; -// 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: -// -// mfft1d Forward(n,-1,M,stride,dist,in,out); -// Forward.fft(in,out); -// -// In-place usage: -// -// mfft1d Forward(n,-1,M,stride,dist); -// Forward.fft(in); -// -// Notes: -// stride is the spacing between the elements of each Complex vector; -// dist is the spacing between the first elements of the vectors. -// -// -class mfft1d : public fftw, public Threadtable<keytype3,keyless3> { - unsigned int nx; +template<class I, class O> +class fftwblock : public virtual fftw { +public: + int nx; unsigned int M; + size_t istride,ostride; + size_t idist,odist; + fftw_plan plan1,plan2; unsigned int T,Q,R; - size_t stride; - size_t dist; - fftw_plan plan1; - fftw_plan plan2; - static Table threadtable; -public: - mfft1d(unsigned int nx, int sign, unsigned int M=1, size_t stride=1, - size_t dist=0, Complex *in=NULL, Complex *out=NULL, - unsigned int Threads=maxthreads) - : fftw(2*((nx-1)*stride+(M-1)*Dist(nx,stride,dist)+1),sign,Threads,nx), - nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)), - plan1(NULL), plan2(NULL) - { + 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) + : fftw(), nx(nx), M(M), istride(istride), ostride(ostride), + idist(Dist(nx,istride,idist)), odist(Dist(nx,ostride,odist)), + plan1(NULL), plan2(NULL) { T=1; Q=M; R=0; + threaddata S1=Setup(in,out); fftw_plan planT1=plan; - T=std::min(M,Threads); - if(T > 1) { - Q=M/T; + 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); @@ -885,7 +781,7 @@ public: plan2=plan1; } - if(ST.mean > S1.mean-S1.stdev) { + if(ST.mean > S1.mean-S1.stdev) { // Use FFTW's multi-threading fftw_destroy_plan(plan); if(R > 0) { fftw_destroy_plan(plan2); @@ -896,63 +792,121 @@ public: R=0; plan=planT1; threads=S1.threads; - } else { + } else { // Do the multi-threading ourselves fftw_destroy_plan(planT1); threads=ST.threads; } } - } + } - unsigned int Threads() {return std::max(T,threads);} + 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); + } - threaddata lookup(bool inplace, unsigned int threads) { - return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace)); + 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); } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data); + + 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) { - int n=(int) nx; if(R > 0) { - plan2=fftw_plan_many_dft(1,&n,Q+1, - (fftw_complex *) in,NULL,stride,dist, - (fftw_complex *) out,NULL,stride,dist, - sign,effort); - if(!plan2) noplan(); + plan2=Plan(Q+1,(I *) in,(O *) out); + if(!plan2) return NULL; if(threads == 1) plan1=plan2; } - return fftw_plan_many_dft(1,&n,Q, - (fftw_complex *) in,NULL,stride,dist, - (fftw_complex *) out,NULL,stride,dist, - sign,effort); + return Plan(Q,(I *) in,(O *) out); } - - ~mfft1d() { - if(plan2) fftw_destroy_plan(plan2); + + void Execute(fftw_plan plan, fftw_complex *in, fftw_complex *out) { + fftw_execute_dft(plan,in,out); } - void Execute(Complex *in, Complex *out, bool=false) { - if(T == 1) { - fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out); - } else { - unsigned int Tdist=T*dist; - unsigned int extra=(T-R)*dist; + void Execute(fftw_plan plan, double *in, fftw_complex *out) { + fftw_execute_dft_r2c(plan,in,out); + } + + void Execute(fftw_plan plan, fftw_complex *in, double *out) { + fftw_execute_dft_c2r(plan,in,out); + } + void Execute(Complex *in, Complex *out, bool=false) { + if(T == 1) + Execute(plan,(I *) in,(O *) out); + else { + unsigned int extra=T-R; #ifndef FFTWPP_SINGLE_THREAD #pragma omp parallel for num_threads(T) #endif - for(unsigned int i=0; i < Tdist; i += dist) { - bool normal=i < extra; - unsigned int offset=normal ? Q*i : Q*i+i-extra; - fftw_execute_dft(normal ? plan : plan2,(fftw_complex *) in+offset, - (fftw_complex *) out+offset); + for(unsigned int i=0; i < T; ++i) { + unsigned int iQ=i*Q; + if(i < extra) + Execute(plan,(I *) in+iQ*idist,(O *) out+iQ*odist); + else { + unsigned int offset=iQ+i-extra; + Execute(plan2,(I *) in+offset*idist,(O *) out+offset*odist); + } } } } - void fftNormalized(Complex *in, Complex *out=NULL) { - fftw::fftNormalized(in,out,nx,M,stride,dist); + 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: +// +// mfft1d Forward(n,-1,M,stride,dist,in,out); +// Forward.fft(in,out); +// +// In-place usage: +// +// mfft1d Forward(n,-1,M,stride,dist); +// Forward.fft(in); +// +// Notes: +// stride is the spacing between the elements of each Complex vector; +// 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; +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) {} + + 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): + 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) {} + + 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); } }; @@ -984,18 +938,7 @@ public: rcfft1d(unsigned int nx, double *in, Complex *out=NULL, unsigned int threads=maxthreads) - : fftw(realsize(nx,in,out),-1,threads,nx), nx(nx) {Setup(in,out);} - -#ifdef __Array_h__ - rcfft1d(unsigned int nx, const Array::array1<Complex>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),-1,threads,nx), nx(nx) {Setup(out,(double*) NULL);} - - rcfft1d(unsigned int nx, const Array::array1<double>& in, - const Array::array1<Complex>& out=Array::NULL1, - unsigned int threads=maxthreads) - : fftw(realsize(nx,in(),out()),-1,threads,nx), nx(nx) {Setup(in,out);} -#endif + : 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)); @@ -1045,20 +988,6 @@ public: unsigned int threads=maxthreads) : fftw(realsize(nx,in,out),1,threads,nx), nx(nx) {Setup(in,out);} -#ifdef __Array_h__ - crfft1d(unsigned int nx, const Array::array1<double>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),1,threads,nx), nx(nx) {Setup(out);} - - crfft1d(unsigned int nx, const Array::array1<Complex>& in, - unsigned int threads=maxthreads) - : fftw(2*in.Size(),1,threads,nx), nx(nx) {Setup(in);} - - crfft1d(unsigned int nx, const Array::array1<Complex>& in, - const Array::array1<double>& out) - : fftw(out.Size(),1,threads,nx), nx(nx) {Setup(in,out);} -#endif - threaddata lookup(bool inplace, unsigned int threads) { return Lookup(threadtable,keytype1(nx,threads,inplace)); } @@ -1074,7 +1003,7 @@ public: fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) 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 @@ -1083,68 +1012,53 @@ public: // // Out-of-place usage: // -// mrcfft1d Forward(n,M,stride,dist,in,out); +// mrcfft1d Forward(n,M,istride,ostride,idist,odist,in,out); // Forward.fft(in,out); // // In-place usage: // -// mrcfft1d Forward(n,M,stride,dist); +// mrcfft1d Forward(n,M,istride,ostride,idist,odist); // Forward.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; +// istride is the spacing between the elements of each real vector; +// ostride is the spacing between the elements of each Complex vector; +// idist is the spacing between the first elements of the real vectors; +// odist is the spacing between the first elements of the Complex vectors; // in contains the n real values stored as a Complex array; // out contains the first n/2+1 Complex Fourier values. // -class mrcfft1d : public fftw, public Threadtable<keytype2,keyless2> { - unsigned int nx; - unsigned int M; - size_t stride; - size_t dist; +class mrcfft1d : public fftwblock<double,fftw_complex>, + public Threadtable<keytype3,keyless3> { static Table threadtable; -public: - mrcfft1d(unsigned int nx, unsigned int M=1, size_t stride=1, - size_t dist=0, Complex *out=NULL, - unsigned int threads=maxthreads) - : fftw(2*(nx/2*stride+(M-1)*Dist(nx,stride,dist)+1),-1,threads,nx), nx(nx), - M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(out);} - - mrcfft1d(unsigned int nx, unsigned int M=1, size_t stride=1, - size_t dist=0, double *in=NULL, Complex *out=NULL, +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(2*(nx/2*stride+(M-1)*Dist(nx,stride,dist)+1),-1,threads,nx), nx(nx), - M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);} + : 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,keytype2(nx,M,threads,inplace)); - } - void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype2(nx,M,data.threads,inplace),data); + return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace)); } - fftw_plan Plan(Complex *in, Complex *out) { - int n=(int) nx; - if(in != out) - return fftw_plan_many_dft_r2c(1,&n,M, - (double *) in,NULL,stride,dist, - (fftw_complex *) out,NULL,stride,dist/2+1, - effort); - return fftw_plan_many_dft_r2c(1,&n,M, - (double *) in,NULL,stride,dist, - (fftw_complex *) out,NULL,stride,dist/2, - effort); + void store(bool inplace, const threaddata& data) { + Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data); } - void Execute(Complex *in, Complex *out, bool=false) { - fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out); + void fftNormalized(double *in, Complex *out=NULL) { + fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out,false); } - void fftNormalized(Complex *in, Complex *out=NULL) { - fftw::fftNormalized(in,out,nx/2+1,M,stride,dist); + void fft0Normalized(double *in, Complex *out=NULL) { + fftw::fftNormalized<double,Complex>(nx/2+1,M,ostride,odist,in,out,true); } }; - + // 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 @@ -1154,12 +1068,12 @@ public: // // Out-of-place usage (input destroyed): // -// mcrfft1d Backward(n,M,stride,dist,in,out); +// mcrfft1d Backward(n,M,istride,ostride,idist,odist,in,out); // Backward.fft(in,out); // // In-place usage: // -// mcrfft1d Backward(n,M,stride,dist); +// mcrfft1d Backward(n,M,istride,ostride,idist,odist); // Backward.fft(out); // // Notes: @@ -1168,61 +1082,29 @@ 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 fftw, public Threadtable<keytype2,keyless2> { - unsigned int nx; - unsigned int M; - size_t stride; - size_t dist; +class mcrfft1d : public fftwblock<fftw_complex,double>, + public Threadtable<keytype3,keyless3> { static Table threadtable; public: - mcrfft1d(unsigned int nx, unsigned int M=1, size_t stride=1, - size_t dist=0, Complex *in=NULL, double *out=NULL, + 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) - : fftw((realsize(nx,in,out)-2)*stride+2*(M-1)*Dist(nx,stride,dist)+2,1, - threads,nx), - nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);} + : 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,keytype2(nx,M,threads,inplace)); + return Lookup(threadtable,keytype3(nx,Q,R,threads,inplace)); } void store(bool inplace, const threaddata& data) { - Store(threadtable,keytype2(nx,M,data.threads,inplace),data); - } - - fftw_plan Plan(Complex *in, Complex *out) { - int n=(int) nx; - if(in != out) - return fftw_plan_many_dft_c2r(1,&n,M, - (fftw_complex *) in,NULL,stride,dist, - (double *) out,NULL,stride,2*(dist-1), - effort); - return fftw_plan_many_dft_c2r(1,&n,M, - (fftw_complex *) in,NULL,stride,dist, - (double *) out,NULL,stride,2*dist, - effort); - } - - void Execute(Complex *in, Complex *out, bool=false) { - fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out); + Store(threadtable,keytype3(nx,Q,R,data.threads,inplace),data); } - void fftNormalized(Complex *in, double *out=NULL) { - // TODO: improve threading (cf: mfft1d). - fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out); - - unsigned int stop=nx*stride; - unsigned int odist=in == (Complex *) out ? 2*dist : 2*(dist-1); - double *outMdist=out+M*odist; - -#ifndef FFTWPP_SINGLE_THREAD -#pragma omp parallel for num_threads(threads) -#endif - for(unsigned int i=0; i < stop; i += stride) { - double *pstop=outMdist+i; - for(double *p=out+i; p < pstop; p += odist) { - *p *= norm; - } - } + 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); } }; @@ -1265,7 +1147,9 @@ public: fft2d(int sign, const Array::array2<Complex>& in, 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);} + : fftw(2*in.Size(),sign,threads), nx(in.Nx()), ny(in.Ny()) { + Setup(in,out); + } #endif threaddata lookup(bool inplace, unsigned int threads) { @@ -1318,24 +1202,10 @@ public: rcfft2d(unsigned int nx, unsigned int ny, double *in, Complex *out=NULL, unsigned int threads=maxthreads) - : fftw(nx*realsize(ny,in,out),-1,threads,nx*ny), nx(nx), ny(ny) { + : fftw(2*nx*(ny/2+1),-1,threads,nx*ny), nx(nx), ny(ny) { Setup(in,out); } -#ifdef __Array_h__ - rcfft2d(unsigned int ny, const Array::array2<Complex>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),-1,threads,out.Nx()*ny), nx(out.Nx()), ny(ny) { - Setup(out); - } - - rcfft2d(unsigned int ny, const Array::array2<double>& in, - const Array::array2<Complex>& out=Array::NULL2, - unsigned int threads=maxthreads) - : fftw(in.Nx()*realsize(ny,in(),out()),-1,threads,in.Nx()*ny), - nx(in.Nx()), ny(ny) {Setup(in,out);} -#endif - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_r2c_2d(nx,ny,(double *) in,(fftw_complex *) out, effort); @@ -1348,6 +1218,23 @@ 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; + if(nx % 2 == 0) +#ifndef FFTWPP_SINGLE_THREAD +#pragma omp parallel for num_threads(threads) +#endif + for(unsigned int 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) + f[(i+1)*nyp-1]=0.0; + } }; // Compute the real two-dimensional inverse Fourier transform of the @@ -1378,31 +1265,16 @@ class crfft2d : public fftw { unsigned int nx; unsigned int ny; public: - crfft2d(unsigned int nx, unsigned int ny, Complex *in=NULL, + 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(in);} + 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, + 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); } -#ifdef __Array_h__ - crfft2d(unsigned int ny, const Array::array2<double>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),1,threads,out.Nx()*ny), nx(out.Nx()), ny(ny) {Setup(out);} - - crfft2d(unsigned int ny, const Array::array2<Complex>& in, - unsigned int threads=maxthreads) - : fftw(2*in.Size(),1,threads,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in);} - - crfft2d(unsigned int ny, const Array::array2<Complex>& in, - const Array::array2<double>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),1,threads,in.Nx()*ny), nx(in.Nx()), ny(ny) {Setup(in,out);} -#endif - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_c2r_2d(nx,ny,(fftw_complex *) in,(double *) out, effort); @@ -1415,6 +1287,23 @@ 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; + if(nx % 2 == 0) +#ifndef FFTWPP_SINGLE_THREAD +#pragma omp parallel for num_threads(threads) +#endif + for(unsigned int 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) + f[(i+1)*nyp-1]=0.0; + } }; // Compute the complex three-dimensional Fourier transform of @@ -1495,27 +1384,17 @@ class rcfft3d : public fftw { unsigned int ny; unsigned int nz; public: - rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *out=NULL) - : fftw(2*nx*ny*(nz/2+1),-1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz) {Setup(out);} + 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) - : fftw(nx*ny*realsize(nz,in,out),-1,threads,nx*ny*nz), + : fftw(2*nx*ny*(nz/2+1),-1,threads,nx*ny*nz), nx(nx), ny(ny), nz(nz) {Setup(in,out);} -#ifdef __Array_h__ - rcfft3d(unsigned int nz, const Array::array3<Complex>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),-1,threads,out.Nx()*out.Ny()*nz), - nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);} - - rcfft3d(unsigned int nz, const Array::array3<double>& in, - const Array::array3<Complex>& out=Array::NULL3, - unsigned int threads=maxthreads) - : fftw(in.Nx()*in.Ny()*realsize(nz,in(),out()),-1,threads,in.Size()), - nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);} -#endif - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_r2c_3d(nx,ny,nz,(double *) in,(fftw_complex *) out, effort); @@ -1528,6 +1407,38 @@ 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; + unsigned int 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) + f[k]=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) { + unsigned int iyz=i*yz; + for(unsigned int k=0; k < nzp; ++k) + f[iyz+k]=0.0; + } + } + + if(nz % 2 == 0) +#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) + f[i*yz+(j+1)*nzp-1]=0.0; + } }; // Compute the real two-dimensional inverse Fourier transform of the @@ -1559,34 +1470,16 @@ class crfft3d : public fftw { unsigned int ny; unsigned int nz; public: - crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in=NULL, + crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, double *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);} + {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);} -#ifdef __Array_h__ - crfft3d(unsigned int nz, const Array::array3<double>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),1,threads,out.Nx()*out.Ny()*nz), - nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);} - - crfft3d(unsigned int nz, const Array::array3<Complex>& in, - unsigned int threads=maxthreads) - : fftw(2*in.Size(),1,threads,in.Nx()*in.Ny()*nz), - nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in);} - - crfft3d(unsigned int nz, const Array::array3<Complex>& in, - const Array::array3<double>& out, - unsigned int threads=maxthreads) - : fftw(out.Size(),1,threads,in.Nx()*in.Ny()*nz), - nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in,out);} -#endif - fftw_plan Plan(Complex *in, Complex *out) { return fftw_plan_dft_c2r_3d(nx,ny,nz,(fftw_complex *) in,(double *) out, effort); @@ -1599,8 +1492,40 @@ 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; + unsigned int 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) + f[k]=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) { + unsigned int iyz=i*yz; + for(unsigned int k=0; k < nzp; ++k) + f[iyz+k]=0.0; + } + } + + if(nz % 2 == 0) +#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) + f[i*yz+(j+1)*nzp-1]=0.0; + } +}; + } #endif |