fix double precision mode

bioem_algorithm.h

@@ -70,7 +70,7 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
 	return(logpro);
 }
 
-__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
+__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, myfloat_t value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
 {
 	// ********************************************************
 	// *********** Calculates the BioEM probability ***********
@@ -379,8 +379,6 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
 	}
 	else
 #endif
-
-		// Summing & Storing total/Orientation Probabilites for each map
 	{
 		update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb);
 	}

bioem_cuda.cu

@@ -19,6 +19,43 @@ using namespace std;
 	} \
 }
 
+static const char *cufftGetErrorStrung(cufftResult error)
+{
+    switch (error)
+    {
+        case CUFFT_SUCCESS:
+            return "CUFFT_SUCCESS";
+
+        case CUFFT_INVALID_PLAN:
+            return "CUFFT_INVALID_PLAN";
+
+        case CUFFT_ALLOC_FAILED:
+            return "CUFFT_ALLOC_FAILED";
+
+        case CUFFT_INVALID_TYPE:
+            return "CUFFT_INVALID_TYPE";
+
+        case CUFFT_INVALID_VALUE:
+            return "CUFFT_INVALID_VALUE";
+
+        case CUFFT_INTERNAL_ERROR:
+            return "CUFFT_INTERNAL_ERROR";
+
+        case CUFFT_EXEC_FAILED:
+            return "CUFFT_EXEC_FAILED";
+
+        case CUFFT_SETUP_FAILED:
+            return "CUFFT_SETUP_FAILED";
+
+        case CUFFT_INVALID_SIZE:
+            return "CUFFT_INVALID_SIZE";
+
+        case CUFFT_UNALIGNED_DATA:
+            return "CUFFT_UNALIGNED_DATA";
+    }
+    return "UNKNOWN";
+}
+
 bioem_cuda::bioem_cuda()
 {
 	deviceInitialized = 0;
@@ -127,9 +164,10 @@ int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const myfloat_t* c
 		{
 			const int num = min(CUDA_FFTS_AT_ONCE, maxRef - i);
 			multComplexMap<<<num, CUDA_THREAD_COUNT, 0, cudaStream>>>(&pConvMapFFT[(iConv & 1) * param.FFTMapSize], pRefMapsFFT, pFFTtmp2, param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D, param.FFTMapSize, num, i);
-			if (mycufftExecC2R(i + CUDA_FFTS_AT_ONCE > maxRef ? plan[1] : plan[0], pFFTtmp2, pFFTtmp) != CUFFT_SUCCESS)
+			cufftResult err = mycufftExecC2R(i + CUDA_FFTS_AT_ONCE > maxRef ? plan[1] : plan[0], pFFTtmp2, pFFTtmp);
+			if (err != CUFFT_SUCCESS)
 			{
-				cout << "Error running CUFFT\n";
+				cout << "Error running CUFFT " << cufftGetErrorStrung(err) << "\n";
 				exit(1);
 			}
 			cuDoRefMapsFFT<<<divup(num, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>>(iProjectionOut, iConv, pFFTtmp, sumC, sumsquareC, pProb_device, param.param_device, *gpumap, num, i);
@@ -301,7 +339,7 @@ int bioem_cuda::deviceStartRun()
 		for (int i = 0; i < 2; i++)
 		{
 			int n[2] = {param.param_device.NumberPixels, param.param_device.NumberPixels};
-			if (cufftPlanMany(&plan[i], 2, n, NULL, 1, 0, NULL, 1, 0, CUFFT_C2R, i ? (maxRef % CUDA_FFTS_AT_ONCE) : CUDA_FFTS_AT_ONCE) != CUFFT_SUCCESS)
+			if (cufftPlanMany(&plan[i], 2, n, NULL, 1, 0, NULL, 1, 0, MY_CUFFT_C2R, i ? (maxRef % CUDA_FFTS_AT_ONCE) : CUDA_FFTS_AT_ONCE) != CUFFT_SUCCESS)
 			{
 				cout << "Error planning CUFFT\n";
 				exit(1);

include/defs.h

 #ifndef BIOEM_DEFS_H
 #define BIOEM_DEFS_H
 
-//#define BIOEM_USE_DOUBLE
+#define BIOEM_USE_DOUBLE
 
 #ifndef BIOEM_USE_DOUBLE
 typedef float myfloat_t;
@@ -16,6 +16,7 @@ typedef float myfloat_t;
 #define myfftw_plan_dft_r2c_2d fftwf_plan_dft_r2c_2d
 #define myfftw_plan_dft_c2r_2d fftwf_plan_dft_c2r_2d
 #define myfftw_plan fftwf_plan
+#define MY_CUFFT_C2R CUFFT_C2R
 #define mycufftExecC2R cufftExecC2R
 #define mycuComplex_t cuComplex
 #else
@@ -33,6 +34,7 @@ typedef double myfloat_t;
 #define myfftw_plan fftw_plan
 #define mycufftExecC2R cufftExecZ2D
 #define mycuComplex_t cuDoubleComplex
+#define MY_CUFFT_C2R CUFFT_Z2D
 #endif
 typedef myfloat_t mycomplex_t[2];
 

include/param.h

@@ -64,7 +64,7 @@ public:
 	myfloat_t gridCTF_amp;
 	// Others
 	//myfloat_t volu;//in device class
-	myfloat3_t angles[MAX_ORIENT];
+	myfloat3_t* angles;
 	int nTotGridAngles;
 	int nTotCTFs;
 	//myfloat_t Ntotpi;//in device class

param.cpp

@@ -33,6 +33,7 @@ bioem_param::bioem_param()
 
 	refCTF = NULL;
 	CtfParam = NULL;
+	angles = NULL;
 }
 
 int bioem_param::readParameters()
@@ -235,6 +236,7 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
 	cos_grid_beta = 2.f / (myfloat_t) angleGridPointsBeta;
 
 	// Euler Angle Array
+	angles = new myfloat3_t[angleGridPointsAlpha * angleGridPointsBeta * angleGridPointsAlpha];
 	for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
 	{
 		for (int ibeta = 0; ibeta < angleGridPointsBeta; ibeta ++)
@@ -344,7 +346,6 @@ int bioem_param::CalculateRefCTF()
 	return(0);
 }
 
-
 bioem_param::~bioem_param()
 {
 	releaseFFTPlans();
@@ -356,6 +357,7 @@ bioem_param::~bioem_param()
 	numberGridPointsCTF_phase = 0;
 	param_device.maxDisplaceCenter = 0;
 	numberGridPointsDisplaceCenter = 0;
-	delete[] refCTF;
-	delete[] CtfParam;
+	if (refCTF) delete[] refCTF;
+	if (CtfParam) delete[] CtfParam;
+	if (angles) delete[] angles;
 }