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