include/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 static const int GMX_FFT_FLAG_NONE = 0;
 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
 188 /*! \brief Setup a 2-dimensional complex-to-complex transform
 189  *
 190  *  \param fft    Pointer to opaque Gromacs FFT datatype
 191  *  \param nx     Length of transform in first dimension
 192  *  \param ny     Length of transform in second dimension
 193  *  \param flags  FFT options
 194  *
 195  *  \return status - 0 or a standard error message.
 196  *
 197  *  \note Since some of the libraries (e.g. MKL) store work array data in their
 198  *        handles this datatype should only be used for one thread at a time,
 199  *        i.e. you should create one copy per thread when executing in parallel.
 200  */
 201 int
 202 gmx_fft_init_2d        (gmx_fft_t *         fft,
 203                         int                 nx,
 204                         int                 ny,
 205                         gmx_fft_flag        flags);
 206
 207
 208 /*! \brief Setup a 2-dimensional real-to-complex transform
 209  *
 210  *  \param fft    Pointer to opaque Gromacs FFT datatype
 211  *  \param nx     Length of transform in first dimension
 212  *  \param ny     Length of transform in second dimension
 213  *  \param flags  FFT options
 214  *
 215  *  The normal space is assumed to be real, while the values in
 216  *  frequency space are complex.
 217  *
 218  *  \return status - 0 or a standard error message.
 219  *
 220  *  \note Since some of the libraries (e.g. MKL) store work array data in their
 221  *        handles this datatype should only be used for one thread at a time,
 222  *        i.e. you should create one copy per thread when executing in parallel.
 223  */
 224 int
 225 gmx_fft_init_2d_real        (gmx_fft_t *         fft,
 226                              int                 nx,
 227                              int                 ny,
 228                              gmx_fft_flag        flags);
 229
 230
 231 /*! \brief Setup a 3-dimensional complex-to-complex transform
 232  *
 233  *  \param fft    Pointer to opaque Gromacs FFT datatype
 234  *  \param nx     Length of transform in first dimension
 235  *  \param ny     Length of transform in second dimension
 236  *  \param nz     Length of transform in third dimension
 237  *  \param flags  FFT options
 238  *
 239  *  \return status - 0 or a standard error message.
 240  *
 241  *  \note Since some of the libraries (e.g. MKL) store work array data in their
 242  *        handles this datatype should only be used for one thread at a time,
 243  *        i.e. you should create one copy per thread when executing in parallel.
 244  */
 245 int
 246 gmx_fft_init_3d        (gmx_fft_t *         fft,
 247                         int                 nx,
 248                         int                 ny,
 249                         int                 nz,
 250                         gmx_fft_flag   flags);
 251
 252
 253 /*! \brief Setup a 3-dimensional real-to-complex transform
 254  *
 255  *  \param fft    Pointer to opaque Gromacs FFT datatype
 256  *  \param nx     Length of transform in first dimension
 257  *  \param ny     Length of transform in second dimension
 258  *  \param nz     Length of transform in third dimension
 259  *  \param flags  FFT options
 260  *
 261  *  The normal space is assumed to be real, while the values in
 262  *  frequency space are complex.
 263  *
 264  *  \return status - 0 or a standard error message.
 265  *
 266  *  \note Since some of the libraries (e.g. MKL) store work array data in their
 267  *        handles this datatype should only be used for one thread at a time,
 268  *        i.e. you should create one copy per thread when executing in parallel.
 269  */
 270 int
 271 gmx_fft_init_3d_real   (gmx_fft_t *         fft,
 272                         int                 nx,
 273                         int                 ny,
 274                         int                 nz,
 275                         gmx_fft_flag   flags);
 276
 277
 278
 279 /*! \brief Perform a 1-dimensional complex-to-complex transform
 280  *
 281  *  Performs an instance of a transform previously initiated.
 282  *
 283  *  \param setup     Setup returned from gmx_fft_init_1d()
 284  *  \param dir       Forward or Backward
 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  * \return 0 on success, or an error code.
 293  *
 294  * \note Data pointers are declared as void, to avoid casting pointers
 295  *       depending on your grid type.
 296  */
 297  int
 298 gmx_fft_1d               (gmx_fft_t                  setup,
 299                           enum gmx_fft_direction     dir,
 300                           void *                     in_data,
 301                           void *                     out_data);
 302
 303
 304 /*! \brief Perform many 1-dimensional complex-to-complex transforms
 305  *
 306  *  Performs many instances of a transform previously initiated.
 307  *
 308  *  \param setup     Setup returned from gmx_fft_init_1d()
 309  *  \param dir       Forward or Backward
 310  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 311  *                   to make it 16-byte aligned for better performance.
 312  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 313  *                   to make it 16-byte aligned for better performance.
 314  *                   You can provide the same pointer for in_data and out_data
 315  *                   to perform an in-place transform.
 316  *
 317  * \return 0 on success, or an error code.
 318  *
 319  * \note Data pointers are declared as void, to avoid casting pointers
 320  *       depending on your grid type.
 321  */
 322  int
 323 gmx_fft_many_1d          (gmx_fft_t                  setup,
 324                           enum gmx_fft_direction     dir,
 325                           void *                     in_data,
 326                           void *                     out_data);
 327
 328
 329 /*! \brief Perform a 1-dimensional real-to-complex transform
 330  *
 331  *  Performs an instance of a transform previously initiated.
 332  *
 333  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 334  *  \param dir       Real-to-complex or complex-to-real
 335  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 336  *                   to make it 16-byte aligned for better performance.
 337  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 338  *                   to make it 16-byte aligned for better performance.
 339  *                   You can provide the same pointer for in_data and out_data
 340  *                   to perform an in-place transform.
 341  *
 342  * If you are doing an in-place transform, the array must be padded up to
 343  * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
 344  * should not be padded (although it doesn't matter in 1d).
 345  *
 346  * \return 0 on success, or an error code.
 347  *
 348  * \note Data pointers are declared as void, to avoid casting pointers
 349  *       depending on transform direction.
 350  */
 351 int
 352 gmx_fft_1d_real          (gmx_fft_t                  setup,
 353                           enum gmx_fft_direction     dir,
 354                           void *                     in_data,
 355                           void *                     out_data);
 356
 357 /*! \brief Perform many 1-dimensional real-to-complex transforms
 358  *
 359  *  Performs many instances of a transform previously initiated.
 360  *
 361  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 362  *  \param dir       Real-to-complex or complex-to-real
 363  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 364  *                   to make it 16-byte aligned for better performance.
 365  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 366  *                   to make it 16-byte aligned for better performance.
 367  *                   You can provide the same pointer for in_data and out_data
 368  *                   to perform an in-place transform.
 369  *
 370  * If you are doing an in-place transform, the array must be padded up to
 371  * an even integer length so n/2 complex numbers can fit. Out-of-place arrays
 372  * should not be padded (although it doesn't matter in 1d).
 373  *
 374  * \return 0 on success, or an error code.
 375  *
 376  * \note Data pointers are declared as void, to avoid casting pointers
 377  *       depending on transform direction.
 378  */
 379 int
 380 gmx_fft_many_1d_real     (gmx_fft_t                  setup,
 381                           enum gmx_fft_direction     dir,
 382                           void *                     in_data,
 383                           void *                     out_data);
 384
 385
 386 /*! \brief Perform a 2-dimensional complex-to-complex transform
 387  *
 388  *  Performs an instance of a transform previously initiated.
 389  *
 390  *  \param setup     Setup returned from gmx_fft_init_1d()
 391  *  \param dir       Forward or Backward
 392  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 393  *                   to make it 16-byte aligned for better performance.
 394  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 395  *                   to make it 16-byte aligned for better performance.
 396  *                   You can provide the same pointer for in_data and out_data
 397  *                   to perform an in-place transform.
 398  *
 399  * \return 0 on success, or an error code.
 400  *
 401  * \note Data pointers are declared as void, to avoid casting pointers
 402  *       depending on your grid type.
 403  */
 404 int
 405 gmx_fft_2d               (gmx_fft_t                  setup,
 406                           enum gmx_fft_direction     dir,
 407                           void *                     in_data,
 408                           void *                     out_data);
 409
 410
 411 /*! \brief Perform a 2-dimensional real-to-complex transform
 412  *
 413  *  Performs an instance of a transform previously initiated.
 414  *
 415  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 416  *  \param dir       Real-to-complex or complex-to-real
 417  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 418  *                   to make it 16-byte aligned for better performance.
 419  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 420  *                   to make it 16-byte aligned for better performance.
 421  *                   You can provide the same pointer for in_data and out_data
 422  *                   to perform an in-place transform.
 423  *
 424  * \return 0 on success, or an error code.
 425  *
 426  * \note If you are doing an in-place transform, the last dimension of the
 427  * array MUST be padded up to an even integer length so n/2 complex numbers can
 428  * fit. Thus, if the real grid e.g. has dimension 5*3, you must allocate it as
 429  * a 5*4 array, where the last element in the second dimension is padding.
 430  * The complex data will be written to the same array, but since that dimension
 431  * is 5*2 it will now fill the entire array. Reverse complex-to-real in-place
 432  * transformation will produce the same sort of padded array.
 433  *
 434  * The padding does NOT apply to out-of-place transformation. In that case the
 435  * input array will simply be 5*3 of real, while the output is 5*2 of complex.
 436  *
 437  * \note Data pointers are declared as void, to avoid casting pointers
 438  *       depending on transform direction.
 439  */
 440 int
 441 gmx_fft_2d_real          (gmx_fft_t                  setup,
 442                           enum gmx_fft_direction     dir,
 443                           void *                     in_data,
 444                           void *                     out_data);
 445
 446
 447 /*! \brief Perform a 3-dimensional complex-to-complex transform
 448  *
 449  *  Performs an instance of a transform previously initiated.
 450  *
 451  *  \param setup     Setup returned from gmx_fft_init_1d()
 452  *  \param dir       Forward or Backward
 453  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 454  *                   to make it 16-byte aligned for better performance.
 455  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 456  *                   to make it 16-byte aligned for better performance.
 457  *                   You can provide the same pointer for in_data and out_data
 458  *                   to perform an in-place transform.
 459  *
 460  * \return 0 on success, or an error code.
 461  *
 462  * \note Data pointers are declared as void, to avoid casting pointers
 463  *       depending on your grid type.
 464  */
 465 int
 466 gmx_fft_3d               (gmx_fft_t                  setup,
 467                           enum gmx_fft_direction     dir,
 468                           void *                     in_data,
 469                           void *                     out_data);
 470
 471
 472 /*! \brief Perform a 3-dimensional real-to-complex transform
 473  *
 474  *  Performs an instance of a transform previously initiated.
 475  *
 476  *  \param setup     Setup returned from gmx_fft_init_1d_real()
 477  *  \param dir       Real-to-complex or complex-to-real
 478  *  \param in_data   Input grid data. This should be allocated with gmx_new()
 479  *                   to make it 16-byte aligned for better performance.
 480  *  \param out_data  Output grid data. This should be allocated with gmx_new()
 481  *                   to make it 16-byte aligned for better performance.
 482  *                   You can provide the same pointer for in_data and out_data
 483  *                   to perform an in-place transform.
 484  *
 485  * \return 0 on success, or an error code.
 486  *
 487  * \note If you are doing an in-place transform, the last dimension of the
 488  * array MUST be padded up to an even integer length so n/2 complex numbers can
 489  * fit. Thus, if the real grid e.g. has dimension 7*5*3, you must allocate it as
 490  * a 7*5*4 array, where the last element in the second dimension is padding.
 491  * The complex data will be written to the same array, but since that dimension
 492  * is 7*5*2 it will now fill the entire array. Reverse complex-to-real in-place
 493  * transformation will produce the same sort of padded array.
 494  *
 495  * The padding does NOT apply to out-of-place transformation. In that case the
 496  * input will simply be 7*5*3 of real, while the output is 7*5*2 of complex.
 497  *
 498  * \note Data pointers are declared as void, to avoid casting pointers
 499  *       depending on transform direction.
 500  */
 501 int
 502 gmx_fft_3d_real          (gmx_fft_t                  setup,
 503                           enum gmx_fft_direction     dir,
 504                           void *                     in_data,
 505                           void *                     out_data);
 506
 507
 508 /*! \brief Release an FFT setup structure
 509  *
 510  *  Destroy setup and release all allocated memory.
 511  *
 512  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 513  *               of the other initializers.
 514  *
 515  */
 516 void
 517 gmx_fft_destroy          (gmx_fft_t                 setup);
 518
 519 /*! \brief Release a many FFT setup structure
 520  *
 521  *  Destroy setup and release all allocated memory.
 522  *
 523  *  \param setup Setup returned from gmx_fft_init_1d(), or one
 524  *               of the other initializers.
 525  *
 526  */
 527 void
 528 gmx_many_fft_destroy          (gmx_fft_t                 setup);
 529
 530
 531 /*! \brief Transpose 2d complex matrix, in-place or out-of-place.
 532  *
 533  * This routines works when the matrix is non-square, i.e. nx!=ny too,
 534  * without allocating an entire matrix of work memory, which is important
 535  * for huge FFT grids.
 536  *
 537  * \param in_data    Input data, to be transposed
 538  * \param out_data   Output, transposed data. If this is identical to
 539  *                   in_data, an in-place transpose is performed.
 540  * \param nx         Number of rows before transpose
 541  * \param ny         Number of columns before transpose
 542  *
 543  * \return GMX_SUCCESS, or an error code from gmx_errno.h
 544  */
 545 int
 546 gmx_fft_transpose_2d   (t_complex *       in_data,
 547                         t_complex *       out_data,
 548                         int               nx,
 549                         int               ny);
 550
 551
 552 /*! \brief Transpose 2d multi-element matrix
 553  *
 554  * This routine is very similar to gmx_fft_transpose_2d(), but it
 555  * supports matrices with more than one data value for each position.
 556  * It is extremely useful when transposing the x/y dimensions of a 3d
 557  * matrix - in that case you just set nelem to nz, and the routine will do
 558  * and x/y transpose where it moves entire columns of z data
 559  *
 560  * This routines works when the matrix is non-square, i.e. nx!=ny too,
 561  * without allocating an entire matrix of work memory, which is important
 562  * for huge FFT grid.
 563  *
 564  * For performance reasons you need to provide a \a small workarray
 565  * with length at least 2*nelem (note that the type is char, not t_complex).
 566  *
 567  * \param in_data    Input data, to be transposed
 568  * \param out_data   Output, transposed data. If this is identical to
 569  *                   in_data, an in-place transpose is performed.
 570  * \param nx         Number of rows before transpose
 571  * \param ny         Number of columns before transpose
 572  * \param nelem      Number of t_complex values in each position. If this
 573  *                   is 1 it is faster to use gmx_fft_transpose_2d() directly.
 574  * \param work       Work array of length 2*nelem, type t_complex.
 575  *
 576  * \return GMX_SUCCESS, or an error code from gmx_errno.h
 577  */
 578 int
 579 gmx_fft_transpose_2d_nelem(t_complex *        in_data,
 580                            t_complex *        out_data,
 581                            int                nx,
 582                            int                ny,
 583                            int                nelem,
 584                            t_complex *        work);
 585
 586
 587
 588
 589 #ifdef __cplusplus
 590 }
 591 #endif
 592
 593 #endif /* _GMX_FFT_H_ */