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