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