1 #include <petsc/private/partitionerimpl.h> /*I "petscpartitioner.h" I*/ 2 3 #if defined(PETSC_HAVE_PARMETIS) 4 #include <parmetis.h> 5 #endif 6 7 PetscBool ParMetisPartitionerCite = PETSC_FALSE; 8 const char ParMetisPartitionerCitation[] = "@article{KarypisKumar98,\n" 9 " author = {George Karypis and Vipin Kumar},\n" 10 " title = {A Parallel Algorithm for Multilevel Graph Partitioning and Sparse Matrix Ordering},\n" 11 " journal = {Journal of Parallel and Distributed Computing},\n" 12 " volume = {48},\n" 13 " pages = {71--85},\n" 14 " year = {1998}\n" 15 " doi = {https://doi.org/10.1006/jpdc.1997.1403}\n" 16 "}\n"; 17 18 typedef struct { 19 MPI_Comm pcomm; 20 PetscInt ptype; 21 PetscReal imbalanceRatio; 22 PetscInt debugFlag; 23 PetscInt randomSeed; 24 } PetscPartitioner_ParMetis; 25 26 static const char *ptypes[] = {"kway", "rb"}; 27 28 static PetscErrorCode PetscPartitionerDestroy_ParMetis(PetscPartitioner part) 29 { 30 PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *)part->data; 31 32 PetscFunctionBegin; 33 PetscCallMPI(MPI_Comm_free(&p->pcomm)); 34 PetscCall(PetscFree(part->data)); 35 PetscFunctionReturn(PETSC_SUCCESS); 36 } 37 38 static PetscErrorCode PetscPartitionerView_ParMetis_ASCII(PetscPartitioner part, PetscViewer viewer) 39 { 40 PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *)part->data; 41 42 PetscFunctionBegin; 43 PetscCall(PetscViewerASCIIPushTab(viewer)); 44 PetscCall(PetscViewerASCIIPrintf(viewer, "ParMetis type: %s\n", ptypes[p->ptype])); 45 PetscCall(PetscViewerASCIIPrintf(viewer, "load imbalance ratio %g\n", (double)p->imbalanceRatio)); 46 PetscCall(PetscViewerASCIIPrintf(viewer, "debug flag %" PetscInt_FMT "\n", p->debugFlag)); 47 PetscCall(PetscViewerASCIIPrintf(viewer, "random seed %" PetscInt_FMT "\n", p->randomSeed)); 48 PetscCall(PetscViewerASCIIPopTab(viewer)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode PetscPartitionerView_ParMetis(PetscPartitioner part, PetscViewer viewer) 53 { 54 PetscBool isascii; 55 56 PetscFunctionBegin; 57 PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1); 58 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 59 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 60 if (isascii) PetscCall(PetscPartitionerView_ParMetis_ASCII(part, viewer)); 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 static PetscErrorCode PetscPartitionerSetFromOptions_ParMetis(PetscPartitioner part, PetscOptionItems PetscOptionsObject) 65 { 66 PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *)part->data; 67 68 PetscFunctionBegin; 69 PetscOptionsHeadBegin(PetscOptionsObject, "PetscPartitioner ParMetis Options"); 70 PetscCall(PetscOptionsEList("-petscpartitioner_parmetis_type", "Partitioning method", "", ptypes, 2, ptypes[p->ptype], &p->ptype, NULL)); 71 PetscCall(PetscOptionsReal("-petscpartitioner_parmetis_imbalance_ratio", "Load imbalance ratio limit", "", p->imbalanceRatio, &p->imbalanceRatio, NULL)); 72 PetscCall(PetscOptionsInt("-petscpartitioner_parmetis_debug", "Debugging flag", "", p->debugFlag, &p->debugFlag, NULL)); 73 PetscCall(PetscOptionsInt("-petscpartitioner_parmetis_seed", "Random seed", "", p->randomSeed, &p->randomSeed, NULL)); 74 PetscOptionsHeadEnd(); 75 PetscFunctionReturn(PETSC_SUCCESS); 76 } 77 78 static PetscErrorCode PetscPartitionerPartition_ParMetis(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection edgeSection, PetscSection targetSection, PetscSection partSection, IS *partition) 79 { 80 #if defined(PETSC_HAVE_PARMETIS) 81 PetscPartitioner_ParMetis *pm = (PetscPartitioner_ParMetis *)part->data; 82 MPI_Comm comm; 83 PetscInt nvtxs = numVertices; /* The number of vertices in full graph */ 84 PetscInt *vtxdist; /* Distribution of vertices across processes */ 85 PetscInt *xadj = start; /* Start of edge list for each vertex */ 86 PetscInt *adjncy = adjacency; /* Edge lists for all vertices */ 87 PetscInt *vwgt = NULL; /* Vertex weights */ 88 PetscInt *adjwgt = NULL; /* Edge weights */ 89 PetscInt wgtflag = 0; /* Indicates which weights are present */ 90 PetscInt numflag = 0; /* Indicates initial offset (0 or 1) */ 91 PetscInt ncon = 1; /* The number of weights per vertex */ 92 PetscInt metis_ptype = pm->ptype; /* kway or recursive bisection */ 93 real_t *tpwgts; /* The fraction of vertex weights assigned to each partition */ 94 real_t *ubvec; /* The balance intolerance for vertex weights */ 95 PetscInt options[64]; /* Options */ 96 PetscInt v, i, *assignment, *points; 97 PetscMPIInt p, size, rank; 98 PetscBool hasempty = PETSC_FALSE; 99 100 PetscFunctionBegin; 101 PetscCall(PetscObjectGetComm((PetscObject)part, &comm)); 102 PetscCallMPI(MPI_Comm_size(comm, &size)); 103 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 104 /* Calculate vertex distribution */ 105 PetscCall(PetscMalloc4(size + 1, &vtxdist, nparts * ncon, &tpwgts, ncon, &ubvec, nvtxs, &assignment)); 106 vtxdist[0] = 0; 107 PetscCallMPI(MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm)); 108 for (p = 2; p <= size; ++p) { 109 hasempty = (PetscBool)(hasempty || !vtxdist[p - 1] || !vtxdist[p]); 110 vtxdist[p] += vtxdist[p - 1]; 111 } 112 /* null graph */ 113 if (vtxdist[size] == 0) { 114 PetscCall(PetscFree4(vtxdist, tpwgts, ubvec, assignment)); 115 PetscCall(ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition)); 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 /* Calculate partition weights */ 119 if (targetSection) { 120 PetscInt p; 121 real_t sumt = 0.0; 122 123 for (p = 0; p < nparts; ++p) { 124 PetscInt tpd; 125 126 PetscCall(PetscSectionGetDof(targetSection, p, &tpd)); 127 sumt += tpd; 128 tpwgts[p] = tpd; 129 } 130 if (sumt) { /* METIS/ParMETIS do not like exactly zero weight */ 131 for (p = 0, sumt = 0.0; p < nparts; ++p) { 132 tpwgts[p] = (real_t)PetscMax(tpwgts[p], PETSC_SMALL); 133 sumt += tpwgts[p]; 134 } 135 for (p = 0; p < nparts; ++p) tpwgts[p] /= sumt; 136 for (p = 0, sumt = 0.0; p < nparts - 1; ++p) sumt += tpwgts[p]; 137 tpwgts[nparts - 1] = (real_t)(1. - sumt); 138 } 139 } else { 140 for (p = 0; p < nparts; ++p) tpwgts[p] = (real_t)(1.0 / nparts); 141 } 142 ubvec[0] = (real_t)pm->imbalanceRatio; 143 144 /* Weight cells */ 145 if (vertSection) { 146 PetscCall(PetscMalloc1(nvtxs, &vwgt)); 147 for (v = 0; v < nvtxs; ++v) PetscCall(PetscSectionGetDof(vertSection, v, &vwgt[v])); 148 wgtflag |= 2; /* have weights on graph vertices */ 149 } 150 // Weight edges 151 if (edgeSection) { 152 PetscCall(PetscMalloc1(xadj[nvtxs], &adjwgt)); 153 for (PetscInt e = 0; e < xadj[nvtxs]; ++e) PetscCall(PetscSectionGetDof(edgeSection, e, &adjwgt[e])); 154 wgtflag |= 1; /* have weights on graph edges */ 155 } 156 157 for (p = 0; !vtxdist[p + 1] && p < size; ++p); 158 if (vtxdist[p + 1] == vtxdist[size]) { 159 if (rank == p) { 160 int err; 161 err = METIS_SetDefaultOptions(options); /* initialize all defaults */ 162 options[METIS_OPTION_DBGLVL] = pm->debugFlag; 163 options[METIS_OPTION_SEED] = pm->randomSeed; 164 PetscCheck(err == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()"); 165 if (metis_ptype == 1) { 166 PetscStackPushExternal("METIS_PartGraphRecursive"); 167 err = METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment); 168 PetscStackPop; 169 PetscCheck(err == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphRecursive()"); 170 } else { 171 /* 172 It would be nice to activate the two options below, but they would need some actual testing. 173 - Turning on these options may exercise path of the METIS code that have bugs and may break production runs. 174 - If CONTIG is set to 1, METIS will exit with error if the graph is disconnected, despite the manual saying the option is ignored in such case. 175 */ 176 /* options[METIS_OPTION_CONTIG] = 1; */ /* try to produce partitions that are contiguous */ 177 /* options[METIS_OPTION_MINCONN] = 1; */ /* minimize the maximum degree of the subdomain graph */ 178 PetscStackPushExternal("METIS_PartGraphKway"); 179 err = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment); 180 PetscStackPop; 181 PetscCheck(err == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()"); 182 } 183 } 184 } else { 185 MPI_Comm pcomm = pm->pcomm; 186 187 options[0] = 1; /*use options */ 188 options[1] = pm->debugFlag; 189 options[2] = (pm->randomSeed == -1) ? 15 : pm->randomSeed; /* default is GLOBAL_SEED=15 from `libparmetis/defs.h` */ 190 191 if (hasempty) { /* parmetis does not support empty graphs on some of the processes */ 192 PetscInt cnt; 193 194 PetscCallMPI(MPI_Comm_split(pm->pcomm, !!nvtxs, rank, &pcomm)); 195 for (p = 0, cnt = 0; p < size; p++) { 196 if (vtxdist[p + 1] != vtxdist[p]) { 197 vtxdist[cnt + 1] = vtxdist[p + 1]; 198 cnt++; 199 } 200 } 201 } 202 if (nvtxs) { 203 int err; 204 PetscStackPushExternal("ParMETIS_V3_PartKway"); 205 err = ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment, &pcomm); 206 PetscStackPop; 207 PetscCheck(err == METIS_OK, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error %d in ParMETIS_V3_PartKway()", err); 208 } 209 if (hasempty) PetscCallMPI(MPI_Comm_free(&pcomm)); 210 } 211 212 /* Convert to PetscSection+IS */ 213 for (v = 0; v < nvtxs; ++v) PetscCall(PetscSectionAddDof(partSection, assignment[v], 1)); 214 PetscCall(PetscMalloc1(nvtxs, &points)); 215 for (p = 0, i = 0; p < nparts; ++p) { 216 for (v = 0; v < nvtxs; ++v) { 217 if (assignment[v] == p) points[i++] = v; 218 } 219 } 220 PetscCheck(i == nvtxs, comm, PETSC_ERR_PLIB, "Number of points %" PetscInt_FMT " should be %" PetscInt_FMT, i, nvtxs); 221 PetscCall(ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition)); 222 PetscCall(PetscFree4(vtxdist, tpwgts, ubvec, assignment)); 223 PetscCall(PetscFree(vwgt)); 224 PetscCall(PetscFree(adjwgt)); 225 PetscFunctionReturn(PETSC_SUCCESS); 226 #else 227 SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis."); 228 #endif 229 } 230 231 static PetscErrorCode PetscPartitionerInitialize_ParMetis(PetscPartitioner part) 232 { 233 PetscFunctionBegin; 234 part->noGraph = PETSC_FALSE; 235 part->ops->view = PetscPartitionerView_ParMetis; 236 part->ops->setfromoptions = PetscPartitionerSetFromOptions_ParMetis; 237 part->ops->destroy = PetscPartitionerDestroy_ParMetis; 238 part->ops->partition = PetscPartitionerPartition_ParMetis; 239 PetscFunctionReturn(PETSC_SUCCESS); 240 } 241 242 /*MC 243 PETSCPARTITIONERPARMETIS = "parmetis" - A PetscPartitioner object using the ParMETIS library 244 245 Level: intermediate 246 247 Options Database Keys: 248 + -petscpartitioner_parmetis_type <string> - ParMETIS partitioning type. Either "kway" or "rb" (recursive bisection) 249 . -petscpartitioner_parmetis_imbalance_ratio <value> - Load imbalance ratio limit 250 . -petscpartitioner_parmetis_debug <int> - Debugging flag passed to ParMETIS/METIS routines 251 - -petscpartitioner_parmetis_seed <int> - Random seed 252 253 Notes: when the graph is on a single process, this partitioner actually calls METIS and not ParMETIS 254 255 .seealso: `PetscPartitionerType`, `PetscPartitionerCreate()`, `PetscPartitionerSetType()` 256 M*/ 257 258 PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMetis(PetscPartitioner part) 259 { 260 PetscPartitioner_ParMetis *p; 261 262 PetscFunctionBegin; 263 PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1); 264 PetscCall(PetscNew(&p)); 265 part->data = p; 266 267 PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)part), &p->pcomm)); 268 p->ptype = 0; 269 p->imbalanceRatio = 1.05; 270 p->debugFlag = 0; 271 p->randomSeed = -1; /* defaults to GLOBAL_SEED=15 from `libparmetis/defs.h` */ 272 273 PetscCall(PetscPartitionerInitialize_ParMetis(part)); 274 PetscCall(PetscCitationsRegister(ParMetisPartitionerCitation, &ParMetisPartitionerCite)); 275 PetscFunctionReturn(PETSC_SUCCESS); 276 } 277