Actual source code: pmetis.c
1: #include <../src/mat/impls/adj/mpi/mpiadj.h>
2: #include <petsc/private/matparmetisimpl.h>
3: #include <parmetis.h>
5: /*
6: The first 5 elements of this structure are the input control array to METIS
7: */
8: typedef struct {
9: PetscInt cuts; /* number of cuts made (output) */
10: PetscInt foldfactor;
11: PetscInt parallel; /* use parallel partitioner for coarse problem */
12: PetscInt indexing; /* 0 indicates C indexing, 1 Fortran */
13: PetscInt printout; /* indicates if one wishes METIS to print info */
14: PetscBool repartition;
15: } MatPartitioning_ParMETIS;
17: static PetscErrorCode MatPartitioningApply_ParMETIS_Private(MatPartitioning part, PetscBool useND, PetscBool isImprove, IS *partitioning)
18: {
19: MatPartitioning_ParMETIS *pmetis = (MatPartitioning_ParMETIS *)part->data;
20: PetscInt *locals = NULL;
21: Mat mat = part->adj, amat, pmat;
22: PetscBool flg;
23: PetscInt bs = 1;
25: PetscFunctionBegin;
27: PetscAssertPointer(partitioning, 4);
28: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIADJ, &flg));
29: if (flg) {
30: amat = mat;
31: PetscCall(PetscObjectReference((PetscObject)amat));
32: } else {
33: /* bs indicates if the converted matrix is "reduced" from the original and hence the
34: resulting partition results need to be stretched to match the original matrix */
35: PetscCall(MatConvert(mat, MATMPIADJ, MAT_INITIAL_MATRIX, &amat));
36: if (amat->rmap->n > 0) bs = mat->rmap->n / amat->rmap->n;
37: }
38: PetscCall(MatMPIAdjCreateNonemptySubcommMat(amat, &pmat));
39: PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)part)));
41: if (pmat) {
42: MPI_Comm pcomm, comm;
43: Mat_MPIAdj *adj = (Mat_MPIAdj *)pmat->data;
44: PetscInt *vtxdist = pmat->rmap->range;
45: PetscInt *xadj = adj->i;
46: PetscInt *adjncy = adj->j;
47: PetscInt *NDorder = NULL;
48: PetscInt itmp = 0, wgtflag = 0, numflag = 0, ncon = part->ncon, nparts = part->n, options[24], i, j;
49: real_t *tpwgts, *ubvec, itr = (real_t)0.1;
51: PetscCall(PetscObjectGetComm((PetscObject)pmat, &pcomm));
52: if (PetscDefined(USE_DEBUG)) {
53: /* check that matrix has no diagonal entries */
54: PetscInt rstart;
55: PetscCall(MatGetOwnershipRange(pmat, &rstart, NULL));
56: for (i = 0; i < pmat->rmap->n; i++) {
57: for (j = xadj[i]; j < xadj[i + 1]; j++) PetscCheck(adjncy[j] != i + rstart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Row %" PetscInt_FMT " has diagonal entry; ParMETIS forbids diagonal entry", i + rstart);
58: }
59: }
61: PetscCall(PetscMalloc1(pmat->rmap->n, &locals));
63: if (isImprove) {
64: PetscInt i;
65: const PetscInt *part_indices;
67: PetscCall(ISGetIndices(*partitioning, &part_indices));
68: for (i = 0; i < pmat->rmap->n; i++) locals[i] = part_indices[i * bs];
69: PetscCall(ISRestoreIndices(*partitioning, &part_indices));
70: PetscCall(ISDestroy(partitioning));
71: }
73: if (adj->values && part->use_edge_weights && !part->vertex_weights) wgtflag = 1;
74: if (part->vertex_weights && !adj->values) wgtflag = 2;
75: if (part->vertex_weights && adj->values && part->use_edge_weights) wgtflag = 3;
77: if (PetscLogPrintInfo) {
78: itmp = pmetis->printout;
79: pmetis->printout = 127;
80: }
81: PetscCall(PetscMalloc1(ncon * nparts, &tpwgts));
82: for (i = 0; i < ncon; i++) {
83: for (j = 0; j < nparts; j++) {
84: if (part->part_weights) {
85: tpwgts[i * nparts + j] = (real_t)part->part_weights[i * nparts + j];
86: } else {
87: tpwgts[i * nparts + j] = (real_t)1. / nparts;
88: }
89: }
90: }
91: PetscCall(PetscMalloc1(ncon, &ubvec));
92: for (i = 0; i < ncon; i++) ubvec[i] = (real_t)1.05;
93: /* This sets the defaults */
94: options[0] = 1;
95: for (i = 1; i < 24; i++) options[i] = -1;
96: options[1] = 0; /* no verbosity */
97: options[2] = 0;
98: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
99: /* Duplicate the communicator to be sure that ParMETIS attribute caching does not interfere with PETSc. */
100: PetscCallMPI(MPI_Comm_dup(pcomm, &comm));
101: if (useND) {
102: PetscInt *sizes, *seps, log2size, subd, *level;
103: PetscMPIInt size;
104: idx_t mtype = PARMETIS_MTYPE_GLOBAL, rtype = PARMETIS_SRTYPE_2PHASE, p_nseps = 1, s_nseps = 1;
105: real_t ubfrac = (real_t)1.05;
107: PetscCallMPI(MPI_Comm_size(comm, &size));
108: PetscCall(PetscMalloc1(pmat->rmap->n, &NDorder));
109: PetscCall(PetscMalloc3(2 * size, &sizes, 4 * size, &seps, size, &level));
110: PetscCallParMETIS(ParMETIS_V32_NodeND, (idx_t *)vtxdist, (idx_t *)xadj, (idx_t *)adjncy, (idx_t *)part->vertex_weights, (idx_t *)&numflag, &mtype, &rtype, &p_nseps, &s_nseps, &ubfrac, NULL /* seed */, NULL /* dbglvl */, (idx_t *)NDorder, (idx_t *)(sizes), &comm);
111: log2size = PetscLog2Real(size);
112: subd = PetscPowInt(2, log2size);
113: PetscCall(MatPartitioningSizesToSep_Private(subd, sizes, seps, level));
114: for (i = 0; i < pmat->rmap->n; i++) {
115: PetscInt loc;
117: PetscCall(PetscFindInt(NDorder[i], 2 * subd, seps, &loc));
118: if (loc < 0) {
119: loc = -(loc + 1);
120: if (loc % 2) { /* part of subdomain */
121: locals[i] = loc / 2;
122: } else {
123: PetscCall(PetscFindInt(NDorder[i], 2 * (subd - 1), seps + 2 * subd, &loc));
124: loc = loc < 0 ? -(loc + 1) / 2 : loc / 2;
125: locals[i] = level[loc];
126: }
127: } else locals[i] = loc / 2;
128: }
129: PetscCall(PetscFree3(sizes, seps, level));
130: } else {
131: if (pmetis->repartition) {
132: PetscCallParMETIS(ParMETIS_V3_AdaptiveRepart, (idx_t *)vtxdist, (idx_t *)xadj, (idx_t *)adjncy, (idx_t *)part->vertex_weights, (idx_t *)part->vertex_weights, (idx_t *)adj->values, (idx_t *)&wgtflag, (idx_t *)&numflag, (idx_t *)&ncon, (idx_t *)&nparts, tpwgts, ubvec, &itr, (idx_t *)options,
133: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm);
134: } else if (isImprove) {
135: PetscCallParMETIS(ParMETIS_V3_RefineKway, (idx_t *)vtxdist, (idx_t *)xadj, (idx_t *)adjncy, (idx_t *)part->vertex_weights, (idx_t *)adj->values, (idx_t *)&wgtflag, (idx_t *)&numflag, (idx_t *)&ncon, (idx_t *)&nparts, tpwgts, ubvec, (idx_t *)options,
136: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm);
137: } else {
138: PetscCallParMETIS(ParMETIS_V3_PartKway, (idx_t *)vtxdist, (idx_t *)xadj, (idx_t *)adjncy, (idx_t *)part->vertex_weights, (idx_t *)adj->values, (idx_t *)&wgtflag, (idx_t *)&numflag, (idx_t *)&ncon, (idx_t *)&nparts, tpwgts, ubvec, (idx_t *)options,
139: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm);
140: }
141: }
142: PetscCallMPI(MPI_Comm_free(&comm));
144: PetscCall(PetscFree(tpwgts));
145: PetscCall(PetscFree(ubvec));
146: if (PetscLogPrintInfo) pmetis->printout = itmp;
148: if (bs > 1) {
149: PetscInt i, j, *newlocals;
150: PetscCall(PetscMalloc1(bs * pmat->rmap->n, &newlocals));
151: for (i = 0; i < pmat->rmap->n; i++) {
152: for (j = 0; j < bs; j++) newlocals[bs * i + j] = locals[i];
153: }
154: PetscCall(PetscFree(locals));
155: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), bs * pmat->rmap->n, newlocals, PETSC_OWN_POINTER, partitioning));
156: } else {
157: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, locals, PETSC_OWN_POINTER, partitioning));
158: }
159: if (useND) {
160: IS ndis;
162: if (bs > 1) {
163: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)part), bs, pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis));
164: } else {
165: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis));
166: }
167: PetscCall(ISSetPermutation(ndis));
168: PetscCall(PetscObjectCompose((PetscObject)*partitioning, "_petsc_matpartitioning_ndorder", (PetscObject)ndis));
169: PetscCall(ISDestroy(&ndis));
170: }
171: } else {
172: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, partitioning));
173: if (useND) {
174: IS ndis;
176: if (bs > 1) {
177: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)part), bs, 0, NULL, PETSC_COPY_VALUES, &ndis));
178: } else {
179: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, &ndis));
180: }
181: PetscCall(ISSetPermutation(ndis));
182: PetscCall(PetscObjectCompose((PetscObject)*partitioning, "_petsc_matpartitioning_ndorder", (PetscObject)ndis));
183: PetscCall(ISDestroy(&ndis));
184: }
185: }
186: PetscCall(MatDestroy(&pmat));
187: PetscCall(MatDestroy(&amat));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*
192: Uses the ParMETIS parallel matrix partitioner to compute a nested dissection ordering of the matrix in parallel
193: */
194: static PetscErrorCode MatPartitioningApplyND_ParMETIS(MatPartitioning part, IS *partitioning)
195: {
196: PetscFunctionBegin;
197: PetscCall(MatPartitioningApply_ParMETIS_Private(part, PETSC_TRUE, PETSC_FALSE, partitioning));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*
202: Uses the ParMETIS parallel matrix partitioner to partition the matrix in parallel
203: */
204: static PetscErrorCode MatPartitioningApply_ParMETIS(MatPartitioning part, IS *partitioning)
205: {
206: PetscFunctionBegin;
207: PetscCall(MatPartitioningApply_ParMETIS_Private(part, PETSC_FALSE, PETSC_FALSE, partitioning));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*
212: Uses the ParMETIS to improve the quality of a partition
213: */
214: static PetscErrorCode MatPartitioningImprove_ParMETIS(MatPartitioning part, IS *partitioning)
215: {
216: PetscFunctionBegin;
217: PetscCall(MatPartitioningApply_ParMETIS_Private(part, PETSC_FALSE, PETSC_TRUE, partitioning));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode MatPartitioningView_ParMETIS(MatPartitioning part, PetscViewer viewer)
222: {
223: MatPartitioning_ParMETIS *pmetis = (MatPartitioning_ParMETIS *)part->data;
224: PetscMPIInt rank;
225: PetscBool isascii;
227: PetscFunctionBegin;
228: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)part), &rank));
229: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
230: if (isascii) {
231: if (pmetis->parallel == 2) {
232: PetscCall(PetscViewerASCIIPrintf(viewer, " Using parallel coarse grid partitioner\n"));
233: } else {
234: PetscCall(PetscViewerASCIIPrintf(viewer, " Using sequential coarse grid partitioner\n"));
235: }
236: PetscCall(PetscViewerASCIIPrintf(viewer, " Using %" PetscInt_FMT " fold factor\n", pmetis->foldfactor));
237: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
238: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d]Number of cuts found %" PetscInt_FMT "\n", rank, pmetis->cuts));
239: PetscCall(PetscViewerFlush(viewer));
240: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
241: }
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*@
246: MatPartitioningParMETISSetCoarseSequential - Use the sequential code to
247: do the partitioning of the coarse grid.
249: Logically Collective
251: Input Parameter:
252: . part - the partitioning context
254: Level: advanced
256: .seealso: `MATPARTITIONINGPARMETIS`
257: @*/
258: PetscErrorCode MatPartitioningParMETISSetCoarseSequential(MatPartitioning part)
259: {
260: MatPartitioning_ParMETIS *pmetis = (MatPartitioning_ParMETIS *)part->data;
262: PetscFunctionBegin;
263: pmetis->parallel = 1;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: /*@
268: MatPartitioningParMETISSetRepartition - Repartition
269: current mesh to rebalance computation.
271: Logically Collective
273: Input Parameter:
274: . part - the partitioning context
276: Level: advanced
278: .seealso: `MATPARTITIONINGPARMETIS`
279: @*/
280: PetscErrorCode MatPartitioningParMETISSetRepartition(MatPartitioning part)
281: {
282: MatPartitioning_ParMETIS *pmetis = (MatPartitioning_ParMETIS *)part->data;
284: PetscFunctionBegin;
285: pmetis->repartition = PETSC_TRUE;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: MatPartitioningParMETISGetEdgeCut - Returns the number of edge cuts in the vertex partition.
292: Input Parameter:
293: . part - the partitioning context
295: Output Parameter:
296: . cut - the edge cut
298: Level: advanced
300: .seealso: `MATPARTITIONINGPARMETIS`
301: @*/
302: PetscErrorCode MatPartitioningParMETISGetEdgeCut(MatPartitioning part, PetscInt *cut)
303: {
304: MatPartitioning_ParMETIS *pmetis = (MatPartitioning_ParMETIS *)part->data;
306: PetscFunctionBegin;
307: *cut = pmetis->cuts;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: static PetscErrorCode MatPartitioningSetFromOptions_ParMETIS(MatPartitioning part, PetscOptionItems PetscOptionsObject)
312: {
313: PetscBool flag = PETSC_FALSE;
315: PetscFunctionBegin;
316: PetscOptionsHeadBegin(PetscOptionsObject, "Set ParMETIS partitioning options");
317: PetscCall(PetscOptionsBool("-mat_partitioning_parmetis_coarse_sequential", "Use sequential coarse partitioner", "MatPartitioningParMETISSetCoarseSequential", flag, &flag, NULL));
318: if (flag) PetscCall(MatPartitioningParMETISSetCoarseSequential(part));
319: PetscCall(PetscOptionsBool("-mat_partitioning_parmetis_repartition", "", "MatPartitioningParMETISSetRepartition", flag, &flag, NULL));
320: if (flag) PetscCall(MatPartitioningParMETISSetRepartition(part));
321: PetscOptionsHeadEnd();
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode MatPartitioningDestroy_ParMETIS(MatPartitioning part)
326: {
327: MatPartitioning_ParMETIS *pmetis = (MatPartitioning_ParMETIS *)part->data;
329: PetscFunctionBegin;
330: PetscCall(PetscFree(pmetis));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /*MC
335: MATPARTITIONINGPARMETIS - Creates a partitioning context via the external package PARMETIS.
337: Collective
339: Input Parameter:
340: . part - the partitioning context
342: Options Database Key:
343: . -mat_partitioning_parmetis_coarse_sequential - use sequential ParMETIS coarse partitioner
345: Level: beginner
347: Note:
348: See https://www-users.cs.umn.edu/~karypis/metis/
350: .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MatPartitioningParMETISSetCoarseSequential()`, `MatPartitioningParMETISSetRepartition()`,
351: `MatPartitioningParMETISGetEdgeCut()`
352: M*/
354: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_ParMETIS(MatPartitioning part)
355: {
356: MatPartitioning_ParMETIS *pmetis;
358: PetscFunctionBegin;
359: PetscCall(PetscNew(&pmetis));
360: part->data = (void *)pmetis;
362: pmetis->cuts = 0; /* output variable */
363: pmetis->foldfactor = 150; /*folding factor */
364: pmetis->parallel = 2; /* use parallel partitioner for coarse grid */
365: pmetis->indexing = 0; /* index numbering starts from 0 */
366: pmetis->printout = 0; /* print no output while running */
367: pmetis->repartition = PETSC_FALSE;
369: part->ops->apply = MatPartitioningApply_ParMETIS;
370: part->ops->applynd = MatPartitioningApplyND_ParMETIS;
371: part->ops->improve = MatPartitioningImprove_ParMETIS;
372: part->ops->view = MatPartitioningView_ParMETIS;
373: part->ops->destroy = MatPartitioningDestroy_ParMETIS;
374: part->ops->setfromoptions = MatPartitioningSetFromOptions_ParMETIS;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*
379: Uses the ParMETIS package to convert a mesh to a cell graph via ParMETIS_V3_Mesh2Dual()
380: */
381: PETSC_EXTERN PetscErrorCode MatMeshToCellGraph_ParMETIS(Mat mesh, PetscInt ncommonnodes, Mat *dual)
382: {
383: PetscInt *newxadj, *newadjncy;
384: PetscInt numflag = 0;
385: Mat_MPIAdj *adj = (Mat_MPIAdj *)mesh->data, *newadj;
386: PetscBool flg;
387: MPI_Comm comm;
389: PetscFunctionBegin;
390: PetscCall(PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg));
391: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Must use MPIAdj matrix type");
393: PetscCall(PetscObjectGetComm((PetscObject)mesh, &comm));
394: PetscCallParMETIS(ParMETIS_V3_Mesh2Dual, (idx_t *)mesh->rmap->range, (idx_t *)adj->i, (idx_t *)adj->j, (idx_t *)&numflag, (idx_t *)&ncommonnodes, (idx_t **)&newxadj, (idx_t **)&newadjncy, &comm);
395: PetscCall(MatCreateMPIAdj(PetscObjectComm((PetscObject)mesh), mesh->rmap->n, mesh->rmap->N, newxadj, newadjncy, NULL, dual));
396: newadj = (Mat_MPIAdj *)(*dual)->data;
398: newadj->freeaijwithfree = PETSC_TRUE; /* signal the matrix should be freed with system free since space was allocated by ParMETIS */
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }