Actual source code: pmetis.c
2: #include <../src/mat/impls/adj/mpi/mpiadj.h>
4: /*
5: Currently using ParMetis-4.0.2
6: */
8: #include <parmetis.h>
10: /*
11: The first 5 elements of this structure are the input control array to Metis
12: */
13: typedef struct {
14: PetscInt cuts; /* number of cuts made (output) */
15: PetscInt foldfactor;
16: PetscInt parallel; /* use parallel partitioner for coarse problem */
17: PetscInt indexing; /* 0 indicates C indexing, 1 Fortran */
18: PetscInt printout; /* indicates if one wishes Metis to print info */
19: PetscBool repartition;
20: } MatPartitioning_Parmetis;
22: #define PetscCallPARMETIS(n, func) \
23: do { \
24: PetscCheck(n != METIS_ERROR_INPUT, PETSC_COMM_SELF, PETSC_ERR_LIB, "ParMETIS error due to wrong inputs and/or options for %s", func); \
25: PetscCheck(n != METIS_ERROR_MEMORY, PETSC_COMM_SELF, PETSC_ERR_LIB, "ParMETIS error due to insufficient memory in %s", func); \
26: PetscCheck(n != METIS_ERROR, PETSC_COMM_SELF, PETSC_ERR_LIB, "ParMETIS general error in %s", func); \
27: } while (0)
29: #define PetscCallParmetis_(name, func, args) \
30: do { \
31: PetscStackPushExternal(name); \
32: int status = func args; \
33: PetscStackPop; \
34: PetscCallPARMETIS(status, name); \
35: } while (0)
37: #define PetscCallParmetis(func, args) PetscCallParmetis_(PetscStringize(func), func, args)
39: static PetscErrorCode MatPartitioningApply_Parmetis_Private(MatPartitioning part, PetscBool useND, PetscBool isImprove, IS *partitioning)
40: {
41: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
42: PetscInt *locals = NULL;
43: Mat mat = part->adj, amat, pmat;
44: PetscBool flg;
45: PetscInt bs = 1;
47: PetscFunctionBegin;
50: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIADJ, &flg));
51: if (flg) {
52: amat = mat;
53: PetscCall(PetscObjectReference((PetscObject)amat));
54: } else {
55: /* bs indicates if the converted matrix is "reduced" from the original and hence the
56: resulting partition results need to be stretched to match the original matrix */
57: PetscCall(MatConvert(mat, MATMPIADJ, MAT_INITIAL_MATRIX, &amat));
58: if (amat->rmap->n > 0) bs = mat->rmap->n / amat->rmap->n;
59: }
60: PetscCall(MatMPIAdjCreateNonemptySubcommMat(amat, &pmat));
61: PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)part)));
63: if (pmat) {
64: MPI_Comm pcomm, comm;
65: Mat_MPIAdj *adj = (Mat_MPIAdj *)pmat->data;
66: PetscInt *vtxdist = pmat->rmap->range;
67: PetscInt *xadj = adj->i;
68: PetscInt *adjncy = adj->j;
69: PetscInt *NDorder = NULL;
70: PetscInt itmp = 0, wgtflag = 0, numflag = 0, ncon = part->ncon, nparts = part->n, options[24], i, j;
71: real_t *tpwgts, *ubvec, itr = 0.1;
73: PetscCall(PetscObjectGetComm((PetscObject)pmat, &pcomm));
74: if (PetscDefined(USE_DEBUG)) {
75: /* check that matrix has no diagonal entries */
76: PetscInt rstart;
77: PetscCall(MatGetOwnershipRange(pmat, &rstart, NULL));
78: for (i = 0; i < pmat->rmap->n; i++) {
79: 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);
80: }
81: }
83: PetscCall(PetscMalloc1(pmat->rmap->n, &locals));
85: if (isImprove) {
86: PetscInt i;
87: const PetscInt *part_indices;
89: PetscCall(ISGetIndices(*partitioning, &part_indices));
90: for (i = 0; i < pmat->rmap->n; i++) locals[i] = part_indices[i * bs];
91: PetscCall(ISRestoreIndices(*partitioning, &part_indices));
92: PetscCall(ISDestroy(partitioning));
93: }
95: if (adj->values && part->use_edge_weights && !part->vertex_weights) wgtflag = 1;
96: if (part->vertex_weights && !adj->values) wgtflag = 2;
97: if (part->vertex_weights && adj->values && part->use_edge_weights) wgtflag = 3;
99: if (PetscLogPrintInfo) {
100: itmp = pmetis->printout;
101: pmetis->printout = 127;
102: }
103: PetscCall(PetscMalloc1(ncon * nparts, &tpwgts));
104: for (i = 0; i < ncon; i++) {
105: for (j = 0; j < nparts; j++) {
106: if (part->part_weights) {
107: tpwgts[i * nparts + j] = part->part_weights[i * nparts + j];
108: } else {
109: tpwgts[i * nparts + j] = 1. / nparts;
110: }
111: }
112: }
113: PetscCall(PetscMalloc1(ncon, &ubvec));
114: for (i = 0; i < ncon; i++) ubvec[i] = 1.05;
115: /* This sets the defaults */
116: options[0] = 1;
117: for (i = 1; i < 24; i++) options[i] = -1;
118: options[1] = 0; /* no verbosity */
119: options[2] = 0;
120: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
121: /* Duplicate the communicator to be sure that ParMETIS attribute caching does not interfere with PETSc. */
122: PetscCallMPI(MPI_Comm_dup(pcomm, &comm));
123: if (useND) {
124: PetscInt *sizes, *seps, log2size, subd, *level;
125: PetscMPIInt size;
126: idx_t mtype = PARMETIS_MTYPE_GLOBAL, rtype = PARMETIS_SRTYPE_2PHASE, p_nseps = 1, s_nseps = 1;
127: real_t ubfrac = 1.05;
129: PetscCallMPI(MPI_Comm_size(comm, &size));
130: PetscCall(PetscMalloc1(pmat->rmap->n, &NDorder));
131: PetscCall(PetscMalloc3(2 * size, &sizes, 4 * size, &seps, size, &level));
132: 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));
133: log2size = PetscLog2Real(size);
134: subd = PetscPowInt(2, log2size);
135: PetscCall(MatPartitioningSizesToSep_Private(subd, sizes, seps, level));
136: for (i = 0; i < pmat->rmap->n; i++) {
137: PetscInt loc;
139: PetscCall(PetscFindInt(NDorder[i], 2 * subd, seps, &loc));
140: if (loc < 0) {
141: loc = -(loc + 1);
142: if (loc % 2) { /* part of subdomain */
143: locals[i] = loc / 2;
144: } else {
145: PetscCall(PetscFindInt(NDorder[i], 2 * (subd - 1), seps + 2 * subd, &loc));
146: loc = loc < 0 ? -(loc + 1) / 2 : loc / 2;
147: locals[i] = level[loc];
148: }
149: } else locals[i] = loc / 2;
150: }
151: PetscCall(PetscFree3(sizes, seps, level));
152: } else {
153: if (pmetis->repartition) {
154: 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,
155: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm));
156: } else if (isImprove) {
157: 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,
158: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm));
159: } else {
160: 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,
161: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm));
162: }
163: }
164: PetscCallMPI(MPI_Comm_free(&comm));
166: PetscCall(PetscFree(tpwgts));
167: PetscCall(PetscFree(ubvec));
168: if (PetscLogPrintInfo) pmetis->printout = itmp;
170: if (bs > 1) {
171: PetscInt i, j, *newlocals;
172: PetscCall(PetscMalloc1(bs * pmat->rmap->n, &newlocals));
173: for (i = 0; i < pmat->rmap->n; i++) {
174: for (j = 0; j < bs; j++) newlocals[bs * i + j] = locals[i];
175: }
176: PetscCall(PetscFree(locals));
177: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), bs * pmat->rmap->n, newlocals, PETSC_OWN_POINTER, partitioning));
178: } else {
179: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, locals, PETSC_OWN_POINTER, partitioning));
180: }
181: if (useND) {
182: IS ndis;
184: if (bs > 1) {
185: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)part), bs, pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis));
186: } else {
187: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis));
188: }
189: PetscCall(ISSetPermutation(ndis));
190: PetscCall(PetscObjectCompose((PetscObject)(*partitioning), "_petsc_matpartitioning_ndorder", (PetscObject)ndis));
191: PetscCall(ISDestroy(&ndis));
192: }
193: } else {
194: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, partitioning));
195: if (useND) {
196: IS ndis;
198: if (bs > 1) {
199: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)part), bs, 0, NULL, PETSC_COPY_VALUES, &ndis));
200: } else {
201: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, &ndis));
202: }
203: PetscCall(ISSetPermutation(ndis));
204: PetscCall(PetscObjectCompose((PetscObject)(*partitioning), "_petsc_matpartitioning_ndorder", (PetscObject)ndis));
205: PetscCall(ISDestroy(&ndis));
206: }
207: }
208: PetscCall(MatDestroy(&pmat));
209: PetscCall(MatDestroy(&amat));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*
214: Uses the ParMETIS parallel matrix partitioner to compute a nested dissection ordering of the matrix in parallel
215: */
216: static PetscErrorCode MatPartitioningApplyND_Parmetis(MatPartitioning part, IS *partitioning)
217: {
218: PetscFunctionBegin;
219: PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_TRUE, PETSC_FALSE, partitioning));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*
224: Uses the ParMETIS parallel matrix partitioner to partition the matrix in parallel
225: */
226: static PetscErrorCode MatPartitioningApply_Parmetis(MatPartitioning part, IS *partitioning)
227: {
228: PetscFunctionBegin;
229: PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_FALSE, partitioning));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*
234: Uses the ParMETIS to improve the quality of a partition
235: */
236: static PetscErrorCode MatPartitioningImprove_Parmetis(MatPartitioning part, IS *partitioning)
237: {
238: PetscFunctionBegin;
239: PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_TRUE, partitioning));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode MatPartitioningView_Parmetis(MatPartitioning part, PetscViewer viewer)
244: {
245: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
246: PetscMPIInt rank;
247: PetscBool iascii;
249: PetscFunctionBegin;
250: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)part), &rank));
251: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
252: if (iascii) {
253: if (pmetis->parallel == 2) {
254: PetscCall(PetscViewerASCIIPrintf(viewer, " Using parallel coarse grid partitioner\n"));
255: } else {
256: PetscCall(PetscViewerASCIIPrintf(viewer, " Using sequential coarse grid partitioner\n"));
257: }
258: PetscCall(PetscViewerASCIIPrintf(viewer, " Using %" PetscInt_FMT " fold factor\n", pmetis->foldfactor));
259: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
260: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d]Number of cuts found %" PetscInt_FMT "\n", rank, pmetis->cuts));
261: PetscCall(PetscViewerFlush(viewer));
262: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
263: }
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: /*@
268: MatPartitioningParmetisSetCoarseSequential - Use the sequential code to
269: do the partitioning of the coarse grid.
271: Logically Collective
273: Input Parameter:
274: . part - the partitioning context
276: Level: advanced
278: .seealso: `MATPARTITIONINGPARMETIS`
279: @*/
280: PetscErrorCode MatPartitioningParmetisSetCoarseSequential(MatPartitioning part)
281: {
282: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
284: PetscFunctionBegin;
285: pmetis->parallel = 1;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: MatPartitioningParmetisSetRepartition - Repartition
291: current mesh to rebalance computation.
293: Logically Collective
295: Input Parameter:
296: . part - the partitioning context
298: Level: advanced
300: .seealso: `MATPARTITIONINGPARMETIS`
301: @*/
302: PetscErrorCode MatPartitioningParmetisSetRepartition(MatPartitioning part)
303: {
304: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
306: PetscFunctionBegin;
307: pmetis->repartition = PETSC_TRUE;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: MatPartitioningParmetisGetEdgeCut - Returns the number of edge cuts in the vertex partition.
314: Input Parameter:
315: . part - the partitioning context
317: Output Parameter:
318: . cut - the edge cut
320: Level: advanced
322: .seealso: `MATPARTITIONINGPARMETIS`
323: @*/
324: PetscErrorCode MatPartitioningParmetisGetEdgeCut(MatPartitioning part, PetscInt *cut)
325: {
326: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
328: PetscFunctionBegin;
329: *cut = pmetis->cuts;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: PetscErrorCode MatPartitioningSetFromOptions_Parmetis(MatPartitioning part, PetscOptionItems *PetscOptionsObject)
334: {
335: PetscBool flag = PETSC_FALSE;
337: PetscFunctionBegin;
338: PetscOptionsHeadBegin(PetscOptionsObject, "Set ParMeTiS partitioning options");
339: PetscCall(PetscOptionsBool("-mat_partitioning_parmetis_coarse_sequential", "Use sequential coarse partitioner", "MatPartitioningParmetisSetCoarseSequential", flag, &flag, NULL));
340: if (flag) PetscCall(MatPartitioningParmetisSetCoarseSequential(part));
341: PetscCall(PetscOptionsBool("-mat_partitioning_parmetis_repartition", "", "MatPartitioningParmetisSetRepartition", flag, &flag, NULL));
342: if (flag) PetscCall(MatPartitioningParmetisSetRepartition(part));
343: PetscOptionsHeadEnd();
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: PetscErrorCode MatPartitioningDestroy_Parmetis(MatPartitioning part)
348: {
349: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
351: PetscFunctionBegin;
352: PetscCall(PetscFree(pmetis));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*MC
357: MATPARTITIONINGPARMETIS - Creates a partitioning context via the external package PARMETIS.
359: Collective
361: Input Parameter:
362: . part - the partitioning context
364: Options Database Key:
365: . -mat_partitioning_parmetis_coarse_sequential - use sequential PARMETIS coarse partitioner
367: Level: beginner
369: Note:
370: See https://www-users.cs.umn.edu/~karypis/metis/
372: .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MatPartitioningParmetisSetCoarseSequential()`, `MatPartitioningParmetisSetRepartition()`,
373: `MatPartitioningParmetisGetEdgeCut()`
374: M*/
376: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Parmetis(MatPartitioning part)
377: {
378: MatPartitioning_Parmetis *pmetis;
380: PetscFunctionBegin;
381: PetscCall(PetscNew(&pmetis));
382: part->data = (void *)pmetis;
384: pmetis->cuts = 0; /* output variable */
385: pmetis->foldfactor = 150; /*folding factor */
386: pmetis->parallel = 2; /* use parallel partitioner for coarse grid */
387: pmetis->indexing = 0; /* index numbering starts from 0 */
388: pmetis->printout = 0; /* print no output while running */
389: pmetis->repartition = PETSC_FALSE;
391: part->ops->apply = MatPartitioningApply_Parmetis;
392: part->ops->applynd = MatPartitioningApplyND_Parmetis;
393: part->ops->improve = MatPartitioningImprove_Parmetis;
394: part->ops->view = MatPartitioningView_Parmetis;
395: part->ops->destroy = MatPartitioningDestroy_Parmetis;
396: part->ops->setfromoptions = MatPartitioningSetFromOptions_Parmetis;
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*@
401: MatMeshToCellGraph - Uses the ParMETIS package to convert a `Mat` that represents coupling of vertices of a mesh to a `Mat` the represents the graph of the coupling
402: between cells (the "dual" graph) and is suitable for partitioning with the `MatPartitioning` object. Use this to partition
403: cells of a mesh.
405: Collective
407: Input Parameters:
408: + mesh - the graph that represents the coupling of the vertices of the mesh
409: - ncommonnodes - mesh elements that share this number of common nodes are considered neighbors, use 2 for triangles and
410: quadrilaterials, 3 for tetrahedrals and 4 for hexahedrals
412: Output Parameter:
413: . dual - the dual graph
415: Level: advanced
417: Notes:
418: Currently requires ParMetis to be installed and uses ParMETIS_V3_Mesh2Dual()
420: Each row of the mesh object represents a single cell in the mesh. For triangles it has 3 entries, quadrilaterials 4 entries,
421: tetrahedrals 4 entries and hexahedrals 8 entries. You can mix triangles and quadrilaterals in the same mesh, but cannot
422: mix tetrahedrals and hexahedrals
423: The columns of each row of the `Mat` mesh are the global vertex numbers of the vertices of that row's cell.
424: The number of rows in mesh is number of cells, the number of columns is the number of vertices.
426: .seealso: `MatCreateMPIAdj()`, `MatPartitioningCreate()`
427: @*/
428: PetscErrorCode MatMeshToCellGraph(Mat mesh, PetscInt ncommonnodes, Mat *dual)
429: {
430: PetscInt *newxadj, *newadjncy;
431: PetscInt numflag = 0;
432: Mat_MPIAdj *adj = (Mat_MPIAdj *)mesh->data, *newadj;
433: PetscBool flg;
434: MPI_Comm comm;
436: PetscFunctionBegin;
437: PetscCall(PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg));
438: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Must use MPIAdj matrix type");
440: PetscCall(PetscObjectGetComm((PetscObject)mesh, &comm));
441: 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));
442: PetscCall(MatCreateMPIAdj(PetscObjectComm((PetscObject)mesh), mesh->rmap->n, mesh->rmap->N, newxadj, newadjncy, NULL, dual));
443: newadj = (Mat_MPIAdj *)(*dual)->data;
445: newadj->freeaijwithfree = PETSC_TRUE; /* signal the matrix should be freed with system free since space was allocated by ParMETIS */
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }