18 #ifdef USE_FACTORY_AUTOREGISTER
20 bool init = FFT_xyz_3dim::register_factory();
24 const std::string FFT_xyz_3dim::class_name =
"FFT_xyz_3dim";
27 void FFT_xyz_3dim::set_parameters(
const Parameters& params)
29 const std::string str_vlevel = params.
get_string(
"verbose_level");
34 string str_fft_direction;
37 err += params.
fetch_string(
"FFT_direction", str_fft_direction);
40 vout.
crucial(m_vl,
"Error at %s: input parameter not found.\n", class_name.c_str());
44 set_parameters(str_fft_direction);
49 void FFT_xyz_3dim::set_parameters(
const string& str_fft_direction)
53 vout.
general(m_vl,
" FFT_direction = %s\n", str_fft_direction.c_str());
58 if (str_fft_direction ==
"Forward") {
60 }
else if (str_fft_direction ==
"Backward") {
63 vout.
crucial(m_vl,
"Error at %s: unsupported FFT direction \"%s\"\n", class_name.c_str(), str_fft_direction.c_str());
70 void FFT_xyz_3dim::init()
78 int threads_ok = fftw_init_threads();
87 if ((NPE_x * NPE_y * NPE_t) != 1) {
88 vout.
crucial(m_vl,
"Error at %s: FFTW supports parallelization only in z-direction.\n",
102 ptrdiff_t fftw_size_p = fftw_mpi_local_size_3d(Lz_p, Ly_p, Lx_p,
104 &m_Nz_p, &m_z_start_p);
106 m_in = fftw_alloc_complex(fftw_size_p);
107 m_out = fftw_alloc_complex(fftw_size_p);
109 if (!m_in || !m_out) {
110 vout.
crucial(m_vl,
"Error at %s: failed to allocate memory %d [Byte].\n",
111 class_name.c_str(), (int)fftw_size_p);
116 const size_t fftw_size =
sizeof(fftw_complex) * Lx * Ly * Lz;
117 m_in = (fftw_complex *)fftw_malloc(fftw_size);
118 m_out = (fftw_complex *)fftw_malloc(fftw_size);
120 if (!m_in || !m_out) {
121 vout.
crucial(m_vl,
"Error at %s: failed to allocate memory %d [Byte].\n",
122 class_name.c_str(), (int)fftw_size);
130 void FFT_xyz_3dim::tidy_up()
132 if (m_in) fftw_free(m_in);
133 if (m_out) fftw_free(m_out);
134 if (m_plan) fftw_destroy_plan(m_plan);
139 void FFT_xyz_3dim::fft(
Field& field)
146 const int Lxyz = Lx * Ly * Lz;
151 const int Nin = field.
nin();
152 const int Nex = field.
nex();
158 fftw_plan_with_nthreads(Nthread);
166 m_plan = fftw_mpi_plan_dft_3d(Lz_p, Ly_p, Lx_p, m_in, m_out,
168 FFTW_FORWARD, FFTW_ESTIMATE);
170 m_plan = fftw_mpi_plan_dft_3d(Lz_p, Ly_p, Lx_p, m_in, m_out,
172 FFTW_BACKWARD, FFTW_ESTIMATE);
176 m_plan = fftw_plan_dft_3d(Lz, Ly, Lx, m_in, m_out,
177 FFTW_FORWARD, FFTW_ESTIMATE);
179 m_plan = fftw_plan_dft_3d(Lz, Ly, Lx, m_in, m_out,
180 FFTW_BACKWARD, FFTW_ESTIMATE);
187 for (
int in2 = 0; in2 < Nin / 2; ++in2) {
188 for (
int t_global = 0; t_global < Lt; t_global++) {
189 for (
int ex = 0; ex < Nex; ++ex) {
191 for (
int z = 0; z < Nz; z++) {
192 for (
int y_global = 0; y_global < Ly; y_global++) {
193 for (
int x_global = 0; x_global < Lx; x_global++) {
194 int isite_xyz_local = x_global + Lx * (y_global + Ly * z);
196 int isite = m_index.site(x_global, y_global, z, t_global);
197 int i_real = 2 * in2;
198 int i_imag = 2 * in2 + 1;
200 m_in[isite_xyz_local][0] = field.
cmp(i_real, isite, ex);
201 m_in[isite_xyz_local][1] = field.
cmp(i_imag, isite, ex);
207 fftw_execute(m_plan);
211 for (
int z = 0; z < Nz; z++) {
212 for (
int y_global = 0; y_global < Ly; y_global++) {
213 for (
int x_global = 0; x_global < Lx; x_global++) {
214 int isite_xyz_local = x_global + Lx * (y_global + Ly * z);
216 int isite = m_index.site(x_global, y_global, z, t_global);
217 int i_real = 2 * in2;
218 int i_imag = 2 * in2 + 1;
220 field.
set(i_real, isite, ex, m_out[isite_xyz_local][0]);
221 field.
set(i_imag, isite, ex, m_out[isite_xyz_local][1]);
232 scal(field, 1.0 / Lxyz);
238 void FFT_xyz_3dim::fft(
Field& field_out,
const Field& field_in)
245 const int Lxyz = Lx * Ly * Lz;
250 const int Nin = field_in.
nin();
251 const int Nex = field_in.
nex();
257 fftw_plan_with_nthreads(Nthread);
265 m_plan = fftw_mpi_plan_dft_3d(Lz_p, Ly_p, Lx_p, m_in, m_out,
267 FFTW_FORWARD, FFTW_ESTIMATE);
269 m_plan = fftw_mpi_plan_dft_3d(Lz_p, Ly_p, Lx_p, m_in, m_out,
271 FFTW_BACKWARD, FFTW_ESTIMATE);
275 m_plan = fftw_plan_dft_3d(Lz, Ly, Lx, m_in, m_out,
276 FFTW_FORWARD, FFTW_ESTIMATE);
278 m_plan = fftw_plan_dft_3d(Lz, Ly, Lx, m_in, m_out,
279 FFTW_BACKWARD, FFTW_ESTIMATE);
286 for (
int in2 = 0; in2 < Nin / 2; ++in2) {
287 for (
int t_global = 0; t_global < Lt; t_global++) {
288 for (
int ex = 0; ex < Nex; ++ex) {
290 for (
int z = 0; z < Nz; z++) {
291 for (
int y_global = 0; y_global < Ly; y_global++) {
292 for (
int x_global = 0; x_global < Lx; x_global++) {
293 int isite_xyz_local = x_global + Lx * (y_global + Ly * z);
295 int isite = m_index.site(x_global, y_global, z, t_global);
296 int i_real = 2 * in2;
297 int i_imag = 2 * in2 + 1;
299 m_in[isite_xyz_local][0] = field_in.
cmp(i_real, isite, ex);
300 m_in[isite_xyz_local][1] = field_in.
cmp(i_imag, isite, ex);
306 fftw_execute(m_plan);
310 for (
int z = 0; z < Nz; z++) {
311 for (
int y_global = 0; y_global < Ly; y_global++) {
312 for (
int x_global = 0; x_global < Lx; x_global++) {
313 int isite_xyz_local = x_global + Lx * (y_global + Ly * z);
315 int isite = m_index.site(x_global, y_global, z, t_global);
316 int i_real = 2 * in2;
317 int i_imag = 2 * in2 + 1;
319 field_out.
set(i_real, isite, ex, m_out[isite_xyz_local][0]);
320 field_out.
set(i_imag, isite, ex, m_out[isite_xyz_local][1]);
331 scal(field_out, 1.0 / Lxyz);
340 bool backup_fwbw = m_is_forward;
343 if (dir == FORWARD) {
345 }
else if (dir == BACKWARD) {
346 m_is_forward =
false;
348 vout.
crucial(m_vl,
"%s: unknown direction %d. failed.\n", class_name.c_str(), dir);
353 fft(field_out, field_in);
356 m_is_forward = backup_fwbw;
void scal(Field &x, const double a)
scal(x, a): x = a * x
static int get_num_threads()
returns available number of threads.
void set(const int jin, const int site, const int jex, double v)
void general(const char *format,...)
Container of Field-type object.
double cmp(const int jin, const int site, const int jex) const
int fetch_string(const string &key, string &value) const
static MPI_Comm & world()
retrieves current communicator.
void crucial(const char *format,...)
string get_string(const string &key) const
static VerboseLevel set_verbose_level(const std::string &str)