2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
46 #include "external/fftpack/fftpack.h"
48 #include "gromacs/fft/fft.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/utility/real.h"
54 * Contents of the FFTPACK fft datatype.
56 * Note that this is one of several possible implementations of gmx_fft_t.
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.
64 struct gmx_fft_fftpack
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 */
76 int gmx_fft_init_1d(gmx_fft_t* pfft, int nx, int gmx_unused flags)
82 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
87 if ((fft = static_cast<struct gmx_fft*>(malloc(sizeof(struct gmx_fft)))) == nullptr)
95 /* Need 4*n storage for 1D complex FFT */
96 if ((fft->work = static_cast<real*>(malloc(sizeof(real) * (4 * nx)))) == nullptr)
104 fftpack_cffti1(nx, fft->work, fft->ifac);
112 int gmx_fft_init_1d_real(gmx_fft_t* pfft, int nx, int gmx_unused flags)
118 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
123 if ((fft = static_cast<struct gmx_fft*>(malloc(sizeof(struct gmx_fft)))) == nullptr)
131 /* Need 2*n storage for 1D real FFT */
132 if ((fft->work = static_cast<real*>(malloc(sizeof(real) * (2 * nx)))) == nullptr)
140 fftpack_rffti1(nx, fft->work, fft->ifac);
147 int gmx_fft_init_2d_real(gmx_fft_t* pfft, int nx, int ny, int flags)
150 int nyc = (ny / 2 + 1);
155 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
160 /* Create the X transform */
161 if ((fft = static_cast<struct gmx_fft*>(malloc(sizeof(struct gmx_fft)))) == nullptr)
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.
171 if ((fft->work = static_cast<real*>(malloc(sizeof(real) * (4 * nx + 2 * nx * nyc)))) == nullptr)
176 fftpack_cffti1(nx, fft->work, fft->ifac);
178 /* Create real Y transform as a link from X */
179 if ((rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
190 int gmx_fft_1d(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
200 p1 = static_cast<real*>(in_data);
201 p2 = static_cast<real*>(out_data);
206 /* FFTPACK only does in-place transforms, so emulate out-of-place
207 * by copying data to the output array first.
209 if (in_data != out_data)
211 p1 = static_cast<real*>(in_data);
212 p2 = static_cast<real*>(out_data);
214 /* n complex = 2*n real elements */
215 for (i = 0; i < 2 * n; i++)
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.
225 if (dir == GMX_FFT_FORWARD)
227 fftpack_cfftf1(n, static_cast<real*>(out_data), fft->work + 2 * n, fft->work, fft->ifac, -1);
229 else if (dir == GMX_FFT_BACKWARD)
231 fftpack_cfftf1(n, static_cast<real*>(out_data), fft->work + 2 * n, fft->work, fft->ifac, 1);
235 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
243 int gmx_fft_1d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
253 p1 = static_cast<real*>(in_data);
254 p2 = static_cast<real*>(out_data);
256 if (dir == GMX_FFT_REAL_TO_COMPLEX)
262 if (dir == GMX_FFT_REAL_TO_COMPLEX)
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.
268 if (in_data != out_data)
270 p1 = static_cast<real*>(in_data);
271 p2 = static_cast<real*>(out_data);
273 for (i = 0; i < 2 * (n / 2 + 1); i++)
279 /* Elements 0 .. n-1 in work are used for ffac values,
280 * Elements n .. 2*n-1 are internal FFTPACK work space.
282 fftpack_rfftf1(n, static_cast<real*>(out_data), fft->work + n, fft->work, fft->ifac);
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.
289 p2 = static_cast<real*>(out_data);
290 for (i = n - 1; i > 0; i--)
294 /* imaginary zero freq. */
303 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
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.
311 p1 = static_cast<real*>(in_data);
312 p2 = static_cast<real*>(out_data);
315 for (i = 1; i < n; i++)
319 fftpack_rfftb1(n, static_cast<real*>(out_data), fft->work + n, fft->work, fft->ifac);
323 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
331 int gmx_fft_2d_real(gmx_fft_t fft, enum gmx_fft_direction dir, void* in_data, void* out_data)
333 int i, j, nx, ny, nyc;
341 /* Number of complex elements in y direction */
344 work = fft->work + 4 * nx;
346 if (dir == GMX_FFT_REAL_TO_COMPLEX)
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.
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.
355 if (in_data != out_data)
357 p1 = static_cast<real*>(in_data);
358 p2 = static_cast<real*>(out_data);
360 for (i = 0; i < nx; i++)
362 for (j = 0; j < ny; j++)
364 p2[i * nyc * 2 + j] = p1[i * ny + j];
368 data = static_cast<t_complex*>(out_data);
370 /* y real-to-complex FFTs */
371 for (i = 0; i < nx; i++)
373 gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data + i * nyc, data + i * nyc);
376 /* Transform to get X data in place */
377 gmx_fft_transpose_2d(data, data, nx, nyc);
379 /* Complex-to-complex X FFTs */
380 for (i = 0; i < nyc; i++)
382 gmx_fft_1d(fft, GMX_FFT_FORWARD, data + i * nx, data + i * nx);
386 gmx_fft_transpose_2d(data, data, nyc, nx);
388 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
390 /* An in-place complex-to-real transform is straightforward,
391 * since the output array must be large enough for the padding to fit.
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.
398 if (in_data != out_data)
400 memcpy(work, in_data, sizeof(t_complex) * nx * nyc);
401 data = reinterpret_cast<t_complex*>(work);
406 data = reinterpret_cast<t_complex*>(out_data);
409 /* Transpose to get X arrays */
410 gmx_fft_transpose_2d(data, data, nx, nyc);
413 for (i = 0; i < nyc; i++)
415 gmx_fft_1d(fft, GMX_FFT_BACKWARD, data + i * nx, data + i * nx);
418 /* Transpose to get Y arrays */
419 gmx_fft_transpose_2d(data, data, nyc, nx);
422 for (i = 0; i < nx; i++)
424 gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data + i * nyc, data + i * nyc);
427 if (in_data != out_data)
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.
432 p1 = reinterpret_cast<real*>(data);
433 p2 = static_cast<real*>(out_data);
435 for (i = 0; i < nx; i++)
437 for (j = 0; j < ny; j++)
439 p2[i * ny + j] = p1[i * nyc * 2 + j];
446 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
453 void gmx_fft_destroy(gmx_fft_t fft)
458 if (fft->next != nullptr)
460 gmx_fft_destroy(fft->next);
466 void gmx_fft_cleanup() {}