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
);