From 659793bcdd387139b284d511b811529727ae829c Mon Sep 17 00:00:00 2001 From: thrud primrose Date: Sat, 21 May 2022 01:51:24 +0200 Subject: [PATCH 1/4] implement routines to generate and get the dual graph representation of the mesh, allow partition() to take the vector fo edge weights as an additional argument --- PartitionMetis.h | 213 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 148 insertions(+), 65 deletions(-) diff --git a/PartitionMetis.h b/PartitionMetis.h index f9cc72f..cb0da9b 100644 --- a/PartitionMetis.h +++ b/PartitionMetis.h @@ -40,13 +40,28 @@ class PartitionMetis const idx_t m_numCells; + int rank = 0; + int procs = 0; + + idx_t numflag = 0; + idx_t ncommonnodes = 3; // TODO adapt for hex + idx_t nparts = 0; // will be set to procs in constructor later + + std::vector vtxdist; + std::vector xadj; + std::vector adjncy; + public: PartitionMetis(const cell_t* cells, unsigned int numCells) : #ifdef USE_MPI m_comm(MPI_COMM_WORLD), #endif // USE_MPI m_cells(cells), m_numCells(numCells) - { } + { + MPI_Comm_rank(m_comm, &rank); + MPI_Comm_size(m_comm, &procs); + nparts = procs; + } #ifdef USE_MPI void setComm(MPI_Comm comm) @@ -55,6 +70,87 @@ class PartitionMetis } #endif // USE_MPI +#ifdef USE_MPI + void generateGraphFromMesh() { + assert((vtxdist.empty() && xadj.empty() && adjncy.empty()) || (!vtxdist.empty() && !xadj.empty() && !adjncy.empty())); + + if (vtxdist.empty() && xadj.empty() && adjncy.empty()) { + std::vector elemdist; + elemdist.resize(procs + 1); + std::fill(elemdist.begin(), elemdist.end(), 0); + + MPI_Allgather(const_cast(&m_numCells), 1, IDX_T, elemdist.data(), 1, IDX_T, m_comm); + + idx_t sum = 0; + for (int i = 0; i < procs; i++) { + idx_t e = elemdist[i]; + elemdist[i] = sum; + sum += e; + } + elemdist[procs] = sum; + + std::vector eptr; + eptr.resize(m_numCells + 1); + std::fill(eptr.begin(), eptr.end(), 0); + + std::vector eind; + eind.resize(m_numCells * internal::Topology::cellvertices()); + std::fill(eind.begin(), eind.end(), 0); + + unsigned long m = 0; + + for (idx_t i = 0; i < m_numCells; i++) { + eptr[i] = i * internal::Topology::cellvertices(); + + for (unsigned int j = 0; j < internal::Topology::cellvertices(); j++) { + m = std::max(m, m_cells[i][j]); + eind[i * internal::Topology::cellvertices() + j] = m_cells[i][j]; + } + } + + eptr[m_numCells] = m_numCells * internal::Topology::cellvertices(); + + idx_t* metis_xadj; + idx_t* metis_adjncy; + + ParMETIS_V3_Mesh2Dual(elemdist.data(), eptr.data(), eind.data(), &numflag, &ncommonnodes, &metis_xadj, + &metis_adjncy, &m_comm); + + vtxdist = std::move(elemdist); + + // the size of xadj is the + // - vtxdist[proc] + vtxdist[proc+1] + // because proc has the index proc to proc +1 elements + + assert(vtxdist.size() == static_cast(procs + 1)); + // the first element is always 0 and on top of that we have n nodes + size_t numElements = vtxdist[rank + 1] - vtxdist[rank] + 1; + xadj.reserve(numElements); + std::copy(metis_xadj, metis_xadj + numElements, std::back_inserter(xadj)); + + // last element of xadj will be the size of adjncy + size_t adjncySize = xadj[numElements - 1]; + adjncy.reserve(adjncySize); + std::copy(metis_adjncy, metis_adjncy + adjncySize, std::back_inserter(adjncy)); + + + + METIS_Free(metis_xadj); + METIS_Free(metis_adjncy); + } + } +#endif + +#ifdef USE_MPI + std::tuple&, const std::vector&, const std::vector&> getGraph() { + if (xadj.empty() && adjncy.empty()) { + generateGraphFromMesh(); + } + + return {vtxdist, xadj, adjncy}; + } +#endif + #ifdef USE_MPI enum Status { Ok, @@ -73,99 +169,86 @@ class PartitionMetis const int* vertexWeights = nullptr, const double* imbalances = nullptr, int nWeightsPerVertex = 1, - const double* nodeWeights = nullptr) { - int rank, procs; - MPI_Comm_rank(m_comm, &rank); - MPI_Comm_size(m_comm, &procs); - - idx_t* elemdist = new idx_t[procs+1]; - MPI_Allgather(const_cast(&m_numCells), 1, IDX_T, elemdist, 1, IDX_T, m_comm); - idx_t sum = 0; - for (int i = 0; i < procs; i++) { - idx_t e = elemdist[i]; - elemdist[i] = sum; - sum += e; - } - elemdist[procs] = sum; + const double* nodeWeights = nullptr, + const int* edgeWeights = nullptr, + size_t edgeCount = 0) + { + generateGraphFromMesh(); - idx_t* eptr = new idx_t[m_numCells+1]; - idx_t* eind = new idx_t[m_numCells * internal::Topology::cellvertices()]; - unsigned long m = 0; - for (idx_t i = 0; i < m_numCells; i++) { - eptr[i] = i * internal::Topology::cellvertices(); + idx_t wgtflag = 0; - for (unsigned int j = 0; j < internal::Topology::cellvertices(); j++) { - m = std::max(m, m_cells[i][j]); - eind[i*internal::Topology::cellvertices() + j] = m_cells[i][j]; - } + // set the flag + if (nodeWeights == nullptr && edgeWeights == nullptr) { + wgtflag = 0; + } else if (nodeWeights != nullptr && edgeWeights != nullptr) { + wgtflag = 3; + } else if (nodeWeights == nullptr && edgeWeights != nullptr) { + wgtflag = 1; + } else { + wgtflag = 2; } - eptr[m_numCells] = m_numCells * internal::Topology::cellvertices(); - - idx_t wgtflag = 0; + idx_t ncon = nWeightsPerVertex; idx_t* elmwgt = nullptr; if (vertexWeights != nullptr) { - wgtflag = 2; - elmwgt = new idx_t[m_numCells * ncon]; - for (idx_t cell = 0; cell < m_numCells; ++cell) { - for (idx_t j = 0; j < ncon; ++j) { - elmwgt[ncon * cell + j] = static_cast(vertexWeights[ncon * cell + j]); - } + elmwgt = new idx_t[m_numCells * ncon]; + for (idx_t cell = 0; cell < m_numCells; ++cell) { + for (idx_t j = 0; j < ncon; ++j) { + elmwgt[ncon * cell + j] = static_cast(vertexWeights[ncon * cell + j]); } } + } + + idx_t* edgewgt = nullptr; + if (edgeWeights != nullptr) { + assert(edgeCount != 0); + edgewgt = new idx_t[edgeCount]; + for (size_t i = 0; i < edgeCount; ++i) { + edgewgt[i] = static_cast(edgeWeights[i]); + } + } idx_t numflag = 0; - idx_t ncommonnodes = 3; // TODO adapt for hex idx_t nparts = procs; real_t* tpwgts = new real_t[nparts * ncon]; if (nodeWeights != nullptr) { - for (idx_t i = 0; i < nparts; i++) { - for (idx_t j = 0; j < ncon; ++j) { - tpwgts[i*ncon + j] = nodeWeights[i]; - } + for (idx_t i = 0; i < nparts; i++) { + for (idx_t j = 0; j < ncon; ++j) { + tpwgts[i * ncon + j] = nodeWeights[i]; } + } } else { - for (idx_t i = 0; i < nparts * ncon; i++) { - tpwgts[i] = static_cast(1.) / nparts; - } + for (idx_t i = 0; i < nparts * ncon; i++) { + tpwgts[i] = static_cast(1.) / nparts; + } } real_t* ubvec = new real_t[ncon]; for (idx_t i = 0; i < ncon; ++i) { - ubvec[i] = imbalances[i]; + ubvec[i] = imbalances[i]; } + idx_t edgecut; idx_t options[3] = {1, 1, METIS_RANDOM_SEED}; idx_t* part = new idx_t[m_numCells]; - auto metisResult = ParMETIS_V3_PartMeshKway(elemdist, - eptr, - eind, - elmwgt, - &wgtflag, - &numflag, - &ncon, - &ncommonnodes, - &nparts, - tpwgts, - ubvec, - options, - &edgecut, - part, - &m_comm); - - delete [] elemdist; - delete [] eptr; - delete [] eind; - delete [] tpwgts; - delete [] ubvec; + assert(xadj.size() == static_cast(vtxdist[rank + 1] - vtxdist[rank] + 1)); + assert(adjncy.size() == static_cast(xadj.back())); + + auto metisResult = ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjncy.data(), elmwgt, edgewgt, &wgtflag, + &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &m_comm); + + delete[] tpwgts; + delete[] ubvec; + delete[] elmwgt; + delete[] edgewgt; for (idx_t i = 0; i < m_numCells; i++) - partition[i] = part[i]; + partition[i] = part[i]; - delete [] part; + delete[] part; return (metisResult == METIS_OK) ? Status::Ok : Status::Error; } From d972ecbd971b3e5e5b5ff31273d3b31cb3d1c46f Mon Sep 17 00:00:00 2001 From: thrud primrose Date: Sat, 21 May 2022 01:54:31 +0200 Subject: [PATCH 2/4] improve formatting --- PartitionMetis.h | 52 +++++++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/PartitionMetis.h b/PartitionMetis.h index cb0da9b..2921eb7 100644 --- a/PartitionMetis.h +++ b/PartitionMetis.h @@ -179,33 +179,34 @@ class PartitionMetis // set the flag if (nodeWeights == nullptr && edgeWeights == nullptr) { - wgtflag = 0; + wgtflag = 0; } else if (nodeWeights != nullptr && edgeWeights != nullptr) { - wgtflag = 3; + wgtflag = 3; } else if (nodeWeights == nullptr && edgeWeights != nullptr) { - wgtflag = 1; + wgtflag = 1; } else { - wgtflag = 2; + wgtflag = 2; } idx_t ncon = nWeightsPerVertex; idx_t* elmwgt = nullptr; + if (vertexWeights != nullptr) { - elmwgt = new idx_t[m_numCells * ncon]; - for (idx_t cell = 0; cell < m_numCells; ++cell) { - for (idx_t j = 0; j < ncon; ++j) { - elmwgt[ncon * cell + j] = static_cast(vertexWeights[ncon * cell + j]); + elmwgt = new idx_t[m_numCells * ncon]; + for (idx_t cell = 0; cell < m_numCells; ++cell) { + for (idx_t j = 0; j < ncon; ++j) { + elmwgt[ncon * cell + j] = static_cast(vertexWeights[ncon * cell + j]); + } } } - } idx_t* edgewgt = nullptr; if (edgeWeights != nullptr) { - assert(edgeCount != 0); - edgewgt = new idx_t[edgeCount]; - for (size_t i = 0; i < edgeCount; ++i) { - edgewgt[i] = static_cast(edgeWeights[i]); - } + assert(edgeCount != 0); + edgewgt = new idx_t[edgeCount]; + for (size_t i = 0; i < edgeCount; ++i) { + edgewgt[i] = static_cast(edgeWeights[i]); + } } idx_t numflag = 0; @@ -213,20 +214,20 @@ class PartitionMetis real_t* tpwgts = new real_t[nparts * ncon]; if (nodeWeights != nullptr) { - for (idx_t i = 0; i < nparts; i++) { - for (idx_t j = 0; j < ncon; ++j) { - tpwgts[i * ncon + j] = nodeWeights[i]; + for (idx_t i = 0; i < nparts; i++) { + for (idx_t j = 0; j < ncon; ++j) { + tpwgts[i * ncon + j] = nodeWeights[i]; + } } - } } else { - for (idx_t i = 0; i < nparts * ncon; i++) { - tpwgts[i] = static_cast(1.) / nparts; - } + for (idx_t i = 0; i < nparts * ncon; i++) { + tpwgts[i] = static_cast(1.) / nparts; + } } real_t* ubvec = new real_t[ncon]; for (idx_t i = 0; i < ncon; ++i) { - ubvec[i] = imbalances[i]; + ubvec[i] = imbalances[i]; } idx_t edgecut; @@ -245,9 +246,10 @@ class PartitionMetis delete[] elmwgt; delete[] edgewgt; - for (idx_t i = 0; i < m_numCells; i++) - partition[i] = part[i]; - + for (idx_t i = 0; i < m_numCells; i++){ + partition[i] = part[i]; + } + delete[] part; return (metisResult == METIS_OK) ? Status::Ok : Status::Error; From 5c40391fc9431b61c3033ed61659b44a1e5cc087 Mon Sep 17 00:00:00 2001 From: thrudprimrose Date: Tue, 12 Jul 2022 12:52:59 +0200 Subject: [PATCH 3/4] use m_ prefix for all private members, remove duplicate member procs, adapt assertions and functions with the new names, and using tabs --- PartitionMetis.h | 154 +++++++++++++++++++++++------------------------ 1 file changed, 76 insertions(+), 78 deletions(-) diff --git a/PartitionMetis.h b/PartitionMetis.h index 2921eb7..4dde495 100644 --- a/PartitionMetis.h +++ b/PartitionMetis.h @@ -40,16 +40,18 @@ class PartitionMetis const idx_t m_numCells; - int rank = 0; - int procs = 0; +#ifdef USE_MPI + int m_rank = 0; + int m_nparts = 0; +#endif - idx_t numflag = 0; - idx_t ncommonnodes = 3; // TODO adapt for hex - idx_t nparts = 0; // will be set to procs in constructor later + // Better the use the same convention as METIS, even though we have the prefix m_ + idx_t m_numflag = 0; // Since we are a C++ library, numflag for Metis will be always 0 + idx_t m_ncommonnodes = 3; // TODO adapt for hex, but until then we can use constexpr for these - std::vector vtxdist; - std::vector xadj; - std::vector adjncy; + std::vector m_vtxdist; + std::vector m_xadj; + std::vector m_adjncy; public: PartitionMetis(const cell_t* cells, unsigned int numCells) : @@ -58,9 +60,10 @@ class PartitionMetis #endif // USE_MPI m_cells(cells), m_numCells(numCells) { - MPI_Comm_rank(m_comm, &rank); - MPI_Comm_size(m_comm, &procs); - nparts = procs; +#ifdef USE_MPI + MPI_Comm_rank(m_comm, &m_rank); + MPI_Comm_size(m_comm, &m_nparts); +#endif } #ifdef USE_MPI @@ -72,82 +75,80 @@ class PartitionMetis #ifdef USE_MPI void generateGraphFromMesh() { - assert((vtxdist.empty() && xadj.empty() && adjncy.empty()) || (!vtxdist.empty() && !xadj.empty() && !adjncy.empty())); - - if (vtxdist.empty() && xadj.empty() && adjncy.empty()) { - std::vector elemdist; - elemdist.resize(procs + 1); - std::fill(elemdist.begin(), elemdist.end(), 0); + assert((m_vtxdist.empty() && m_xadj.empty() && m_adjncy.empty()) || (!m_vtxdist.empty() && !m_xadj.empty() && !m_adjncy.empty())); - MPI_Allgather(const_cast(&m_numCells), 1, IDX_T, elemdist.data(), 1, IDX_T, m_comm); + if (m_vtxdist.empty() && m_xadj.empty() && m_adjncy.empty()) { + std::vector elemdist; + elemdist.resize(m_nparts + 1); + std::fill(elemdist.begin(), elemdist.end(), 0); - idx_t sum = 0; - for (int i = 0; i < procs; i++) { - idx_t e = elemdist[i]; - elemdist[i] = sum; - sum += e; - } - elemdist[procs] = sum; + MPI_Allgather(const_cast(&m_numCells), 1, IDX_T, elemdist.data(), 1, IDX_T, m_comm); - std::vector eptr; - eptr.resize(m_numCells + 1); - std::fill(eptr.begin(), eptr.end(), 0); + idx_t sum = 0; + for (int i = 0; i < m_nparts; i++) { + idx_t e = elemdist[i]; + elemdist[i] = sum; + sum += e; + } + elemdist[m_nparts] = sum; - std::vector eind; - eind.resize(m_numCells * internal::Topology::cellvertices()); - std::fill(eind.begin(), eind.end(), 0); + std::vector eptr; + eptr.resize(m_numCells + 1); + std::fill(eptr.begin(), eptr.end(), 0); - unsigned long m = 0; + std::vector eind; + eind.resize(m_numCells * internal::Topology::cellvertices()); + std::fill(eind.begin(), eind.end(), 0); - for (idx_t i = 0; i < m_numCells; i++) { - eptr[i] = i * internal::Topology::cellvertices(); + unsigned long m = 0; - for (unsigned int j = 0; j < internal::Topology::cellvertices(); j++) { - m = std::max(m, m_cells[i][j]); - eind[i * internal::Topology::cellvertices() + j] = m_cells[i][j]; - } - } + for (idx_t i = 0; i < m_numCells; i++) { + eptr[i] = i * internal::Topology::cellvertices(); - eptr[m_numCells] = m_numCells * internal::Topology::cellvertices(); - - idx_t* metis_xadj; - idx_t* metis_adjncy; + for (unsigned int j = 0; j < internal::Topology::cellvertices(); j++) { + m = std::max(m, m_cells[i][j]); + eind[i * internal::Topology::cellvertices() + j] = m_cells[i][j]; + } + } - ParMETIS_V3_Mesh2Dual(elemdist.data(), eptr.data(), eind.data(), &numflag, &ncommonnodes, &metis_xadj, - &metis_adjncy, &m_comm); + eptr[m_numCells] = m_numCells * internal::Topology::cellvertices(); - vtxdist = std::move(elemdist); + idx_t* metis_xadj; + idx_t* metis_adjncy; - // the size of xadj is the - // - vtxdist[proc] + vtxdist[proc+1] - // because proc has the index proc to proc +1 elements + ParMETIS_V3_Mesh2Dual(elemdist.data(), eptr.data(), eind.data(), &m_numflag, &m_ncommonnodes, &metis_xadj, + &metis_adjncy, &m_comm); - assert(vtxdist.size() == static_cast(procs + 1)); - // the first element is always 0 and on top of that we have n nodes - size_t numElements = vtxdist[rank + 1] - vtxdist[rank] + 1; - xadj.reserve(numElements); - std::copy(metis_xadj, metis_xadj + numElements, std::back_inserter(xadj)); + m_vtxdist = std::move(elemdist); - // last element of xadj will be the size of adjncy - size_t adjncySize = xadj[numElements - 1]; - adjncy.reserve(adjncySize); - std::copy(metis_adjncy, metis_adjncy + adjncySize, std::back_inserter(adjncy)); + // the size of xadj is the + // - vtxdist[proc] + vtxdist[proc+1] + // because proc has the index proc to proc +1 elements + assert(m_vtxdist.size() == static_cast(m_nparts + 1)); + // the first element is always 0 and on top of that we have n nodes + size_t numElements = m_vtxdist[m_rank + 1] - m_vtxdist[m_rank] + 1; + m_xadj.reserve(numElements); + std::copy(metis_xadj, metis_xadj + numElements, std::back_inserter(m_xadj)); + // last element of xadj will be the size of adjncy + size_t adjncySize = m_xadj[numElements - 1]; + m_adjncy.reserve(adjncySize); + std::copy(metis_adjncy, metis_adjncy + adjncySize, std::back_inserter(m_adjncy)); - METIS_Free(metis_xadj); - METIS_Free(metis_adjncy); - } + METIS_Free(metis_xadj); + METIS_Free(metis_adjncy); + } } #endif #ifdef USE_MPI std::tuple&, const std::vector&, const std::vector&> getGraph() { - if (xadj.empty() && adjncy.empty()) { - generateGraphFromMesh(); - } + if (m_xadj.empty() && m_adjncy.empty()) { + generateGraphFromMesh(); + } - return {vtxdist, xadj, adjncy}; + return {m_vtxdist, m_xadj, m_adjncy}; } #endif @@ -171,7 +172,7 @@ class PartitionMetis int nWeightsPerVertex = 1, const double* nodeWeights = nullptr, const int* edgeWeights = nullptr, - size_t edgeCount = 0) + size_t edgeCount = 0) { generateGraphFromMesh(); @@ -209,19 +210,16 @@ class PartitionMetis } } - idx_t numflag = 0; - idx_t nparts = procs; - - real_t* tpwgts = new real_t[nparts * ncon]; + real_t* tpwgts = new real_t[m_nparts* ncon]; if (nodeWeights != nullptr) { - for (idx_t i = 0; i < nparts; i++) { + for (idx_t i = 0; i < m_nparts; i++) { for (idx_t j = 0; j < ncon; ++j) { tpwgts[i * ncon + j] = nodeWeights[i]; } } } else { - for (idx_t i = 0; i < nparts * ncon; i++) { - tpwgts[i] = static_cast(1.) / nparts; + for (idx_t i = 0; i < m_nparts * ncon; i++) { + tpwgts[i] = static_cast(1.) / m_nparts; } } @@ -235,11 +233,11 @@ class PartitionMetis idx_t* part = new idx_t[m_numCells]; - assert(xadj.size() == static_cast(vtxdist[rank + 1] - vtxdist[rank] + 1)); - assert(adjncy.size() == static_cast(xadj.back())); + assert(m_xadj.size() == static_cast(m_vtxdist[m_rank + 1] - m_vtxdist[m_rank] + 1)); + assert(m_adjncy.size() == static_cast(m_xadj.back())); - auto metisResult = ParMETIS_V3_PartKway(vtxdist.data(), xadj.data(), adjncy.data(), elmwgt, edgewgt, &wgtflag, - &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &m_comm); + auto metisResult = ParMETIS_V3_PartKway(m_vtxdist.data(), m_xadj.data(), m_adjncy.data(), elmwgt, edgewgt, &wgtflag, + &m_numflag, &ncon, &m_nparts, tpwgts, ubvec, options, &edgecut, part, &m_comm); delete[] tpwgts; delete[] ubvec; @@ -247,7 +245,7 @@ class PartitionMetis delete[] edgewgt; for (idx_t i = 0; i < m_numCells; i++){ - partition[i] = part[i]; + partition[i] = static_cast(part[i]); } delete[] part; From d69a89670b21e00248f4aa077dc302efeb0685fc Mon Sep 17 00:00:00 2001 From: thrudprimrose Date: Tue, 12 Jul 2022 12:54:59 +0200 Subject: [PATCH 4/4] remove unused parameter m --- PartitionMetis.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/PartitionMetis.h b/PartitionMetis.h index 4dde495..fb8b122 100644 --- a/PartitionMetis.h +++ b/PartitionMetis.h @@ -100,13 +100,10 @@ class PartitionMetis eind.resize(m_numCells * internal::Topology::cellvertices()); std::fill(eind.begin(), eind.end(), 0); - unsigned long m = 0; - for (idx_t i = 0; i < m_numCells; i++) { eptr[i] = i * internal::Topology::cellvertices(); for (unsigned int j = 0; j < internal::Topology::cellvertices(); j++) { - m = std::max(m, m_cells[i][j]); eind[i * internal::Topology::cellvertices() + j] = m_cells[i][j]; } }