Bridge++  Ver. 2.0.2
fft_3d_local.cpp
Go to the documentation of this file.
1 
13 #ifdef USE_FFTWLIB
14 
15 #include "fft_3d_local.h"
16 #include "Field/index_lex.h"
17 #include <cstring>
18 
19 #ifdef USE_OPENMP
21 #endif
22 
23 const std::string FFT_3d_local::class_name = "FFT_3d_local";
24 
25 #ifdef USE_FACTORY_AUTOREGISTER
26 namespace {
27  bool init = FFT_3d_local::register_factory();
28 }
29 #endif
30 
31 //====================================================================
32 FFT_3d_local::FFT_3d_local()
33  : m_vl(CommonParameters::Vlevel())
34  , m_ndim(0)
35  , m_vol(0)
36  , m_nv(0)
37  , m_buf_in(NULL)
38  , m_buf_out(NULL)
39  , m_plan_fw(NULL)
40  , m_plan_bw(NULL)
41  , m_direction(UNDEF)
42 {
43  if (check_ok()) {
44  initialize();
45  }
46 }
47 
48 
49 //====================================================================
50 FFT_3d_local::FFT_3d_local(const Parameters& params)
51  : m_vl(CommonParameters::Vlevel())
52  , m_ndim(0)
53  , m_vol(0)
54  , m_nv(0)
55  , m_buf_in(NULL)
56  , m_buf_out(NULL)
57  , m_plan_fw(NULL)
58  , m_plan_bw(NULL)
59  , m_direction(UNDEF)
60 {
61  if (check_ok()) {
62  initialize();
63  }
64  set_parameters(params);
65 }
66 
67 
68 //====================================================================
69 FFT_3d_local::~FFT_3d_local()
70 {
71  finalize();
72 }
73 
74 
75 //====================================================================
76 void FFT_3d_local::set_parameters(const Parameters& params)
77 {
78  std::string vlevel;
79  if (!params.fetch_string("verbose_level", vlevel)) {
80  m_vl = vout.set_verbose_level(vlevel);
81  }
82 
83  std::string direction;
84  if (!params.fetch_string("FFT_direction", direction)) {
85  set_parameters(direction);
86  }
87 }
88 
89 
90 //====================================================================
91 void FFT_3d_local::get_parameters(Parameters& params) const
92 {
93  if (m_direction == FORWARD) {
94  params.set_string("FFT_direction", "Forward");
95  } else if (m_direction == BACKWARD) {
96  params.set_string("FFT_direction", "Backward");
97  } else {
98  params.set_string("FFT_direction", "None");
99  }
100 
101  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
102 }
103 
104 
105 //====================================================================
106 void FFT_3d_local::set_parameters(const std::string& direction)
107 {
108  if (direction == "Forward") {
109  m_direction = FORWARD;
110  } else if (direction == "Backward") {
111  m_direction = BACKWARD;
112  } else {
113  m_direction = UNDEF;
114 
115  vout.crucial(m_vl, "Error at %s: unsupported FFT direction \"%s\"\n", class_name.c_str(), direction.c_str());
116  exit(EXIT_FAILURE);
117  }
118 }
119 
120 
121 //====================================================================
122 bool FFT_3d_local::check_ok()
123 {
124  int npe_xyz = Communicator::npe(0) * Communicator::npe(1) * Communicator::npe(2);
125 
126  if (npe_xyz > 1) {
127  vout.crucial(m_vl, "%s: incompatible with xyz parallelization.\n", class_name.c_str());
128  exit(EXIT_FAILURE);
129  }
130 
131  return true;
132 }
133 
134 
135 //====================================================================
136 void FFT_3d_local::initialize()
137 {
138 #ifdef USE_OPENMP
139  int thread_ok = fftw_init_threads();
140 
141  if (thread_ok) {
142  fftw_plan_with_nthreads(ThreadManager::get_num_threads());
143  }
144 #endif
145 }
146 
147 
148 //====================================================================
149 void FFT_3d_local::initialize_plan(const Field& src)
150 {
151  if ((m_nv == src.nin() / 2) && (m_vol == CommonParameters::Lvol() / CommonParameters::Lt())) {
152  vout.detailed(m_vl, "%s: plan recycled.\n", class_name.c_str());
153  return;
154  } else {
155  vout.detailed(m_vl, "%s: create plan.\n", class_name.c_str());
156  }
157 
158  // first, clear pre-existing plan, if any.
159  clear_plan();
160 
161  // assume 3dimensional
162  m_ndim = 3;
163  // row-major fortran order
164  m_nsize[0] = CommonParameters::Lz();
165  m_nsize[1] = CommonParameters::Ly();
166  m_nsize[2] = CommonParameters::Lx();
167 
168  // local volume
169  m_vol = 1;
170  for (int i = 0; i < m_ndim; ++i) {
171  m_vol *= m_nsize[i];
172  }
173 
174  // number of complex elements. run nin at a time.
175  m_nv = src.nin() / 2;
176 
177  m_buf_in = fftw_alloc_complex(m_nv * m_vol);
178  m_buf_out = fftw_alloc_complex(m_nv * m_vol);
179 
180  if ((!m_buf_in) || (!m_buf_out)) {
181  vout.crucial(m_vl, "%s: memory allocation failed.\n", class_name.c_str());
182  exit(EXIT_FAILURE);
183  }
184 
185  m_plan_fw = fftw_plan_many_dft(m_ndim, m_nsize, m_nv,
186  m_buf_in, m_nsize, m_nv, 1,
187  m_buf_out, m_nsize, m_nv, 1,
188  FFTW_FORWARD, FFTW_ESTIMATE);
189 
190  m_plan_bw = fftw_plan_many_dft(m_ndim, m_nsize, m_nv,
191  m_buf_in, m_nsize, m_nv, 1,
192  m_buf_out, m_nsize, m_nv, 1,
193  FFTW_BACKWARD, FFTW_ESTIMATE);
194 
195  if ((!m_plan_fw) || (!m_plan_bw)) {
196  vout.crucial(m_vl, "%s: create plan failed.\n", class_name.c_str());
197  exit(EXIT_FAILURE);
198  }
199 }
200 
201 
202 //====================================================================
203 void FFT_3d_local::clear_plan()
204 {
205  if (m_buf_in) fftw_free(m_buf_in);
206  if (m_buf_out) fftw_free(m_buf_out);
207 
208  if (m_plan_fw) fftw_destroy_plan(m_plan_fw);
209  if (m_plan_bw) fftw_destroy_plan(m_plan_bw);
210 }
211 
212 
213 //====================================================================
214 void FFT_3d_local::finalize()
215 {
216  clear_plan();
217 }
218 
219 
220 //====================================================================
221 void FFT_3d_local::fft(Field& dst, const Field& src, enum Direction dir)
222 {
223  if (not ((dir == FORWARD) || (dir == BACKWARD))) {
224  vout.crucial(m_vl, "%s: unsupported direction. %d\n", class_name.c_str(), dir);
225  exit(EXIT_FAILURE);
226  }
227 
228  initialize_plan(src);
229 
230  int nex = src.nex();
231  int nt = CommonParameters::Nt();
232 
233  Index_lex index;
234 
235  size_t count = m_nv * m_vol; // count in complex numbers
236 
237  for (int iex = 0; iex < nex; ++iex) {
238  for (int it = 0; it < nt; ++it) {
239  memcpy(m_buf_in, src.ptr(0, index.site(0, 0, 0, it), iex), sizeof(double) * count * 2);
240 
241  if (dir == FORWARD) {
242  fftw_execute(m_plan_fw);
243  } else if (dir == BACKWARD) {
244  fftw_execute(m_plan_bw);
245  } else {
246  vout.crucial(m_vl, "%s: unsupported direction. %d\n", class_name.c_str(), dir);
247  exit(EXIT_FAILURE);
248  }
249 
250  memcpy(dst.ptr(0, index.site(0, 0, 0, it), iex), m_buf_out, sizeof(double) * count * 2);
251  }
252  }
253 
254  if (dir == BACKWARD) {
255  scal(dst, 1.0 / m_vol);
256  }
257 }
258 
259 
260 //====================================================================
261 void FFT_3d_local::fft(Field& dst, const Field& src)
262 {
263  return fft(dst, src, m_direction);
264 }
265 
266 
267 //====================================================================
268 void FFT_3d_local::fft(Field& field)
269 {
270  // return fft(field, field, m_direction);
271  vout.crucial(m_vl, "Error at %s: fft on-the-fly unsupported.\n", class_name.c_str());
272  exit(EXIT_FAILURE);
273 }
274 
275 
276 //====================================================================
277 #endif /* USE_FFTWLIB */
278 
279 //====================================================================
280 //============================================================END=====
CommonParameters::Lvol
static long_t Lvol()
Definition: commonParameters.h:95
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
CommonParameters
Common parameter class: provides parameters as singleton.
Definition: commonParameters.h:42
Index_lex
Lexical site index.
Definition: index_lex.h:34
ThreadManager::get_num_threads
static int get_num_threads()
returns available number of threads.
Definition: threadManager.cpp:246
Parameters
Class for parameters.
Definition: parameters.h:46
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
Field::nex
int nex() const
Definition: field.h:128
CommonParameters::Ly
static int Ly()
Definition: commonParameters.h:92
Direction
Direction
Definition: bridge_defs.h:24
Field::nin
int nin() const
Definition: field.h:126
CommonParameters::Lt
static int Lt()
Definition: commonParameters.h:94
CommonParameters::Lx
static int Lx()
Definition: commonParameters.h:91
CommonParameters::Lz
static int Lz()
Definition: commonParameters.h:93
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Communicator::npe
static int npe(const int dir)
logical grid extent
Definition: communicator.cpp:112
threadManager.h
index_lex.h
Index_lex::site
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:55
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Field
Container of Field-type object.
Definition: field.h:46
fft_3d_local.h
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154