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"
34 * Contents of the FFTPACK fft datatype.
36 * Note that this is one of several possible implementations of gmx_fft_t.
38 * FFTPACK only does 1d transforms, so we use a pointers to another fft for
39 * the transform in the next dimension.
40 * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
41 * a pointer to a 1d. The 1d structure has next==NULL.
44 struct gmx_fft_fftpack
49 int ndim; /**< Dimensions, including our subdimensions. */
50 int n; /**< Number of points in this dimension. */
51 int ifac[15]; /**< 15 bytes needed for cfft and rfft */
52 struct gmx_fft *next; /**< Pointer to next dimension, or NULL. */
53 real * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
61 gmx_fft_init_1d(gmx_fft_t * pfft,
69 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
74 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
82 /* Need 4*n storage for 1D complex FFT */
83 if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
91 fftpack_cffti1(nx, fft->work, fft->ifac);
101 gmx_fft_init_1d_real(gmx_fft_t * pfft,
109 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
114 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
122 /* Need 2*n storage for 1D real FFT */
123 if ((fft->work = (real *)malloc(sizeof(real)*(2*nx))) == NULL)
131 fftpack_rffti1(nx, fft->work, fft->ifac);
139 gmx_fft_init_2d_real(gmx_fft_t * pfft,
145 int nyc = (ny/2 + 1);
150 gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
155 /* Create the X transform */
156 if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
163 /* Need 4*nx storage for 1D complex FFT, and another
164 * 2*nx*nyc elements for complex-to-real storage in our high-level routine.
166 if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
171 fftpack_cffti1(nx, fft->work, fft->ifac);
173 /* Create real Y transform as a link from X */
174 if ( (rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
186 gmx_fft_1d (gmx_fft_t fft,
187 enum gmx_fft_direction dir,
199 p1 = (real *)in_data;
200 p2 = (real *)out_data;
205 /* FFTPACK only does in-place transforms, so emulate out-of-place
206 * by copying data to the output array first.
208 if (in_data != out_data)
210 p1 = (real *)in_data;
211 p2 = (real *)out_data;
213 /* n complex = 2*n real elements */
214 for (i = 0; i < 2*n; i++)
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.
224 if (dir == GMX_FFT_FORWARD)
226 fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, -1);
228 else if (dir == GMX_FFT_BACKWARD)
230 fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, 1);
234 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
244 gmx_fft_1d_real (gmx_fft_t fft,
245 enum gmx_fft_direction dir,
257 p1 = (real *)in_data;
258 p2 = (real *)out_data;
260 if (dir == GMX_FFT_REAL_TO_COMPLEX)
266 if (dir == GMX_FFT_REAL_TO_COMPLEX)
268 /* FFTPACK only does in-place transforms, so emulate out-of-place
269 * by copying data to the output array first. This works fine, since
270 * the complex array must be larger than the real.
272 if (in_data != out_data)
274 p1 = (real *)in_data;
275 p2 = (real *)out_data;
277 for (i = 0; i < 2*(n/2+1); i++)
283 /* Elements 0 .. n-1 in work are used for ffac values,
284 * Elements n .. 2*n-1 are internal FFTPACK work space.
286 fftpack_rfftf1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
289 * FFTPACK has a slightly more compact storage than we, time to
290 * convert it: ove most of the array one step up to make room for
291 * zero imaginary parts.
293 p2 = (real *)out_data;
294 for (i = n-1; i > 0; i--)
298 /* imaginary zero freq. */
308 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
310 /* FFTPACK only does in-place transforms, and we cannot just copy
311 * input to output first here since our real array is smaller than
312 * the complex one. However, since the FFTPACK complex storage format
313 * is more compact than ours (2 reals) it will fit, so compact it
314 * and copy on-the-fly to the output array.
316 p1 = (real *) in_data;
317 p2 = (real *)out_data;
320 for (i = 1; i < n; i++)
324 fftpack_rfftb1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
328 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
337 gmx_fft_2d_real (gmx_fft_t fft,
338 enum gmx_fft_direction dir,
342 int i, j, nx, ny, nyc;
350 /* Number of complex elements in y direction */
353 work = fft->work+4*nx;
355 if (dir == GMX_FFT_REAL_TO_COMPLEX)
357 /* If we are doing an in-place transform the 2D array is already
358 * properly padded by the user, and we are all set.
360 * For out-of-place there is no array padding, but FFTPACK only
361 * does in-place FFTs internally, so we need to start by copying
362 * data from the input to the padded (larger) output array.
364 if (in_data != out_data)
366 p1 = (real *)in_data;
367 p2 = (real *)out_data;
369 for (i = 0; i < nx; i++)
371 for (j = 0; j < ny; j++)
373 p2[i*nyc*2+j] = p1[i*ny+j];
377 data = (t_complex *)out_data;
379 /* y real-to-complex FFTs */
380 for (i = 0; i < nx; i++)
382 gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data+i*nyc, data+i*nyc);
385 /* Transform to get X data in place */
386 gmx_fft_transpose_2d(data, data, nx, nyc);
388 /* Complex-to-complex X FFTs */
389 for (i = 0; i < nyc; i++)
391 gmx_fft_1d(fft, GMX_FFT_FORWARD, data+i*nx, data+i*nx);
395 gmx_fft_transpose_2d(data, data, nyc, nx);
398 else if (dir == GMX_FFT_COMPLEX_TO_REAL)
400 /* An in-place complex-to-real transform is straightforward,
401 * since the output array must be large enough for the padding to fit.
403 * For out-of-place complex-to-real transforms we cannot just copy
404 * data to the output array, since it is smaller than the input.
405 * In this case there's nothing to do but employing temporary work data,
406 * starting at work+4*nx and using nx*nyc*2 elements.
408 if (in_data != out_data)
410 memcpy(work, in_data, sizeof(t_complex)*nx*nyc);
411 data = (t_complex *)work;
416 data = (t_complex *)out_data;
419 /* Transpose to get X arrays */
420 gmx_fft_transpose_2d(data, data, nx, nyc);
423 for (i = 0; i < nyc; i++)
425 gmx_fft_1d(fft, GMX_FFT_BACKWARD, data+i*nx, data+i*nx);
428 /* Transpose to get Y arrays */
429 gmx_fft_transpose_2d(data, data, nyc, nx);
432 for (i = 0; i < nx; i++)
434 gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data+i*nyc, data+i*nyc);
437 if (in_data != out_data)
439 /* Output (pointed to by data) is now in padded format.
440 * Pack it into out_data if we were doing an out-of-place transform.
443 p2 = (real *)out_data;
445 for (i = 0; i < nx; i++)
447 for (j = 0; j < ny; j++)
449 p2[i*ny+j] = p1[i*nyc*2+j];
456 gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
464 gmx_fft_destroy(gmx_fft_t fft)
469 if (fft->next != NULL)
471 gmx_fft_destroy(fft->next);
477 void gmx_fft_cleanup()
481 const char *gmx_fft_get_version_info()
483 return "fftpack (built-in)";