2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2003 David van der Spoel, Erik Lindahl, University of Groningen.
5 * Copyright (c) 2013,2014,2015,2017,2019, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
37 * \brief Fast Fourier Transforms.
39 * This file provides an abstract Gromacs interface to Fourier transforms,
40 * including multi-dimensional and real-to-complex transforms.
42 * Internally it is implemented as wrappers to external libraries such
43 * as FFTW or the Intel Math Kernel Library, but we also have a built-in
44 * version of FFTPACK in case the faster alternatives are unavailable.
46 * We also provide our own multi-dimensional transform setups even when
47 * the underlying library does not support it directly.
57 #include "gromacs/math/gmxcomplex.h"
58 #include "gromacs/utility/real.h"
60 /*! \brief Datatype for FFT setup
62 * The gmx_fft_t type contains all the setup information, e.g. twiddle
63 * factors, necessary to perform an FFT. Internally it is mapped to
64 * whatever FFT library we are using, or the built-in FFTPACK if no fast
65 * external library is available.
67 * Since some of the libraries (e.g. MKL) store work array data in their
68 * handles this datatype should only be used for one thread at a time, i.e.
69 * they should allocate one instance each when executing in parallel.
71 typedef struct gmx_fft* gmx_fft_t;
74 /*! \brief Specifier for FFT direction.
76 * The definition of the 1D forward transform from input x[] to output y[] is
78 * y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{-i 2 \pi j k /N}
81 * while the corresponding backward transform is
84 * y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{i 2 \pi j k /N}
87 * A forward-backward transform pair will this result in data scaled by N.
89 * For complex-to-complex transforms you can only use one of
90 * GMX_FFT_FORWARD or GMX_FFT_BACKWARD, and for real-complex transforms you
91 * can only use GMX_FFT_REAL_TO_COMPLEX or GMX_FFT_COMPLEX_TO_REAL.
93 enum gmx_fft_direction
95 GMX_FFT_FORWARD, /**< Forward complex-to-complex transform */
96 GMX_FFT_BACKWARD, /**< Backward complex-to-complex transform */
97 GMX_FFT_REAL_TO_COMPLEX, /**< Real-to-complex valued FFT */
98 GMX_FFT_COMPLEX_TO_REAL /**< Complex-to-real valued FFT */
101 /*! \brief Specifier for FFT flags.
103 * Some FFT libraries (FFTW, in particular) can do timings and other
104 * tricks to try and optimize the FFT for the current architecture. However,
105 * this can also lead to results that differ between consecutive runs with
107 * To avoid this, the conservative flag will attempt to disable such
108 * optimization, but there are no guarantees since we cannot control what
109 * the FFT libraries do internally.
112 typedef int gmx_fft_flag;
113 /** Macro to indicate no special flags for FFT routines. */
114 static const int GMX_FFT_FLAG_NONE = 0;
115 /** Flag to disable FFT optimizations based on timings, see ::gmx_fft_flag. */
116 static const int GMX_FFT_FLAG_CONSERVATIVE = (1 << 0);
118 /*! \brief Setup a 1-dimensional complex-to-complex transform
120 * \param fft Pointer to opaque Gromacs FFT datatype
121 * \param nx Length of transform
122 * \param flags FFT options
124 * \return status - 0 or a standard error message.
126 * \note Since some of the libraries (e.g. MKL) store work array data in their
127 * handles this datatype should only be used for one thread at a time,
128 * i.e. you should create one copy per thread when executing in parallel.
130 int gmx_fft_init_1d(gmx_fft_t* fft, int nx, gmx_fft_flag flags);
133 /*! \brief Setup multiple 1-dimensional complex-to-complex transform
135 * \param fft Pointer to opaque Gromacs FFT datatype
136 * \param nx Length of transform
137 * \param howmany Howmany 1D FFT
138 * \param flags FFT options
140 * \return status - 0 or a standard error message.
142 * \note Since some of the libraries (e.g. MKL) store work array data in their
143 * handles this datatype should only be used for one thread at a time,
144 * i.e. you should create one copy per thread when executing in parallel.
146 int gmx_fft_init_many_1d(gmx_fft_t* fft, int nx, int howmany, gmx_fft_flag flags);
149 /*! \brief Setup a 1-dimensional real-to-complex transform
151 * \param fft Pointer to opaque Gromacs FFT datatype
152 * \param nx Length of transform in real space
153 * \param flags FFT options
155 * \return status - 0 or a standard error message.
157 * \note Since some of the libraries (e.g. MKL) store work array data in their
158 * handles this datatype should only be used for one thread at a time,
159 * i.e. you should create one copy per thread when executing in parallel.
161 int gmx_fft_init_1d_real(gmx_fft_t* fft, int nx, gmx_fft_flag flags);
164 /*! \brief Setup multiple 1-dimensional real-to-complex transform
166 * \param fft Pointer to opaque Gromacs FFT datatype
167 * \param nx Length of transform in real space
168 * \param howmany Homany 1D FFTs
169 * \param flags FFT options
171 * \return status - 0 or a standard error message.
173 * \note Since some of the libraries (e.g. MKL) store work array data in their
174 * handles this datatype should only be used for one thread at a time,
175 * i.e. you should create one copy per thread when executing in parallel.
177 int gmx_fft_init_many_1d_real(gmx_fft_t* fft, int nx, int howmany, gmx_fft_flag flags);
180 /*! \brief Setup a 2-dimensional real-to-complex transform
182 * \param fft Pointer to opaque Gromacs FFT datatype
183 * \param nx Length of transform in first dimension
184 * \param ny Length of transform in second dimension
185 * \param flags FFT options
187 * The normal space is assumed to be real, while the values in
188 * frequency space are complex.
190 * \return status - 0 or a standard error message.
192 * \note Since some of the libraries (e.g. MKL) store work array data in their
193 * handles this datatype should only be used for one thread at a time,
194 * i.e. you should create one copy per thread when executing in parallel.
196 int gmx_fft_init_2d_real(gmx_fft_t* fft, int nx, int ny, gmx_fft_flag flags);
199 /*! \brief Perform a 1-dimensional complex-to-complex transform
201 * Performs an instance of a transform previously initiated.
203 * \param setup Setup returned from gmx_fft_init_1d()
204 * \param dir Forward or Backward
205 * \param in_data Input grid data. This should be allocated with gmx_new()
206 * to make it 16-byte aligned for better performance.
207 * \param out_data Output grid data. This should be allocated with gmx_new()
208 * to make it 16-byte aligned for better performance.
209 * You can provide the same pointer for in_data and out_data
210 * to perform an in-place transform.
212 * \return 0 on success, or an error code.
214 * \note Data pointers are declared as void, to avoid casting pointers
215 * depending on your grid type.
217 int gmx_fft_1d(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
220 /*! \brief Perform many 1-dimensional complex-to-complex transforms
222 * Performs many instances of a transform previously initiated.
224 * \param setup Setup returned from gmx_fft_init_1d()
225 * \param dir Forward or Backward
226 * \param in_data Input grid data. This should be allocated with gmx_new()
227 * to make it 16-byte aligned for better performance.
228 * \param out_data Output grid data. This should be allocated with gmx_new()
229 * to make it 16-byte aligned for better performance.
230 * You can provide the same pointer for in_data and out_data
231 * to perform an in-place transform.
233 * \return 0 on success, or an error code.
235 * \note Data pointers are declared as void, to avoid casting pointers
236 * depending on your grid type.
238 int gmx_fft_many_1d(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
241 /*! \brief Perform a 1-dimensional real-to-complex transform
243 * Performs an instance of a transform previously initiated.
245 * \param setup Setup returned from gmx_fft_init_1d_real()
246 * \param dir Real-to-complex or complex-to-real
247 * \param in_data Input grid data. This should be allocated with gmx_new()
248 * to make it 16-byte aligned for better performance.
249 * \param out_data Output grid data. This should be allocated with gmx_new()
250 * to make it 16-byte aligned for better performance.
251 * You can provide the same pointer for in_data and out_data
252 * to perform an in-place transform.
254 * If you are doing an in-place transform, the array must be padded up to
255 * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
256 * should not be padded (although it doesn't matter in 1d).
258 * \return 0 on success, or an error code.
260 * \note Data pointers are declared as void, to avoid casting pointers
261 * depending on transform direction.
263 int gmx_fft_1d_real(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
265 /*! \brief Perform many 1-dimensional real-to-complex transforms
267 * Performs many instances of a transform previously initiated.
269 * \param setup Setup returned from gmx_fft_init_1d_real()
270 * \param dir Real-to-complex or complex-to-real
271 * \param in_data Input grid data. This should be allocated with gmx_new()
272 * to make it 16-byte aligned for better performance.
273 * \param out_data Output grid data. This should be allocated with gmx_new()
274 * to make it 16-byte aligned for better performance.
275 * You can provide the same pointer for in_data and out_data
276 * to perform an in-place transform.
278 * If you are doing an in-place transform, the array must be padded up to
279 * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
280 * should not be padded (although it doesn't matter in 1d).
282 * \return 0 on success, or an error code.
284 * \note Data pointers are declared as void, to avoid casting pointers
285 * depending on transform direction.
287 int gmx_fft_many_1d_real(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
289 /*! \brief Perform a 2-dimensional real-to-complex transform
291 * Performs an instance of a transform previously initiated.
293 * \param setup Setup returned from gmx_fft_init_1d_real()
294 * \param dir Real-to-complex or complex-to-real
295 * \param in_data Input grid data. This should be allocated with gmx_new()
296 * to make it 16-byte aligned for better performance.
297 * \param out_data Output grid data. This should be allocated with gmx_new()
298 * to make it 16-byte aligned for better performance.
299 * You can provide the same pointer for in_data and out_data
300 * to perform an in-place transform.
302 * \return 0 on success, or an error code.
304 * \note If you are doing an in-place transform, the last dimension of the
305 * array MUST be padded up to an even integer length so n/2 complex numbers can
306 * fit. Thus, if the real grid e.g. has dimension 5*3, you must allocate it as
307 * a 5*4 array, where the last element in the second dimension is padding.
308 * The complex data will be written to the same array, but since that dimension
309 * is 5*2 it will now fill the entire array. Reverse complex-to-real in-place
310 * transformation will produce the same sort of padded array.
312 * The padding does NOT apply to out-of-place transformation. In that case the
313 * input array will simply be 5*3 of real, while the output is 5*2 of complex.
315 * \note Data pointers are declared as void, to avoid casting pointers
316 * depending on transform direction.
318 int gmx_fft_2d_real(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
320 /*! \brief Release an FFT setup structure
322 * Destroy setup and release all allocated memory.
324 * \param setup Setup returned from gmx_fft_init_1d(), or one
325 * of the other initializers.
328 void gmx_fft_destroy(gmx_fft_t setup);
330 /*! \brief Release a many FFT setup structure
332 * Destroy setup and release all allocated memory.
334 * \param setup Setup returned from gmx_fft_init_1d(), or one
335 * of the other initializers.
338 void gmx_many_fft_destroy(gmx_fft_t setup);
341 /*! \brief Transpose 2d complex matrix, in-place or out-of-place.
343 * This routines works when the matrix is non-square, i.e. nx!=ny too,
344 * without allocating an entire matrix of work memory, which is important
345 * for huge FFT grids.
347 * \param in_data Input data, to be transposed
348 * \param out_data Output, transposed data. If this is identical to
349 * in_data, an in-place transpose is performed.
350 * \param nx Number of rows before transpose
351 * \param ny Number of columns before transpose
353 * \return GMX_SUCCESS, or an error code from gmx_errno.h
355 int gmx_fft_transpose_2d(t_complex* in_data, t_complex* out_data, int nx, int ny);
357 /*! \brief Cleanup global data of FFT
359 * Any plans are invalid after this function. Should be called
360 * after all plans have been destroyed.
362 void gmx_fft_cleanup();