Skip to content

GitLab

B
BioEM
Commit 25d2e7dc authored by David Rohr's avatar David Rohr
Browse files
Options

finish first mpi implementation

parent 5fd06c57
Hide whitespace changes
Inline Side-by-side
Showing with 193 additions and 77 deletions
bioem.cpp
View file @ 25d2e7dc
......@@ -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);
}
......
bioem_algorithm.h
View file @ 25d2e7dc
......@@ -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();
......
bioem_cuda.cu
View file @ 25d2e7dc
......@@ -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);
}
......
include/defs.h
View file @ 25d2e7dc
......@@ -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];
......
include/map.h
View file @ 25d2e7dc
......@@ -5,6 +5,7 @@
#include "param.h"
#include <complex>
#include <math.h>
#include <boost/concept_check.hpp>
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
......
main.cpp
View file @ 25d2e7dc
......@@ -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;
......
param.cpp
View file @ 25d2e7dc
......@@ -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);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment