11#include " fft.h"
2- #include < iostream>
3-
42namespace ModulePW
53{
64
@@ -34,7 +32,7 @@ FFT::~FFT()
3432#endif
3533}
3634
37- void FFT:: initfft(int nx_in, int bigny_in, int nz_in, int lix_in , int rix_in , int ns_in, int nplane_in, int nproc_in, bool gamma_only_in, bool mpifft_in)
35+ void FFT:: initfft(int nx_in, int bigny_in, int nz_in, int liy_in , int riy_in , int ns_in, int nplane_in, int nproc_in, bool gamma_only_in, bool mpifft_in)
3836{
3937 this ->gamma_only = gamma_only_in;
4038 this ->nx = nx_in;
@@ -43,8 +41,8 @@ void FFT:: initfft(int nx_in, int bigny_in, int nz_in, int lix_in, int rix_in, i
4341 else this ->ny = this ->bigny ;
4442 this ->nz = nz_in;
4543 this ->ns = ns_in;
46- this ->lix = lix_in ;
47- this ->rix = rix_in ;
44+ this ->liy = liy_in ;
45+ this ->riy = riy_in ;
4846 this ->nplane = nplane_in;
4947 this ->nproc = nproc_in;
5048 this ->mpifft = mpifft_in;
@@ -136,11 +134,11 @@ void FFT :: initplan()
136134 // (fftw_complex*) aux2, rankc, this->nplane, 1,
137135 // r_rspace, rankd, this->nplane, 1, FFTW_MEASURE);
138136
139- int npy = this ->nplane * this ->ny ;
137+ // int npy = this->nplane * this->ny;
140138 int bignpy = this ->nplane * this ->bigny ;
141- this ->planxfor = fftw_plan_many_dft ( 1 , &this ->nx , npy , (fftw_complex *)aux2, embed, bignpy, 1 ,
139+ this ->planxfor = fftw_plan_many_dft ( 1 , &this ->nx , this -> nplane , (fftw_complex *)aux2, embed, bignpy, 1 ,
142140 (fftw_complex *)aux2, embed, bignpy, 1 , FFTW_FORWARD, FFTW_MEASURE );
143- this ->planxbac = fftw_plan_many_dft ( 1 , &this ->nx , npy , (fftw_complex *)aux2, embed, bignpy, 1 ,
141+ this ->planxbac = fftw_plan_many_dft ( 1 , &this ->nx , this -> nplane , (fftw_complex *)aux2, embed, bignpy, 1 ,
144142 (fftw_complex *)aux2, embed, bignpy, 1 , FFTW_BACKWARD, FFTW_MEASURE );
145143 this ->planyr2c = fftw_plan_many_dft_r2c ( 1 , &this ->bigny , this ->nplane , r_rspace , embed, this ->nplane , 1 ,
146144 (fftw_complex*)aux1, embed, this ->nplane , 1 , FFTW_MEASURE );
@@ -169,14 +167,14 @@ void FFT :: initplan()
169167 // (fftw_complex*) aux2, embed, this->nplane, 1, FFTW_BACKWARD, FFTW_MEASURE);
170168
171169 // It is better to use in-place for stride > 1
172- int npy = this ->nplane * this ->ny ;
173- this ->planxfor = fftw_plan_many_dft ( 1 , &this ->nx , npy , (fftw_complex *)aux2, embed, npy , 1 ,
174- (fftw_complex *)aux2, embed, npy , 1 , FFTW_FORWARD, FFTW_MEASURE );
175- this ->planxbac = fftw_plan_many_dft ( 1 , &this ->nx , npy , (fftw_complex *)aux2, embed, npy , 1 ,
176- (fftw_complex *)aux2, embed, npy , 1 , FFTW_BACKWARD, FFTW_MEASURE );
177- this ->planyfor = fftw_plan_many_dft ( 1 , &this ->ny , this ->nplane , (fftw_complex*)aux2 , embed, this ->nplane , 1 ,
170+ int bignpy = this ->nplane * this ->bigny ;
171+ this ->planxfor = fftw_plan_many_dft ( 1 , &this ->nx , this -> nplane , (fftw_complex *)aux2, embed, bignpy , 1 ,
172+ (fftw_complex *)aux2, embed, bignpy , 1 , FFTW_FORWARD, FFTW_MEASURE );
173+ this ->planxbac = fftw_plan_many_dft ( 1 , &this ->nx , this -> nplane , (fftw_complex *)aux2, embed, bignpy , 1 ,
174+ (fftw_complex *)aux2, embed, bignpy , 1 , FFTW_BACKWARD, FFTW_MEASURE );
175+ this ->planyfor = fftw_plan_many_dft ( 1 , &this ->bigny , this ->nplane , (fftw_complex*)aux2 , embed, this ->nplane , 1 ,
178176 (fftw_complex*)aux2, embed, this ->nplane , 1 , FFTW_FORWARD, FFTW_MEASURE );
179- this ->planybac = fftw_plan_many_dft ( 1 , &this ->ny , this ->nplane , (fftw_complex*)aux2 , embed, this ->nplane , 1 ,
177+ this ->planybac = fftw_plan_many_dft ( 1 , &this ->bigny , this ->nplane , (fftw_complex*)aux2 , embed, this ->nplane , 1 ,
180178 (fftw_complex*)aux2, embed, this ->nplane , 1 , FFTW_BACKWARD, FFTW_MEASURE );
181179 }
182180
@@ -308,61 +306,80 @@ void FFT::fftzbac(std::complex<double>* & in, std::complex<double>* & out)
308306
309307void FFT::fftxyfor (std::complex <double >* & in, std::complex <double >* & out)
310308{
311- int npy = this ->nplane * this -> ny;
312- fftw_execute_dft ( this ->planxfor , (fftw_complex *)in, (fftw_complex *)out);
313- for (int i=0 ; i<=this ->lix ;++i)
309+ int bignpy = this ->nplane * this -> bigny;
310+ for (int i=0 ; i<this ->nx ;++i)
314311 {
315- fftw_execute_dft ( this ->planyfor , (fftw_complex*)&in[i*npy ], (fftw_complex*)&out[i*npy] );
312+ fftw_execute_dft ( this ->planyfor , (fftw_complex *)&in[i*bignpy ], (fftw_complex *)&out[i*bignpy] );
316313 }
317- for (int i=this ->rix ; i<this ->nx ;++i)
314+
315+ for (int i=0 ; i<=this ->liy ;++i)
318316 {
319- fftw_execute_dft ( this ->planyfor , (fftw_complex*)&in[i*npy], (fftw_complex*)&out[i*npy] );
317+ fftw_execute_dft ( this ->planxfor , (fftw_complex *)&in[i*nplane], (fftw_complex *)&out[i*nplane]);
318+ }
319+ for (int i=this ->riy ; i<this ->ny ;++i)
320+ {
321+ fftw_execute_dft ( this ->planxfor , (fftw_complex *)&in[i*nplane], (fftw_complex *)&out[i*nplane]);
320322 }
321323 return ;
322324}
323325
324326void FFT::fftxybac (std::complex <double >* & in, std::complex <double >* & out)
325327{
326- int npy = this ->nplane * this -> ny;
327-
328- for (int i=0 ; i<=this ->lix ;++i)
328+ int bignpy = this ->nplane * this -> bigny;
329+ // x-direction
330+ for (int i=0 ; i<=this ->liy ;++i)
331+ {
332+ fftw_execute_dft ( this ->planxbac , (fftw_complex *)&in[i*nplane], (fftw_complex *)&out[i*nplane]);
333+ }
334+ for (int i=this ->riy ; i<this ->ny ;++i)
329335 {
330- fftw_execute_dft ( this ->planybac , (fftw_complex*)&in[i*npy ], (fftw_complex*)&out[i*npy] );
336+ fftw_execute_dft ( this ->planxbac , (fftw_complex *)&in[i*nplane ], (fftw_complex *)&out[i*nplane] );
331337 }
332- for (int i=this ->rix ; i<this ->nx ;++i)
338+
339+ // //y-direction
340+ for (int i=0 ; i<this ->nx ;++i)
333341 {
334- fftw_execute_dft ( this ->planybac , (fftw_complex*)&in[i*npy ], (fftw_complex*)&out[i*npy ] );
342+ fftw_execute_dft ( this ->planybac , (fftw_complex*)&in[i*bignpy ], (fftw_complex*)&out[i*bignpy ] );
335343 }
336- fftw_execute_dft ( this ->planxbac , (fftw_complex *)in, (fftw_complex *)out);
337344 return ;
338345}
339346
340347void FFT::fftxyr2c (double * &in, std::complex <double >* & out)
341348{
342- // int npy = this->nplane * this-> ny;
343349 int bignpy = this ->nplane * this -> bigny;
344- // int padnpy = this->nplane * this-> ny * 2;
350+
345351 for (int i=0 ; i<this ->nx ;++i)
346352 {
347353 fftw_execute_dft_r2c ( this ->planyr2c , &in[i*bignpy*2 ], (fftw_complex*)&out[i*bignpy] );
348- if ((double *)&out[i*bignpy] != &in[i*bignpy*2 ] ) std::cout<<i<<" wrond\n " ;
349- // fftw_execute_dft_r2c( this->planyfor, &r_rspace[4*i*padnpy], (fftw_complex*)&aux2[i*padnpy] );
350354 }
351- fftw_execute_dft ( this ->planxfor , (fftw_complex *)out, (fftw_complex *)out);
355+
356+ for (int i=0 ; i<=this ->liy ;++i)
357+ {
358+ fftw_execute_dft ( this ->planxfor , (fftw_complex *)&out[i*nplane], (fftw_complex *)&out[i*nplane]);
359+ }
360+ for (int i=this ->riy ; i<this ->ny ;++i)
361+ {
362+ fftw_execute_dft ( this ->planxfor , (fftw_complex *)&out[i*nplane], (fftw_complex *)&out[i*nplane]);
363+ }
352364 return ;
353365}
354366
355367
356368void FFT::fftxyc2r (std::complex <double >* & in, double * & out)
357369{
358- // int npy = this->nplane * this-> ny;
359370 int bignpy = this ->nplane * this -> bigny;
360- // int padnpy = this->nplane * this-> ny * 2;
361- fftw_execute_dft ( this ->planxbac , (fftw_complex *)in, (fftw_complex *)in);
371+ for (int i=0 ; i<=this ->liy ;++i)
372+ {
373+ fftw_execute_dft ( this ->planxbac , (fftw_complex *)&in[i*nplane], (fftw_complex *)&in[i*nplane]);
374+ }
375+ for (int i=this ->riy ; i<this ->ny ;++i)
376+ {
377+ fftw_execute_dft ( this ->planxbac , (fftw_complex *)&in[i*nplane], (fftw_complex *)&in[i*nplane]);
378+ }
379+
362380 for (int i=0 ; i<this ->nx ;++i)
363381 {
364382 fftw_execute_dft_c2r ( this ->planyc2r , (fftw_complex*)&in[i*bignpy], &out[i*bignpy*2 ] );
365- // fftw_execute_dft_c2r( this->planybac, (fftw_complex*)&aux2[i*padnpy], &r_rspace[4*i*padnpy] );
366383 }
367384 return ;
368385}
0 commit comments