3 * Gromacs 4.0 Copyright (c) 1991-2003
4 * David van der Spoel, Erik Lindahl, University of Groningen.
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU General Public License
8 * as published by the Free Software Foundation; either version 2
9 * of the License, or (at your option) any later version.
11 * To help us fund GROMACS development, we humbly ask that you cite
12 * the research papers on the package. Check out http://www.gromacs.org
15 * Gnomes, ROck Monsters And Chili Sauce
29 #include "gmx_fatal.h"
33 #define FFTWPREFIX(name) fftw_ ## name
35 #define FFTWPREFIX(name) fftwf_ ## name
39 #include "thread_mpi/threads.h"
44 /* none of the fftw3 calls, except execute(), are thread-safe, so
45 we need to serialize them with this mutex. */
46 static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
47 static gmx_bool gmx_fft_threads_initialized=FALSE;
48 #define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex)
49 #define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex)
50 #else /* GMX_THREAD_MPI */
53 #endif /* GMX_THREAD_MPI */
55 /* 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.
57 - It is not allowed to use these FFT plans from memory which doesn't have a starting address as a multiple of 8 bytes.
58 This is OK as long as the memory directly comes from malloc and is not some subarray within alloated memory.
59 - This has to be fixed if any future architecute requires memory to be aligned to multiples of 32 bytes.
64 /* Three alternatives (unaligned/aligned, out-of-place/in-place, forward/backward)
65 * results in 8 different FFTW plans. Keep track of them with 3 array indices:
66 * first index: 0=unaligned, 1=aligned
67 * second index: 0=out-of-place, 1=in-place
68 * third index: 0=backward, 1=forward
70 FFTWPREFIX(plan) plan[2][2][2];
71 /* Catch user mistakes */
79 gmx_fft_init_1d(gmx_fft_t * pfft,
83 return gmx_fft_init_many_1d(pfft,nx,1,flags);
88 gmx_fft_init_many_1d(gmx_fft_t * pfft,
94 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
99 #ifdef GMX_DISABLE_FFTW_MEASURE
100 flags |= GMX_FFT_FLAG_CONSERVATIVE;
103 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
107 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
113 if( (fft = (gmx_fft_t)FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
119 /* allocate aligned, and extra memory to make it unaligned */
120 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
123 FFTWPREFIX(free)(fft);
128 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
131 FFTWPREFIX(free)(p1);
132 FFTWPREFIX(free)(fft);
137 /* make unaligned pointers.
138 * In double precision the actual complex datatype will be 16 bytes,
139 * so go to a char pointer and force an offset of 8 bytes instead.
143 up1 = (FFTWPREFIX(complex) *)pc;
147 up2 = (FFTWPREFIX(complex) *)pc;
149 /* int rank, const int *n, int howmany,
150 fftw_complex *in, const int *inembed,
151 int istride, int idist,
152 fftw_complex *out, const int *onembed,
153 int ostride, int odist,
154 int sign, unsigned flags */
155 fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
156 fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
157 fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
158 fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
159 fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
160 fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
161 fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
162 fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
170 if(fft->plan[i][j][k] == NULL)
172 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
174 gmx_fft_destroy(fft);
176 FFTWPREFIX(free)(p1);
177 FFTWPREFIX(free)(p2);
185 FFTWPREFIX(free)(p1);
186 FFTWPREFIX(free)(p2);
188 fft->real_transform = 0;
198 gmx_fft_init_1d_real(gmx_fft_t * pfft,
202 return gmx_fft_init_many_1d_real(pfft, nx, 1, flags);
206 gmx_fft_init_many_1d_real(gmx_fft_t * pfft,
212 real *p1,*p2,*up1,*up2;
217 #ifdef GMX_DISABLE_FFTW_MEASURE
218 flags |= GMX_FFT_FLAG_CONSERVATIVE;
221 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
225 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
231 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
237 /* allocate aligned, and extra memory to make it unaligned */
238 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
241 FFTWPREFIX(free)(fft);
246 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
249 FFTWPREFIX(free)(p1);
250 FFTWPREFIX(free)(fft);
255 /* make unaligned pointers.
256 * In double precision the actual complex datatype will be 16 bytes,
257 * so go to a char pointer and force an offset of 8 bytes instead.
267 /* int rank, const int *n, int howmany,
268 double *in, const int *inembed,
269 int istride, int idist,
270 fftw_complex *out, const int *onembed,
271 int ostride, int odist,
273 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);
274 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);
275 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);
276 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);
278 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);
279 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);
280 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);
281 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);
289 if(fft->plan[i][j][k] == NULL)
291 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
293 gmx_fft_destroy(fft);
295 FFTWPREFIX(free)(p1);
296 FFTWPREFIX(free)(p2);
304 FFTWPREFIX(free)(p1);
305 FFTWPREFIX(free)(p2);
307 fft->real_transform = 1;
318 gmx_fft_init_2d(gmx_fft_t * pfft,
324 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
329 #ifdef GMX_DISABLE_FFTW_MEASURE
330 flags |= GMX_FFT_FLAG_CONSERVATIVE;
333 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
337 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
343 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
349 /* allocate aligned, and extra memory to make it unaligned */
350 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
353 FFTWPREFIX(free)(fft);
358 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
361 FFTWPREFIX(free)(p1);
362 FFTWPREFIX(free)(fft);
367 /* make unaligned pointers.
368 * In double precision the actual complex datatype will be 16 bytes,
369 * so go to a char pointer and force an offset of 8 bytes instead.
373 up1 = (FFTWPREFIX(complex) *)pc;
377 up2 = (FFTWPREFIX(complex) *)pc;
380 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_BACKWARD,fftw_flags);
381 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_FORWARD,fftw_flags);
382 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_BACKWARD,fftw_flags);
383 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_FORWARD,fftw_flags);
385 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_BACKWARD,fftw_flags);
386 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_FORWARD,fftw_flags);
387 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_BACKWARD,fftw_flags);
388 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_FORWARD,fftw_flags);
397 if(fft->plan[i][j][k] == NULL)
399 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
401 gmx_fft_destroy(fft);
403 FFTWPREFIX(free)(p1);
404 FFTWPREFIX(free)(p2);
412 FFTWPREFIX(free)(p1);
413 FFTWPREFIX(free)(p2);
415 fft->real_transform = 0;
426 gmx_fft_init_2d_real(gmx_fft_t * pfft,
432 real *p1,*p2,*up1,*up2;
437 #ifdef GMX_DISABLE_FFTW_MEASURE
438 flags |= GMX_FFT_FLAG_CONSERVATIVE;
441 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
445 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
451 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
457 /* allocate aligned, and extra memory to make it unaligned */
458 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
461 FFTWPREFIX(free)(fft);
466 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
469 FFTWPREFIX(free)(p1);
470 FFTWPREFIX(free)(fft);
475 /* make unaligned pointers.
476 * In double precision the actual complex datatype will be 16 bytes,
477 * so go to a char pointer and force an offset of 8 bytes instead.
488 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
489 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
490 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
491 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
493 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
494 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
495 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
496 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
505 if(fft->plan[i][j][k] == NULL)
507 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
509 gmx_fft_destroy(fft);
511 FFTWPREFIX(free)(p1);
512 FFTWPREFIX(free)(p2);
520 FFTWPREFIX(free)(p1);
521 FFTWPREFIX(free)(p2);
523 fft->real_transform = 1;
534 gmx_fft_init_3d(gmx_fft_t * pfft,
541 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
546 #ifdef GMX_DISABLE_FFTW_MEASURE
547 flags |= GMX_FFT_FLAG_CONSERVATIVE;
550 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
554 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
560 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
566 /* allocate aligned, and extra memory to make it unaligned */
567 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
570 FFTWPREFIX(free)(fft);
575 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
578 FFTWPREFIX(free)(p1);
579 FFTWPREFIX(free)(fft);
584 /* make unaligned pointers.
585 * In double precision the actual complex datatype will be 16 bytes,
586 * so go to a char pointer and force an offset of 8 bytes instead.
590 up1 = (FFTWPREFIX(complex) *)pc;
594 up2 = (FFTWPREFIX(complex) *)pc;
597 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_BACKWARD,fftw_flags);
598 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_FORWARD,fftw_flags);
599 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_BACKWARD,fftw_flags);
600 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_FORWARD,fftw_flags);
602 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_BACKWARD,fftw_flags);
603 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_FORWARD,fftw_flags);
604 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_BACKWARD,fftw_flags);
605 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_FORWARD,fftw_flags);
614 if(fft->plan[i][j][k] == NULL)
616 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
618 gmx_fft_destroy(fft);
620 FFTWPREFIX(free)(p1);
621 FFTWPREFIX(free)(p2);
629 FFTWPREFIX(free)(p1);
630 FFTWPREFIX(free)(p2);
632 fft->real_transform = 0;
643 gmx_fft_init_3d_real(gmx_fft_t * pfft,
650 real *p1,*p2,*up1,*up2;
655 #ifdef GMX_DISABLE_FFTW_MEASURE
656 flags |= GMX_FFT_FLAG_CONSERVATIVE;
659 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
663 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
669 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
675 /* allocate aligned, and extra memory to make it unaligned */
676 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
679 FFTWPREFIX(free)(fft);
684 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
687 FFTWPREFIX(free)(p1);
688 FFTWPREFIX(free)(fft);
693 /* make unaligned pointers.
694 * In double precision the actual complex datatype will be 16 bytes,
695 * so go to a void pointer and force an offset of 8 bytes instead.
706 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
707 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
708 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
709 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
711 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
712 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
713 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
714 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
723 if(fft->plan[i][j][k] == NULL)
725 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
727 gmx_fft_destroy(fft);
729 FFTWPREFIX(free)(p1);
730 FFTWPREFIX(free)(p2);
738 FFTWPREFIX(free)(p1);
739 FFTWPREFIX(free)(p2);
741 fft->real_transform = 1;
751 gmx_fft_1d (gmx_fft_t fft,
752 enum gmx_fft_direction dir,
756 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
757 int inplace = (in_data == out_data);
758 int isforward = (dir == GMX_FFT_FORWARD);
761 if( (fft->real_transform == 1) || (fft->ndim != 1) ||
762 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
764 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
768 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
769 (FFTWPREFIX(complex) *)in_data,
770 (FFTWPREFIX(complex) *)out_data);
776 gmx_fft_many_1d (gmx_fft_t fft,
777 enum gmx_fft_direction dir,
781 return gmx_fft_1d(fft,dir,in_data,out_data);
785 gmx_fft_1d_real (gmx_fft_t fft,
786 enum gmx_fft_direction dir,
790 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
791 int inplace = (in_data == out_data);
792 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
795 if( (fft->real_transform != 1) || (fft->ndim != 1) ||
796 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
798 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
804 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
805 (real *)in_data,(FFTWPREFIX(complex) *)out_data);
809 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
810 (FFTWPREFIX(complex) *)in_data,(real *)out_data);
817 gmx_fft_many_1d_real (gmx_fft_t fft,
818 enum gmx_fft_direction dir,
822 return gmx_fft_1d_real(fft,dir,in_data,out_data);
826 gmx_fft_2d (gmx_fft_t fft,
827 enum gmx_fft_direction dir,
831 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
832 int inplace = (in_data == out_data);
833 int isforward = (dir == GMX_FFT_FORWARD);
836 if( (fft->real_transform == 1) || (fft->ndim != 2) ||
837 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
839 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
843 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
844 (FFTWPREFIX(complex) *)in_data,
845 (FFTWPREFIX(complex) *)out_data);
852 gmx_fft_2d_real (gmx_fft_t fft,
853 enum gmx_fft_direction dir,
857 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
858 int inplace = (in_data == out_data);
859 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
862 if( (fft->real_transform != 1) || (fft->ndim != 2) ||
863 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
865 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
871 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
873 (FFTWPREFIX(complex) *)out_data);
877 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
878 (FFTWPREFIX(complex) *)in_data,
888 gmx_fft_3d (gmx_fft_t fft,
889 enum gmx_fft_direction dir,
893 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
894 int inplace = (in_data == out_data);
895 int isforward = (dir == GMX_FFT_FORWARD);
898 if( (fft->real_transform == 1) || (fft->ndim != 3) ||
899 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
901 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
905 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
906 (FFTWPREFIX(complex) *)in_data,
907 (FFTWPREFIX(complex) *)out_data);
914 gmx_fft_3d_real (gmx_fft_t fft,
915 enum gmx_fft_direction dir,
919 int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
920 int inplace = (in_data == out_data);
921 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
924 if( (fft->real_transform != 1) || (fft->ndim != 3) ||
925 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
927 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
933 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
935 (FFTWPREFIX(complex) *)out_data);
939 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
940 (FFTWPREFIX(complex) *)in_data,
950 gmx_fft_destroy(gmx_fft_t fft)
962 if(fft->plan[i][j][k] != NULL)
965 FFTWPREFIX(destroy_plan)(fft->plan[i][j][k]);
967 fft->plan[i][j][k] = NULL;
973 FFTWPREFIX(free)(fft);
980 gmx_many_fft_destroy(gmx_fft_t fft)
982 gmx_fft_destroy(fft);
988 #endif /* GMX_FFT_FFTW3 */