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