src/gromacs/fft/fft_fftpack.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,2016,2017 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 <errno.h>
  40 #include <stdio.h>
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include <cmath>
  45
  46 #include "external/fftpack/fftpack.h"
  47
  48 #include "gromacs/fft/fft.h"
  49 #include "gromacs/utility/fatalerror.h"
  50 #include "gromacs/utility/real.h"
  51
  52 /*! \internal
  53  * \brief
  54  * Contents of the FFTPACK fft datatype.
  55  *
  56  * Note that this is one of several possible implementations of gmx_fft_t.
  57  *
  58  *  FFTPACK only does 1d transforms, so we use a pointers to another fft for
  59  *  the transform in the next dimension.
  60  * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
  61  * a pointer to a 1d. The 1d structure has next==NULL.
  62  */
  63 #ifdef DOXYGEN
  64 struct gmx_fft_fftpack
  65 #else
  66 struct gmx_fft
  67 #endif
  68 {
  69     int             ndim;     /**< Dimensions, including our subdimensions.  */
  70     int             n;        /**< Number of points in this dimension.       */
  71     int             ifac[15]; /**< 15 bytes needed for cfft and rfft         */
  72     struct gmx_fft* next;     /**< Pointer to next dimension, or NULL.       */
  73     real*           work;     /**< 1st 4n reserved for cfft, 1st 2n for rfft */
  74 };
  75
  76 int gmx_fft_init_1d(gmx_fft_t* pfft, int nx, int gmx_unused flags)
  77 {
  78     gmx_fft_t fft;
  79
  80     if (pfft == nullptr)
  81     {
  82         gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
  83         return EINVAL;
  84     }
  85     *pfft = nullptr;
  86
  87     if ((fft = static_cast<struct gmx_fft*>(malloc(sizeof(struct gmx_fft)))) == nullptr)
  88     {
  89         return ENOMEM;
  90     }
  91
  92     fft->next = nullptr;
  93     fft->n    = nx;
  94
  95     /* Need 4*n storage for 1D complex FFT */
  96     if ((fft->work = static_cast<real*>(malloc(sizeof(real) * (4 * nx)))) == nullptr)
  97     {
  98         free(fft);
  99         return ENOMEM;
 100     }
 101
 102     if (fft->n > 1)
 103     {
 104         fftpack_cffti1(nx, fft->work, fft->ifac);
 105     }
 106
 107     *pfft = fft;
 108     return 0;
 109 };
 110
 111
 112 int gmx_fft_init_1d_real(gmx_fft_t* pfft, int nx, int gmx_unused flags)
 113 {
 114     gmx_fft_t fft;
 115
 116     if (pfft == nullptr)
 117     {
 118         gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
 119         return EINVAL;
 120     }
 121     *pfft = nullptr;
 122
 123     if ((fft = static_cast<struct gmx_fft*>(malloc(sizeof(struct gmx_fft)))) == nullptr)
 124     {
 125         return ENOMEM;
 126     }
 127
 128     fft->next = nullptr;
 129     fft->n    = nx;
 130
 131     /* Need 2*n storage for 1D real FFT */
 132     if ((fft->work = static_cast<real*>(malloc(sizeof(real) * (2 * nx)))) == nullptr)
 133     {
 134         free(fft);
 135         return ENOMEM;
 136     }
 137
 138     if (fft->n > 1)
 139     {
 140         fftpack_rffti1(nx, fft->work, fft->ifac);
 141     }
 142
 143     *pfft = fft;
 144     return 0;
 145 }
 146
 147 int gmx_fft_init_2d_real(gmx_fft_t* pfft, int nx, int ny, int flags)
 148 {
 149     gmx_fft_t fft;
 150     int       nyc = (ny / 2 + 1);
 151     int       rc;
 152
 153     if (pfft == nullptr)
 154     {
 155         gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
 156         return EINVAL;
 157     }
 158     *pfft = nullptr;
 159
 160     /* Create the X transform */
 161     if ((fft = static_cast<struct gmx_fft*>(malloc(sizeof(struct gmx_fft)))) == nullptr)
 162     {
 163         return ENOMEM;
 164     }
 165
 166     fft->n = nx;
 167
 168     /* Need 4*nx storage for 1D complex FFT, and another
 169      * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
 170      */
 171     if ((fft->work = static_cast<real*>(malloc(sizeof(real) * (4 * nx + 2 * nx * nyc)))) == nullptr)
 172     {
 173         free(fft);
 174         return ENOMEM;
 175     }
 176     fftpack_cffti1(nx, fft->work, fft->ifac);
 177
 178     /* Create real Y transform as a link from X */
 179     if ((rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
 180     {
 181         free(fft);
 182         return rc;
 183     }
 184
 185     *pfft = fft;
 186     return 0;
 187 }
 188
 189
 190 int gmx_fft_1d(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 191 {
 192     int   i, n;
 193     real* p1;
 194     real* p2;
 195
 196     n = fft->n;
 197
 198     if (n == 1)
 199     {
 200         p1    = static_cast<real*>(in_data);
 201         p2    = static_cast<real*>(out_data);
 202         p2[0] = p1[0];
 203         p2[1] = p1[1];
 204     }
 205
 206     /* FFTPACK only does in-place transforms, so emulate out-of-place
 207      * by copying data to the output array first.
 208      */
 209     if (in_data != out_data)
 210     {
 211         p1 = static_cast<real*>(in_data);
 212         p2 = static_cast<real*>(out_data);
 213
 214         /* n complex = 2*n real elements */
 215         for (i = 0; i < 2 * n; i++)
 216         {
 217             p2[i] = p1[i];
 218         }
 219     }
 220
 221     /* Elements 0   .. 2*n-1 in work are used for ffac values,
 222      * Elements 2*n .. 4*n-1 are internal FFTPACK work space.
 223      */
 224
 225     if (dir == GMX_FFT_FORWARD)
 226     {
 227         fftpack_cfftf1(n, static_cast<real*>(out_data), fft->work + 2 * n, fft->work, fft->ifac, -1);
 228     }
 229     else if (dir == GMX_FFT_BACKWARD)
 230     {
 231         fftpack_cfftf1(n, static_cast<real*>(out_data), fft->work + 2 * n, fft->work, fft->ifac, 1);
 232     }
 233     else
 234     {
 235         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 236         return EINVAL;
 237     }
 238
 239     return 0;
 240 }
 241
 242
 243 int gmx_fft_1d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 244 {
 245     int   i, n;
 246     real* p1;
 247     real* p2;
 248
 249     n = fft->n;
 250
 251     if (n == 1)
 252     {
 253         p1    = static_cast<real*>(in_data);
 254         p2    = static_cast<real*>(out_data);
 255         p2[0] = p1[0];
 256         if (dir == GMX_FFT_REAL_TO_COMPLEX)
 257         {
 258             p2[1] = 0.0;
 259         }
 260     }
 261
 262     if (dir == GMX_FFT_REAL_TO_COMPLEX)
 263     {
 264         /* FFTPACK only does in-place transforms, so emulate out-of-place
 265          * by copying data to the output array first. This works fine, since
 266          * the complex array must be larger than the real.
 267          */
 268         if (in_data != out_data)
 269         {
 270             p1 = static_cast<real*>(in_data);
 271             p2 = static_cast<real*>(out_data);
 272
 273             for (i = 0; i < 2 * (n / 2 + 1); i++)
 274             {
 275                 p2[i] = p1[i];
 276             }
 277         }
 278
 279         /* Elements 0 ..   n-1 in work are used for ffac values,
 280          * Elements n .. 2*n-1 are internal FFTPACK work space.
 281          */
 282         fftpack_rfftf1(n, static_cast<real*>(out_data), fft->work + n, fft->work, fft->ifac);
 283
 284         /*
 285          * FFTPACK has a slightly more compact storage than we, time to
 286          * convert it: ove most of the array one step up to make room for
 287          * zero imaginary parts.
 288          */
 289         p2 = static_cast<real*>(out_data);
 290         for (i = n - 1; i > 0; i--)
 291         {
 292             p2[i + 1] = p2[i];
 293         }
 294         /* imaginary zero freq. */
 295         p2[1] = 0;
 296
 297         /* Is n even? */
 298         if ((n & 0x1) == 0)
 299         {
 300             p2[n + 1] = 0;
 301         }
 302     }
 303     else if (dir == GMX_FFT_COMPLEX_TO_REAL)
 304     {
 305         /* FFTPACK only does in-place transforms, and we cannot just copy
 306          * input to output first here since our real array is smaller than
 307          * the complex one. However, since the FFTPACK complex storage format
 308          * is more compact than ours (2 reals) it will fit, so compact it
 309          * and copy on-the-fly to the output array.
 310          */
 311         p1 = static_cast<real*>(in_data);
 312         p2 = static_cast<real*>(out_data);
 313
 314         p2[0] = p1[0];
 315         for (i = 1; i < n; i++)
 316         {
 317             p2[i] = p1[i + 1];
 318         }
 319         fftpack_rfftb1(n, static_cast<real*>(out_data), fft->work + n, fft->work, fft->ifac);
 320     }
 321     else
 322     {
 323         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 324         return EINVAL;
 325     }
 326
 327     return 0;
 328 }
 329
 330
 331 int gmx_fft_2d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
 332 {
 333     int        i, j, nx, ny, nyc;
 334     t_complex* data;
 335     real*      work;
 336     real*      p1;
 337     real*      p2;
 338
 339     nx = fft->n;
 340     ny = fft->next->n;
 341     /* Number of complex elements in y direction */
 342     nyc = (ny / 2 + 1);
 343
 344     work = fft->work + 4 * nx;
 345
 346     if (dir == GMX_FFT_REAL_TO_COMPLEX)
 347     {
 348         /* If we are doing an in-place transform the 2D array is already
 349          * properly padded by the user, and we are all set.
 350          *
 351          * For out-of-place there is no array padding, but FFTPACK only
 352          * does in-place FFTs internally, so we need to start by copying
 353          * data from the input to the padded (larger) output array.
 354          */
 355         if (in_data != out_data)
 356         {
 357             p1 = static_cast<real*>(in_data);
 358             p2 = static_cast<real*>(out_data);
 359
 360             for (i = 0; i < nx; i++)
 361             {
 362                 for (j = 0; j < ny; j++)
 363                 {
 364                     p2[i * nyc * 2 + j] = p1[i * ny + j];
 365                 }
 366             }
 367         }
 368         data = static_cast<t_complex*>(out_data);
 369
 370         /* y real-to-complex FFTs */
 371         for (i = 0; i < nx; i++)
 372         {
 373             gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data + i * nyc, data + i * nyc);
 374         }
 375
 376         /* Transform to get X data in place */
 377         gmx_fft_transpose_2d(data, data, nx, nyc);
 378
 379         /* Complex-to-complex X FFTs */
 380         for (i = 0; i < nyc; i++)
 381         {
 382             gmx_fft_1d(fft, GMX_FFT_FORWARD, data + i * nx, data + i * nx);
 383         }
 384
 385         /* Transpose back */
 386         gmx_fft_transpose_2d(data, data, nyc, nx);
 387     }
 388     else if (dir == GMX_FFT_COMPLEX_TO_REAL)
 389     {
 390         /* An in-place complex-to-real transform is straightforward,
 391          * since the output array must be large enough for the padding to fit.
 392          *
 393          * For out-of-place complex-to-real transforms we cannot just copy
 394          * data to the output array, since it is smaller than the input.
 395          * In this case there's nothing to do but employing temporary work data,
 396          * starting at work+4*nx and using nx*nyc*2 elements.
 397          */
 398         if (in_data != out_data)
 399         {
 400             memcpy(work, in_data, sizeof(t_complex) * nx * nyc);
 401             data = reinterpret_cast<t_complex*>(work);
 402         }
 403         else
 404         {
 405             /* in-place */
 406             data = reinterpret_cast<t_complex*>(out_data);
 407         }
 408
 409         /* Transpose to get X arrays */
 410         gmx_fft_transpose_2d(data, data, nx, nyc);
 411
 412         /* Do X iFFTs */
 413         for (i = 0; i < nyc; i++)
 414         {
 415             gmx_fft_1d(fft, GMX_FFT_BACKWARD, data + i * nx, data + i * nx);
 416         }
 417
 418         /* Transpose to get Y arrays */
 419         gmx_fft_transpose_2d(data, data, nyc, nx);
 420
 421         /* Do Y iFFTs */
 422         for (i = 0; i < nx; i++)
 423         {
 424             gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data + i * nyc, data + i * nyc);
 425         }
 426
 427         if (in_data != out_data)
 428         {
 429             /* Output (pointed to by data) is now in padded format.
 430              * Pack it into out_data if we were doing an out-of-place transform.
 431              */
 432             p1 = reinterpret_cast<real*>(data);
 433             p2 = static_cast<real*>(out_data);
 434
 435             for (i = 0; i < nx; i++)
 436             {
 437                 for (j = 0; j < ny; j++)
 438                 {
 439                     p2[i * ny + j] = p1[i * nyc * 2 + j];
 440                 }
 441             }
 442         }
 443     }
 444     else
 445     {
 446         gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
 447         return EINVAL;
 448     }
 449
 450     return 0;
 451 }
 452
 453 void gmx_fft_destroy(gmx_fft_t fft)
 454 {
 455     if (fft != nullptr)
 456     {
 457         free(fft->work);
 458         if (fft->next != nullptr)
 459         {
 460             gmx_fft_destroy(fft->next);
 461         }
 462         free(fft);
 463     }
 464 }
 465
 466 void gmx_fft_cleanup() {}