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