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