#include <config.h>
#endif
+#include <algorithm>
#include <stdio.h>
#include <stdlib.h>
#ifdef GMX_FFT_FFTW3
-#ifdef GMX_THREAD_MPI
+#include "thread_mpi/mutex.h"
+#include "gromacs/utility/exceptions.h"
/* none of the fftw3 calls, except execute(), are thread-safe, so
we need to serialize them with this mutex. */
-static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
-
-#define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex)
-#define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex)
-#else /* GMX_THREAD_MPI */
-#define FFTW_LOCK
-#define FFTW_UNLOCK
-#endif /* GMX_THREAD_MPI */
+static tMPI::mutex big_fftw_mutex;
+#define FFTW_LOCK try { big_fftw_mutex.lock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
+#define FFTW_UNLOCK try { big_fftw_mutex.unlock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
#endif /* GMX_FFT_FFTW3 */
-static double fft5d_fmax(double a, double b){
- return (a>b)?a:b;
-}
-
/* largest factor smaller than sqrt */
static int lfactor(int z) {
int i;
- for (i=sqrt(z);;i--)
+ for (i=static_cast<int>(sqrt(static_cast<double>(z)));;i--)
if (z%i==0) return i;
return 1;
}
static int lpfactor(int z) {
int f = l2factor(z);
if (f==1) return z;
- return fft5d_fmax(lpfactor(f),lpfactor(z/f));
+ return std::max(lpfactor(f),lpfactor(z/f));
}
#ifndef GMX_MPI
*/
/* int lsize = fmax(N[0]*M[0]*K[0]*nP[0],N[1]*M[1]*K[1]*nP[1]); */
- lsize = fft5d_fmax(N[0]*M[0]*K[0]*nP[0],fft5d_fmax(N[1]*M[1]*K[1]*nP[1],C[2]*M[2]*K[2]));
+ lsize = std::max(N[0]*M[0]*K[0]*nP[0],std::max(N[1]*M[1]*K[1]*nP[1],C[2]*M[2]*K[2]));
/* int lsize = fmax(C[0]*M[0]*K[0],fmax(C[1]*M[1]*K[1],C[2]*M[2]*K[2])); */
if (!(flags&FFT5D_NOMALLOC)) {
snew_aligned(lin, lsize, 32);
}
static void print_localdata(const t_complex* lin, const char* txt, int s, fft5d_plan plan) {
- int x,y,z,l,t;
+ int x,y,z,l;
int *coor = plan->coor;
int xs[3],xl[3],xc[3],NG[3];
int ll=(plan->flags&FFT5D_REALCOMPLEX)?1:2;
#ifdef GMX_MPI
MPI_Comm *cart=plan->cart;
#endif
-
+#ifdef NOGMX
double time_fft=0,time_local=0,time_mpi[2]={0},time=0;
+#endif
int *N=plan->N,*M=plan->M,*K=plan->K,*pN=plan->pN,*pM=plan->pM,*pK=plan->pK,
*C=plan->C,*P=plan->P,**iNin=plan->iNin,**oNin=plan->oNin,**iNout=plan->iNout,**oNout=plan->oNout;
int s=0,tstart,tend,bParallelDim;
{
FFTW(destroy_plan)(plan->mpip[s]);
}
+#endif /* FFT5D_MPI_TRANSPOS */
if (plan->p3d)
{
FFTW(destroy_plan)(plan->p3d);
}
-#endif /* FFT5D_MPI_TRANSPOS */
+ FFTW_UNLOCK;
#endif /* GMX_FFT_FFTW3 */
if (!(plan->flags&FFT5D_NOMALLOC))
#include "gmx_fft.h"
#include "gmx_fatal.h"
-
#ifdef GMX_DOUBLE
#define FFTWPREFIX(name) fftw_ ## name
#else
#define FFTWPREFIX(name) fftwf_ ## name
#endif
-#ifdef GMX_THREAD_MPI
-#include "thread_mpi/threads.h"
-#endif
-
+#include "thread_mpi/mutex.h"
+#include "gromacs/utility/exceptions.h"
-#ifdef GMX_THREAD_MPI
/* none of the fftw3 calls, except execute(), are thread-safe, so
we need to serialize them with this mutex. */
-static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
-static gmx_bool gmx_fft_threads_initialized=FALSE;
-#define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex)
-#define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex)
-#else /* GMX_THREAD_MPI */
-#define FFTW_LOCK
-#define FFTW_UNLOCK
-#endif /* GMX_THREAD_MPI */
+static tMPI::mutex big_fftw_mutex;
+#define FFTW_LOCK try { big_fftw_mutex.lock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
+#define FFTW_UNLOCK try { big_fftw_mutex.unlock(); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
/* 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.
Consequesences:
int
gmx_fft_init_1d(gmx_fft_t * pfft,
int nx,
- gmx_fft_flag flags)
+ gmx_fft_flag flags)
{
return gmx_fft_init_many_1d(pfft,nx,1,flags);
}
int
gmx_fft_init_many_1d(gmx_fft_t * pfft,
- int nx,
- int howmany,
- gmx_fft_flag flags)
+ int nx,
+ int howmany,
+ gmx_fft_flag flags)
{
gmx_fft_t fft;
FFTWPREFIX(complex) *p1,*p2,*up1,*up2;
size_t pc;
int i,j,k;
int fftw_flags;
-
+
#ifdef GMX_DISABLE_FFTW_MEASURE
flags |= GMX_FFT_FLAG_CONSERVATIVE;
#endif
-
- fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
+
+ fftw_flags = (flags & GMX_FFT_FLAG_CONSERVATIVE) ? FFTW_ESTIMATE : FFTW_MEASURE;
if(pfft==NULL)
{
return EINVAL;
}
*pfft = NULL;
-
+
FFTW_LOCK;
if( (fft = (gmx_fft_t)FFTWPREFIX(malloc)(sizeof(struct gmx_fft))) == NULL)
{
FFTW_UNLOCK;
return ENOMEM;
- }
-
+ }
+
/* allocate aligned, and extra memory to make it unaligned */
p1 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
if(p1==NULL)
FFTW_UNLOCK;
return ENOMEM;
}
-
+
p2 = (FFTWPREFIX(complex) *) FFTWPREFIX(malloc)(sizeof(FFTWPREFIX(complex))*(nx+2)*howmany);
if(p2==NULL)
{
FFTW_UNLOCK;
return ENOMEM;
}
-
- /* make unaligned pointers.
+
+ /* make unaligned pointers.
* In double precision the actual complex datatype will be 16 bytes,
* so go to a char pointer and force an offset of 8 bytes instead.
*/
pc = (size_t)p1;
- pc += 8;
+ pc += 8;
up1 = (FFTWPREFIX(complex) *)pc;
-
+
pc = (size_t)p2;
- pc += 8;
+ pc += 8;
up2 = (FFTWPREFIX(complex) *)pc;
-
+
/* int rank, const int *n, int howmany,
fftw_complex *in, const int *inembed,
int istride, int idist,
fftw_complex *out, const int *onembed,
int ostride, int odist,
int sign, unsigned flags */
- fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
- fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
- fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
- fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
- fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
- fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
- fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
- fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
+ fft->plan[0][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
+ fft->plan[0][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
+ fft->plan[0][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
+ fft->plan[0][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,up1,&nx,1,nx,up1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
+ fft->plan[1][0][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
+ fft->plan[1][0][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p2,&nx,1,nx,FFTW_FORWARD,fftw_flags);
+ fft->plan[1][1][0] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_BACKWARD,fftw_flags);
+ fft->plan[1][1][1] = FFTWPREFIX(plan_many_dft)(1,&nx,howmany,p1,&nx,1,nx,p1,&nx,1,nx,FFTW_FORWARD,fftw_flags);
for(i=0;i<2;i++)
{
}
}
}
- }
-
+ }
+
FFTWPREFIX(free)(p1);
FFTWPREFIX(free)(p2);
-
+
fft->real_transform = 0;
fft->ndim = 1;
-
+
*pfft = fft;
FFTW_UNLOCK;
return 0;
}
-
int
gmx_fft_init_1d_real(gmx_fft_t * pfft,
int nx,
int aligned = ((((size_t)in_data | (size_t)out_data) & 0xf)==0);
int inplace = (in_data == out_data);
int isforward = (dir == GMX_FFT_FORWARD);
-
+
/* Some checks */
if( (fft->real_transform == 1) || (fft->ndim != 1) ||
((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
{
gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
return EINVAL;
- }
+ }
FFTWPREFIX(execute_dft)(fft->plan[aligned][inplace][isforward],
(FFTWPREFIX(complex) *)in_data,
(FFTWPREFIX(complex) *)out_data);
-
+
return 0;
}
int
gmx_fft_many_1d (gmx_fft_t fft,
- enum gmx_fft_direction dir,
- void * in_data,
- void * out_data)
+ enum gmx_fft_direction dir,
+ void * in_data,
+ void * out_data)
{
return gmx_fft_1d(fft,dir,in_data,out_data);
}
-int
+int
gmx_fft_1d_real (gmx_fft_t fft,
enum gmx_fft_direction dir,
void * in_data,
gmx_fft_destroy(fft);
}
-#else
-int
-gmx_fft_fftw3_empty;
+void gmx_fft_cleanup()
+{
+ FFTWPREFIX(cleanup)();
+}
#endif /* GMX_FFT_FFTW3 */
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff