src/gromacs/fft/fft_fftw3.cpp

   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 #include "config.h"
  37
  38 #include <errno.h>
  39 #include <stdlib.h>
  40
  41 #include <fftw3.h>
  42
  43 #include "gromacs/fft/fft.h"
  44 #include "gromacs/utility/fatalerror.h"
  45
  46 #ifdef GMX_DOUBLE
  47 #define FFTWPREFIX(name) fftw_ ## name
  48 #else
  49 #define FFTWPREFIX(name) fftwf_ ## name
  50 #endif
  51
  52 #include "thread_mpi/mutex.h"
  53 #include "gromacs/utility/exceptions.h"
  54
  55 /* none of the fftw3 calls, except execute(), are thread-safe, so
  56    we need to serialize them with this mutex. */
  57 static tMPI::mutex big_fftw_mutex;
  58 #define FFTW_LOCK try { big_fftw_mutex.lock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
  59 #define FFTW_UNLOCK try { big_fftw_mutex.unlock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
  60
  61 /* We assume here that aligned memory starts at multiple of 16 bytes and unaligned memory starts at multiple of 8 bytes. The later is guranteed for all malloc implementation.
  62    Consequesences:
  63    - It is not allowed to use these FFT plans from memory which doesn't have a starting address as a multiple of 8 bytes.
  64      This is OK as long as the memory directly comes from malloc and is not some subarray within alloated memory.
  65    - This has to be fixed if any future architecute requires memory to be aligned to multiples of 32 bytes.
  66  */
  67 /*! \internal
  68  * \brief
  69  * Contents of the FFTW3 fft datatype.
  70  *
  71  * Note that this is one of several possible implementations of gmx_fft_t.
  72  */
  73 #ifdef DOXYGEN
  74 struct gmx_fft_fftw3
  75 #else
  76 struct gmx_fft
  77 #endif
  78 {
  79     /*! \brief
  80      * FFTW plans.
  81      *
  82      * Three alternatives (unaligned/aligned, out-of-place/in-place, forward/backward)
  83      * results in 8 different FFTW plans. Keep track of them with 3 array indices:
  84      * first index:   0=unaligned, 1=aligned
  85      * second index:  0=out-of-place, 1=in-place
  86      * third index:   0=backward, 1=forward
  87      */
  88     FFTWPREFIX(plan)         plan[2][2][2];
  89     /** Used to catch user mistakes */
  90     int                      real_transform;
  91     /** Number of dimensions in the FFT */
  92     int                      ndim;
  93 };
  94
  95 int
  96 gmx_fft_init_1d(gmx_fft_t *        pfft,
  97                 int                nx,
  98                 gmx_fft_flag       flags)
  99 {
 100     return gmx_fft_init_many_1d(pfft, nx, 1, flags);
 101 }
 102
 103
 104 int
 105 gmx_fft_init_many_1d(gmx_fft_t *        pfft,
 106                      int                nx,
 107                      int                howmany,
 108                      gmx_fft_flag       flags)
 109 {
 110     gmx_fft_t              fft;
 111     FFTWPREFIX(complex)   *p1, *p2, *up1, *up2;
 112     size_t                 pc;
 113     int                    i, j, k;
 114     int                    fftw_flags;
 115
 116 #ifdef GMX_DISABLE_FFTW_MEASURE
 117     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 118 #endif
 119
 120     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 121
 122     if (pfft == NULL)
 123     {
 124         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 125         return EINVAL;
 126     }
 127     *pfft = NULL;
 128
 129     FFTW_LOCK;
 130     if ( (fft = (gmx_fft_t)FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 131     {
 132         FFTW_UNLOCK;
 133         return ENOMEM;
 134     }
 135
 136     /* allocate aligned, and extra memory to make it unaligned */
 137     p1  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
 138     if (p1 == NULL)
 139     {
 140         FFTWPREFIX(free)(fft);
 141         FFTW_UNLOCK;
 142         return ENOMEM;
 143     }
 144
 145     p2  = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
 146     if (p2 == NULL)
 147     {
 148         FFTWPREFIX(free)(p1);
 149         FFTWPREFIX(free)(fft);
 150         FFTW_UNLOCK;
 151         return ENOMEM;
 152     }
 153
 154     /* make unaligned pointers.
 155      * In double precision the actual complex datatype will be 16 bytes,
 156      * so go to a char pointer and force an offset of 8 bytes instead.
 157      */
 158     pc  = (size_t)p1;
 159     pc += 8;
 160     up1 = (FFTWPREFIX(complex) *)pc;
 161
 162     pc  = (size_t)p2;
 163     pc += 8;
 164     up2 = (FFTWPREFIX(complex) *)pc;
 165
 166     /*                            int rank, const int *n, int howmany,
 167                                   fftw_complex *in, const int *inembed,
 168                                   int istride, int idist,
 169                                   fftw_complex *out, const int *onembed,
 170                                   int ostride, int odist,
 171                                   int sign, unsigned flags */
 172     fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, up1, &nx, 1, nx, up2, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 173     fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, up1, &nx, 1, nx, up2, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 174     fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, up1, &nx, 1, nx, up1, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 175     fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, up1, &nx, 1, nx, up1, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 176     fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, p1, &nx, 1, nx, p2, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 177     fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, p1, &nx, 1, nx, p2, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 178     fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, p1, &nx, 1, nx, p1, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 179     fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1, &nx, howmany, p1, &nx, 1, nx, p1, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 180
 181     for (i = 0; i < 2; i++)
 182     {
 183         for (j = 0; j < 2; j++)
 184         {
 185             for (k = 0; k < 2; k++)
 186             {
 187                 if (fft->plan[i][j][k] == NULL)
 188                 {
 189                     gmx_fatal(FARGS, "Error initializing FFTW3 plan.");
 190                     FFTW_UNLOCK;
 191                     gmx_fft_destroy(fft);
 192                     FFTW_LOCK;
 193                     FFTWPREFIX(free)(p1);
 194                     FFTWPREFIX(free)(p2);
 195                     FFTW_UNLOCK;
 196                     return -1;
 197                 }
 198             }
 199         }
 200     }
 201
 202     FFTWPREFIX(free)(p1);
 203     FFTWPREFIX(free)(p2);
 204
 205     fft->real_transform = 0;
 206     fft->ndim           = 1;
 207
 208     *pfft = fft;
 209     FFTW_UNLOCK;
 210     return 0;
 211 }
 212
 213 int
 214 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 215                      int                nx,
 216                      gmx_fft_flag       flags)
 217 {
 218     return gmx_fft_init_many_1d_real(pfft, nx, 1, flags);
 219 }
 220
 221 int
 222 gmx_fft_init_many_1d_real(gmx_fft_t *        pfft,
 223                           int                nx,
 224                           int                howmany,
 225                           gmx_fft_flag       flags)
 226 {
 227     gmx_fft_t              fft;
 228     real                  *p1, *p2, *up1, *up2;
 229     size_t                 pc;
 230     int                    i, j, k;
 231     int                    fftw_flags;
 232
 233 #ifdef GMX_DISABLE_FFTW_MEASURE
 234     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 235 #endif
 236
 237     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 238
 239     if (pfft == NULL)
 240     {
 241         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 242         return EINVAL;
 243     }
 244     *pfft = NULL;
 245
 246     FFTW_LOCK;
 247     if ( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 248     {
 249         FFTW_UNLOCK;
 250         return ENOMEM;
 251     }
 252
 253     /* allocate aligned, and extra memory to make it unaligned */
 254     p1  = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
 255     if (p1 == NULL)
 256     {
 257         FFTWPREFIX(free)(fft);
 258         FFTW_UNLOCK;
 259         return ENOMEM;
 260     }
 261
 262     p2  = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
 263     if (p2 == NULL)
 264     {
 265         FFTWPREFIX(free)(p1);
 266         FFTWPREFIX(free)(fft);
 267         FFTW_UNLOCK;
 268         return ENOMEM;
 269     }
 270
 271     /* make unaligned pointers.
 272      * In double precision the actual complex datatype will be 16 bytes,
 273      * so go to a char pointer and force an offset of 8 bytes instead.
 274      */
 275     pc  = (size_t)p1;
 276     pc += 8;
 277     up1 = (real *)pc;
 278
 279     pc  = (size_t)p2;
 280     pc += 8;
 281     up2 = (real *)pc;
 282
 283     /*                                int rank, const int *n, int howmany,
 284                                       double *in, const int *inembed,
 285                                       int istride, int idist,
 286                                       fftw_complex *out, const int *onembed,
 287                                       int ostride, int odist,
 288                                       unsigned flag    */
 289     fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1, &nx, howmany, up1, 0, 1, (nx/2+1) *2, (FFTWPREFIX(complex) *) up2, 0, 1, (nx/2+1), fftw_flags);
 290     fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1, &nx, howmany, up1, 0, 1, (nx/2+1) *2, (FFTWPREFIX(complex) *) up1, 0, 1, (nx/2+1), fftw_flags);
 291     fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1, &nx, howmany, p1, 0, 1, (nx/2+1) *2, (FFTWPREFIX(complex) *) p2, 0, 1, (nx/2+1), fftw_flags);
 292     fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1, &nx, howmany, p1, 0, 1, (nx/2+1) *2, (FFTWPREFIX(complex) *) p1, 0, 1, (nx/2+1), fftw_flags);
 293
 294     fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1, &nx, howmany, (FFTWPREFIX(complex) *) up1, 0, 1, (nx/2+1), up2, 0, 1, (nx/2+1) *2, fftw_flags);
 295     fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1, &nx, howmany, (FFTWPREFIX(complex) *) up1, 0, 1, (nx/2+1), up1, 0, 1, (nx/2+1) *2, fftw_flags);
 296     fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1, &nx, howmany, (FFTWPREFIX(complex) *) p1, 0, 1, (nx/2+1), p2, 0, 1, (nx/2+1) *2, fftw_flags);
 297     fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1, &nx, howmany, (FFTWPREFIX(complex) *) p1, 0, 1, (nx/2+1), p1, 0, 1, (nx/2+1) *2, fftw_flags);
 298
 299     for (i = 0; i < 2; i++)
 300     {
 301         for (j = 0; j < 2; j++)
 302         {
 303             for (k = 0; k < 2; k++)
 304             {
 305                 if (fft->plan[i][j][k] == NULL)
 306                 {
 307                     gmx_fatal(FARGS, "Error initializing FFTW3 plan.");
 308                     FFTW_UNLOCK;
 309                     gmx_fft_destroy(fft);
 310                     FFTW_LOCK;
 311                     FFTWPREFIX(free)(p1);
 312                     FFTWPREFIX(free)(p2);
 313                     FFTW_UNLOCK;
 314                     return -1;
 315                 }
 316             }
 317         }
 318     }
 319
 320     FFTWPREFIX(free)(p1);
 321     FFTWPREFIX(free)(p2);
 322
 323     fft->real_transform = 1;
 324     fft->ndim           = 1;
 325
 326     *pfft = fft;
 327     FFTW_UNLOCK;
 328     return 0;
 329 }
 330
 331
 332 int
 333 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 334                      int                nx,
 335                      int                ny,
 336                      gmx_fft_flag       flags)
 337 {
 338     gmx_fft_t              fft;
 339     real                  *p1, *p2, *up1, *up2;
 340     size_t                 pc;
 341     int                    i, j, k;
 342     int                    fftw_flags;
 343
 344 #ifdef GMX_DISABLE_FFTW_MEASURE
 345     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 346 #endif
 347
 348     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 349
 350     if (pfft == NULL)
 351     {
 352         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 353         return EINVAL;
 354     }
 355     *pfft = NULL;
 356
 357     FFTW_LOCK;
 358     if ( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
 359     {
 360         FFTW_UNLOCK;
 361         return ENOMEM;
 362     }
 363
 364     /* allocate aligned, and extra memory to make it unaligned */
 365     p1  = (real *) FFTWPREFIX(malloc)(sizeof(real) *( nx*(ny/2+1)*2 + 2) );
 366     if (p1 == NULL)
 367     {
 368         FFTWPREFIX(free)(fft);
 369         FFTW_UNLOCK;
 370         return ENOMEM;
 371     }
 372
 373     p2  = (real *) FFTWPREFIX(malloc)(sizeof(real) *( nx*(ny/2+1)*2 + 2) );
 374     if (p2 == NULL)
 375     {
 376         FFTWPREFIX(free)(p1);
 377         FFTWPREFIX(free)(fft);
 378         FFTW_UNLOCK;
 379         return ENOMEM;
 380     }
 381
 382     /* make unaligned pointers.
 383      * In double precision the actual complex datatype will be 16 bytes,
 384      * so go to a char pointer and force an offset of 8 bytes instead.
 385      */
 386     pc  = (size_t)p1;
 387     pc += 8;
 388     up1 = (real *)pc;
 389
 390     pc  = (size_t)p2;
 391     pc += 8;
 392     up2 = (real *)pc;
 393
 394
 395     fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx, ny, (FFTWPREFIX(complex) *) up1, up2, fftw_flags);
 396     fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx, ny, up1, (FFTWPREFIX(complex) *) up2, fftw_flags);
 397     fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx, ny, (FFTWPREFIX(complex) *) up1, up1, fftw_flags);
 398     fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx, ny, up1, (FFTWPREFIX(complex) *) up1, fftw_flags);
 399
 400     fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx, ny, (FFTWPREFIX(complex) *) p1, p2, fftw_flags);
 401     fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx, ny, p1, (FFTWPREFIX(complex) *) p2, fftw_flags);
 402     fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx, ny, (FFTWPREFIX(complex) *) p1, p1, fftw_flags);
 403     fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx, ny, p1, (FFTWPREFIX(complex) *) p1, fftw_flags);
 404
 405
 406     for (i = 0; i < 2; i++)
 407     {
 408         for (j = 0; j < 2; j++)
 409         {
 410             for (k = 0; k < 2; k++)
 411             {
 412                 if (fft->plan[i][j][k] == NULL)
 413                 {
 414                     gmx_fatal(FARGS, "Error initializing FFTW3 plan.");
 415                     FFTW_UNLOCK;
 416                     gmx_fft_destroy(fft);
 417                     FFTW_LOCK;
 418                     FFTWPREFIX(free)(p1);
 419                     FFTWPREFIX(free)(p2);
 420                     FFTW_UNLOCK;
 421                     return -1;
 422                 }
 423             }
 424         }
 425     }
 426
 427     FFTWPREFIX(free)(p1);
 428     FFTWPREFIX(free)(p2);
 429
 430     fft->real_transform = 1;
 431     fft->ndim           = 2;
 432
 433     *pfft = fft;
 434     FFTW_UNLOCK;
 435     return 0;
 436 }
 437
 438 int
 439 gmx_fft_1d               (gmx_fft_t                  fft,
 440                           enum gmx_fft_direction     dir,
 441                           void *                     in_data,
 442                           void *                     out_data)
 443 {
 444     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf) == 0);
 445     int           inplace   = (in_data == out_data);
 446     int           isforward = (dir == GMX_FFT_FORWARD);
 447
 448     /* Some checks */
 449     if ( (fft->real_transform == 1) || (fft->ndim != 1) ||
 450          ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 451     {
 452         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 453         return EINVAL;
 454     }
 455
 456     FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
 457                             (FFTWPREFIX(complex) *) in_data,
 458                             (FFTWPREFIX(complex) *) out_data);
 459
 460     return 0;
 461 }
 462
 463 int
 464 gmx_fft_many_1d               (gmx_fft_t                  fft,
 465                                enum gmx_fft_direction     dir,
 466                                void *                     in_data,
 467                                void *                     out_data)
 468 {
 469     return gmx_fft_1d(fft, dir, in_data, out_data);
 470 }
 471
 472 int
 473 gmx_fft_1d_real          (gmx_fft_t                  fft,
 474                           enum gmx_fft_direction     dir,
 475                           void *                     in_data,
 476                           void *                     out_data)
 477 {
 478     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf) == 0);
 479     int           inplace   = (in_data == out_data);
 480     int           isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 481
 482     /* Some checks */
 483     if ( (fft->real_transform != 1) || (fft->ndim != 1) ||
 484          ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 485     {
 486         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 487         return EINVAL;
 488     }
 489
 490     if (isforward)
 491     {
 492         FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
 493                                     (real *)in_data, (FFTWPREFIX(complex) *) out_data);
 494     }
 495     else
 496     {
 497         FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
 498                                     (FFTWPREFIX(complex) *) in_data, (real *)out_data);
 499     }
 500
 501     return 0;
 502 }
 503
 504 int
 505 gmx_fft_many_1d_real     (gmx_fft_t                  fft,
 506                           enum gmx_fft_direction     dir,
 507                           void *                     in_data,
 508                           void *                     out_data)
 509 {
 510     return gmx_fft_1d_real(fft, dir, in_data, out_data);
 511 }
 512
 513 int
 514 gmx_fft_2d_real          (gmx_fft_t                  fft,
 515                           enum gmx_fft_direction     dir,
 516                           void *                     in_data,
 517                           void *                     out_data)
 518 {
 519     int           aligned   = ((((size_t)in_data | (size_t)out_data) & 0xf) == 0);
 520     int           inplace   = (in_data == out_data);
 521     int           isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 522
 523     /* Some checks */
 524     if ( (fft->real_transform != 1) || (fft->ndim != 2) ||
 525          ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 526     {
 527         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 528         return EINVAL;
 529     }
 530
 531     if (isforward)
 532     {
 533         FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
 534                                     (real *)in_data,
 535                                     (FFTWPREFIX(complex) *) out_data);
 536     }
 537     else
 538     {
 539         FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
 540                                     (FFTWPREFIX(complex) *) in_data,
 541                                     (real *)out_data);
 542     }
 543
 544
 545     return 0;
 546 }
 547
 548 void
 549 gmx_fft_destroy(gmx_fft_t      fft)
 550 {
 551     int                   i, j, k;
 552
 553     if (fft != NULL)
 554     {
 555         for (i = 0; i < 2; i++)
 556         {
 557             for (j = 0; j < 2; j++)
 558             {
 559                 for (k = 0; k < 2; k++)
 560                 {
 561                     if (fft->plan[i][j][k] != NULL)
 562                     {
 563                         FFTW_LOCK;
 564                         FFTWPREFIX(destroy_plan)(fft->plan[i][j][k]);
 565                         FFTW_UNLOCK;
 566                         fft->plan[i][j][k] = NULL;
 567                     }
 568                 }
 569             }
 570         }
 571         FFTW_LOCK;
 572         FFTWPREFIX(free)(fft);
 573         FFTW_UNLOCK;
 574     }
 575
 576 }
 577
 578 void
 579 gmx_many_fft_destroy(gmx_fft_t    fft)
 580 {
 581     gmx_fft_destroy(fft);
 582 }
 583
 584 void gmx_fft_cleanup()
 585 {
 586     FFTWPREFIX(cleanup)();
 587 }
 588
 589 const char *gmx_fft_get_version_info()
 590 {
 591 #ifdef GMX_NATIVE_WINDOWS
 592     return "fftw3";
 593 #else
 594     return FFTWPREFIX(version);
 595 #endif
 596 }