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.
38 #include "gromacs/fft/fft.h"
47 #include "gromacs/math/gmxcomplex.h"
48 #include "gromacs/utility/fatalerror.h"
49 #include "gromacs/utility/real.h"
51 /* This file contains common fft utility functions, but not
52 * the actual transform implementations. Check the
53 * files like fft_fftw3.c or fft_mkl.c for that.
64 typedef struct gmx_many_fft* gmx_many_fft_t;
67 gmx_fft_init_many_1d(gmx_fft_t * pfft,
75 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
80 if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
85 gmx_fft_init_1d(&fft->fft, nx, flags);
86 fft->howmany = howmany;
89 *pfft = (gmx_fft_t)fft;
94 gmx_fft_init_many_1d_real(gmx_fft_t * pfft,
102 gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
107 if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
112 gmx_fft_init_1d_real(&fft->fft, nx, flags);
113 fft->howmany = howmany;
114 fft->dist = 2*(nx/2+1);
116 *pfft = (gmx_fft_t)fft;
121 gmx_fft_many_1d (gmx_fft_t fft,
122 enum gmx_fft_direction dir,
126 gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
128 for (i = 0; i < mfft->howmany; i++)
130 ret = gmx_fft_1d(mfft->fft, dir, in_data, out_data);
135 in_data = (real*)in_data+mfft->dist;
136 out_data = (real*)out_data+mfft->dist;
142 gmx_fft_many_1d_real (gmx_fft_t fft,
143 enum gmx_fft_direction dir,
147 gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
149 for (i = 0; i < mfft->howmany; i++)
151 ret = gmx_fft_1d_real(mfft->fft, dir, in_data, out_data);
156 in_data = (real*)in_data+mfft->dist;
157 out_data = (real*)out_data+mfft->dist;
164 gmx_many_fft_destroy(gmx_fft_t fft)
166 gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
169 if (mfft->fft != NULL)
171 gmx_fft_destroy(mfft->fft);
177 #endif //not GMX_FFT_FFTW3
179 int gmx_fft_transpose_2d(t_complex * in_data,
180 t_complex * out_data,
184 int i, j, k, im, n, ncount, done1, done2;
185 int i1, i1c, i2, i2c, kmi, max;
187 t_complex tmp1, tmp2, tmp3;
190 /* Use 500 bytes on stack to indicate moves.
191 * This is just for optimization, it does not limit any dimensions.
196 if (nx < 2 || ny < 2)
198 if (in_data != out_data)
200 memcpy(out_data, in_data, sizeof(t_complex)*nx*ny);
205 /* Out-of-place transposes are easy */
206 if (in_data != out_data)
208 for (i = 0; i < nx; i++)
210 for (j = 0; j < ny; j++)
212 out_data[j*nx+i].re = in_data[i*ny+j].re;
213 out_data[j*nx+i].im = in_data[i*ny+j].im;
219 /* In-place transform. in_data=out_data=data */
224 /* trivial case, just swap elements */
225 for (i = 0; i < nx; i++)
227 for (j = i+1; j < nx; j++)
229 tmp1.re = data[i*nx+j].re;
230 tmp1.im = data[i*nx+j].im;
231 data[i*nx+j].re = data[j*nx+i].re;
232 data[i*nx+j].im = data[j*nx+i].im;
233 data[j*nx+i].re = tmp1.re;
234 data[j*nx+i].im = tmp1.im;
240 for (i = 0; i < nmove; i++)
247 if (nx > 2 && ny > 2)
271 tmp1.re = data[i1].re;
272 tmp1.im = data[i1].im;
274 tmp2.re = data[i1c].re;
275 tmp2.im = data[i1c].im;
280 i2 = ny*i1-k*(i1/nx);
307 data[i1].re = data[i2].re;
308 data[i1].im = data[i2].im;
309 data[i1c].re = data[i2c].re;
310 data[i1c].im = data[i2c].im;
317 data[i1].re = tmp1.re;
318 data[i1].im = tmp1.im;
319 data[i1c].re = tmp2.re;
320 data[i1c].im = tmp2.im;
343 while (i2 > i && i2 < max)
346 i2 = ny*i1-k*(i1/nx);