1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * Gromacs 4.0 Copyright (c) 1991-2003
5 * David van der Spoel, Erik Lindahl, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
19 * \brief Fast Fourier Transforms.
21 * This file provides an abstract Gromacs interface to Fourier transforms,
22 * including multi-dimensional and real-to-complex transforms.
24 * Internally it is implemented as wrappers to external libraries such
25 * as FFTW or the Intel Math Kernel Library, but we also have a built-in
26 * version of FFTPACK in case the faster alternatives are unavailable.
28 * We also provide our own multi-dimensional transform setups even when
29 * the underlying library does not support it directly.
38 #include "../legacyheaders/types/simple.h"
39 #include "../legacyheaders/gmxcomplex.h"
45 } /* fixes auto-indentation problems */
50 /*! \brief Datatype for FFT setup
52 * The gmx_fft_t type contains all the setup information, e.g. twiddle
53 * factors, necessary to perform an FFT. Internally it is mapped to
54 * whatever FFT library we are using, or the built-in FFTPACK if no fast
55 * external library is available.
57 * Since some of the libraries (e.g. MKL) store work array data in their
58 * handles this datatype should only be used for one thread at a time, i.e.
59 * they should allocate one instance each when executing in parallel.
61 typedef struct gmx_fft *
67 /*! \brief Specifier for FFT direction.
69 * The definition of the 1D forward transform from input x[] to output y[] is
71 * y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{-i 2 \pi j k /N}
74 * while the corresponding backward transform is
77 * y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{i 2 \pi j k /N}
80 * A forward-backward transform pair will this result in data scaled by N.
82 * For complex-to-complex transforms you can only use one of
83 * GMX_FFT_FORWARD or GMX_FFT_BACKWARD, and for real-complex transforms you
84 * can only use GMX_FFT_REAL_TO_COMPLEX or GMX_FFT_COMPLEX_TO_REAL.
86 typedef enum gmx_fft_direction
88 GMX_FFT_FORWARD, /*!< Forward complex-to-complex transform */
89 GMX_FFT_BACKWARD, /*!< Backward complex-to-complex transform */
90 GMX_FFT_REAL_TO_COMPLEX, /*!< Real-to-complex valued fft */
91 GMX_FFT_COMPLEX_TO_REAL /*!< Complex-to-real valued fft */
94 /*! \brief Specifier for FFT flags.
96 * Some FFT libraries (FFTW, in particular) can do timings and other
97 * tricks to try and optimize the FFT for the current architecture. However,
98 * this can also lead to results that differ between consecutive runs with
100 * To avoid this, the conservative flag will attempt to disable such
101 * optimization, but there are no guarantees since we cannot control what
102 * the FFT libraries do internally.
105 typedef int gmx_fft_flag;
106 /** Macro to indicate no special flags for FFT routines. */
107 static const int GMX_FFT_FLAG_NONE = 0;
108 /** Flag to disable FFT optimizations based on timings, see ::gmx_fft_flag. */
109 static const int GMX_FFT_FLAG_CONSERVATIVE = (1<<0);
111 /*! \brief Setup a 1-dimensional complex-to-complex transform
113 * \param fft Pointer to opaque Gromacs FFT datatype
114 * \param nx Length of transform
115 * \param flags FFT options
117 * \return status - 0 or a standard error message.
119 * \note Since some of the libraries (e.g. MKL) store work array data in their
120 * handles this datatype should only be used for one thread at a time,
121 * i.e. you should create one copy per thread when executing in parallel.
124 gmx_fft_init_1d (gmx_fft_t * fft,
129 /*! \brief Setup multiple 1-dimensional complex-to-complex transform
131 * \param fft Pointer to opaque Gromacs FFT datatype
132 * \param nx Length of transform
133 * \param howmany Howmany 1D FFT
134 * \param flags FFT options
136 * \return status - 0 or a standard error message.
138 * \note Since some of the libraries (e.g. MKL) store work array data in their
139 * handles this datatype should only be used for one thread at a time,
140 * i.e. you should create one copy per thread when executing in parallel.
143 gmx_fft_init_many_1d (gmx_fft_t * fft,
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.
162 gmx_fft_init_1d_real (gmx_fft_t * fft,
167 /*! \brief Setup multiple 1-dimensional real-to-complex transform
169 * \param fft Pointer to opaque Gromacs FFT datatype
170 * \param nx Length of transform in real space
171 * \param howmany Homany 1D FFTs
172 * \param flags FFT options
174 * \return status - 0 or a standard error message.
176 * \note Since some of the libraries (e.g. MKL) store work array data in their
177 * handles this datatype should only be used for one thread at a time,
178 * i.e. you should create one copy per thread when executing in parallel.
181 gmx_fft_init_many_1d_real (gmx_fft_t * fft,
187 /*! \brief Setup a 2-dimensional real-to-complex transform
189 * \param fft Pointer to opaque Gromacs FFT datatype
190 * \param nx Length of transform in first dimension
191 * \param ny Length of transform in second dimension
192 * \param flags FFT options
194 * The normal space is assumed to be real, while the values in
195 * frequency space are complex.
197 * \return status - 0 or a standard error message.
199 * \note Since some of the libraries (e.g. MKL) store work array data in their
200 * handles this datatype should only be used for one thread at a time,
201 * i.e. you should create one copy per thread when executing in parallel.
204 gmx_fft_init_2d_real (gmx_fft_t * fft,
210 /*! \brief Perform a 1-dimensional complex-to-complex transform
212 * Performs an instance of a transform previously initiated.
214 * \param setup Setup returned from gmx_fft_init_1d()
215 * \param dir Forward or Backward
216 * \param in_data Input grid data. This should be allocated with gmx_new()
217 * to make it 16-byte aligned for better performance.
218 * \param out_data Output grid data. This should be allocated with gmx_new()
219 * to make it 16-byte aligned for better performance.
220 * You can provide the same pointer for in_data and out_data
221 * to perform an in-place transform.
223 * \return 0 on success, or an error code.
225 * \note Data pointers are declared as void, to avoid casting pointers
226 * depending on your grid type.
229 gmx_fft_1d (gmx_fft_t setup,
230 enum gmx_fft_direction dir,
235 /*! \brief Perform many 1-dimensional complex-to-complex transforms
237 * Performs many instances of a transform previously initiated.
239 * \param setup Setup returned from gmx_fft_init_1d()
240 * \param dir Forward or Backward
241 * \param in_data Input grid data. This should be allocated with gmx_new()
242 * to make it 16-byte aligned for better performance.
243 * \param out_data Output grid data. This should be allocated with gmx_new()
244 * to make it 16-byte aligned for better performance.
245 * You can provide the same pointer for in_data and out_data
246 * to perform an in-place transform.
248 * \return 0 on success, or an error code.
250 * \note Data pointers are declared as void, to avoid casting pointers
251 * depending on your grid type.
254 gmx_fft_many_1d (gmx_fft_t setup,
255 enum gmx_fft_direction dir,
260 /*! \brief Perform a 1-dimensional real-to-complex transform
262 * Performs an instance of a transform previously initiated.
264 * \param setup Setup returned from gmx_fft_init_1d_real()
265 * \param dir Real-to-complex or complex-to-real
266 * \param in_data Input grid data. This should be allocated with gmx_new()
267 * to make it 16-byte aligned for better performance.
268 * \param out_data Output grid data. This should be allocated with gmx_new()
269 * to make it 16-byte aligned for better performance.
270 * You can provide the same pointer for in_data and out_data
271 * to perform an in-place transform.
273 * If you are doing an in-place transform, the array must be padded up to
274 * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
275 * should not be padded (although it doesn't matter in 1d).
277 * \return 0 on success, or an error code.
279 * \note Data pointers are declared as void, to avoid casting pointers
280 * depending on transform direction.
283 gmx_fft_1d_real (gmx_fft_t setup,
284 enum gmx_fft_direction dir,
288 /*! \brief Perform many 1-dimensional real-to-complex transforms
290 * Performs many instances of a transform previously initiated.
292 * \param setup Setup returned from gmx_fft_init_1d_real()
293 * \param dir Real-to-complex or complex-to-real
294 * \param in_data Input grid data. This should be allocated with gmx_new()
295 * to make it 16-byte aligned for better performance.
296 * \param out_data Output grid data. This should be allocated with gmx_new()
297 * to make it 16-byte aligned for better performance.
298 * You can provide the same pointer for in_data and out_data
299 * to perform an in-place transform.
301 * If you are doing an in-place transform, the array must be padded up to
302 * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
303 * should not be padded (although it doesn't matter in 1d).
305 * \return 0 on success, or an error code.
307 * \note Data pointers are declared as void, to avoid casting pointers
308 * depending on transform direction.
311 gmx_fft_many_1d_real (gmx_fft_t setup,
312 enum gmx_fft_direction dir,
316 /*! \brief Perform a 2-dimensional real-to-complex transform
318 * Performs an instance of a transform previously initiated.
320 * \param setup Setup returned from gmx_fft_init_1d_real()
321 * \param dir Real-to-complex or complex-to-real
322 * \param in_data Input grid data. This should be allocated with gmx_new()
323 * to make it 16-byte aligned for better performance.
324 * \param out_data Output grid data. This should be allocated with gmx_new()
325 * to make it 16-byte aligned for better performance.
326 * You can provide the same pointer for in_data and out_data
327 * to perform an in-place transform.
329 * \return 0 on success, or an error code.
331 * \note If you are doing an in-place transform, the last dimension of the
332 * array MUST be padded up to an even integer length so n/2 complex numbers can
333 * fit. Thus, if the real grid e.g. has dimension 5*3, you must allocate it as
334 * a 5*4 array, where the last element in the second dimension is padding.
335 * The complex data will be written to the same array, but since that dimension
336 * is 5*2 it will now fill the entire array. Reverse complex-to-real in-place
337 * transformation will produce the same sort of padded array.
339 * The padding does NOT apply to out-of-place transformation. In that case the
340 * input array will simply be 5*3 of real, while the output is 5*2 of complex.
342 * \note Data pointers are declared as void, to avoid casting pointers
343 * depending on transform direction.
346 gmx_fft_2d_real (gmx_fft_t setup,
347 enum gmx_fft_direction dir,
351 /*! \brief Release an FFT setup structure
353 * Destroy setup and release all allocated memory.
355 * \param setup Setup returned from gmx_fft_init_1d(), or one
356 * of the other initializers.
360 gmx_fft_destroy (gmx_fft_t setup);
362 /*! \brief Release a many FFT setup structure
364 * Destroy setup and release all allocated memory.
366 * \param setup Setup returned from gmx_fft_init_1d(), or one
367 * of the other initializers.
371 gmx_many_fft_destroy (gmx_fft_t setup);
374 /*! \brief Transpose 2d complex matrix, in-place or out-of-place.
376 * This routines works when the matrix is non-square, i.e. nx!=ny too,
377 * without allocating an entire matrix of work memory, which is important
378 * for huge FFT grids.
380 * \param in_data Input data, to be transposed
381 * \param out_data Output, transposed data. If this is identical to
382 * in_data, an in-place transpose is performed.
383 * \param nx Number of rows before transpose
384 * \param ny Number of columns before transpose
386 * \return GMX_SUCCESS, or an error code from gmx_errno.h
389 gmx_fft_transpose_2d (t_complex * in_data,
390 t_complex * out_data,
394 /*! \brief Cleanup global data of FFT
396 * Any plans are invalid after this function. Should be called
397 * after all plans have been destroyed.
399 void gmx_fft_cleanup();
401 /*! \brief Return string describing the underlying FFT implementation.
403 * Used to print out information about the used FFT library where needed.
405 const char *gmx_fft_get_version_info();