2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2003 Erik Lindahl, David van der Spoel, 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.
36 #include "gromacs/fft/fft.h"
45 #include "gromacs/math/gmxcomplex.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/real.h"
49 /* This file contains common fft utility functions, but not
50 * the actual transform implementations. Check the
51 * files like fft_fftw3.c or fft_mkl.c for that.
62 typedef struct gmx_many_fft* gmx_many_fft_t;
65 gmx_fft_init_many_1d(gmx_fft_t * pfft,
73 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
78 if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
83 gmx_fft_init_1d(&fft->fft, nx, flags);
84 fft->howmany = howmany;
87 *pfft = (gmx_fft_t)fft;
92 gmx_fft_init_many_1d_real(gmx_fft_t * pfft,
100 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
105 if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
110 gmx_fft_init_1d_real(&fft->fft, nx, flags);
111 fft->howmany = howmany;
112 fft->dist = 2*(nx/2+1);
114 *pfft = (gmx_fft_t)fft;
119 gmx_fft_many_1d (gmx_fft_t fft,
120 enum gmx_fft_direction dir,
124 gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
126 for (i = 0; i < mfft->howmany; i++)
128 ret = gmx_fft_1d(mfft->fft, dir, in_data, out_data);
133 in_data = (real*)in_data+mfft->dist;
134 out_data = (real*)out_data+mfft->dist;
140 gmx_fft_many_1d_real (gmx_fft_t fft,
141 enum gmx_fft_direction dir,
145 gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
147 for (i = 0; i < mfft->howmany; i++)
149 ret = gmx_fft_1d_real(mfft->fft, dir, in_data, out_data);
154 in_data = (real*)in_data+mfft->dist;
155 out_data = (real*)out_data+mfft->dist;
162 gmx_many_fft_destroy(gmx_fft_t fft)
164 gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
167 if (mfft->fft != NULL)
169 gmx_fft_destroy(mfft->fft);
175 #endif //not GMX_FFT_FFTW3
177 int gmx_fft_transpose_2d(t_complex * in_data,
178 t_complex * out_data,
182 int i, j, k, im, n, ncount, done1, done2;
183 int i1, i1c, i2, i2c, kmi, max;
185 t_complex tmp1, tmp2, tmp3;
188 /* Use 500 bytes on stack to indicate moves.
189 * This is just for optimization, it does not limit any dimensions.
194 if (nx < 2 || ny < 2)
196 if (in_data != out_data)
198 memcpy(out_data, in_data, sizeof(t_complex)*nx*ny);
203 /* Out-of-place transposes are easy */
204 if (in_data != out_data)
206 for (i = 0; i < nx; i++)
208 for (j = 0; j < ny; j++)
210 out_data[j*nx+i].re = in_data[i*ny+j].re;
211 out_data[j*nx+i].im = in_data[i*ny+j].im;
217 /* In-place transform. in_data=out_data=data */
222 /* trivial case, just swap elements */
223 for (i = 0; i < nx; i++)
225 for (j = i+1; j < nx; j++)
227 tmp1.re = data[i*nx+j].re;
228 tmp1.im = data[i*nx+j].im;
229 data[i*nx+j].re = data[j*nx+i].re;
230 data[i*nx+j].im = data[j*nx+i].im;
231 data[j*nx+i].re = tmp1.re;
232 data[j*nx+i].im = tmp1.im;
238 for (i = 0; i < nmove; i++)
245 if (nx > 2 && ny > 2)
269 tmp1.re = data[i1].re;
270 tmp1.im = data[i1].im;
272 tmp2.re = data[i1c].re;
273 tmp2.im = data[i1c].im;
278 i2 = ny*i1-k*(i1/nx);
305 data[i1].re = data[i2].re;
306 data[i1].im = data[i2].im;
307 data[i1c].re = data[i2c].re;
308 data[i1c].im = data[i2c].im;
315 data[i1].re = tmp1.re;
316 data[i1].im = tmp1.im;
317 data[i1c].re = tmp2.re;
318 data[i1c].im = tmp2.im;
341 while (i2 > i && i2 < max)
344 i2 = ny*i1-k*(i1/nx);