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
PetscPartitionerDestroy_ParMetis(PetscPartitioner part)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
PetscPartitionerView_ParMetis_ASCII(PetscPartitioner part,PetscViewer viewer)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
PetscPartitionerView_ParMetis(PetscPartitioner part,PetscViewer viewer)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
PetscPartitionerSetFromOptions_ParMetis(PetscPartitioner part,PetscOptionItems PetscOptionsObject)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
PetscPartitionerPartition_ParMetis(PetscPartitioner part,PetscInt nparts,PetscInt numVertices,PetscInt start[],PetscInt adjacency[],PetscSection vertSection,PetscSection edgeSection,PetscSection targetSection,PetscSection partSection,IS * partition)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
PetscPartitionerInitialize_ParMetis(PetscPartitioner part)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
PetscPartitionerCreate_ParMetis(PetscPartitioner part)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