Commit 566b2ae0 by Francois Gygi

restored ESSL FFT functionality

git-svn-id: http://qboxcode.org/svn/qb/branches/fftw3@1531 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 19e75157
......@@ -585,9 +585,7 @@ void FourierTransform::bwd(complex<double>* val)
dcft_(&initflag,&zvec_[0],&inc1,&inc2,&zvec_[0],&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
#endif
#if USE_FFTW2
#elif USE_FFTW2
/*
* void fftw(fftw_plan plan, int howmany,
* FFTW_COMPLEX *in, int istride, int idist,
......@@ -804,7 +802,7 @@ void FourierTransform::bwd(complex<double>* val)
}
#endif // USE_FFTW3_2D
#elif USE_ESSL_FFT // USE_FFTW3
#elif USE_ESSL_FFT
for ( int k = 0; k < np2_loc_[myproc_]; k++ )
{
// transform along x for non-zero vectors only
......@@ -981,7 +979,7 @@ void FourierTransform::fwd(complex<double>* val)
for ( int k = 0; k < np2_loc_[myproc_]; k++ )
fftw_execute_dft ( fwplan2d, (fftw_complex*)&val[k*np0_*np1_],
(fftw_complex*)&val[k*np0_*np1_] );
#endif
#endif // USE_FFTW3_THREADS
#else // USE_FFTW3_2D
for ( int k = 0; k < np2_loc_[myproc_]; k++ )
{
......@@ -1032,7 +1030,7 @@ void FourierTransform::fwd(complex<double>* val)
tm_f_x.stop();
#endif
}
#endif // use FFTW3_2D
#endif
#elif USE_ESSL_FFT
for ( int k = 0; k < np2_loc_[myproc_]; k++ )
{
......@@ -1122,7 +1120,7 @@ void FourierTransform::fwd(complex<double>* val)
(FFTW_COMPLEX*)0);
}
}
#else
#else // _OPENMP
int inc1, inc2, istart;
// transform along y for all values of x
......@@ -1146,8 +1144,8 @@ void FourierTransform::fwd(complex<double>* val)
istart = np0_ * ( (np1_-ntrans) + k * np1_ );
fftw(fwplan0,ntrans,(FFTW_COMPLEX*)&val[istart],inc1,inc2,
(FFTW_COMPLEX*)0,0,0);
#endif // _OPENMP
} // k
#endif
#else
// No library
for ( int k = 0; k < np2_loc_[myproc_]; k++ )
......
......@@ -19,31 +19,55 @@
#include "sinft.h"
#include <math.h>
#include <assert.h>
#include <vector>
#include <complex>
using namespace std;
#if USE_FFTW2
#if USE_DFFTW
#include "dfftw.h"
#else
#include "fftw.h"
#endif
#endif
#if USE_FFTW3
#elif USE_FFTW3
#include "fftw3.h"
#elif USE_ESSL_FFT
extern "C" {
void dcft_(int *initflag, std::complex<double> *x, int *inc2x, int *inc3x,
std::complex<double> *y, int *inc2y, int *inc3y,
int *length, int *ntrans, int *isign,
double *scale, double *aux1, int *naux1,
double *aux2, int *naux2);
}
#else
// no FFT library
void cfft ( int idir, complex<double> *z1, complex<double> *z2, int n,
int *inzee );
void fftstp ( int idir, complex<double> *zin, int after,
int now, int before, complex<double> *zout );
#endif
#include <vector>
#include <complex>
using namespace std;
void sinft(int n, double *y)
{
vector<complex<double> > zin(2*n), zout(2*n);
#if defined(USE_FFTW2) || defined(USE_FFTW3)
fftw_plan fwplan;
#endif
#if USE_FFTW2
fwplan = fftw_create_plan(2*n,FFTW_FORWARD,FFTW_ESTIMATE);
#elif USE_FFTW3
fwplan = fftw_plan_dft_1d(2*n,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0],
FFTW_FORWARD, FFTW_ESTIMATE);
#elif USE_ESSL_FFT
int np = 2 * n;
int naux1 = (int) (30000 + 2.28 * np);
std::vector<double> aux1(naux1);
int ntrans = 1;
int naux2 = (int) (20000 + 2.28 * np + (256 + 2*np)*min(64,ntrans));
std::vector<double> aux2(naux2);
#else
#error
// no FFT library
// no initialization needed
#endif
zin[0] = 0.0;
......@@ -57,28 +81,56 @@ void sinft(int n, double *y)
fftw_one(fwplan,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0]);
#elif USE_FFTW3
fftw_execute(fwplan);
#elif USE_ESSL_FFT
// initialize forward transform
int initflag = 1;
int inc1 = 1, inc2 = np;
int isign = 1;
double scale = 1.0;
complex<double> *p = 0;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np,&ntrans,
&isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);
// call transform
initflag = 0;
dcft_(&initflag,&zin[0],&inc1,&inc2,&zout[0],&inc1,&inc2,&np,&ntrans,
&isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);
#else
#error
// no FFT library
int idir = 1, inzee = 1, np = 2*n;
cfft ( idir, &zin[0],&zout[0],np, &inzee );
#endif
for ( int i = 0; i < n; i++ )
{
y[i] = -0.5 * imag(zout[i]);
}
#if defined(USE_FFTW2) || defined(USE_FFTW3)
fftw_destroy_plan(fwplan);
#endif
}
void cosft1(int n, double *y)
{
/* Note: the array y contains n+1 elements */
vector<complex<double> > zin(2*n), zout(2*n);
#if defined(USE_FFTW2) || defined(USE_FFTW3)
fftw_plan fwplan;
#endif
#if USE_FFTW2
fwplan = fftw_create_plan(2*n,FFTW_FORWARD,FFTW_ESTIMATE);
#elif USE_FFTW3
fwplan = fftw_plan_dft_1d(2*n,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0],
FFTW_FORWARD, FFTW_ESTIMATE);
#elif USE_ESSL_FFT
int np = 2 * n;
int naux1 = (int) (30000 + 2.28 * np);
std::vector<double> aux1(naux1);
int ntrans = 1;
int naux2 = (int) (20000 + 2.28 * np + (256 + 2*np)*min(64,ntrans));
std::vector<double> aux2(naux2);
#else
#error
// no FFT library
// no initialization needed
#endif
zin[0] = y[0];
......@@ -92,13 +144,31 @@ void cosft1(int n, double *y)
fftw_one(fwplan,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0]);
#elif USE_FFTW3
fftw_execute(fwplan);
#elif USE_ESSL_FFT
// initialize forward transform
int initflag = 1;
int inc1 = 1, inc2 = np;
int isign = 1;
double scale = 1.0;
complex<double> *p = 0;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np,&ntrans,
&isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);
// call transform
initflag = 0;
dcft_(&initflag,&zin[0],&inc1,&inc2,&zout[0],&inc1,&inc2,&np,&ntrans,
&isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);
#else
#error
// no FFT library
int idir = 1, inzee = 1, np = 2*n;
cfft ( idir, &zin[0],&zout[0],np, &inzee );
#endif
y[0] = 0.5 * real(zout[0]);
for ( int i = 1; i < n; i++ )
{
y[i] = 0.5 * real(zout[i]);
}
#if defined(USE_FFTW2) || defined(USE_FFTW3)
fftw_destroy_plan(fwplan);
#endif
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment