2 * This file is part of the GROMACS molecular simulation package.
4 * Gromacs 4.0 Copyright (c) 1991-2003
5 * David van der Spoel, Erik Lindahl, University of Groningen.
6 * Copyright (c) 2012,2013, by the GROMACS development team, led by
7 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
8 * others, as listed in the AUTHORS file in the top-level source
9 * directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gmx_fatal.h"
53 #define FFTWPREFIX(name) fftw_ ## name
55 #define FFTWPREFIX(name) fftwf_ ## name
59 #include "thread_mpi/threads.h"
64 /* none of the fftw3 calls, except execute(), are thread-safe, so
65 we need to serialize them with this mutex. */
66 static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
67 static gmx_bool gmx_fft_threads_initialized=FALSE;
68 #define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex)
69 #define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex)
70 #else /* GMX_THREAD_MPI */
73 #endif /* GMX_THREAD_MPI */
75 /* We assume here that aligned memory starts at multiple of 16 bytes and unaligned memory starts at multiple of 8 bytes. The later is guranteed for all malloc implementation.
77 - It is not allowed to use these FFT plans from memory which doesn't have a starting address as a multiple of 8 bytes.
78 This is OK as long as the memory directly comes from malloc and is not some subarray within alloated memory.
79 - This has to be fixed if any future architecute requires memory to be aligned to multiples of 32 bytes.
84 /* Three alternatives (unaligned/aligned, out-of-place/in-place, forward/backward)
85 * results in 8 different FFTW plans. Keep track of them with 3 array indices:
86 * first index: 0=unaligned, 1=aligned
87 * second index: 0=out-of-place, 1=in-place
88 * third index: 0=backward, 1=forward
90 FFTWPREFIX(plan) plan[2][2][2];
91 /* Catch user mistakes */
99 gmx_fft_init_1d(gmx_fft_t * pfft,
103 return gmx_fft_init_many_1d(pfft,nx,1,flags);
108 gmx_fft_init_many_1d(gmx_fft_t * pfft,
114 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
119 #ifdef GMX_DISABLE_FFTW_MEASURE
120 flags |= GMX_FFT_FLAG_CONSERVATIVE;
123 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
127 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
133 if( (fft = (gmx_fft_t)FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
139 /* allocate aligned, and extra memory to make it unaligned */
140 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
143 FFTWPREFIX(free)(fft);
148 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
151 FFTWPREFIX(free)(p1);
152 FFTWPREFIX(free)(fft);
157 /* make unaligned pointers.
158 * In double precision the actual complex datatype will be 16 bytes,
159 * so go to a char pointer and force an offset of 8 bytes instead.
163 up1 = (FFTWPREFIX(complex) *)pc;
167 up2 = (FFTWPREFIX(complex) *)pc;
169 /* int rank, const int *n, int howmany,
170 fftw_complex *in, const int *inembed,
171 int istride, int idist,
172 fftw_complex *out, const int *onembed,
173 int ostride, int odist,
174 int sign, unsigned flags */
175 fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
176 fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
177 fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
178 fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
179 fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
180 fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
181 fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
182 fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
190 if(fft->plan[i][j][k] == NULL)
192 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
194 gmx_fft_destroy(fft);
196 FFTWPREFIX(free)(p1);
197 FFTWPREFIX(free)(p2);
205 FFTWPREFIX(free)(p1);
206 FFTWPREFIX(free)(p2);
208 fft->real_transform = 0;
218 gmx_fft_init_1d_real(gmx_fft_t * pfft,
222 return gmx_fft_init_many_1d_real(pfft, nx, 1, flags);
226 gmx_fft_init_many_1d_real(gmx_fft_t * pfft,
232 real *p1,*p2,*up1,*up2;
237 #ifdef GMX_DISABLE_FFTW_MEASURE
238 flags |= GMX_FFT_FLAG_CONSERVATIVE;
241 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
245 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
251 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
257 /* allocate aligned, and extra memory to make it unaligned */
258 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
261 FFTWPREFIX(free)(fft);
266 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
269 FFTWPREFIX(free)(p1);
270 FFTWPREFIX(free)(fft);
275 /* make unaligned pointers.
276 * In double precision the actual complex datatype will be 16 bytes,
277 * so go to a char pointer and force an offset of 8 bytes instead.
287 /* int rank, const int *n, int howmany,
288 double *in, const int *inembed,
289 int istride, int idist,
290 fftw_complex *out, const int *onembed,
291 int ostride, int odist,
293 fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany,up1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)up2,0,1,(nx/2+1),fftw_flags);
294 fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany,up1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)up1,0,1,(nx/2+1),fftw_flags);
295 fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany, p1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)p2 ,0,1,(nx/2+1),fftw_flags);
296 fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft_r2c)(1,&nx,howmany, p1,0,1,(nx/2+1)*2,(FFTWPREFIX(complex) *)p1 ,0,1,(nx/2+1),fftw_flags);
298 fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *)up1,0,1,(nx/2+1),up2,0,1,(nx/2+1)*2,fftw_flags);
299 fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *)up1,0,1,(nx/2+1),up1,0,1,(nx/2+1)*2,fftw_flags);
300 fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *) p1,0,1,(nx/2+1), p2,0,1,(nx/2+1)*2,fftw_flags);
301 fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft_c2r)(1,&nx,howmany,(FFTWPREFIX(complex) *) p1,0,1,(nx/2+1), p1,0,1,(nx/2+1)*2,fftw_flags);
309 if(fft->plan[i][j][k] == NULL)
311 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
313 gmx_fft_destroy(fft);
315 FFTWPREFIX(free)(p1);
316 FFTWPREFIX(free)(p2);
324 FFTWPREFIX(free)(p1);
325 FFTWPREFIX(free)(p2);
327 fft->real_transform = 1;
338 gmx_fft_init_2d(gmx_fft_t * pfft,
344 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
349 #ifdef GMX_DISABLE_FFTW_MEASURE
350 flags |= GMX_FFT_FLAG_CONSERVATIVE;
353 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
357 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
363 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
369 /* allocate aligned, and extra memory to make it unaligned */
370 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
373 FFTWPREFIX(free)(fft);
378 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
381 FFTWPREFIX(free)(p1);
382 FFTWPREFIX(free)(fft);
387 /* make unaligned pointers.
388 * In double precision the actual complex datatype will be 16 bytes,
389 * so go to a char pointer and force an offset of 8 bytes instead.
393 up1 = (FFTWPREFIX(complex) *)pc;
397 up2 = (FFTWPREFIX(complex) *)pc;
400 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_BACKWARD,fftw_flags);
401 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_FORWARD,fftw_flags);
402 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_BACKWARD,fftw_flags);
403 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_FORWARD,fftw_flags);
405 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_BACKWARD,fftw_flags);
406 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_FORWARD,fftw_flags);
407 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_BACKWARD,fftw_flags);
408 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_FORWARD,fftw_flags);
417 if(fft->plan[i][j][k] == NULL)
419 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
421 gmx_fft_destroy(fft);
423 FFTWPREFIX(free)(p1);
424 FFTWPREFIX(free)(p2);
432 FFTWPREFIX(free)(p1);
433 FFTWPREFIX(free)(p2);
435 fft->real_transform = 0;
446 gmx_fft_init_2d_real(gmx_fft_t * pfft,
452 real *p1,*p2,*up1,*up2;
457 #ifdef GMX_DISABLE_FFTW_MEASURE
458 flags |= GMX_FFT_FLAG_CONSERVATIVE;
461 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
465 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
471 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
477 /* allocate aligned, and extra memory to make it unaligned */
478 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
481 FFTWPREFIX(free)(fft);
486 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
489 FFTWPREFIX(free)(p1);
490 FFTWPREFIX(free)(fft);
495 /* make unaligned pointers.
496 * In double precision the actual complex datatype will be 16 bytes,
497 * so go to a char pointer and force an offset of 8 bytes instead.
508 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
509 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
510 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
511 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
513 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
514 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
515 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
516 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
525 if(fft->plan[i][j][k] == NULL)
527 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
529 gmx_fft_destroy(fft);
531 FFTWPREFIX(free)(p1);
532 FFTWPREFIX(free)(p2);
540 FFTWPREFIX(free)(p1);
541 FFTWPREFIX(free)(p2);
543 fft->real_transform = 1;
554 gmx_fft_init_3d(gmx_fft_t * pfft,
561 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
566 #ifdef GMX_DISABLE_FFTW_MEASURE
567 flags |= GMX_FFT_FLAG_CONSERVATIVE;
570 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
574 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
580 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
586 /* allocate aligned, and extra memory to make it unaligned */
587 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
590 FFTWPREFIX(free)(fft);
595 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
598 FFTWPREFIX(free)(p1);
599 FFTWPREFIX(free)(fft);
604 /* make unaligned pointers.
605 * In double precision the actual complex datatype will be 16 bytes,
606 * so go to a char pointer and force an offset of 8 bytes instead.
610 up1 = (FFTWPREFIX(complex) *)pc;
614 up2 = (FFTWPREFIX(complex) *)pc;
617 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_BACKWARD,fftw_flags);
618 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_FORWARD,fftw_flags);
619 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_BACKWARD,fftw_flags);
620 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_FORWARD,fftw_flags);
622 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_BACKWARD,fftw_flags);
623 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_FORWARD,fftw_flags);
624 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_BACKWARD,fftw_flags);
625 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_FORWARD,fftw_flags);
634 if(fft->plan[i][j][k] == NULL)
636 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
638 gmx_fft_destroy(fft);
640 FFTWPREFIX(free)(p1);
641 FFTWPREFIX(free)(p2);
649 FFTWPREFIX(free)(p1);
650 FFTWPREFIX(free)(p2);
652 fft->real_transform = 0;
663 gmx_fft_init_3d_real(gmx_fft_t * pfft,
670 real *p1,*p2,*up1,*up2;
675 #ifdef GMX_DISABLE_FFTW_MEASURE
676 flags |= GMX_FFT_FLAG_CONSERVATIVE;
679 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
683 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
689 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
695 /* allocate aligned, and extra memory to make it unaligned */
696 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
699 FFTWPREFIX(free)(fft);
704 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
707 FFTWPREFIX(free)(p1);
708 FFTWPREFIX(free)(fft);
713 /* make unaligned pointers.
714 * In double precision the actual complex datatype will be 16 bytes,
715 * so go to a void pointer and force an offset of 8 bytes instead.
726 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
727 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
728 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
729 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
731 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
732 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
733 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
734 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
743 if(fft->plan[i][j][k] == NULL)
745 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
747 gmx_fft_destroy(fft);
749 FFTWPREFIX(free)(p1);
750 FFTWPREFIX(free)(p2);
758 FFTWPREFIX(free)(p1);
759 FFTWPREFIX(free)(p2);
761 fft->real_transform = 1;
771 gmx_fft_1d (gmx_fft_t fft,
772 enum gmx_fft_direction dir,
776 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
777 int inplace = (in_data == out_data);
778 int isforward = (dir == GMX_FFT_FORWARD);
781 if( (fft->real_transform == 1) || (fft->ndim != 1) ||
782 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
784 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
788 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
789 (FFTWPREFIX(complex) *)in_data,
790 (FFTWPREFIX(complex) *)out_data);
796 gmx_fft_many_1d (gmx_fft_t fft,
797 enum gmx_fft_direction dir,
801 return gmx_fft_1d(fft,dir,in_data,out_data);
805 gmx_fft_1d_real (gmx_fft_t fft,
806 enum gmx_fft_direction dir,
810 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
811 int inplace = (in_data == out_data);
812 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
815 if( (fft->real_transform != 1) || (fft->ndim != 1) ||
816 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
818 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
824 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
825 (real *)in_data,(FFTWPREFIX(complex) *)out_data);
829 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
830 (FFTWPREFIX(complex) *)in_data,(real *)out_data);
837 gmx_fft_many_1d_real (gmx_fft_t fft,
838 enum gmx_fft_direction dir,
842 return gmx_fft_1d_real(fft,dir,in_data,out_data);
846 gmx_fft_2d (gmx_fft_t fft,
847 enum gmx_fft_direction dir,
851 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
852 int inplace = (in_data == out_data);
853 int isforward = (dir == GMX_FFT_FORWARD);
856 if( (fft->real_transform == 1) || (fft->ndim != 2) ||
857 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
859 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
863 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
864 (FFTWPREFIX(complex) *)in_data,
865 (FFTWPREFIX(complex) *)out_data);
872 gmx_fft_2d_real (gmx_fft_t fft,
873 enum gmx_fft_direction dir,
877 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
878 int inplace = (in_data == out_data);
879 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
882 if( (fft->real_transform != 1) || (fft->ndim != 2) ||
883 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
885 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
891 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
893 (FFTWPREFIX(complex) *)out_data);
897 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
898 (FFTWPREFIX(complex) *)in_data,
908 gmx_fft_3d (gmx_fft_t fft,
909 enum gmx_fft_direction dir,
913 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
914 int inplace = (in_data == out_data);
915 int isforward = (dir == GMX_FFT_FORWARD);
918 if( (fft->real_transform == 1) || (fft->ndim != 3) ||
919 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
921 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
925 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
926 (FFTWPREFIX(complex) *)in_data,
927 (FFTWPREFIX(complex) *)out_data);
934 gmx_fft_3d_real (gmx_fft_t fft,
935 enum gmx_fft_direction dir,
939 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
940 int inplace = (in_data == out_data);
941 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
944 if( (fft->real_transform != 1) || (fft->ndim != 3) ||
945 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
947 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
953 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
955 (FFTWPREFIX(complex) *)out_data);
959 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
960 (FFTWPREFIX(complex) *)in_data,
970 gmx_fft_destroy(gmx_fft_t fft)
982 if(fft->plan[i][j][k] != NULL)
985 FFTWPREFIX(destroy_plan)(fft->plan[i][j][k]);
987 fft->plan[i][j][k] = NULL;
993 FFTWPREFIX(free)(fft);
1000 gmx_many_fft_destroy(gmx_fft_t fft)
1002 gmx_fft_destroy(fft);
1007 gmx_fft_fftw3_empty;
1008 #endif /* GMX_FFT_FFTW3 */