diff --git a/bioem.cpp b/bioem.cpp index e48f1ec60bf34a2a8a81cbe15714933d6221ac70..3e62798d8b8e2a5b001c0e4d0b105509228f69aa 100644 --- a/bioem.cpp +++ b/bioem.cpp @@ -227,12 +227,13 @@ int bioem::configure(int ac, char* av[]) //refCtf, CtfParam, angles automatically filled by precalculare function below MPI_Bcast(&Model, sizeof(Model), MPI_BYTE, 0, MPI_COMM_WORLD); - Model.points = (bioem_model::bioem_model_point*) mallocchk(sizeof(bioem_model::bioem_model_point) * Model.nPointsModel); + if (mpi_rank != 0) Model.points = (bioem_model::bioem_model_point*) mallocchk(sizeof(bioem_model::bioem_model_point) * Model.nPointsModel); MPI_Bcast(Model.points, sizeof(bioem_model::bioem_model_point) * Model.nPointsModel, MPI_BYTE, 0, MPI_COMM_WORLD); - + MPI_Bcast(&RefMap, sizeof(RefMap), MPI_BYTE, 0, MPI_COMM_WORLD); + if (mpi_rank != 0) RefMap.maps = (myfloat_t*) mallocchk(RefMap.refMapSize * sizeof(myfloat_t) * RefMap.ntotRefMap); + MPI_Bcast(RefMap.maps, RefMap.refMapSize * sizeof(myfloat_t) * RefMap.ntotRefMap, MPI_BYTE, 0, MPI_COMM_WORLD); #endif - if (mpi_rank != 0) return(1); // ****************** Precalculating Necessary Stuff ********************* param.PrepareFFTs(); @@ -286,7 +287,7 @@ int bioem::run() // **** If we want to control the number of threads -> omp_set_num_threads(XX); ****** // ****************** Declarying class of Probability Pointer ************************* - printf("\tInitializing Probabilities\n"); + if (mpi_rank == 0) printf("\tInitializing Probabilities\n"); // Inizialzing Probabilites to zero and constant to -Infinity for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++) { @@ -294,7 +295,6 @@ int bioem::run() pProbMap.Total = 0.0; pProbMap.Constoadd = -9999999; - pProbMap.max_prob = -9999999; if (param.param_device.writeAngles) { for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++) @@ -333,9 +333,13 @@ int bioem::run() HighResTimer timer; - if (DebugOutput >= 1) printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %d²\n", param.nTotGridAngles, param.nTotCTFs, RefMap.ntotRefMap, 2 * param.param_device.maxDisplaceCenter + param.param_device.GridSpaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels); -// printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1))); - for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++) + if (DebugOutput >= 1 && mpi_rank == 0) printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %d²\n", param.nTotGridAngles, param.nTotCTFs, RefMap.ntotRefMap, 2 * param.param_device.maxDisplaceCenter + param.param_device.GridSpaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels); + + const int iOrientStart = (int) ((long long int) mpi_rank * param.nTotGridAngles / mpi_size); + int iOrientEnd = (int) ((long long int) (mpi_rank + 1) * param.nTotGridAngles / mpi_size); + if (iOrientEnd > param.nTotGridAngles) iOrientEnd = param.nTotGridAngles; + + for (int iOrient = iOrientStart; iOrient < iOrientEnd; iOrient++) { // *************************************************************************************** // ***** Creating Projection for given orientation and transforming to Fourier space ***** @@ -347,7 +351,6 @@ int bioem::run() // ***** **** Internal Loop over convolutions **** ***** for (int iConv = 0; iConv < param.nTotCTFs; iConv++) { - if (DebugOutput >= 2) printf("\t\tConvolution %d %d\n", iOrient, iConv); // *** Calculating convolutions of projection map and crosscorrelations *** if (DebugOutput >= 2) timer.ResetStart(); @@ -393,54 +396,163 @@ int bioem::run() // ************* Writing Out Probabilities *************** // *** Angular Probability *** + +#ifdef WITH_MPI + if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart(); + //Reduce Constant and summarize probabilities + { + myfloat_t* tmp1 = new myfloat_t[RefMap.ntotRefMap]; + myfloat_t* tmp2 = new myfloat_t[RefMap.ntotRefMap]; + myfloat_t* tmp3 = new myfloat_t[RefMap.ntotRefMap]; + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + tmp1[i] = pProb.getProbMap(i).Constoadd; + } + MPI_Allreduce(tmp1, tmp2, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD); + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + bioem_Probability_map& pProbMap = pProb.getProbMap(i); + tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]); + } + MPI_Reduce(tmp1, tmp3, RefMap.ntotRefMap, MY_MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); - ofstream angProbfile; - if(param.param_device.writeAngles) + //Find MaxProb + { + int* tmpi1 = new int[RefMap.ntotRefMap]; + int* tmpi2 = new int[RefMap.ntotRefMap]; + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + bioem_Probability_map& pProbMap = pProb.getProbMap(i); + tmpi1[i] = tmp2[i] <= pProbMap.Constoadd ? mpi_rank : -1; + } + MPI_Allreduce(tmpi1, tmpi2, RefMap.ntotRefMap, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + if (tmpi2[i] == -1) + { + if (mpi_rank == 0) printf("Error: Could not find highest probability\n"); + } + else if (tmpi2[i] != 0) //Skip if rank 0 already has highest probability + { + if (mpi_rank == 0) + { + MPI_Recv(&pProb.getProbMap(i).max, sizeof(pProb.getProbMap(i).max), MPI_BYTE, tmpi2[i], i, MPI_COMM_WORLD, NULL); + } + else if (mpi_rank == tmpi2[i]) + { + MPI_Send(&pProb.getProbMap(i).max, sizeof(pProb.getProbMap(i).max), MPI_BYTE, 0, i, MPI_COMM_WORLD); + } + } + } + delete[] tmpi1; + delete[] tmpi2; + } + + if (mpi_rank == 0) + { + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + bioem_Probability_map& pProbMap = pProb.getProbMap(i); + pProbMap.Total = tmp3[i]; + pProbMap.Constoadd = tmp2[i]; + } + } + + delete[] tmp1; + delete[] tmp2; + delete[] tmp3; + if (DebugOutput >= 2 && mpi_rank == 0) printf("Time MPI Reduction: %f\n", timer.GetCurrentElapsedTime()); + } + + //Angle Reduction and Probability summation for individual angles + if (param.param_device.writeAngles) { - angProbfile.open ("ANG_PROB"); + const int count = RefMap.ntotRefMap * param.nTotGridAngles; + myfloat_t* tmp1 = new myfloat_t[count]; + myfloat_t* tmp2 = new myfloat_t[count]; + myfloat_t* tmp3 = new myfloat_t[count]; + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + tmp1[i] = pProb.getProbMap(i).Constoadd; + } + MPI_Allreduce(tmp1, tmp2, count, MY_MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD); + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + for (int j = 0;j < param.nTotGridAngles;j++) + { + bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j); + tmp1[i * param.nTotGridAngles + j] = pProbAngle.forAngles * exp(pProbAngle.ConstAngle - tmp2[i * param.nTotGridAngles + j]); + } + } + MPI_Reduce(tmp1, tmp3, count, MY_MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD); + if (mpi_rank == 0) + { + for (int i = 0;i < RefMap.ntotRefMap;i++) + { + for (int j = 0;j < param.nTotGridAngles;j++) + { + bioem_Probability_angle& pProbAngle = pProb.getProbAngle(i, j); + pProbAngle.forAngles = tmp3[i * param.nTotGridAngles + j]; + pProbAngle.ConstAngle = tmp2[i * param.nTotGridAngles + j]; + } + } + } + delete[] tmp1; + delete[] tmp2; + delete[] tmp3; } +#endif - ofstream outputProbFile; - outputProbFile.open ("Output_Probabilities"); - for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++) + if (mpi_rank == 0) { - // **** Total Probability *** - bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap); + ofstream angProbfile; + if(param.param_device.writeAngles) + { + angProbfile.open ("ANG_PROB"); + } + + ofstream outputProbFile; + outputProbFile.open ("Output_Probabilities"); + for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++) + { + // **** Total Probability *** + bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap); - outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProbMap.Constoadd << "\n"; + outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProbMap.Constoadd << "\n"; - outputProbFile << "RefMap " << iRefMap << " Maximizing Param: "; + outputProbFile << "RefMap " << iRefMap << " Maximizing Param: "; - // *** Param that maximize probability**** - outputProbFile << (pProbMap.max_prob + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " "; - outputProbFile << param.angles[pProbMap.max_prob_orient].pos[0] << " "; - outputProbFile << param.angles[pProbMap.max_prob_orient].pos[1] << " "; - outputProbFile << param.angles[pProbMap.max_prob_orient].pos[2] << " "; - outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[0] << " "; - outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[1] << " "; - outputProbFile << param.CtfParam[pProbMap.max_prob_conv].pos[2] << " "; - outputProbFile << pProbMap.max_prob_cent_x << " "; - outputProbFile << pProbMap.max_prob_cent_y; - outputProbFile << "\n"; + // *** Param that maximize probability**** + outputProbFile << (pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " "; + outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[0] << " "; + outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[1] << " "; + outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[2] << " "; + outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " "; + outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " "; + outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " "; + outputProbFile << pProbMap.max.max_prob_cent_x << " "; + outputProbFile << pProbMap.max.max_prob_cent_y; + outputProbFile << "\n"; - // *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap); + // *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap); - if(param.param_device.writeAngles) - { - for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++) + if(param.param_device.writeAngles) { - bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); + for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++) + { + bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); - angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n"; + angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n"; + } } } - } - if(param.param_device.writeAngles) - { - angProbfile.close(); + if(param.param_device.writeAngles) + { + angProbfile.close(); + } + outputProbFile.close(); } - outputProbFile.close(); return(0); } diff --git a/bioem_algorithm.h b/bioem_algorithm.h index d8ab2a48c588e1cd6b8ad647c0d93c1798085933..5b8d357cfe29c8970331779e9801f42063d585c4 100644 --- a/bioem_algorithm.h +++ b/bioem_algorithm.h @@ -30,6 +30,20 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef { pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd); pProbMap.Constoadd = logpro; + + // ********** Getting parameters that maximize the probability *********** + if (GPUAlgo == 2) + { + bufint[0] = 1; + buf3[1] = logpro; + } + else + { + pProbMap.max.max_prob_cent_x = cent_x; + pProbMap.max.max_prob_cent_y = cent_y; + } + pProbMap.max.max_prob_orient = iOrient; + pProbMap.max.max_prob_conv = iConv; } if (GPUAlgo != 2) pProbMap.Total += exp(logpro - pProbMap.Constoadd); @@ -46,25 +60,6 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef if (GPUAlgo != 2) pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); } - - // ********** Getting parameters that maximize the probability *********** - if(pProbMap.max_prob < logpro) - { - pProbMap.max_prob = logpro; - - if (GPUAlgo == 2) - { - bufint[0] = 1; - buf3[1] = logpro; - } - else - { - pProbMap.max_prob_cent_x = cent_x; - pProbMap.max_prob_cent_y = cent_y; - } - pProbMap.max_prob_orient = iOrient; - pProbMap.max_prob_conv = iConv; - } } __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param, const myfloat_t sum, const myfloat_t sumsquare, const myfloat_t crossproMapConv, const myfloat_t sumref, const myfloat_t sumsquareref) @@ -97,6 +92,12 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo { pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd); pProbMap.Constoadd = logpro; + + // ********** Getting parameters that maximize the probability *********** + pProbMap.max.max_prob_cent_x = disx; + pProbMap.max.max_prob_cent_y = disy; + pProbMap.max.max_prob_orient = iOrient; + pProbMap.max.max_prob_conv = iConv; } pProbMap.Total += exp(logpro - pProbMap.Constoadd); @@ -111,15 +112,6 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo } pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); } - - if(pProbMap.max_prob < logpro) - { - pProbMap.max_prob = logpro; - pProbMap.max_prob_cent_x = disx; - pProbMap.max_prob_cent_y = disy; - pProbMap.max_prob_orient = iOrient; - pProbMap.max_prob_conv = iConv; - } } __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, const int iConv, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability& pProb, const bioem_param_device& param, const bioem_RefMap& RefMap) @@ -305,8 +297,8 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient if (bufint[0] == 1 && buf3[1] == logpro && iRefMap < RefMap.ntotRefMap && atomicAdd(&bufint[0], 1) == 1) { - pProbMap.max_prob_cent_x = cent_x; - pProbMap.max_prob_cent_y = cent_y; + pProbMap.max.max_prob_cent_x = cent_x; + pProbMap.max.max_prob_cent_y = cent_y; } __syncthreads(); diff --git a/bioem_cuda.cu b/bioem_cuda.cu index 09e1db5504d7a0926b88d4ab24a5ecd405cdcf79..e5ccf87bc6ec27f12d3c5a86cde00e84f59f3084 100644 --- a/bioem_cuda.cu +++ b/bioem_cuda.cu @@ -336,6 +336,11 @@ int bioem_cuda::selectCudaDevice() printf("textureAlignment = %lld", (unsigned long long int) deviceProp.textureAlignment); } + if (DebugOutput >= 1) + { + printf("BioEM for CUDA initialized, %d GPUs found, using GPU %d\n", count, bestDevice); + } + return(0); } diff --git a/include/defs.h b/include/defs.h index 0f31d6e863efa7196ee0df4bb0a19b80a84caace..82d00f70a568dad7043334b38955a6a4764a6794 100644 --- a/include/defs.h +++ b/include/defs.h @@ -20,6 +20,7 @@ typedef float myfloat_t; #define MY_CUFFT_C2R CUFFT_C2R #define mycufftExecC2R cufftExecC2R #define mycuComplex_t cuComplex +#define MY_MPI_FLOAT MPI_FLOAT #else typedef double myfloat_t; #define myfftw_malloc fftw_malloc @@ -37,6 +38,7 @@ typedef double myfloat_t; #define mycufftExecC2R cufftExecZ2D #define mycuComplex_t cuDoubleComplex #define MY_CUFFT_C2R CUFFT_Z2D +#define MY_MPI_FLOAT MPI_DOUBLE #endif typedef myfloat_t mycomplex_t[2]; diff --git a/include/map.h b/include/map.h index 1495ec4fde77c6aa99edbe9d0409feeec971b5bd..907a833d595b41e62948040e865467ddca0cfff3 100644 --- a/include/map.h +++ b/include/map.h @@ -5,6 +5,7 @@ #include "param.h" #include #include +#include class bioem_param; class bioem; @@ -85,8 +86,12 @@ class bioem_Probability_map public: myfloat_t Total; myfloat_t Constoadd; - myfloat_t max_prob; - int max_prob_cent_x, max_prob_cent_y, max_prob_orient, max_prob_conv; + + class bioem_Probability_map_max + { + public: + int max_prob_cent_x, max_prob_cent_y, max_prob_orient, max_prob_conv; + } max; }; class bioem_Probability_angle diff --git a/main.cpp b/main.cpp index 11a39322802a8af73d01d1c88123bd4ebf162883..bfd5f7b10a218cf6052a6f3ec248b33fcba4d68d 100644 --- a/main.cpp +++ b/main.cpp @@ -70,18 +70,18 @@ int main(int argc, char* argv[]) } // ************ Configuration and Pre-calculating necessary objects ***************** - printf("Configuring\n"); + if (mpi_rank == 0) printf("Configuring\n"); if (bio->configure(argc, argv) == 0) { // ******************************* Run BioEM routine ****************************** - printf("Running\n"); + if (mpi_rank == 0) printf("Running\n"); timer.Start(); bio->run(); timer.Stop(); // ************************************ End ********************************** - printf ("The code ran for %f seconds.\n", timer.GetElapsedTime()); + printf ("The code ran for %f seconds (rank %d).\n", timer.GetElapsedTime(), mpi_rank); bio->cleanup(); } delete bio; diff --git a/param.cpp b/param.cpp index e1853c101ddf6b674d218779015a2cf58030f6c4..db86fd03fb8e157298891e61c70f9f158c786a19 100644 --- a/param.cpp +++ b/param.cpp @@ -273,7 +273,7 @@ int bioem_param::readParameters(const char* fileinput) void bioem_param::PrepareFFTs() { - cout << "Preparing FFTs\n"; + if (mpi_rank) cout << "Preparing FFTs\n"; releaseFFTPlans(); mycomplex_t *tmp_map, *tmp_map2; tmp_map = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels);