Commit a2aaef9e authored by Kirill Terekhov's avatar Kirill Terekhov

parallel K-means clustering grid tool

Add algorithm to balance number of points among processors
parent 18dd51e0
......@@ -54,7 +54,8 @@ int main(int argc, char ** argv)
std::cout << "Time to fix normals: " << Timer() - tt << std::endl;
std::cout << "Total face normals fixed: " << fixed << std::endl;
*/
int total_points = 0, total_global_points = 0;
int total_points = 0;
#if defined(USE_OMP)
#pragma omp parallel for reduction(+:total_points)
......@@ -64,30 +65,17 @@ int main(int argc, char ** argv)
Cell n = m.CellByLocalID(q);
if( n->GetStatus() != Element::Ghost ) total_points++;
}
#if defined(USE_MPI)
MPI_Allreduce(&total_global_points,&total,points,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD)
#else
total_global_points = total_points;
#endif
std::vector< double > points_center(total_points*3);
std::vector< int > points_node(total_points);
std::vector< double > points_center(total_points*3);
std::vector< int > points_cluster(total_points,-1);
std::vector< double > cluster_center(K*3);
std::vector< int > cluster_npoints(K);
#if defined(USE_MPI)
std::vector< double > cluster_center_tmp(K*3);
std::vector< int > cluster_npoints_tmp(K);
#endif
int k = 0;
for(int q = 0; q < m.CellLastLocalID(); ++q) if( m.isValidCell(q) )
{
Cell n = m.CellByLocalID(q);
if( n->GetStatus() != Element::Ghost ) points_node[k++] = n->LocalID();
}
#if defined(USE_OMP)
#pragma omp paralell for
#endif
......@@ -101,7 +89,70 @@ int main(int argc, char ** argv)
points_center[q*3+2] = cnt[2];
}
std::cout << m.GetProcessorRank() << " init clusters" << std::endl;
#if defined(USE_MPI)
std::vector<int> displs(m.GetProcessorsNumber());
std::vector<int> counts(m.GetProcessorsNumber());
std::vector<int> npoints(m.GetProcessorsNumber());
int total_global_points = 0;
int total_local_points = 0;
bool balance = false;
MPI_Allgather(&total_points,1,MPI_INT,&npoints[0],1,MPI_INT,MPI_COMM_WORLD);
double imbalance = 1;
for(int k = 0; k < (int) m.GetProcessorsNumber(); ++k)
{
for(int l = k+1; l < (int) m.GetProcessorsNumber(); ++l)
{
imbalance = std::max(imbalance,npoints[k]/(double)npoints[l]);
imbalance = std::max(imbalance,npoints[l]/(double)npoints[k]);
}
}
if( m.GetProcessorRank() == 0 )
std::cout << "Imbalance is " << imbalance << std::endl;
if( imbalance > 1.2 )
balance = true;
if( balance )//redistribute points
{
for(int k = 0; k < (int) m.GetProcessorsNumber(); ++k)
total_global_points += npoints[k];
total_local_points = (int)ceil((double)total_global_points/(double)m.GetProcessorsNumber());
std::vector<double> points_center_global(m.GetProcessorRank() == 0 ? total_global_points*3 : 1);
displs[0] = 0;
counts[0] = npoints[0]*3;
for(int k = 1; k < (int) m.GetProcessorsNumber(); ++k)
{
displs[k] = displs[k-1] + counts[k-1];
counts[k] = npoints[k]*3;
}
MPI_Gatherv(&points_center[0],total_points*3,MPI_DOUBLE,&points_center_global[0],&counts[0],&displs[0],MPI_DOUBLE,0,MPI_COMM_WORLD);
displs[0] = 0;
counts[0] = total_local_points*3;
for(int k = 1; k < (int) m.GetProcessorsNumber(); ++k)
{
displs[k] = displs[k-1] + counts[k-1];
counts[k] = total_local_points*3;
}
counts.back() = total_global_points*3 - displs.back();
total_points = counts[m.GetProcessorRank()]/3;
points_center.resize(total_points*3);
MPI_Scatterv(&points_center_global[0],&counts[0],&displs[0],MPI_DOUBLE,&points_center[0],total_points*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
points_cluster.resize(total_points,-1);
}
#endif
std::vector< double > cluster_center(K*3);
std::vector< int > cluster_npoints(K);
#if defined(USE_MPI)
std::vector< double > cluster_center_tmp(K*3);
std::vector< int > cluster_npoints_tmp(K);
#endif
if( m.GetProcessorRank() == 0 )
std::cout << "Init clusters" << std::endl;
// choose K distinct values for the centers of the clusters
{
......@@ -131,32 +182,28 @@ int main(int argc, char ** argv)
//guter cluster centers over all processes
#if defined(USE_MPI)
{
std::vector<int> recvcounts(m.GetProcessorsNumber()), displs(m.GetProcessorsNumber());
//std::vector<int> counts(m.GetProcessorsNumber()), displs(m.GetProcessorsNumber());
displs[0] = 0;
recvcounts[0] = Kpart*3;
counts[0] = Kpart*3;
for(int i = 1; i < (int)m.GetProcessorsNumber(); ++i)
{
displs[i] = displs[i-1] + recvcounts[i-1];
recvcounts[i] = Kpart*3;
}
recvcounts.back() = K*3 - displs.back();
std::cout << m.GetProcessorRank() << " " << Kstart << " " << Kend << std::endl;
for(int i = 0; i < (int)m.GetProcessorsNumber(); ++i)
{
std::cout << m.GetProcessorRank() << " " << i << " " << displs[i] << " " << recvcounts[i] << std::endl;
displs[i] = displs[i-1] + counts[i-1];
counts[i] = Kpart*3;
}
counts.back() = K*3 - displs.back();
for(int i = Kstart; i < Kend; ++i)
{
cluster_center_tmp[i*3+0] = cluster_center[i*3+0];
cluster_center_tmp[i*3+1] = cluster_center[i*3+1];
cluster_center_tmp[i*3+2] = cluster_center[i*3+2];
}
MPI_Allgatherv(&cluster_center_tmp[Kstart*3],(Kend-Kstart)*3,MPI_DOUBLE,&cluster_center[0],&recvcounts[0],&displs[0],MPI_DOUBLE,MPI_COMM_WORLD);
MPI_Allgatherv(&cluster_center_tmp[Kstart*3],(Kend-Kstart)*3,MPI_DOUBLE,&cluster_center[0],&counts[0],&displs[0],MPI_DOUBLE,MPI_COMM_WORLD);
}
#endif
}
std::cout << m.GetProcessorRank() << " start clustering" << std::endl;
if( m.GetProcessorRank() == 0 )
std::cout << "Start clustering" << std::endl;
int iter = 1;
double t = Timer();
......@@ -203,7 +250,8 @@ int main(int argc, char ** argv)
if(changed == 0 || iter >= max_iterations)
{
std::cout << m.GetProcessorRank() << " break in iteration " << iter << std::endl;
if( m.GetProcessorRank() == 0 )
std::cout << "Break in iteration " << iter << std::endl;
break;
}
......@@ -254,10 +302,36 @@ int main(int argc, char ** argv)
cluster_center[i*3+2] /= (double) cluster_npoints[i];
}
std::cout << "Iteration " << iter << std::endl;
if( m.GetProcessorRank() == 0 )
std::cout << "Iteration " << iter << std::endl;
iter++;
}
std::cout << "Clustering in " << Timer() - t << " secs " << std::endl;
#if defined(USE_MPI)
if( balance )
{
std::vector<int> points_cluster_global(total_global_points);
displs[0] = 0;
counts[0] = total_local_points;
for(int k = 1; k < (int) m.GetProcessorsNumber(); ++k)
{
displs[k] = displs[k-1] + counts[k-1];
counts[k] = total_local_points;
}
counts.back() = total_global_points - displs.back();
MPI_Gatherv(&points_cluster[0],total_points,MPI_INT,&points_cluster_global[0],&counts[0],&displs[0],MPI_INT,0,MPI_COMM_WORLD);
displs[0] = 0;
counts[0] = npoints[0];
for(int k = 1; k < (int) m.GetProcessorsNumber(); ++k)
{
displs[k] = displs[k-1] + counts[k-1];
counts[k] = npoints[k];
}
total_points = counts[m.GetProcessorRank()];
points_cluster.resize(total_points);
MPI_Scatterv(&points_cluster_global[0],&counts[0],&displs[0],MPI_INT,&points_cluster[0],total_points,MPI_INT,0,MPI_COMM_WORLD);
}
#endif
// shows elements of clusters
TagInteger mat = m.CreateTag("MATERIAL",DATA_INTEGER,CELL,NONE,1);
#if defined(USE_OMP)
......
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