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
41 /* none of the fftw3 calls, except execute(), are thread-safe, so
42 we need to serialize them with this mutex. */
43 static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
44 static bool gmx_fft_threads_initialized=FALSE;
45 #define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex);
46 #define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex);
47 #else /* GMX_THREADS */
50 #endif /* GMX_THREADS */
54 /* Three alternatives (unaligned/aligned, out-of-place/in-place, forward/backward)
55 * results in 8 different FFTW plans. Keep track of them with 3 array indices:
56 * first index: 0=unaligned, 1=aligned
57 * second index: 0=out-of-place, 1=in-place
58 * third index: 0=backward, 1=forward
60 FFTWPREFIX(plan) plan[2][2][2];
61 /* Catch user mistakes */
69 gmx_fft_init_1d(gmx_fft_t * pfft,
73 return gmx_fft_init_many_1d(pfft,nx,1,flags);
78 gmx_fft_init_many_1d(gmx_fft_t * pfft,
84 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
89 #ifdef GMX_DISABLE_FFTW_MEASURE
90 flags |= GMX_FFT_FLAG_CONSERVATIVE;
93 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
97 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
103 if( (fft = (gmx_fft_t)FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
109 /* allocate aligned, and extra memory to make it unaligned */
110 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
113 FFTWPREFIX(free)(fft);
118 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
121 FFTWPREFIX(free)(p1);
122 FFTWPREFIX(free)(fft);
127 /* make unaligned pointers.
128 * In double precision the actual complex datatype will be 16 bytes,
129 * so go to a char pointer and force an offset of 8 bytes instead.
133 up1 = (FFTWPREFIX(complex) *)pc;
137 up2 = (FFTWPREFIX(complex) *)pc;
139 /* int rank, const int *n, int howmany,
140 fftw_complex *in, const int *inembed,
141 int istride, int idist,
142 fftw_complex *out, const int *onembed,
143 int ostride, int odist,
144 int sign, unsigned flags */
145 fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
146 fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
147 fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
148 fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
149 fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
150 fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
151 fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
152 fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
160 if(fft->plan[i][j][k] == NULL)
162 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
164 gmx_fft_destroy(fft);
166 FFTWPREFIX(free)(p1);
167 FFTWPREFIX(free)(p2);
175 FFTWPREFIX(free)(p1);
176 FFTWPREFIX(free)(p2);
178 fft->real_transform = 0;
188 gmx_fft_init_1d_real(gmx_fft_t * pfft,
192 return gmx_fft_init_many_1d_real(pfft, nx, 1, flags);
196 gmx_fft_init_many_1d_real(gmx_fft_t * pfft,
202 real *p1,*p2,*up1,*up2;
207 #ifdef GMX_DISABLE_FFTW_MEASURE
208 flags |= GMX_FFT_FLAG_CONSERVATIVE;
211 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
215 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
221 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
227 /* allocate aligned, and extra memory to make it unaligned */
228 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
231 FFTWPREFIX(free)(fft);
236 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*(nx/2+1)*2*howmany + 8);
239 FFTWPREFIX(free)(p1);
240 FFTWPREFIX(free)(fft);
245 /* make unaligned pointers.
246 * In double precision the actual complex datatype will be 16 bytes,
247 * so go to a char pointer and force an offset of 8 bytes instead.
257 /* int rank, const int *n, int howmany,
258 double *in, const int *inembed,
259 int istride, int idist,
260 fftw_complex *out, const int *onembed,
261 int ostride, int odist,
263 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);
264 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);
265 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);
266 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);
268 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);
269 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);
270 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);
271 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);
279 if(fft->plan[i][j][k] == NULL)
281 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
283 gmx_fft_destroy(fft);
285 FFTWPREFIX(free)(p1);
286 FFTWPREFIX(free)(p2);
294 FFTWPREFIX(free)(p1);
295 FFTWPREFIX(free)(p2);
297 fft->real_transform = 1;
308 gmx_fft_init_2d(gmx_fft_t * pfft,
314 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
319 #ifdef GMX_DISABLE_FFTW_MEASURE
320 flags |= GMX_FFT_FLAG_CONSERVATIVE;
323 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
327 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
333 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
339 /* allocate aligned, and extra memory to make it unaligned */
340 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
343 FFTWPREFIX(free)(fft);
348 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny+2));
351 FFTWPREFIX(free)(p1);
352 FFTWPREFIX(free)(fft);
357 /* make unaligned pointers.
358 * In double precision the actual complex datatype will be 16 bytes,
359 * so go to a char pointer and force an offset of 8 bytes instead.
363 up1 = (FFTWPREFIX(complex) *)pc;
367 up2 = (FFTWPREFIX(complex) *)pc;
370 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_BACKWARD,fftw_flags);
371 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up2,FFTW_FORWARD,fftw_flags);
372 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_BACKWARD,fftw_flags);
373 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,up1,up1,FFTW_FORWARD,fftw_flags);
375 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_BACKWARD,fftw_flags);
376 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p2,FFTW_FORWARD,fftw_flags);
377 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_BACKWARD,fftw_flags);
378 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_2d)(nx,ny,p1,p1,FFTW_FORWARD,fftw_flags);
387 if(fft->plan[i][j][k] == NULL)
389 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
391 gmx_fft_destroy(fft);
393 FFTWPREFIX(free)(p1);
394 FFTWPREFIX(free)(p2);
402 FFTWPREFIX(free)(p1);
403 FFTWPREFIX(free)(p2);
405 fft->real_transform = 0;
416 gmx_fft_init_2d_real(gmx_fft_t * pfft,
422 real *p1,*p2,*up1,*up2;
427 #ifdef GMX_DISABLE_FFTW_MEASURE
428 flags |= GMX_FFT_FLAG_CONSERVATIVE;
431 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
435 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
441 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
447 /* allocate aligned, and extra memory to make it unaligned */
448 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
451 FFTWPREFIX(free)(fft);
456 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*(ny/2+1)*2 + 2) );
459 FFTWPREFIX(free)(p1);
460 FFTWPREFIX(free)(fft);
465 /* make unaligned pointers.
466 * In double precision the actual complex datatype will be 16 bytes,
467 * so go to a char pointer and force an offset of 8 bytes instead.
478 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
479 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
480 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
481 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
483 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
484 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
485 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_2d)(nx,ny,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
486 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_2d)(nx,ny,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
495 if(fft->plan[i][j][k] == NULL)
497 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
499 gmx_fft_destroy(fft);
501 FFTWPREFIX(free)(p1);
502 FFTWPREFIX(free)(p2);
510 FFTWPREFIX(free)(p1);
511 FFTWPREFIX(free)(p2);
513 fft->real_transform = 1;
524 gmx_fft_init_3d(gmx_fft_t * pfft,
531 FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
536 #ifdef GMX_DISABLE_FFTW_MEASURE
537 flags |= GMX_FFT_FLAG_CONSERVATIVE;
540 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
544 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
550 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
556 /* allocate aligned, and extra memory to make it unaligned */
557 p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
560 FFTWPREFIX(free)(fft);
565 p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx*ny*nz+2));
568 FFTWPREFIX(free)(p1);
569 FFTWPREFIX(free)(fft);
574 /* make unaligned pointers.
575 * In double precision the actual complex datatype will be 16 bytes,
576 * so go to a char pointer and force an offset of 8 bytes instead.
580 up1 = (FFTWPREFIX(complex) *)pc;
584 up2 = (FFTWPREFIX(complex) *)pc;
587 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_BACKWARD,fftw_flags);
588 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up2,FFTW_FORWARD,fftw_flags);
589 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_BACKWARD,fftw_flags);
590 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,up1,up1,FFTW_FORWARD,fftw_flags);
592 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_BACKWARD,fftw_flags);
593 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p2,FFTW_FORWARD,fftw_flags);
594 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_BACKWARD,fftw_flags);
595 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_3d)(nx,ny,nz,p1,p1,FFTW_FORWARD,fftw_flags);
604 if(fft->plan[i][j][k] == NULL)
606 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
608 gmx_fft_destroy(fft);
610 FFTWPREFIX(free)(p1);
611 FFTWPREFIX(free)(p2);
619 FFTWPREFIX(free)(p1);
620 FFTWPREFIX(free)(p2);
622 fft->real_transform = 0;
633 gmx_fft_init_3d_real(gmx_fft_t * pfft,
640 real *p1,*p2,*up1,*up2;
645 #ifdef GMX_DISABLE_FFTW_MEASURE
646 flags |= GMX_FFT_FLAG_CONSERVATIVE;
649 fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
653 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
659 if( (fft = (gmx_fft_t) FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
665 /* allocate aligned, and extra memory to make it unaligned */
666 p1 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
669 FFTWPREFIX(free)(fft);
674 p2 = (real *) FFTWPREFIX(malloc)(sizeof(real)*( nx*ny*(nz/2+1)*2 + 2) );
677 FFTWPREFIX(free)(p1);
678 FFTWPREFIX(free)(fft);
683 /* make unaligned pointers.
684 * In double precision the actual complex datatype will be 16 bytes,
685 * so go to a void pointer and force an offset of 8 bytes instead.
696 fft->plan[0][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up2,fftw_flags);
697 fft->plan[0][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up2,fftw_flags);
698 fft->plan[0][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)up1,up1,fftw_flags);
699 fft->plan[0][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,up1,(FFTWPREFIX(complex) *)up1,fftw_flags);
701 fft->plan[1][0][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p2,fftw_flags);
702 fft->plan[1][0][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p2,fftw_flags);
703 fft->plan[1][1][0] = FFTWPREFIX(plan_dft_c2r_3d)(nx,ny,nz,(FFTWPREFIX(complex) *)p1,p1,fftw_flags);
704 fft->plan[1][1][1] = FFTWPREFIX(plan_dft_r2c_3d)(nx,ny,nz,p1,(FFTWPREFIX(complex) *)p1,fftw_flags);
713 if(fft->plan[i][j][k] == NULL)
715 gmx_fatal(FARGS,"Error initializing FFTW3 plan.");
717 gmx_fft_destroy(fft);
719 FFTWPREFIX(free)(p1);
720 FFTWPREFIX(free)(p2);
728 FFTWPREFIX(free)(p1);
729 FFTWPREFIX(free)(p2);
731 fft->real_transform = 1;
741 gmx_fft_1d (gmx_fft_t fft,
742 enum gmx_fft_direction dir,
746 int aligned = (((size_t)in_data & (size_t)out_data & 0xf)==0);
747 int inplace = (in_data == out_data);
748 int isforward = (dir == GMX_FFT_FORWARD);
751 if( (fft->real_transform == 1) || (fft->ndim != 1) ||
752 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
754 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
758 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
759 (FFTWPREFIX(complex) *)in_data,
760 (FFTWPREFIX(complex) *)out_data);
766 gmx_fft_many_1d (gmx_fft_t fft,
767 enum gmx_fft_direction dir,
771 return gmx_fft_1d(fft,dir,in_data,out_data);
775 gmx_fft_1d_real (gmx_fft_t fft,
776 enum gmx_fft_direction dir,
780 int aligned = (((size_t)in_data & (size_t)out_data & 0xf)==0);
781 int inplace = (in_data == out_data);
782 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
785 if( (fft->real_transform != 1) || (fft->ndim != 1) ||
786 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
788 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
794 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
795 (real *)in_data,(FFTWPREFIX(complex) *)out_data);
799 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
800 (FFTWPREFIX(complex) *)in_data,(real *)out_data);
807 gmx_fft_many_1d_real (gmx_fft_t fft,
808 enum gmx_fft_direction dir,
812 return gmx_fft_1d_real(fft,dir,in_data,out_data);
816 gmx_fft_2d (gmx_fft_t fft,
817 enum gmx_fft_direction dir,
821 int aligned = (((size_t)in_data & (size_t)out_data & 0xf)==0);
822 int inplace = (in_data == out_data);
823 int isforward = (dir == GMX_FFT_FORWARD);
826 if( (fft->real_transform == 1) || (fft->ndim != 2) ||
827 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
829 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
833 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
834 (FFTWPREFIX(complex) *)in_data,
835 (FFTWPREFIX(complex) *)out_data);
842 gmx_fft_2d_real (gmx_fft_t fft,
843 enum gmx_fft_direction dir,
847 int aligned = (((size_t)in_data & (size_t)out_data & 0xf)==0);
848 int inplace = (in_data == out_data);
849 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
852 if( (fft->real_transform != 1) || (fft->ndim != 2) ||
853 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
855 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
861 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
863 (FFTWPREFIX(complex) *)out_data);
867 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
868 (FFTWPREFIX(complex) *)in_data,
878 gmx_fft_3d (gmx_fft_t fft,
879 enum gmx_fft_direction dir,
883 int aligned = (((size_t)in_data & (size_t)out_data & 0xf)==0);
884 int inplace = (in_data == out_data);
885 int isforward = (dir == GMX_FFT_FORWARD);
888 if( (fft->real_transform == 1) || (fft->ndim != 3) ||
889 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
891 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
895 FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
896 (FFTWPREFIX(complex) *)in_data,
897 (FFTWPREFIX(complex) *)out_data);
904 gmx_fft_3d_real (gmx_fft_t fft,
905 enum gmx_fft_direction dir,
909 int aligned = (((size_t)in_data & (size_t)out_data & 0xf)==0);
910 int inplace = (in_data == out_data);
911 int isforward = (dir == GMX_FFT_REAL_TO_COMPLEX);
914 if( (fft->real_transform != 1) || (fft->ndim != 3) ||
915 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
917 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
923 FFTWPREFIX(execute_dft_r2c)(fft->plan[aligned][inplace][isforward],
925 (FFTWPREFIX(complex) *)out_data);
929 FFTWPREFIX(execute_dft_c2r)(fft->plan[aligned][inplace][isforward],
930 (FFTWPREFIX(complex) *)in_data,
940 gmx_fft_destroy(gmx_fft_t fft)
952 if(fft->plan[i][j][k] != NULL)
955 FFTWPREFIX(destroy_plan)(fft->plan[i][j][k]);
957 fft->plan[i][j][k] = NULL;
963 FFTWPREFIX(free)(fft);
970 gmx_many_fft_destroy(gmx_fft_t fft)
972 gmx_fft_destroy(fft);
978 #endif /* GMX_FFT_FFTW2 */