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, 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.
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.
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.
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.
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.
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.
38 #include "gromacs/fft/fft.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/utility/real.h"
50 #include "external/fftpack/fftpack.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 */
81 gmx_fft_init_1d(gmx_fft_t * pfft,
89 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
94 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
102 /* Need 4*n storage for 1D complex FFT */
103 if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
111 fftpack_cffti1(nx, fft->work, fft->ifac);
121 gmx_fft_init_1d_real(gmx_fft_t * pfft,
123 int gmx_unused flags)
129 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
134 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
142 /* Need 2*n storage for 1D real FFT */
143 if ((fft->work = (real *)malloc(sizeof(real)*(2*nx))) == NULL)
151 fftpack_rffti1(nx, fft->work, fft->ifac);
159 gmx_fft_init_2d_real(gmx_fft_t * pfft,
165 int nyc = (ny/2 + 1);
170 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
175 /* Create the X transform */
176 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
183 /* Need 4*nx storage for 1D complex FFT, and another
184 * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
186 if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
191 fftpack_cffti1(nx, fft->work, fft->ifac);
193 /* Create real Y transform as a link from X */
194 if ( (rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
206 gmx_fft_1d (gmx_fft_t fft,
207 enum gmx_fft_direction dir,
219 p1 = (real *)in_data;
220 p2 = (real *)out_data;
225 /* FFTPACK only does in-place transforms, so emulate out-of-place
226 * by copying data to the output array first.
228 if (in_data != out_data)
230 p1 = (real *)in_data;
231 p2 = (real *)out_data;
233 /* n complex = 2*n real elements */
234 for (i = 0; i < 2*n; i++)
240 /* Elements 0 .. 2*n-1 in work are used for ffac values,
241 * Elements 2*n .. 4*n-1 are internal FFTPACK work space.
244 if (dir == GMX_FFT_FORWARD)
246 fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, -1);
248 else if (dir == GMX_FFT_BACKWARD)
250 fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, 1);
254 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
264 gmx_fft_1d_real (gmx_fft_t fft,
265 enum gmx_fft_direction dir,
277 p1 = (real *)in_data;
278 p2 = (real *)out_data;
280 if (dir == GMX_FFT_REAL_TO_COMPLEX)
286 if (dir == GMX_FFT_REAL_TO_COMPLEX)
288 /* FFTPACK only does in-place transforms, so emulate out-of-place
289 * by copying data to the output array first. This works fine, since
290 * the complex array must be larger than the real.
292 if (in_data != out_data)
294 p1 = (real *)in_data;
295 p2 = (real *)out_data;
297 for (i = 0; i < 2*(n/2+1); i++)
303 /* Elements 0 .. n-1 in work are used for ffac values,
304 * Elements n .. 2*n-1 are internal FFTPACK work space.
306 fftpack_rfftf1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
309 * FFTPACK has a slightly more compact storage than we, time to
310 * convert it: ove most of the array one step up to make room for
311 * zero imaginary parts.
313 p2 = (real *)out_data;
314 for (i = n-1; i > 0; i--)
318 /* imaginary zero freq. */
328 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
330 /* FFTPACK only does in-place transforms, and we cannot just copy
331 * input to output first here since our real array is smaller than
332 * the complex one. However, since the FFTPACK complex storage format
333 * is more compact than ours (2 reals) it will fit, so compact it
334 * and copy on-the-fly to the output array.
336 p1 = (real *) in_data;
337 p2 = (real *)out_data;
340 for (i = 1; i < n; i++)
344 fftpack_rfftb1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
348 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
357 gmx_fft_2d_real (gmx_fft_t fft,
358 enum gmx_fft_direction dir,
362 int i, j, nx, ny, nyc;
370 /* Number of complex elements in y direction */
373 work = fft->work+4*nx;
375 if (dir == GMX_FFT_REAL_TO_COMPLEX)
377 /* If we are doing an in-place transform the 2D array is already
378 * properly padded by the user, and we are all set.
380 * For out-of-place there is no array padding, but FFTPACK only
381 * does in-place FFTs internally, so we need to start by copying
382 * data from the input to the padded (larger) output array.
384 if (in_data != out_data)
386 p1 = (real *)in_data;
387 p2 = (real *)out_data;
389 for (i = 0; i < nx; i++)
391 for (j = 0; j < ny; j++)
393 p2[i*nyc*2+j] = p1[i*ny+j];
397 data = (t_complex *)out_data;
399 /* y real-to-complex FFTs */
400 for (i = 0; i < nx; i++)
402 gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data+i*nyc, data+i*nyc);
405 /* Transform to get X data in place */
406 gmx_fft_transpose_2d(data, data, nx, nyc);
408 /* Complex-to-complex X FFTs */
409 for (i = 0; i < nyc; i++)
411 gmx_fft_1d(fft, GMX_FFT_FORWARD, data+i*nx, data+i*nx);
415 gmx_fft_transpose_2d(data, data, nyc, nx);
418 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
420 /* An in-place complex-to-real transform is straightforward,
421 * since the output array must be large enough for the padding to fit.
423 * For out-of-place complex-to-real transforms we cannot just copy
424 * data to the output array, since it is smaller than the input.
425 * In this case there's nothing to do but employing temporary work data,
426 * starting at work+4*nx and using nx*nyc*2 elements.
428 if (in_data != out_data)
430 memcpy(work, in_data, sizeof(t_complex)*nx*nyc);
431 data = (t_complex *)work;
436 data = (t_complex *)out_data;
439 /* Transpose to get X arrays */
440 gmx_fft_transpose_2d(data, data, nx, nyc);
443 for (i = 0; i < nyc; i++)
445 gmx_fft_1d(fft, GMX_FFT_BACKWARD, data+i*nx, data+i*nx);
448 /* Transpose to get Y arrays */
449 gmx_fft_transpose_2d(data, data, nyc, nx);
452 for (i = 0; i < nx; i++)
454 gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data+i*nyc, data+i*nyc);
457 if (in_data != out_data)
459 /* Output (pointed to by data) is now in padded format.
460 * Pack it into out_data if we were doing an out-of-place transform.
463 p2 = (real *)out_data;
465 for (i = 0; i < nx; i++)
467 for (j = 0; j < ny; j++)
469 p2[i*ny+j] = p1[i*nyc*2+j];
476 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
484 gmx_fft_destroy(gmx_fft_t fft)
489 if (fft->next != NULL)
491 gmx_fft_destroy(fft->next);
497 void gmx_fft_cleanup()
501 const char *gmx_fft_get_version_info()
503 return "fftpack (built-in)";