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