Actual source code: partparmetis.c
1: #include <petsc/private/partitionerimpl.h>
2: #include <petsc/private/matmetisimpl.h>
3: #if PetscDefined(HAVE_PARMETIS)
4: #include <petsc/private/matparmetisimpl.h>
5: #include <parmetis.h>
6: #endif
8: PetscBool ParMETISPartitionerCite = PETSC_FALSE;
9: const char ParMETISPartitionerCitation[] = "@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";
19: typedef struct {
20: MPI_Comm pcomm;
21: PetscInt ptype;
22: PetscReal imbalanceRatio;
23: PetscInt debugFlag;
24: PetscInt randomSeed;
25: } PetscPartitioner_ParMETIS;
27: static const char *ptypes[] = {"kway", "rb"};
29: static PetscErrorCode PetscPartitionerDestroy_ParMETIS(PetscPartitioner part)
30: {
31: PetscPartitioner_ParMETIS *p = (PetscPartitioner_ParMETIS *)part->data;
33: PetscFunctionBegin;
34: PetscCallMPI(MPI_Comm_free(&p->pcomm));
35: PetscCall(PetscFree(part->data));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode PetscPartitionerView_ParMETIS_ASCII(PetscPartitioner part, PetscViewer viewer)
40: {
41: PetscPartitioner_ParMETIS *p = (PetscPartitioner_ParMETIS *)part->data;
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(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PetscPartitionerView_ParMETIS(PetscPartitioner part, PetscViewer viewer)
54: {
55: PetscBool isascii;
57: PetscFunctionBegin;
60: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
61: if (isascii) PetscCall(PetscPartitionerView_ParMETIS_ASCII(part, viewer));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode PetscPartitionerSetFromOptions_ParMETIS(PetscPartitioner part, PetscOptionItems PetscOptionsObject)
66: {
67: PetscPartitioner_ParMETIS *p = (PetscPartitioner_ParMETIS *)part->data;
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(PETSC_SUCCESS);
77: }
79: static PetscErrorCode PetscPartitionerPartition_ParMETIS(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection edgeSection, PetscSection targetSection, PetscSection partSection, IS *partition)
80: {
81: #if PetscDefined(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;
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(PETSC_SUCCESS);
118: }
119: /* Calculate partition weights */
120: if (targetSection) {
121: PetscInt p;
122: real_t sumt = 0.0;
124: for (p = 0; p < nparts; ++p) {
125: PetscInt tpd;
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] = (real_t)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] = (real_t)(1. - sumt);
139: }
140: } else {
141: for (p = 0; p < nparts; ++p) tpwgts[p] = (real_t)(1.0 / nparts);
142: }
143: ubvec[0] = (real_t)pm->imbalanceRatio;
145: /* Weight cells */
146: if (vertSection) {
147: PetscCall(PetscMalloc1(nvtxs, &vwgt));
148: for (v = 0; v < nvtxs; ++v) PetscCall(PetscSectionGetDof(vertSection, v, &vwgt[v]));
149: wgtflag |= 2; /* have weights on graph vertices */
150: }
151: // Weight edges
152: if (edgeSection) {
153: PetscCall(PetscMalloc1(xadj[nvtxs], &adjwgt));
154: for (PetscInt e = 0; e < xadj[nvtxs]; ++e) PetscCall(PetscSectionGetDof(edgeSection, e, &adjwgt[e]));
155: wgtflag |= 1; /* have weights on graph edges */
156: }
158: for (p = 0; !vtxdist[p + 1] && p < size; ++p);
159: if (vtxdist[p + 1] == vtxdist[size]) {
160: if (rank == p) {
161: PetscCallMETIS(METIS_SetDefaultOptions, options);
162: options[METIS_OPTION_DBGLVL] = pm->debugFlag;
163: options[METIS_OPTION_SEED] = pm->randomSeed;
164: /*
165: It would be nice to activate the two METIS_PartGraphKway() options below, but they would need some actual testing.
166: - Turning on these options may exercise path of the METIS code that have bugs and may break production runs.
167: - 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.
168: */
169: /* options[METIS_OPTION_CONTIG] = 1; */ /* try to produce partitions that are contiguous */
170: /* options[METIS_OPTION_MINCONN] = 1; */ /* minimize the maximum degree of the subdomain graph */
171: if (metis_ptype == 1) PetscCallMETIS(METIS_PartGraphRecursive, &nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
172: else PetscCallMETIS(METIS_PartGraphKway, &nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
173: }
174: } else {
175: MPI_Comm pcomm = pm->pcomm;
177: options[0] = 1; /*use options */
178: options[1] = pm->debugFlag;
179: options[2] = (pm->randomSeed == -1) ? 15 : pm->randomSeed; /* default is GLOBAL_SEED=15 from `libparmetis/defs.h` */
181: if (hasempty) { /* ParMETIS does not support empty graphs on some of the processes */
182: PetscInt cnt;
184: PetscCallMPI(MPI_Comm_split(pm->pcomm, !!nvtxs, rank, &pcomm));
185: for (p = 0, cnt = 0; p < size; p++) {
186: if (vtxdist[p + 1] != vtxdist[p]) {
187: vtxdist[cnt + 1] = vtxdist[p + 1];
188: cnt++;
189: }
190: }
191: }
192: if (nvtxs) PetscCallParMETIS(ParMETIS_V3_PartKway, vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment, &pcomm);
193: if (hasempty) PetscCallMPI(MPI_Comm_free(&pcomm));
194: }
196: /* Convert to PetscSection+IS */
197: for (v = 0; v < nvtxs; ++v) PetscCall(PetscSectionAddDof(partSection, assignment[v], 1));
198: PetscCall(PetscMalloc1(nvtxs, &points));
199: for (p = 0, i = 0; p < nparts; ++p) {
200: for (v = 0; v < nvtxs; ++v) {
201: if (assignment[v] == p) points[i++] = v;
202: }
203: }
204: PetscCheck(i == nvtxs, comm, PETSC_ERR_PLIB, "Number of points %" PetscInt_FMT " should be %" PetscInt_FMT, i, nvtxs);
205: PetscCall(ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition));
206: PetscCall(PetscFree4(vtxdist, tpwgts, ubvec, assignment));
207: PetscCall(PetscFree(vwgt));
208: PetscCall(PetscFree(adjwgt));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: #else
211: SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
212: #endif
213: }
215: static PetscErrorCode PetscPartitionerInitialize_ParMETIS(PetscPartitioner part)
216: {
217: PetscFunctionBegin;
218: part->noGraph = PETSC_FALSE;
219: part->ops->view = PetscPartitionerView_ParMETIS;
220: part->ops->setfromoptions = PetscPartitionerSetFromOptions_ParMETIS;
221: part->ops->destroy = PetscPartitionerDestroy_ParMETIS;
222: part->ops->partition = PetscPartitionerPartition_ParMETIS;
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*MC
227: PETSCPARTITIONERPARMETIS = "parmetis" - A PetscPartitioner object using the ParMETIS library
229: Level: intermediate
231: Options Database Keys:
232: + -petscpartitioner_parmetis_type (kway|rb) - ParMETIS partitioning type.
233: . -petscpartitioner_parmetis_imbalance_ratio value - Load imbalance ratio limit
234: . -petscpartitioner_parmetis_debug level - Debugging flag passed to ParMETIS/METIS routines
235: - -petscpartitioner_parmetis_seed seed - Random seed
237: Notes: when the graph is on a single process, this partitioner actually calls METIS and not ParMETIS
239: .seealso: `PetscPartitionerType`, `PetscPartitionerCreate()`, `PetscPartitionerSetType()`
240: M*/
242: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMETIS(PetscPartitioner part)
243: {
244: PetscPartitioner_ParMETIS *p;
246: PetscFunctionBegin;
248: PetscCall(PetscNew(&p));
249: part->data = p;
251: PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)part), &p->pcomm));
252: p->ptype = 0;
253: p->imbalanceRatio = 1.05;
254: p->debugFlag = 0;
255: p->randomSeed = -1; /* defaults to GLOBAL_SEED=15 from `libparmetis/defs.h` */
257: PetscCall(PetscPartitionerInitialize_ParMETIS(part));
258: PetscCall(PetscCitationsRegister(ParMETISPartitionerCitation, &ParMETISPartitionerCite));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }