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 by the GROMACS development team.
   6  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <cerrno>
  42 #include <cstdlib>
  43
  44 #include <fftw3.h>
  45
  46 #include "gromacs/fft/fft.h"
  47 #include "gromacs/utility/exceptions.h"
  48 #include "gromacs/utility/fatalerror.h"
  49 #include "gromacs/utility/mutex.h"
  50
  51 #if 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 gmx::Mutex big_fftw_mutex;
  60 #define FFTW_LOCK              \
  61     try                        \
  62     {                          \
  63         big_fftw_mutex.lock(); \
  64     }                          \
  65     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
  66 #define FFTW_UNLOCK              \
  67     try                          \
  68     {                            \
  69         big_fftw_mutex.unlock(); \
  70     }                            \
  71     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
  72
  73 /* 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.
  74    Consequesences:
  75    - It is not allowed to use these FFT plans from memory which doesn't have a starting address as a multiple of 8 bytes.
  76      This is OK as long as the memory directly comes from malloc and is not some subarray within alloated memory.
  77    - This has to be fixed if any future architecute requires memory to be aligned to multiples of 32 bytes.
  78  */
  79 /*! \internal
  80  * \brief
  81  * Contents of the FFTW3 fft datatype.
  82  *
  83  * Note that this is one of several possible implementations of gmx_fft_t.
  84  */
  85 #ifdef DOXYGEN
  86 struct gmx_fft_fftw3
  87 #else
  88 struct gmx_fft
  89 #endif
  90 {
  91     /*! \brief
  92      * FFTW plans.
  93      *
  94      * Three alternatives (unaligned/aligned, out-of-place/in-place, forward/backward)
  95      * results in 8 different FFTW plans. Keep track of them with 3 array indices:
  96      * first index:   0=unaligned, 1=aligned
  97      * second index:  0=out-of-place, 1=in-place
  98      * third index:   0=backward, 1=forward
  99      */
 100     FFTWPREFIX(plan) plan[2][2][2];
 101     /** Used to catch user mistakes */
 102     int real_transform;
 103     /** Number of dimensions in the FFT */
 104     int ndim;
 105 };
 106
 107 int gmx_fft_init_1d(gmx_fft_t* pfft, int nx, gmx_fft_flag flags)
 108 {
 109     return gmx_fft_init_many_1d(pfft, nx, 1, flags);
 110 }
 111
 112
 113 int gmx_fft_init_many_1d(gmx_fft_t* pfft, int nx, int howmany, gmx_fft_flag flags)
 114 {
 115     gmx_fft_t fft;
 116     FFTWPREFIX(complex) * p1, *p2, *up1, *up2;
 117     char* pc;
 118     int   i, j, k;
 119     int   fftw_flags;
 120
 121 #if GMX_DISABLE_FFTW_MEASURE
 122     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 123 #endif
 124
 125     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 126
 127     if (pfft == nullptr)
 128     {
 129         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 130         return EINVAL;
 131     }
 132     *pfft = nullptr;
 133
 134     FFTW_LOCK
 135     if ((fft = static_cast<gmx_fft_t>(FFTWPREFIX(malloc)(sizeof(struct gmx_fft)))) == nullptr)
 136     {
 137         FFTW_UNLOCK
 138         return ENOMEM;
 139     }
 140
 141     /* allocate aligned, and extra memory to make it unaligned */
 142     p1 = static_cast<FFTWPREFIX(complex)*>(
 143             FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex)) * (nx + 2) * howmany));
 144     if (p1 == nullptr)
 145     {
 146         FFTWPREFIX(free)(fft);
 147         FFTW_UNLOCK
 148         return ENOMEM;
 149     }
 150
 151     p2 = static_cast<FFTWPREFIX(complex)*>(
 152             FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex)) * (nx + 2) * howmany));
 153     if (p2 == nullptr)
 154     {
 155         FFTWPREFIX(free)(p1);
 156         FFTWPREFIX(free)(fft);
 157         FFTW_UNLOCK
 158         return ENOMEM;
 159     }
 160
 161     /* make unaligned pointers.
 162      * In double precision the actual complex datatype will be 16 bytes,
 163      * so go to a char pointer and force an offset of 8 bytes instead.
 164      */
 165     pc = reinterpret_cast<char*>(p1);
 166     pc += 8;
 167     up1 = reinterpret_cast<FFTWPREFIX(complex)*>(pc);
 168
 169     pc = reinterpret_cast<char*>(p2);
 170     pc += 8;
 171     up2 = reinterpret_cast<FFTWPREFIX(complex)*>(pc);
 172
 173     /*                            int rank, const int *n, int howmany,
 174                                   fftw_complex *in, const int *inembed,
 175                                   int istride, int idist,
 176                                   fftw_complex *out, const int *onembed,
 177                                   int ostride, int odist,
 178                                   int sign, unsigned flags */
 179     fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(
 180             1, &nx, howmany, up1, &nx, 1, nx, up2, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 181     fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(
 182             1, &nx, howmany, up1, &nx, 1, nx, up2, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 183     fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(
 184             1, &nx, howmany, up1, &nx, 1, nx, up1, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 185     fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(
 186             1, &nx, howmany, up1, &nx, 1, nx, up1, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 187     fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(
 188             1, &nx, howmany, p1, &nx, 1, nx, p2, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 189     fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(
 190             1, &nx, howmany, p1, &nx, 1, nx, p2, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 191     fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(
 192             1, &nx, howmany, p1, &nx, 1, nx, p1, &nx, 1, nx, FFTW_BACKWARD, fftw_flags);
 193     fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(
 194             1, &nx, howmany, p1, &nx, 1, nx, p1, &nx, 1, nx, FFTW_FORWARD, fftw_flags);
 195
 196     for (i = 0; i < 2; i++)
 197     {
 198         for (j = 0; j < 2; j++)
 199         {
 200             for (k = 0; k < 2; k++)
 201             {
 202                 if (fft->plan[i][j][k] == nullptr)
 203                 {
 204                     gmx_fatal(FARGS, "Error initializing FFTW3 plan.");
 205                     FFTW_UNLOCK
 206                     gmx_fft_destroy(fft);
 207                     FFTW_LOCK
 208                     FFTWPREFIX(free)(p1);
 209                     FFTWPREFIX(free)(p2);
 210                     FFTW_UNLOCK
 211                     return -1;
 212                 }
 213             }
 214         }
 215     }
 216
 217     FFTWPREFIX(free)(p1);
 218     FFTWPREFIX(free)(p2);
 219
 220     fft->real_transform = 0;
 221     fft->ndim           = 1;
 222
 223     *pfft = fft;
 224     FFTW_UNLOCK
 225     return 0;
 226 }
 227
 228 int gmx_fft_init_1d_real(gmx_fft_t* pfft, int nx, gmx_fft_flag flags)
 229 {
 230     return gmx_fft_init_many_1d_real(pfft, nx, 1, flags);
 231 }
 232
 233 int gmx_fft_init_many_1d_real(gmx_fft_t* pfft, int nx, int howmany, gmx_fft_flag flags)
 234 {
 235     gmx_fft_t fft;
 236     real *    p1, *p2, *up1, *up2;
 237     char*     pc;
 238     int       i, j, k;
 239     int       fftw_flags;
 240
 241 #if GMX_DISABLE_FFTW_MEASURE
 242     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 243 #endif
 244
 245     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 246
 247     if (pfft == nullptr)
 248     {
 249         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 250         return EINVAL;
 251     }
 252     *pfft = nullptr;
 253
 254     FFTW_LOCK
 255     if ((fft = static_cast<gmx_fft_t>(FFTWPREFIX(malloc)(sizeof(struct gmx_fft)))) == nullptr)
 256     {
 257         FFTW_UNLOCK
 258         return ENOMEM;
 259     }
 260
 261     /* allocate aligned, and extra memory to make it unaligned */
 262     p1 = static_cast<real*>(FFTWPREFIX(malloc)(sizeof(real) * (nx / 2 + 1) * 2 * howmany + 8));
 263     if (p1 == nullptr)
 264     {
 265         FFTWPREFIX(free)(fft);
 266         FFTW_UNLOCK
 267         return ENOMEM;
 268     }
 269
 270     p2 = static_cast<real*>(FFTWPREFIX(malloc)(sizeof(real) * (nx / 2 + 1) * 2 * howmany + 8));
 271     if (p2 == nullptr)
 272     {
 273         FFTWPREFIX(free)(p1);
 274         FFTWPREFIX(free)(fft);
 275         FFTW_UNLOCK
 276         return ENOMEM;
 277     }
 278
 279     /* make unaligned pointers.
 280      * In double precision the actual complex datatype will be 16 bytes,
 281      * so go to a char pointer and force an offset of 8 bytes instead.
 282      */
 283     pc = reinterpret_cast<char*>(p1);
 284     pc += 8;
 285     up1 = reinterpret_cast<real*>(pc);
 286
 287     pc = reinterpret_cast<char*>(p2);
 288     pc += 8;
 289     up2 = reinterpret_cast<real*>(pc);
 290
 291     /*                                int rank, const int *n, int howmany,
 292                                       double *in, const int *inembed,
 293                                       int istride, int idist,
 294                                       fftw_complex *out, const int *onembed,
 295                                       int ostride, int odist,
 296                                       unsigned flag    */
 297     fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1,
 298                                                        &nx,
 299                                                        howmany,
 300                                                        up1,
 301                                                        nullptr,
 302                                                        1,
 303                                                        (nx / 2 + 1) * 2,
 304                                                        reinterpret_cast<FFTWPREFIX(complex)*>(up2),
 305                                                        nullptr,
 306                                                        1,
 307                                                        (nx / 2 + 1),
 308                                                        fftw_flags);
 309     fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1,
 310                                                        &nx,
 311                                                        howmany,
 312                                                        up1,
 313                                                        nullptr,
 314                                                        1,
 315                                                        (nx / 2 + 1) * 2,
 316                                                        reinterpret_cast<FFTWPREFIX(complex)*>(up1),
 317                                                        nullptr,
 318                                                        1,
 319                                                        (nx / 2 + 1),
 320                                                        fftw_flags);
 321     fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1,
 322                                                        &nx,
 323                                                        howmany,
 324                                                        p1,
 325                                                        nullptr,
 326                                                        1,
 327                                                        (nx / 2 + 1) * 2,
 328                                                        reinterpret_cast<FFTWPREFIX(complex)*>(p2),
 329                                                        nullptr,
 330                                                        1,
 331                                                        (nx / 2 + 1),
 332                                                        fftw_flags);
 333     fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1,
 334                                                        &nx,
 335                                                        howmany,
 336                                                        p1,
 337                                                        nullptr,
 338                                                        1,
 339                                                        (nx / 2 + 1) * 2,
 340                                                        reinterpret_cast<FFTWPREFIX(complex)*>(p1),
 341                                                        nullptr,
 342                                                        1,
 343                                                        (nx / 2 + 1),
 344                                                        fftw_flags);
 345
 346     fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1,
 347                                                        &nx,
 348                                                        howmany,
 349                                                        reinterpret_cast<FFTWPREFIX(complex)*>(up1),
 350                                                        nullptr,
 351                                                        1,
 352                                                        (nx / 2 + 1),
 353                                                        up2,
 354                                                        nullptr,
 355                                                        1,
 356                                                        (nx / 2 + 1) * 2,
 357                                                        fftw_flags);
 358     fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1,
 359                                                        &nx,
 360                                                        howmany,
 361                                                        reinterpret_cast<FFTWPREFIX(complex)*>(up1),
 362                                                        nullptr,
 363                                                        1,
 364                                                        (nx / 2 + 1),
 365                                                        up1,
 366                                                        nullptr,
 367                                                        1,
 368                                                        (nx / 2 + 1) * 2,
 369                                                        fftw_flags);
 370     fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1,
 371                                                        &nx,
 372                                                        howmany,
 373                                                        reinterpret_cast<FFTWPREFIX(complex)*>(p1),
 374                                                        nullptr,
 375                                                        1,
 376                                                        (nx / 2 + 1),
 377                                                        p2,
 378                                                        nullptr,
 379                                                        1,
 380                                                        (nx / 2 + 1) * 2,
 381                                                        fftw_flags);
 382     fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1,
 383                                                        &nx,
 384                                                        howmany,
 385                                                        reinterpret_cast<FFTWPREFIX(complex)*>(p1),
 386                                                        nullptr,
 387                                                        1,
 388                                                        (nx / 2 + 1),
 389                                                        p1,
 390                                                        nullptr,
 391                                                        1,
 392                                                        (nx / 2 + 1) * 2,
 393                                                        fftw_flags);
 394
 395     for (i = 0; i < 2; i++)
 396     {
 397         for (j = 0; j < 2; j++)
 398         {
 399             for (k = 0; k < 2; k++)
 400             {
 401                 if (fft->plan[i][j][k] == nullptr)
 402                 {
 403                     gmx_fatal(FARGS, "Error initializing FFTW3 plan.");
 404                     FFTW_UNLOCK
 405                     gmx_fft_destroy(fft);
 406                     FFTW_LOCK
 407                     FFTWPREFIX(free)(p1);
 408                     FFTWPREFIX(free)(p2);
 409                     FFTW_UNLOCK
 410                     return -1;
 411                 }
 412             }
 413         }
 414     }
 415
 416     FFTWPREFIX(free)(p1);
 417     FFTWPREFIX(free)(p2);
 418
 419     fft->real_transform = 1;
 420     fft->ndim           = 1;
 421
 422     *pfft = fft;
 423     FFTW_UNLOCK
 424     return 0;
 425 }
 426
 427
 428 int gmx_fft_init_2d_real(gmx_fft_t* pfft, int nx, int ny, gmx_fft_flag flags)
 429 {
 430     gmx_fft_t fft;
 431     real *    p1, *p2, *up1, *up2;
 432     char*     pc;
 433     int       i, j, k;
 434     int       fftw_flags;
 435
 436 #if GMX_DISABLE_FFTW_MEASURE
 437     flags |= GMX_FFT_FLAG_CONSERVATIVE;
 438 #endif
 439
 440     fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
 441
 442     if (pfft == nullptr)
 443     {
 444         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
 445         return EINVAL;
 446     }
 447     *pfft = nullptr;
 448
 449     FFTW_LOCK
 450     if ((fft = static_cast<gmx_fft_t>(FFTWPREFIX(malloc)(sizeof(struct gmx_fft)))) == nullptr)
 451     {
 452         FFTW_UNLOCK
 453         return ENOMEM;
 454     }
 455
 456     /* allocate aligned, and extra memory to make it unaligned */
 457     p1 = static_cast<real*>(FFTWPREFIX(malloc)(sizeof(real) * (nx * (ny / 2 + 1) * 2 + 2)));
 458     if (p1 == nullptr)
 459     {
 460         FFTWPREFIX(free)(fft);
 461         FFTW_UNLOCK
 462         return ENOMEM;
 463     }
 464
 465     p2 = static_cast<real*>(FFTWPREFIX(malloc)(sizeof(real) * (nx * (ny / 2 + 1) * 2 + 2)));
 466     if (p2 == nullptr)
 467     {
 468         FFTWPREFIX(free)(p1);
 469         FFTWPREFIX(free)(fft);
 470         FFTW_UNLOCK
 471         return ENOMEM;
 472     }
 473
 474     /* make unaligned pointers.
 475      * In double precision the actual complex datatype will be 16 bytes,
 476      * so go to a char pointer and force an offset of 8 bytes instead.
 477      */
 478     pc = reinterpret_cast<char*>(p1);
 479     pc += 8;
 480     up1 = reinterpret_cast<real*>(pc);
 481
 482     pc = reinterpret_cast<char*>(p2);
 483     pc += 8;
 484     up2 = reinterpret_cast<real*>(pc);
 485
 486
 487     fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(
 488             nx, ny, reinterpret_cast<FFTWPREFIX(complex)*>(up1), up2, fftw_flags);
 489     fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(
 490             nx, ny, up1, reinterpret_cast<FFTWPREFIX(complex)*>(up2), fftw_flags);
 491     fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(
 492             nx, ny, reinterpret_cast<FFTWPREFIX(complex)*>(up1), up1, fftw_flags);
 493     fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(
 494             nx, ny, up1, reinterpret_cast<FFTWPREFIX(complex)*>(up1), fftw_flags);
 495
 496     fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(
 497             nx, ny, reinterpret_cast<FFTWPREFIX(complex)*>(p1), p2, fftw_flags);
 498     fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(
 499             nx, ny, p1, reinterpret_cast<FFTWPREFIX(complex)*>(p2), fftw_flags);
 500     fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(
 501             nx, ny, reinterpret_cast<FFTWPREFIX(complex)*>(p1), p1, fftw_flags);
 502     fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(
 503             nx, ny, p1, reinterpret_cast<FFTWPREFIX(complex)*>(p1), fftw_flags);
 504
 505
 506     for (i = 0; i < 2; i++)
 507     {
 508         for (j = 0; j < 2; j++)
 509         {
 510             for (k = 0; k < 2; k++)
 511             {
 512                 if (fft->plan[i][j][k] == nullptr)
 513                 {
 514                     gmx_fatal(FARGS, "Error initializing FFTW3 plan.");
 515                     FFTW_UNLOCK
 516                     gmx_fft_destroy(fft);
 517                     FFTW_LOCK
 518                     FFTWPREFIX(free)(p1);
 519                     FFTWPREFIX(free)(p2);
 520                     FFTW_UNLOCK
 521                     return -1;
 522                 }
 523             }
 524         }
 525     }
 526
 527     FFTWPREFIX(free)(p1);
 528     FFTWPREFIX(free)(p2);
 529
 530     fft->real_transform = 1;
 531     fft->ndim           = 2;
 532
 533     *pfft = fft;
 534     FFTW_UNLOCK
 535     return 0;
 536 }
 537
 538 int gmx_fft_1d(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 539 {
 540     bool aligned   = (((size_t(in_data) | size_t(out_data)) & 0xf) == 0);
 541     bool inplace   = (in_data == out_data);
 542     bool isforward = (dir == GMX_FFT_FORWARD);
 543
 544     /* Some checks */
 545     if ((fft->real_transform == 1) || (fft->ndim != 1)
 546         || ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)))
 547     {
 548         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 549         return EINVAL;
 550     }
 551
 552     FFTWPREFIX(execute_dft)
 553     (fft->plan[aligned][inplace][isforward],
 554      static_cast<FFTWPREFIX(complex)*>(in_data),
 555      static_cast<FFTWPREFIX(complex)*>(out_data));
 556
 557     return 0;
 558 }
 559
 560 int gmx_fft_many_1d(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 561 {
 562     return gmx_fft_1d(fft, dir, in_data, out_data);
 563 }
 564
 565 int gmx_fft_1d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 566 {
 567     bool aligned   = (((size_t(in_data) | size_t(out_data)) & 0xf) == 0);
 568     bool inplace   = (in_data == out_data);
 569     bool isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 570
 571     /* Some checks */
 572     if ((fft->real_transform != 1) || (fft->ndim != 1)
 573         || ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
 574     {
 575         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 576         return EINVAL;
 577     }
 578
 579     if (isforward)
 580     {
 581         FFTWPREFIX(execute_dft_r2c)
 582         (fft->plan[aligned][inplace][isforward],
 583          static_cast<real*>(in_data),
 584          static_cast<FFTWPREFIX(complex)*>(out_data));
 585     }
 586     else
 587     {
 588         FFTWPREFIX(execute_dft_c2r)
 589         (fft->plan[aligned][inplace][isforward],
 590          static_cast<FFTWPREFIX(complex)*>(in_data),
 591          static_cast<real*>(out_data));
 592     }
 593
 594     return 0;
 595 }
 596
 597 int gmx_fft_many_1d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 598 {
 599     return gmx_fft_1d_real(fft, dir, in_data, out_data);
 600 }
 601
 602 int gmx_fft_2d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 603 {
 604     bool aligned   = (((size_t(in_data) | size_t(out_data)) & 0xf) == 0);
 605     bool inplace   = (in_data == out_data);
 606     bool isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
 607
 608     /* Some checks */
 609     if ((fft->real_transform != 1) || (fft->ndim != 2)
 610         || ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)))
 611     {
 612         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 613         return EINVAL;
 614     }
 615
 616     if (isforward)
 617     {
 618         FFTWPREFIX(execute_dft_r2c)
 619         (fft->plan[aligned][inplace][isforward],
 620          static_cast<real*>(in_data),
 621          static_cast<FFTWPREFIX(complex)*>(out_data));
 622     }
 623     else
 624     {
 625         FFTWPREFIX(execute_dft_c2r)
 626         (fft->plan[aligned][inplace][isforward],
 627          static_cast<FFTWPREFIX(complex)*>(in_data),
 628          static_cast<real*>(out_data));
 629     }
 630
 631
 632     return 0;
 633 }
 634
 635 void gmx_fft_destroy(gmx_fft_t fft)
 636 {
 637     int i, j, k;
 638
 639     if (fft != nullptr)
 640     {
 641         for (i = 0; i < 2; i++)
 642         {
 643             for (j = 0; j < 2; j++)
 644             {
 645                 for (k = 0; k < 2; k++)
 646                 {
 647                     if (fft->plan[i][j][k] != nullptr)
 648                     {
 649                         FFTW_LOCK
 650                         FFTWPREFIX(destroy_plan)(fft->plan[i][j][k]);
 651                         FFTW_UNLOCK
 652                         fft->plan[i][j][k] = nullptr;
 653                     }
 654                 }
 655             }
 656         }
 657         FFTW_LOCK
 658         FFTWPREFIX(free)(fft);
 659         FFTW_UNLOCK
 660     }
 661 }
 662
 663 void gmx_many_fft_destroy(gmx_fft_t fft)
 664 {
 665     gmx_fft_destroy(fft);
 666 }
 667
 668 void gmx_fft_cleanup()
 669 {
 670     FFTWPREFIX(cleanup)();
 671 }