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, 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 #ifdef __cplusplus
  61 extern "C" {
  62 #endif
  63 #if 0
  64 } /* fixes auto-indentation problems */
  65 #endif
  66
  67
  68
  69 /*! \brief Datatype for FFT setup
  70  *
  71  *  The gmx_fft_t type contains all the setup information, e.g. twiddle
  72  *  factors, necessary to perform an FFT. Internally it is mapped to
  73  *  whatever FFT library we are using, or the built-in FFTPACK if no fast
  74  *  external library is available.
  75  *
  76  *  Since some of the libraries (e.g. MKL) store work array data in their
  77  *  handles this datatype should only be used for one thread at a time, i.e.
  78  *  they should allocate one instance each when executing in parallel.
  79  */
  80 typedef struct gmx_fft *
  81     gmx_fft_t;
  82
  83
  84
  85
  86 /*! \brief Specifier for FFT direction.
  87  *
  88  *  The definition of the 1D forward transform from input x[] to output y[] is
  89  *  \f[
  90  *  y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{-i 2 \pi j k /N}
  91  *  \f]
  92  *
  93  *  while the corresponding backward transform is
  94  *
  95  *  \f[
  96  *  y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{i 2 \pi j k /N}
  97  *  \f]
  98  *
  99  *  A forward-backward transform pair will this result in data scaled by N.
 100  *
 101  *  For complex-to-complex transforms you can only use one of
 102  *  GMX_FFT_FORWARD or GMX_FFT_BACKWARD, and for real-complex transforms you
 103  *  can only use GMX_FFT_REAL_TO_COMPLEX or GMX_FFT_COMPLEX_TO_REAL.
 104  */
 105 typedef enum gmx_fft_direction
 106 {
 107     GMX_FFT_FORWARD,         /**< Forward complex-to-complex transform  */
 108     GMX_FFT_BACKWARD,        /**< Backward complex-to-complex transform */
 109     GMX_FFT_REAL_TO_COMPLEX, /**< Real-to-complex valued FFT            */
 110     GMX_FFT_COMPLEX_TO_REAL  /**< Complex-to-real valued FFT            */
 111 } gmx_fft_direction;
 112
 113 /*! \brief Specifier for FFT flags.
 114  *
 115  *  Some FFT libraries (FFTW, in particular) can do timings and other
 116  *  tricks to try and optimize the FFT for the current architecture. However,
 117  *  this can also lead to results that differ between consecutive runs with
 118  *  identical input.
 119  *  To avoid this, the conservative flag will attempt to disable such
 120  *  optimization, but there are no guarantees since we cannot control what
 121  *  the FFT libraries do internally.
 122  */
 123
 124 typedef int gmx_fft_flag;
 125 /** Macro to indicate no special flags for FFT routines. */
 126 static const int GMX_FFT_FLAG_NONE = 0;
 127 /** Flag to disable FFT optimizations based on timings, see ::gmx_fft_flag. */
 128 static const int GMX_FFT_FLAG_CONSERVATIVE = (1<<0);
 129
 130 /*! \brief Setup a 1-dimensional complex-to-complex transform
 131  *
 132  *  \param fft    Pointer to opaque Gromacs FFT datatype
 133  *  \param nx     Length of transform
 134  *  \param flags  FFT options
 135  *
 136  *  \return status - 0 or a standard error message.
 137  *
 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.
 141  */
 142 int
 143 gmx_fft_init_1d        (gmx_fft_t *       fft,
 144                         int               nx,
 145                         gmx_fft_flag      flags);
 146
 147
 148 /*! \brief Setup multiple 1-dimensional complex-to-complex transform
 149  *
 150  *  \param fft    Pointer to opaque Gromacs FFT datatype
 151  *  \param nx     Length of transform
 152  *  \param howmany Howmany 1D FFT
 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
 162 gmx_fft_init_many_1d        (gmx_fft_t *       fft,
 163                              int               nx,
 164                              int               howmany,
 165                              gmx_fft_flag      flags);
 166
 167
 168 /*! \brief Setup a 1-dimensional real-to-complex transform
 169  *
 170  *  \param fft    Pointer to opaque Gromacs FFT datatype
 171  *  \param nx     Length of transform in real space
 172  *  \param flags  FFT options
 173  *
 174  *  \return status - 0 or a standard error message.
 175  *
 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.
 179  */
 180 int
 181 gmx_fft_init_1d_real        (gmx_fft_t *       fft,
 182                              int               nx,
 183                              gmx_fft_flag      flags);
 184
 185
 186 /*! \brief Setup multiple 1-dimensional real-to-complex transform
 187  *
 188  *  \param fft    Pointer to opaque Gromacs FFT datatype
 189  *  \param nx     Length of transform in real space
 190  *  \param howmany Homany 1D FFTs
 191  *  \param flags  FFT options
 192  *
 193  *  \return status - 0 or a standard error message.
 194  *
 195  *  \note Since some of the libraries (e.g. MKL) store work array data in their
 196  *        handles this datatype should only be used for one thread at a time,
 197  *        i.e. you should create one copy per thread when executing in parallel.
 198  */
 199 int
 200 gmx_fft_init_many_1d_real        (gmx_fft_t *       fft,
 201                                   int               nx,
 202                                   int               howmany,
 203                                   gmx_fft_flag      flags);
 204
 205
 206 /*! \brief Setup a 2-dimensional real-to-complex transform
 207  *
 208  *  \param fft    Pointer to opaque Gromacs FFT datatype
 209  *  \param nx     Length of transform in first dimension
 210  *  \param ny     Length of transform in second dimension
 211  *  \param flags  FFT options
 212  *
 213  *  The normal space is assumed to be real, while the values in
 214  *  frequency space are complex.
 215  *
 216  *  \return status - 0 or a standard error message.
 217  *
 218  *  \note Since some of the libraries (e.g. MKL) store work array data in their
 219  *        handles this datatype should only be used for one thread at a time,
 220  *        i.e. you should create one copy per thread when executing in parallel.
 221  */
 222 int
 223 gmx_fft_init_2d_real        (gmx_fft_t *         fft,
 224                              int                 nx,
 225                              int                 ny,
 226                              gmx_fft_flag        flags);
 227
 228
 229 /*! \brief Perform a 1-dimensional complex-to-complex transform
 230  *
 231  *  Performs an instance of a transform previously initiated.
 232  *
 233  *  \param setup     Setup returned from gmx_fft_init_1d()
 234  *  \param dir       Forward or Backward
 235  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 236  *                   to make it 16-byte aligned for better performance.
 237  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 238  *                   to make it 16-byte aligned for better performance.
 239  *                   You can provide the same pointer for in_data and out_data
 240  *                   to perform an in-place transform.
 241  *
 242  * \return 0 on success, or an error code.
 243  *
 244  * \note Data pointers are declared as void, to avoid casting pointers
 245  *       depending on your grid type.
 246  */
 247 int
 248 gmx_fft_1d               (gmx_fft_t                  setup,
 249                           enum gmx_fft_direction     dir,
 250                           void *                     in_data,
 251                           void *                     out_data);
 252
 253
 254 /*! \brief Perform many 1-dimensional complex-to-complex transforms
 255  *
 256  *  Performs many instances of a transform previously initiated.
 257  *
 258  *  \param setup     Setup returned from gmx_fft_init_1d()
 259  *  \param dir       Forward or Backward
 260  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 261  *                   to make it 16-byte aligned for better performance.
 262  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 263  *                   to make it 16-byte aligned for better performance.
 264  *                   You can provide the same pointer for in_data and out_data
 265  *                   to perform an in-place transform.
 266  *
 267  * \return 0 on success, or an error code.
 268  *
 269  * \note Data pointers are declared as void, to avoid casting pointers
 270  *       depending on your grid type.
 271  */
 272 int
 273 gmx_fft_many_1d          (gmx_fft_t                  setup,
 274                           enum gmx_fft_direction     dir,
 275                           void *                     in_data,
 276                           void *                     out_data);
 277
 278
 279 /*! \brief Perform a 1-dimensional real-to-complex transform
 280  *
 281  *  Performs an instance of a transform previously initiated.
 282  *
 283  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 284  *  \param dir       Real-to-complex or complex-to-real
 285  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 286  *                   to make it 16-byte aligned for better performance.
 287  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 288  *                   to make it 16-byte aligned for better performance.
 289  *                   You can provide the same pointer for in_data and out_data
 290  *                   to perform an in-place transform.
 291  *
 292  * If you are doing an in-place transform, the array must be padded up to
 293  * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
 294  * should not be padded (although it doesn't matter in 1d).
 295  *
 296  * \return 0 on success, or an error code.
 297  *
 298  * \note Data pointers are declared as void, to avoid casting pointers
 299  *       depending on transform direction.
 300  */
 301 int
 302 gmx_fft_1d_real          (gmx_fft_t                  setup,
 303                           enum gmx_fft_direction     dir,
 304                           void *                     in_data,
 305                           void *                     out_data);
 306
 307 /*! \brief Perform many 1-dimensional real-to-complex transforms
 308  *
 309  *  Performs many instances of a transform previously initiated.
 310  *
 311  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 312  *  \param dir       Real-to-complex or complex-to-real
 313  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 314  *                   to make it 16-byte aligned for better performance.
 315  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 316  *                   to make it 16-byte aligned for better performance.
 317  *                   You can provide the same pointer for in_data and out_data
 318  *                   to perform an in-place transform.
 319  *
 320  * If you are doing an in-place transform, the array must be padded up to
 321  * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
 322  * should not be padded (although it doesn't matter in 1d).
 323  *
 324  * \return 0 on success, or an error code.
 325  *
 326  * \note Data pointers are declared as void, to avoid casting pointers
 327  *       depending on transform direction.
 328  */
 329 int
 330 gmx_fft_many_1d_real     (gmx_fft_t                  setup,
 331                           enum gmx_fft_direction     dir,
 332                           void *                     in_data,
 333                           void *                     out_data);
 334
 335 /*! \brief Perform a 2-dimensional real-to-complex transform
 336  *
 337  *  Performs an instance of a transform previously initiated.
 338  *
 339  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 340  *  \param dir       Real-to-complex or complex-to-real
 341  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 342  *                   to make it 16-byte aligned for better performance.
 343  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 344  *                   to make it 16-byte aligned for better performance.
 345  *                   You can provide the same pointer for in_data and out_data
 346  *                   to perform an in-place transform.
 347  *
 348  * \return 0 on success, or an error code.
 349  *
 350  * \note If you are doing an in-place transform, the last dimension of the
 351  * array MUST be padded up to an even integer length so n/2 complex numbers can
 352  * fit. Thus, if the real grid e.g. has dimension 5*3, you must allocate it as
 353  * a 5*4 array, where the last element in the second dimension is padding.
 354  * The complex data will be written to the same array, but since that dimension
 355  * is 5*2 it will now fill the entire array. Reverse complex-to-real in-place
 356  * transformation will produce the same sort of padded array.
 357  *
 358  * The padding does NOT apply to out-of-place transformation. In that case the
 359  * input array will simply be 5*3 of real, while the output is 5*2 of complex.
 360  *
 361  * \note Data pointers are declared as void, to avoid casting pointers
 362  *       depending on transform direction.
 363  */
 364 int
 365 gmx_fft_2d_real          (gmx_fft_t                  setup,
 366                           enum gmx_fft_direction     dir,
 367                           void *                     in_data,
 368                           void *                     out_data);
 369
 370 /*! \brief Release an FFT setup structure
 371  *
 372  *  Destroy setup and release all allocated memory.
 373  *
 374  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 375  *               of the other initializers.
 376  *
 377  */
 378 void
 379 gmx_fft_destroy          (gmx_fft_t                 setup);
 380
 381 /*! \brief Release a many FFT setup structure
 382  *
 383  *  Destroy setup and release all allocated memory.
 384  *
 385  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 386  *               of the other initializers.
 387  *
 388  */
 389 void
 390 gmx_many_fft_destroy          (gmx_fft_t                 setup);
 391
 392
 393 /*! \brief Transpose 2d complex matrix, in-place or out-of-place.
 394  *
 395  * This routines works when the matrix is non-square, i.e. nx!=ny too,
 396  * without allocating an entire matrix of work memory, which is important
 397  * for huge FFT grids.
 398  *
 399  * \param in_data    Input data, to be transposed
 400  * \param out_data   Output, transposed data. If this is identical to
 401  *                   in_data, an in-place transpose is performed.
 402  * \param nx         Number of rows before transpose
 403  * \param ny         Number of columns before transpose
 404  *
 405  * \return GMX_SUCCESS, or an error code from gmx_errno.h
 406  */
 407 int
 408 gmx_fft_transpose_2d   (t_complex *       in_data,
 409                         t_complex *       out_data,
 410                         int               nx,
 411                         int               ny);
 412
 413 /*! \brief Cleanup global data of FFT
 414  *
 415  *  Any plans are invalid after this function. Should be called
 416  *  after all plans have been destroyed.
 417  */
 418 void gmx_fft_cleanup();
 419
 420 /*! \brief Return string describing the underlying FFT implementation.
 421  *
 422  * Used to print out information about the used FFT library where needed.
 423  */
 424 const char *gmx_fft_get_version_info();
 425
 426 #ifdef __cplusplus
 427 }
 428 #endif
 429
 430 #endif