@@ -87,69 +87,69 @@ void FFT_CPU<double>::setupFFT()
8787 default :
8888 break ;
8989 }
90- z_auxg = (std::complex <double >*)fftw_malloc (sizeof (fftw_complex) * this ->maxgrids );
91- z_auxr = (std::complex <double >*)fftw_malloc (sizeof (fftw_complex) * this ->maxgrids );
92- d_rspace = (double *)z_auxg;
93- this ->planzfor = fftw_plan_many_dft (1 , &this ->nz , this ->ns , (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz ,
94- (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz , FFTW_FORWARD, flag);
90+ z_auxg = (std::complex <double >*)fftw_malloc (sizeof (fftw_complex) * this ->maxgrids );
91+ z_auxr = (std::complex <double >*)fftw_malloc (sizeof (fftw_complex) * this ->maxgrids );
92+ d_rspace = (double *)z_auxg;
93+ this ->planzfor = fftw_plan_many_dft (1 , &this ->nz , this ->ns , (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz ,
94+ (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz , FFTW_FORWARD, flag);
9595
96- this ->planzbac = fftw_plan_many_dft (1 , &this ->nz , this ->ns , (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz ,
97- (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz , FFTW_BACKWARD, flag);
96+ this ->planzbac = fftw_plan_many_dft (1 , &this ->nz , this ->ns , (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz ,
97+ (fftw_complex*)z_auxg, &this ->nz , 1 , this ->nz , FFTW_BACKWARD, flag);
9898
99- // ---------------------------------------------------------
100- // 2 D - XY
101- // ---------------------------------------------------------
102- // 1D+1D is much faster than 2D FFT!
103- // in-place fft is better for c2c and out-of-place fft is better for c2r
104- int * embed = nullptr ;
105- int npy = this ->nplane * this ->ny ;
106- if (this ->xprime )
99+ // ---------------------------------------------------------
100+ // 2 D - XY
101+ // ---------------------------------------------------------
102+ // 1D+1D is much faster than 2D FFT!
103+ // in-place fft is better for c2c and out-of-place fft is better for c2r
104+ int * embed = nullptr ;
105+ int npy = this ->nplane * this ->ny ;
106+ if (this ->xprime )
107+ {
108+ this ->planyfor = fftw_plan_many_dft (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed,this ->nplane , 1 ,
109+ (fftw_complex*)z_auxr, embed,this ->nplane , 1 , FFTW_FORWARD, flag);
110+ this ->planybac = fftw_plan_many_dft (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed,this ->nplane , 1 ,
111+ (fftw_complex*)z_auxr, embed,this ->nplane , 1 , FFTW_BACKWARD, flag);
112+ if (this ->gamma_only )
113+ {
114+ this ->planxr2c = fftw_plan_many_dft_r2c (1 , &this ->nx , npy, d_rspace, embed, npy, 1 , (fftw_complex*)z_auxr,
115+ embed, npy, 1 , flag);
116+ this ->planxc2r = fftw_plan_many_dft_c2r (1 , &this ->nx , npy, (fftw_complex*)z_auxr, embed, npy, 1 , d_rspace,
117+ embed, npy, 1 , flag);
118+ }
119+ else
107120 {
108- this ->planyfor = fftw_plan_many_dft (1 , &this ->ny , this -> nplane , (fftw_complex*)z_auxr, embed,this -> nplane , 1 ,
109- (fftw_complex*)z_auxr, embed,this -> nplane , 1 , FFTW_FORWARD, flag);
110- this ->planybac = fftw_plan_many_dft (1 , &this ->ny , this -> nplane , (fftw_complex*)z_auxr, embed,this -> nplane , 1 ,
111- (fftw_complex*)z_auxr, embed,this -> nplane , 1 , FFTW_BACKWARD, flag);
112- if ( this -> gamma_only )
113- {
114- this -> planxr2c = fftw_plan_many_dft_r2c ( 1 , & this -> nx , npy, d_rspace, embed, npy, 1 , (fftw_complex*)z_auxr,
115- embed, npy, 1 , flag);
116- this ->planxc2r = fftw_plan_many_dft_c2r (1 , &this ->nx , npy , (fftw_complex*)z_auxr, embed, npy, 1 , d_rspace ,
117- embed, npy, 1 , flag);
118- }
119- else
120- {
121- this -> planxfor1 = fftw_plan_many_dft ( 1 , & this -> nx , npy, (fftw_complex*)z_auxr, embed, npy, 1 ,
122- (fftw_complex*)z_auxr , embed, npy , 1 , FFTW_FORWARD, flag);
123- this -> planxbac1 = fftw_plan_many_dft ( 1 , & this -> nx , npy, (fftw_complex*)z_auxr, embed, npy , 1 ,
124- (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_BACKWARD, flag);
125- }
121+ this ->planxfor1 = fftw_plan_many_dft (1 , &this ->nx , npy , (fftw_complex*)z_auxr, embed, npy , 1 ,
122+ (fftw_complex*)z_auxr, embed, npy , 1 , FFTW_FORWARD, flag);
123+ this ->planxbac1 = fftw_plan_many_dft (1 , &this ->nx , npy , (fftw_complex*)z_auxr, embed, npy , 1 ,
124+ (fftw_complex*)z_auxr, embed, npy , 1 , FFTW_BACKWARD, flag);
125+ }
126+ }
127+ else
128+ {
129+ this ->planxfor1 = fftw_plan_many_dft (1 , &this ->nx , this -> nplane * ( this -> lixy + 1 ) , (fftw_complex*)z_auxr, embed, npy,
130+ 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_FORWARD , flag);
131+ this -> planxbac1 = fftw_plan_many_dft ( 1 , & this -> nx , this -> nplane * ( this -> lixy + 1 ), (fftw_complex*)z_auxr, embed, npy,
132+ 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_BACKWARD, flag);
133+ if ( this -> gamma_only )
134+ {
135+ this -> planyr2c = fftw_plan_many_dft_r2c ( 1 , & this -> ny , this -> nplane , d_rspace , embed, this -> nplane , 1 ,
136+ (fftw_complex*)z_auxr, embed, this -> nplane , 1 , flag);
137+ this -> planyc2r = fftw_plan_many_dft_c2r ( 1 , & this -> ny , this -> nplane , (fftw_complex*)z_auxr, embed,
138+ this -> nplane , 1 , d_rspace, embed, this -> nplane , 1 , flag);
126139 }
127140 else
128141 {
129- this ->planxfor1 = fftw_plan_many_dft (1 , &this ->nx , this ->nplane * (this ->lixy + 1 ), (fftw_complex*)z_auxr, embed, npy,
130- 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_FORWARD, flag);
131- this ->planxbac1 = fftw_plan_many_dft (1 , &this ->nx , this ->nplane * (this ->lixy + 1 ), (fftw_complex*)z_auxr, embed, npy,
132- 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_BACKWARD, flag);
133- if (this ->gamma_only )
134- {
135- this ->planyr2c = fftw_plan_many_dft_r2c (1 , &this ->ny , this ->nplane , d_rspace, embed, this ->nplane , 1 ,
136- (fftw_complex*)z_auxr, embed, this ->nplane , 1 , flag);
137- this ->planyc2r = fftw_plan_many_dft_c2r (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed,
138- this ->nplane , 1 , d_rspace, embed, this ->nplane , 1 , flag);
139- }
140- else
141- {
142142
143- this ->planxfor2 = fftw_plan_many_dft (1 , &this ->nx , this ->nplane * (this ->ny - this ->rixy ), (fftw_complex*)z_auxr, embed,
144- npy, 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_FORWARD, flag);
145- this ->planxbac2 = fftw_plan_many_dft (1 , &this ->nx , this ->nplane * (this ->ny - this ->rixy ), (fftw_complex*)z_auxr, embed,
146- npy, 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_BACKWARD, flag);
147- this ->planyfor = fftw_plan_many_dft (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed, this ->nplane ,
148- 1 , (fftw_complex*)z_auxr, embed, this ->nplane , 1 , FFTW_FORWARD, flag);
149- this ->planybac = fftw_plan_many_dft (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed, this ->nplane ,
150- 1 , (fftw_complex*)z_auxr, embed, this ->nplane , 1 , FFTW_BACKWARD, flag);
151- }
143+ this ->planxfor2 = fftw_plan_many_dft (1 , &this ->nx , this ->nplane * (this ->ny - this ->rixy ), (fftw_complex*)z_auxr, embed,
144+ npy, 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_FORWARD, flag);
145+ this ->planxbac2 = fftw_plan_many_dft (1 , &this ->nx , this ->nplane * (this ->ny - this ->rixy ), (fftw_complex*)z_auxr, embed,
146+ npy, 1 , (fftw_complex*)z_auxr, embed, npy, 1 , FFTW_BACKWARD, flag);
147+ this ->planyfor = fftw_plan_many_dft (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed, this ->nplane ,
148+ 1 , (fftw_complex*)z_auxr, embed, this ->nplane , 1 , FFTW_FORWARD, flag);
149+ this ->planybac = fftw_plan_many_dft (1 , &this ->ny , this ->nplane , (fftw_complex*)z_auxr, embed, this ->nplane ,
150+ 1 , (fftw_complex*)z_auxr, embed, this ->nplane , 1 , FFTW_BACKWARD, flag);
152151 }
152+ }
153153 return ;
154154}
155155
0 commit comments