src/gromacs/fft/fft.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   9  *
  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.
  14  *
  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.
  19  *
  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.
  24  *
  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.
  32  *
  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.
  35  */
  36 /*! \file
  37  *  \brief Fast Fourier Transforms.
  38  *
  39  *  This file provides an abstract Gromacs interface to Fourier transforms,
  40  *  including multi-dimensional and real-to-complex transforms.
  41  *
  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.
  45  *
  46  *  We also provide our own multi-dimensional transform setups even when
  47  *  the underlying library does not support it directly.
  48  *
  49  * \inpublicapi
  50  * \ingroup module_fft
  51  */
  52 #ifndef GMX_FFT_FFT_H
  53 #define GMX_FFT_FFT_H
  54
  55 #include <stdio.h>
  56
  57 #include "gromacs/math/gmxcomplex.h"
  58 #include "gromacs/utility/real.h"
  59
  60 /*! \brief Datatype for FFT setup
  61  *
  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.
  66  *
  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.
  70  */
  71 typedef struct gmx_fft* gmx_fft_t;
  72
  73
  74 /*! \brief Specifier for FFT direction.
  75  *
  76  *  The definition of the 1D forward transform from input x[] to output y[] is
  77  *  \f[
  78  *  y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{-i 2 \pi j k /N}
  79  *  \f]
  80  *
  81  *  while the corresponding backward transform is
  82  *
  83  *  \f[
  84  *  y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{i 2 \pi j k /N}
  85  *  \f]
  86  *
  87  *  A forward-backward transform pair will this result in data scaled by N.
  88  *
  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.
  92  */
  93 enum gmx_fft_direction
  94 {
  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            */
  99 };
 100
 101 /*! \brief Specifier for FFT flags.
 102  *
 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
 106  *  identical input.
 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.
 110  */
 111
 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);
 117
 118 /*! \brief Setup a 1-dimensional complex-to-complex transform
 119  *
 120  *  \param fft    Pointer to opaque Gromacs FFT datatype
 121  *  \param nx     Length of transform
 122  *  \param flags  FFT options
 123  *
 124  *  \return status - 0 or a standard error message.
 125  *
 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.
 129  */
 130 int gmx_fft_init_1d(gmx_fft_t* fft, int nx, gmx_fft_flag flags);
 131
 132
 133 /*! \brief Setup multiple 1-dimensional complex-to-complex transform
 134  *
 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
 139  *
 140  *  \return status - 0 or a standard error message.
 141  *
 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.
 145  */
 146 int gmx_fft_init_many_1d(gmx_fft_t* fft, int nx, int howmany, gmx_fft_flag flags);
 147
 148
 149 /*! \brief Setup a 1-dimensional real-to-complex transform
 150  *
 151  *  \param fft    Pointer to opaque Gromacs FFT datatype
 152  *  \param nx     Length of transform in real space
 153  *  \param flags  FFT options
 154  *
 155  *  \return status - 0 or a standard error message.
 156  *
 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.
 160  */
 161 int gmx_fft_init_1d_real(gmx_fft_t* fft, int nx, gmx_fft_flag flags);
 162
 163
 164 /*! \brief Setup multiple 1-dimensional real-to-complex transform
 165  *
 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
 170  *
 171  *  \return status - 0 or a standard error message.
 172  *
 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.
 176  */
 177 int gmx_fft_init_many_1d_real(gmx_fft_t* fft, int nx, int howmany, gmx_fft_flag flags);
 178
 179
 180 /*! \brief Setup a 2-dimensional real-to-complex transform
 181  *
 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
 186  *
 187  *  The normal space is assumed to be real, while the values in
 188  *  frequency space are complex.
 189  *
 190  *  \return status - 0 or a standard error message.
 191  *
 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.
 195  */
 196 int gmx_fft_init_2d_real(gmx_fft_t* fft, int nx, int ny, gmx_fft_flag flags);
 197
 198
 199 /*! \brief Perform a 1-dimensional complex-to-complex transform
 200  *
 201  *  Performs an instance of a transform previously initiated.
 202  *
 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.
 211  *
 212  * \return 0 on success, or an error code.
 213  *
 214  * \note Data pointers are declared as void, to avoid casting pointers
 215  *       depending on your grid type.
 216  */
 217 int gmx_fft_1d(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
 218
 219
 220 /*! \brief Perform many 1-dimensional complex-to-complex transforms
 221  *
 222  *  Performs many instances of a transform previously initiated.
 223  *
 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.
 232  *
 233  * \return 0 on success, or an error code.
 234  *
 235  * \note Data pointers are declared as void, to avoid casting pointers
 236  *       depending on your grid type.
 237  */
 238 int gmx_fft_many_1d(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
 239
 240
 241 /*! \brief Perform a 1-dimensional real-to-complex transform
 242  *
 243  *  Performs an instance of a transform previously initiated.
 244  *
 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.
 253  *
 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).
 257  *
 258  * \return 0 on success, or an error code.
 259  *
 260  * \note Data pointers are declared as void, to avoid casting pointers
 261  *       depending on transform direction.
 262  */
 263 int gmx_fft_1d_real(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
 264
 265 /*! \brief Perform many 1-dimensional real-to-complex transforms
 266  *
 267  *  Performs many instances of a transform previously initiated.
 268  *
 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.
 277  *
 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).
 281  *
 282  * \return 0 on success, or an error code.
 283  *
 284  * \note Data pointers are declared as void, to avoid casting pointers
 285  *       depending on transform direction.
 286  */
 287 int gmx_fft_many_1d_real(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
 288
 289 /*! \brief Perform a 2-dimensional real-to-complex transform
 290  *
 291  *  Performs an instance of a transform previously initiated.
 292  *
 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.
 301  *
 302  * \return 0 on success, or an error code.
 303  *
 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.
 311  *
 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.
 314  *
 315  * \note Data pointers are declared as void, to avoid casting pointers
 316  *       depending on transform direction.
 317  */
 318 int gmx_fft_2d_real(gmx_fft_t setup, enum gmx_fft_direction dir, void* in_data, void* out_data);
 319
 320 /*! \brief Release an FFT setup structure
 321  *
 322  *  Destroy setup and release all allocated memory.
 323  *
 324  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 325  *               of the other initializers.
 326  *
 327  */
 328 void gmx_fft_destroy(gmx_fft_t setup);
 329
 330 /*! \brief Release a many FFT setup structure
 331  *
 332  *  Destroy setup and release all allocated memory.
 333  *
 334  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 335  *               of the other initializers.
 336  *
 337  */
 338 void gmx_many_fft_destroy(gmx_fft_t setup);
 339
 340
 341 /*! \brief Transpose 2d complex matrix, in-place or out-of-place.
 342  *
 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.
 346  *
 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
 352  *
 353  * \return GMX_SUCCESS, or an error code from gmx_errno.h
 354  */
 355 int gmx_fft_transpose_2d(t_complex* in_data, t_complex* out_data, int nx, int ny);
 356
 357 /*! \brief Cleanup global data of FFT
 358  *
 359  *  Any plans are invalid after this function. Should be called
 360  *  after all plans have been destroyed.
 361  */
 362 void gmx_fft_cleanup();
 363
 364 #endif