1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * Gromacs 4.0 Copyright (c) 1991-2003
5 * David van der Spoel, Erik Lindahl, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
27 #include "gromacs/fft/fft.h"
28 #include "gmx_fatal.h"
30 #include "external/fftpack/fftpack.h"
33 * Contents of the FFTPACK fft datatype.
35 * Note that this is one of several possible implementations of gmx_fft_t.
37 * FFTPACK only does 1d transforms, so we use a pointers to another fft for
38 * the transform in the next dimension.
39 * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
40 * a pointer to a 1d. The 1d structure has next==NULL.
43 struct gmx_fft_fftpack
48 int ndim; /**< Dimensions, including our subdimensions. */
49 int n; /**< Number of points in this dimension. */
50 int ifac[15]; /**< 15 bytes needed for cfft and rfft */
51 struct gmx_fft *next; /**< Pointer to next dimension, or NULL. */
52 real * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
60 gmx_fft_init_1d(gmx_fft_t * pfft,
68 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
73 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
81 /* Need 4*n storage for 1D complex FFT */
82 if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
90 fftpack_cffti1(nx, fft->work, fft->ifac);
100 gmx_fft_init_1d_real(gmx_fft_t * pfft,
108 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
113 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
121 /* Need 2*n storage for 1D real FFT */
122 if ((fft->work = (real *)malloc(sizeof(real)*(2*nx))) == NULL)
130 fftpack_rffti1(nx, fft->work, fft->ifac);
138 gmx_fft_init_2d_real(gmx_fft_t * pfft,
144 int nyc = (ny/2 + 1);
149 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
154 /* Create the X transform */
155 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
162 /* Need 4*nx storage for 1D complex FFT, and another
163 * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
165 if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
170 fftpack_cffti1(nx, fft->work, fft->ifac);
172 /* Create real Y transform as a link from X */
173 if ( (rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
185 gmx_fft_1d (gmx_fft_t fft,
186 enum gmx_fft_direction dir,
198 p1 = (real *)in_data;
199 p2 = (real *)out_data;
204 /* FFTPACK only does in-place transforms, so emulate out-of-place
205 * by copying data to the output array first.
207 if (in_data != out_data)
209 p1 = (real *)in_data;
210 p2 = (real *)out_data;
212 /* n complex = 2*n real elements */
213 for (i = 0; i < 2*n; i++)
219 /* Elements 0 .. 2*n-1 in work are used for ffac values,
220 * Elements 2*n .. 4*n-1 are internal FFTPACK work space.
223 if (dir == GMX_FFT_FORWARD)
225 fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, -1);
227 else if (dir == GMX_FFT_BACKWARD)
229 fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, 1);
233 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
243 gmx_fft_1d_real (gmx_fft_t fft,
244 enum gmx_fft_direction dir,
256 p1 = (real *)in_data;
257 p2 = (real *)out_data;
259 if (dir == GMX_FFT_REAL_TO_COMPLEX)
265 if (dir == GMX_FFT_REAL_TO_COMPLEX)
267 /* FFTPACK only does in-place transforms, so emulate out-of-place
268 * by copying data to the output array first. This works fine, since
269 * the complex array must be larger than the real.
271 if (in_data != out_data)
273 p1 = (real *)in_data;
274 p2 = (real *)out_data;
276 for (i = 0; i < 2*(n/2+1); i++)
282 /* Elements 0 .. n-1 in work are used for ffac values,
283 * Elements n .. 2*n-1 are internal FFTPACK work space.
285 fftpack_rfftf1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
288 * FFTPACK has a slightly more compact storage than we, time to
289 * convert it: ove most of the array one step up to make room for
290 * zero imaginary parts.
292 p2 = (real *)out_data;
293 for (i = n-1; i > 0; i--)
297 /* imaginary zero freq. */
307 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
309 /* FFTPACK only does in-place transforms, and we cannot just copy
310 * input to output first here since our real array is smaller than
311 * the complex one. However, since the FFTPACK complex storage format
312 * is more compact than ours (2 reals) it will fit, so compact it
313 * and copy on-the-fly to the output array.
315 p1 = (real *) in_data;
316 p2 = (real *)out_data;
319 for (i = 1; i < n; i++)
323 fftpack_rfftb1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
327 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
336 gmx_fft_2d_real (gmx_fft_t fft,
337 enum gmx_fft_direction dir,
341 int i, j, nx, ny, nyc;
349 /* Number of complex elements in y direction */
352 work = fft->work+4*nx;
354 if (dir == GMX_FFT_REAL_TO_COMPLEX)
356 /* If we are doing an in-place transform the 2D array is already
357 * properly padded by the user, and we are all set.
359 * For out-of-place there is no array padding, but FFTPACK only
360 * does in-place FFTs internally, so we need to start by copying
361 * data from the input to the padded (larger) output array.
363 if (in_data != out_data)
365 p1 = (real *)in_data;
366 p2 = (real *)out_data;
368 for (i = 0; i < nx; i++)
370 for (j = 0; j < ny; j++)
372 p2[i*nyc*2+j] = p1[i*ny+j];
376 data = (t_complex *)out_data;
378 /* y real-to-complex FFTs */
379 for (i = 0; i < nx; i++)
381 gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data+i*nyc, data+i*nyc);
384 /* Transform to get X data in place */
385 gmx_fft_transpose_2d(data, data, nx, nyc);
387 /* Complex-to-complex X FFTs */
388 for (i = 0; i < nyc; i++)
390 gmx_fft_1d(fft, GMX_FFT_FORWARD, data+i*nx, data+i*nx);
394 gmx_fft_transpose_2d(data, data, nyc, nx);
397 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
399 /* An in-place complex-to-real transform is straightforward,
400 * since the output array must be large enough for the padding to fit.
402 * For out-of-place complex-to-real transforms we cannot just copy
403 * data to the output array, since it is smaller than the input.
404 * In this case there's nothing to do but employing temporary work data,
405 * starting at work+4*nx and using nx*nyc*2 elements.
407 if (in_data != out_data)
409 memcpy(work, in_data, sizeof(t_complex)*nx*nyc);
410 data = (t_complex *)work;
415 data = (t_complex *)out_data;
418 /* Transpose to get X arrays */
419 gmx_fft_transpose_2d(data, data, nx, nyc);
422 for (i = 0; i < nyc; i++)
424 gmx_fft_1d(fft, GMX_FFT_BACKWARD, data+i*nx, data+i*nx);
427 /* Transpose to get Y arrays */
428 gmx_fft_transpose_2d(data, data, nyc, nx);
431 for (i = 0; i < nx; i++)
433 gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data+i*nyc, data+i*nyc);
436 if (in_data != out_data)
438 /* Output (pointed to by data) is now in padded format.
439 * Pack it into out_data if we were doing an out-of-place transform.
442 p2 = (real *)out_data;
444 for (i = 0; i < nx; i++)
446 for (j = 0; j < ny; j++)
448 p2[i*ny+j] = p1[i*nyc*2+j];
455 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
463 gmx_fft_destroy(gmx_fft_t fft)
468 if (fft->next != NULL)
470 gmx_fft_destroy(fft->next);
476 void gmx_fft_cleanup()
480 const char *gmx_fft_get_version_info()
482 return "fftpack (built-in)";