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