Commit e6d44f5b authored by David Rohr's avatar David Rohr
Browse files

Make FFT version work multi-threaded

parent a33835d9
...@@ -214,7 +214,6 @@ int bioem::run() ...@@ -214,7 +214,6 @@ int bioem::run()
/**** If we want to control the number of threads -> omp_set_num_threads(XX); ******/ /**** If we want to control the number of threads -> omp_set_num_threads(XX); ******/
/****************** Declarying class of Probability Pointer *************************/ /****************** Declarying class of Probability Pointer *************************/
pProb = new bioem_Probability[RefMap.ntotRefMap]; pProb = new bioem_Probability[RefMap.ntotRefMap];
crossCor = new bioem_crossCor[RefMap.ntotRefMap];
printf("\tInitializing\n"); printf("\tInitializing\n");
// Inizialzing Probabilites to zero and constant to -Infinity // Inizialzing Probabilites to zero and constant to -Infinity
...@@ -249,7 +248,6 @@ int bioem::run() ...@@ -249,7 +248,6 @@ int bioem::run()
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("\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))); 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)));
//#pragma omp parallel for
for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++) for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
{ {
/***************************************************************************************/ /***************************************************************************************/
...@@ -358,12 +356,6 @@ int bioem::run() ...@@ -358,12 +356,6 @@ int bioem::run()
param.refCTF =NULL; param.refCTF =NULL;
} }
if (crossCor)
{
delete[] crossCor;
crossCor = NULL;
}
if(RefMap.RefMapFFT) if(RefMap.RefMapFFT)
{ {
delete[] RefMap.RefMapFFT; delete[] RefMap.RefMapFFT;
...@@ -384,42 +376,31 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m ...@@ -384,42 +376,31 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myfloat_t sumC,myfloat_t sumsquareC) int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myfloat_t sumC,myfloat_t sumsquareC)
{ {
#pragma omp parallel for #pragma omp parallel
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{ {
mycomplex_t *localCCT, *lCC;
mycomplex_t* localCCT;
localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
mycomplex_t* lCC;
lCC= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); lCC= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
//setting crossCor value to zero for each projection const int num_threads = omp_get_num_threads();
for(int n=0; n < param.param_device.NtotDist ; n++) const int thread_id = omp_get_thread_num();
const int mapsPerThread = (RefMap.ntotRefMap + num_threads - 1) / num_threads;
const int iStart = thread_id * mapsPerThread;
const int iEnd = min(RefMap.ntotRefMap, (thread_id + 1) * mapsPerThread);
for (int iRefMap = iStart; iRefMap < iEnd; iRefMap ++)
{ {
crossCor[iRefMap].value[n]=0.0; calculateCCFFT(iRefMap,iOrient, iConv, sumC,sumsquareC, localConvFFT, localCCT,lCC);
crossCor[iRefMap].disx[n]=-99999;
crossCor[iRefMap].disy[n]=-99999;
} }
calculateCCFFT(iRefMap,iConv, localConvFFT, localCCT,lCC);
myfftw_free(localCCT); myfftw_free(localCCT);
myfftw_free(lCC); myfftw_free(lCC);
} }
//Not in openMP loop SUM OVER PROBABILITIES
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
calProb(iRefMap,iOrient,iConv,sumC,sumsquareC);
}
return(0); return(0);
} }
/////////////NEW ROUTINE //////////////// /////////////NEW ROUTINE ////////////////
int bioem::calculateCCFFT(int iRefMap, int iConv, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC) inline int bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC)
{ {
for(int i=0; i < param.param_device.NumberPixels ; i++ ) for(int i=0; i < param.param_device.NumberPixels ; i++ )
{ {
...@@ -433,53 +414,33 @@ int bioem::calculateCCFFT(int iRefMap, int iConv, mycomplex_t* localConvFFT,myco ...@@ -433,53 +414,33 @@ int bioem::calculateCCFFT(int iRefMap, int iConv, mycomplex_t* localConvFFT,myco
myfftw_execute_dft(param.fft_plan_c2c_backward,localCCT,lCC); myfftw_execute_dft(param.fft_plan_c2c_backward,localCCT,lCC);
// Storing CORRELATIONS FOR CORRESPONDING DISPLACEMENTS & Normalizing after Backward FFT // Storing CORRELATIONS FOR CORRESPONDING DISPLACEMENTS & Normalizing after Backward FFT
int n=0;
for (int cent_x = 0; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter) for (int cent_x = 0; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
{ {
for (int cent_y = 0; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = 0; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels * param.param_device.NumberPixels), cent_x, cent_y);
crossCor[iRefMap].disx[n]=cent_x;
crossCor[iRefMap].disy[n]=cent_y;
n++;
} }
for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y < param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y < param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels), cent_x, param.param_device.NumberPixels-cent_y);
crossCor[iRefMap].disx[n]=cent_x;
crossCor[iRefMap].disy[n]=param.param_device.NumberPixels-cent_y;
n++;
} }
} }
for (int cent_x = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_x < param.param_device.NumberPixels; cent_x=cent_x+param.param_device.GridSpaceCenter) for (int cent_x = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_x < param.param_device.NumberPixels; cent_x=cent_x+param.param_device.GridSpaceCenter)
{ {
for (int cent_y = 0; cent_y < param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = 0; cent_y < param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels), param.param_device.NumberPixels-cent_x, cent_y);
crossCor[iRefMap].disx[n]=param.param_device.NumberPixels-cent_x;
crossCor[iRefMap].disy[n]=cent_y;
n++;
} }
for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y <= param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y <= param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
crossCor[iRefMap].value[n]=float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, float(lCC[cent_x*param.param_device.NumberPixels+cent_y][0])/float(param.param_device.NumberPixels*param.param_device.NumberPixels), param.param_device.NumberPixels-cent_x, param.param_device.NumberPixels-cent_y);
crossCor[iRefMap].disx[n]=param.param_device.NumberPixels-cent_x;
crossCor[iRefMap].disy[n]=param.param_device.NumberPixels-cent_y;
n++;
} }
} }
/* Controll but slows down the parallelisim
if(n> MAX_DISPLACE && n> param.param_device.NtotDist)
{
cout << "Problem with displace center and Max allowed" ;
exit(1);
}*/
//
return (0); return (0);
} }
int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC) inline int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC, int value, int disx, int disy)
{ {
/********************************************************/ /********************************************************/
...@@ -488,12 +449,9 @@ int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t s ...@@ -488,12 +449,9 @@ int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t s
const myfloat_t ForLogProb = (sumsquareC * param.param_device.Ntotpi - sumC * sumC); const myfloat_t ForLogProb = (sumsquareC * param.param_device.Ntotpi - sumC * sumC);
// Loop again over displacements
for(int n=0; n < param.param_device.NtotDist ; n++)
{
// Products of different cross-correlations (first element in formula) // Products of different cross-correlations (first element in formula)
const myfloat_t firstele = param.param_device.Ntotpi * (RefMap.sumsquare_RefMap[iRefMap] * sumsquareC - crossCor[iRefMap].value[n]* crossCor[iRefMap].value[n]) + const myfloat_t firstele = param.param_device.Ntotpi * (RefMap.sumsquare_RefMap[iRefMap] * sumsquareC - value * value) +
2 * RefMap.sum_RefMap[iRefMap] * sumC * crossCor[iRefMap].value[n] - RefMap.sumsquare_RefMap[iRefMap] * sumC * sumC - RefMap.sum_RefMap[iRefMap] * RefMap.sum_RefMap[iRefMap] * sumsquareC; 2 * RefMap.sum_RefMap[iRefMap] * sumC * value - RefMap.sumsquare_RefMap[iRefMap] * sumC * sumC - RefMap.sum_RefMap[iRefMap] * RefMap.sum_RefMap[iRefMap] * sumsquareC;
//******* Calculating log of Prob*********/ //******* Calculating log of Prob*********/
// As in fortran code: logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1); // As in fortran code: logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1);
...@@ -521,13 +479,12 @@ int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t s ...@@ -521,13 +479,12 @@ int bioem::calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t s
if(pProb[iRefMap].max_prob < logpro) if(pProb[iRefMap].max_prob < logpro)
{ {
pProb[iRefMap].max_prob = logpro; pProb[iRefMap].max_prob = logpro;
pProb[iRefMap].max_prob_cent_x = crossCor[iRefMap].disx[n]; pProb[iRefMap].max_prob_cent_x = disx;
pProb[iRefMap].max_prob_cent_y = crossCor[iRefMap].disy[n]; pProb[iRefMap].max_prob_cent_y = disy;
pProb[iRefMap].max_prob_orient = iOrient; pProb[iRefMap].max_prob_orient = iOrient;
pProb[iRefMap].max_prob_conv = iConv; pProb[iRefMap].max_prob_conv = iConv;
} }
} }
}
return (0); return (0);
} }
...@@ -594,8 +551,6 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT) ...@@ -594,8 +551,6 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+float(param.param_device.NumberPixels)/2.0+0.5); j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+float(param.param_device.NumberPixels)/2.0+0.5);
localproj[i*param.param_device.NumberPixels+j][0]+=Model.densityPointsModel[n]/Model.NormDen; localproj[i*param.param_device.NumberPixels+j][0]+=Model.densityPointsModel[n]/Model.NormDen;
} }
/**** Output Just to check****/ /**** Output Just to check****/
...@@ -648,11 +603,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b ...@@ -648,11 +603,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
} }
/**** Bringing convoluted Map to real Space ****/ /**** Bringing convoluted Map to real Space ****/
//#pragma omp critical myfftw_execute_dft(param.fft_plan_c2c_backward,localmultFFT,localconvFFT);
{
myfftw_execute_dft(param.fft_plan_c2c_backward,localmultFFT,localconvFFT);
}
/****Asigning convolution fftw_complex to bioem_map ****/ /****Asigning convolution fftw_complex to bioem_map ****/
for(int i=0; i < param.param_device.NumberPixels ; i++ ) for(int i=0; i < param.param_device.NumberPixels ; i++ )
......
...@@ -29,11 +29,10 @@ public: ...@@ -29,11 +29,10 @@ public:
int createProjection(int iMap, mycomplex_t* map); int createProjection(int iMap, mycomplex_t* map);
int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare); int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare);
int calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC); int calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC, int value, int disx, int disy);
int calculateCCFFT(int iMap, int iConv, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC); int calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC);
bioem_Probability* pProb; bioem_Probability* pProb;
bioem_crossCor* crossCor;
protected: protected:
virtual int deviceInit(); virtual int deviceInit();
......
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