src/gromacs/fft/fft.h

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * Gromacs 4.0                         Copyright (c) 1991-2003
   5  * David van der Spoel, Erik Lindahl, University of Groningen.
   6  *
   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.
  11  *
  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
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18 /*! \file
  19  *  \brief Fast Fourier Transforms.
  20  *
  21  *  This file provides an abstract Gromacs interface to Fourier transforms,
  22  *  including multi-dimensional and real-to-complex transforms.
  23  *
  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.
  27  *
  28  *  We also provide our own multi-dimensional transform setups even when
  29  *  the underlying library does not support it directly.
  30  *
  31  * \ingroup module_fft
  32  */
  33 #ifndef GMX_FFT_FFT_H
  34 #define GMX_FFT_FFT_H
  35
  36 #include <stdio.h>
  37
  38 #include "../legacyheaders/types/simple.h"
  39 #include "../legacyheaders/gmxcomplex.h"
  40
  41 #ifdef __cplusplus
  42 extern "C" {
  43 #endif
  44 #if 0
  45 } /* fixes auto-indentation problems */
  46 #endif
  47
  48
  49
  50 /*! \brief Datatype for FFT setup
  51  *
  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.
  56  *
  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.
  60  */
  61 typedef struct gmx_fft *
  62     gmx_fft_t;
  63
  64
  65
  66
  67 /*! \brief Specifier for FFT direction.
  68  *
  69  *  The definition of the 1D forward transform from input x[] to output y[] is
  70  *  \f[
  71  *  y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{-i 2 \pi j k /N}
  72  *  \f]
  73  *
  74  *  while the corresponding backward transform is
  75  *
  76  *  \f[
  77  *  y_{k} = \sum_{j=0}^{N-1} x_{j} \exp{i 2 \pi j k /N}
  78  *  \f]
  79  *
  80  *  A forward-backward transform pair will this result in data scaled by N.
  81  *
  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.
  85  */
  86 typedef enum gmx_fft_direction
  87 {
  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            */
  92 } gmx_fft_direction;
  93
  94 /*! \brief Specifier for FFT flags.
  95  *
  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
  99  *  identical input.
 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.
 103  */
 104
 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);
 110
 111 /*! \brief Setup a 1-dimensional complex-to-complex transform
 112  *
 113  *  \param fft    Pointer to opaque Gromacs FFT datatype
 114  *  \param nx     Length of transform
 115  *  \param flags  FFT options
 116  *
 117  *  \return status - 0 or a standard error message.
 118  *
 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.
 122  */
 123 int
 124 gmx_fft_init_1d        (gmx_fft_t *       fft,
 125                         int               nx,
 126                         gmx_fft_flag      flags);
 127
 128
 129 /*! \brief Setup multiple 1-dimensional complex-to-complex transform
 130  *
 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
 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_many_1d        (gmx_fft_t *       fft,
 144                              int               nx,
 145                              int               howmany,
 146                              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
 162 gmx_fft_init_1d_real        (gmx_fft_t *       fft,
 163                              int               nx,
 164                              gmx_fft_flag      flags);
 165
 166
 167 /*! \brief Setup multiple 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 howmany Homany 1D FFTs
 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_many_1d_real        (gmx_fft_t *       fft,
 182                                   int               nx,
 183                                   int               howmany,
 184                                   gmx_fft_flag      flags);
 185
 186
 187 /*! \brief Setup a 2-dimensional real-to-complex transform
 188  *
 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
 193  *
 194  *  The normal space is assumed to be real, while the values in
 195  *  frequency space are complex.
 196  *
 197  *  \return status - 0 or a standard error message.
 198  *
 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.
 202  */
 203 int
 204 gmx_fft_init_2d_real        (gmx_fft_t *         fft,
 205                              int                 nx,
 206                              int                 ny,
 207                              gmx_fft_flag        flags);
 208
 209
 210 /*! \brief Perform a 1-dimensional complex-to-complex transform
 211  *
 212  *  Performs an instance of a transform previously initiated.
 213  *
 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.
 222  *
 223  * \return 0 on success, or an error code.
 224  *
 225  * \note Data pointers are declared as void, to avoid casting pointers
 226  *       depending on your grid type.
 227  */
 228 int
 229 gmx_fft_1d               (gmx_fft_t                  setup,
 230                           enum gmx_fft_direction     dir,
 231                           void *                     in_data,
 232                           void *                     out_data);
 233
 234
 235 /*! \brief Perform many 1-dimensional complex-to-complex transforms
 236  *
 237  *  Performs many instances of a transform previously initiated.
 238  *
 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.
 247  *
 248  * \return 0 on success, or an error code.
 249  *
 250  * \note Data pointers are declared as void, to avoid casting pointers
 251  *       depending on your grid type.
 252  */
 253 int
 254 gmx_fft_many_1d          (gmx_fft_t                  setup,
 255                           enum gmx_fft_direction     dir,
 256                           void *                     in_data,
 257                           void *                     out_data);
 258
 259
 260 /*! \brief Perform a 1-dimensional real-to-complex transform
 261  *
 262  *  Performs an instance of a transform previously initiated.
 263  *
 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.
 272  *
 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).
 276  *
 277  * \return 0 on success, or an error code.
 278  *
 279  * \note Data pointers are declared as void, to avoid casting pointers
 280  *       depending on transform direction.
 281  */
 282 int
 283 gmx_fft_1d_real          (gmx_fft_t                  setup,
 284                           enum gmx_fft_direction     dir,
 285                           void *                     in_data,
 286                           void *                     out_data);
 287
 288 /*! \brief Perform many 1-dimensional real-to-complex transforms
 289  *
 290  *  Performs many instances of a transform previously initiated.
 291  *
 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.
 300  *
 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).
 304  *
 305  * \return 0 on success, or an error code.
 306  *
 307  * \note Data pointers are declared as void, to avoid casting pointers
 308  *       depending on transform direction.
 309  */
 310 int
 311 gmx_fft_many_1d_real     (gmx_fft_t                  setup,
 312                           enum gmx_fft_direction     dir,
 313                           void *                     in_data,
 314                           void *                     out_data);
 315
 316 /*! \brief Perform a 2-dimensional real-to-complex transform
 317  *
 318  *  Performs an instance of a transform previously initiated.
 319  *
 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.
 328  *
 329  * \return 0 on success, or an error code.
 330  *
 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.
 338  *
 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.
 341  *
 342  * \note Data pointers are declared as void, to avoid casting pointers
 343  *       depending on transform direction.
 344  */
 345 int
 346 gmx_fft_2d_real          (gmx_fft_t                  setup,
 347                           enum gmx_fft_direction     dir,
 348                           void *                     in_data,
 349                           void *                     out_data);
 350
 351 /*! \brief Release an FFT setup structure
 352  *
 353  *  Destroy setup and release all allocated memory.
 354  *
 355  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 356  *               of the other initializers.
 357  *
 358  */
 359 void
 360 gmx_fft_destroy          (gmx_fft_t                 setup);
 361
 362 /*! \brief Release a many FFT setup structure
 363  *
 364  *  Destroy setup and release all allocated memory.
 365  *
 366  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 367  *               of the other initializers.
 368  *
 369  */
 370 void
 371 gmx_many_fft_destroy          (gmx_fft_t                 setup);
 372
 373
 374 /*! \brief Transpose 2d complex matrix, in-place or out-of-place.
 375  *
 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.
 379  *
 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
 385  *
 386  * \return GMX_SUCCESS, or an error code from gmx_errno.h
 387  */
 388 int
 389 gmx_fft_transpose_2d   (t_complex *       in_data,
 390                         t_complex *       out_data,
 391                         int               nx,
 392                         int               ny);
 393
 394 /*! \brief Cleanup global data of FFT
 395  *
 396  *  Any plans are invalid after this function. Should be called
 397  *  after all plans have been destroyed.
 398  */
 399 void gmx_fft_cleanup();
 400
 401 /*! \brief Return string describing the underlying FFT implementation.
 402  *
 403  * Used to print out information about the used FFT library where needed.
 404  */
 405 const char *gmx_fft_get_version_info();
 406
 407 #ifdef __cplusplus
 408 }
 409 #endif
 410
 411 #endif