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 { \
27: } while (0)
29: #define PetscCallParmetis_(name, func, args) \
30: do { \
31: PetscStackPushExternal(name); \
32: int status = func args; \
33: PetscStackPop; \
34: 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;
49: PetscObjectTypeCompare((PetscObject)mat, MATMPIADJ, &flg);
50: if (flg) {
51: amat = mat;
52: PetscObjectReference((PetscObject)amat);
53: } else {
54: /* bs indicates if the converted matrix is "reduced" from the original and hence the
55: resulting partition results need to be stretched to match the original matrix */
56: MatConvert(mat, MATMPIADJ, MAT_INITIAL_MATRIX, &amat);
57: if (amat->rmap->n > 0) bs = mat->rmap->n / amat->rmap->n;
58: }
59: MatMPIAdjCreateNonemptySubcommMat(amat, &pmat);
60: MPI_Barrier(PetscObjectComm((PetscObject)part));
62: if (pmat) {
63: MPI_Comm pcomm, comm;
64: Mat_MPIAdj *adj = (Mat_MPIAdj *)pmat->data;
65: PetscInt *vtxdist = pmat->rmap->range;
66: PetscInt *xadj = adj->i;
67: PetscInt *adjncy = adj->j;
68: PetscInt *NDorder = NULL;
69: PetscInt itmp = 0, wgtflag = 0, numflag = 0, ncon = part->ncon, nparts = part->n, options[24], i, j;
70: real_t *tpwgts, *ubvec, itr = 0.1;
72: PetscObjectGetComm((PetscObject)pmat, &pcomm);
73: if (PetscDefined(USE_DEBUG)) {
74: /* check that matrix has no diagonal entries */
75: PetscInt rstart;
76: MatGetOwnershipRange(pmat, &rstart, NULL);
77: for (i = 0; i < pmat->rmap->n; i++) {
79: }
80: }
82: PetscMalloc1(pmat->rmap->n, &locals);
84: if (isImprove) {
85: PetscInt i;
86: const PetscInt *part_indices;
88: ISGetIndices(*partitioning, &part_indices);
89: for (i = 0; i < pmat->rmap->n; i++) locals[i] = part_indices[i * bs];
90: ISRestoreIndices(*partitioning, &part_indices);
91: ISDestroy(partitioning);
92: }
94: if (adj->values && part->use_edge_weights && !part->vertex_weights) wgtflag = 1;
95: if (part->vertex_weights && !adj->values) wgtflag = 2;
96: if (part->vertex_weights && adj->values && part->use_edge_weights) wgtflag = 3;
98: if (PetscLogPrintInfo) {
99: itmp = pmetis->printout;
100: pmetis->printout = 127;
101: }
102: PetscMalloc1(ncon * nparts, &tpwgts);
103: for (i = 0; i < ncon; i++) {
104: for (j = 0; j < nparts; j++) {
105: if (part->part_weights) {
106: tpwgts[i * nparts + j] = part->part_weights[i * nparts + j];
107: } else {
108: tpwgts[i * nparts + j] = 1. / nparts;
109: }
110: }
111: }
112: PetscMalloc1(ncon, &ubvec);
113: for (i = 0; i < ncon; i++) ubvec[i] = 1.05;
114: /* This sets the defaults */
115: options[0] = 1;
116: for (i = 1; i < 24; i++) options[i] = -1;
117: options[1] = 0; /* no verbosity */
118: options[2] = 0;
119: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
120: /* Duplicate the communicator to be sure that ParMETIS attribute caching does not interfere with PETSc. */
121: MPI_Comm_dup(pcomm, &comm);
122: if (useND) {
123: PetscInt *sizes, *seps, log2size, subd, *level;
124: PetscMPIInt size;
125: idx_t mtype = PARMETIS_MTYPE_GLOBAL, rtype = PARMETIS_SRTYPE_2PHASE, p_nseps = 1, s_nseps = 1;
126: real_t ubfrac = 1.05;
128: MPI_Comm_size(comm, &size);
129: PetscMalloc1(pmat->rmap->n, &NDorder);
130: PetscMalloc3(2 * size, &sizes, 4 * size, &seps, size, &level);
131: 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));
132: log2size = PetscLog2Real(size);
133: subd = PetscPowInt(2, log2size);
134: MatPartitioningSizesToSep_Private(subd, sizes, seps, level);
135: for (i = 0; i < pmat->rmap->n; i++) {
136: PetscInt loc;
138: PetscFindInt(NDorder[i], 2 * subd, seps, &loc);
139: if (loc < 0) {
140: loc = -(loc + 1);
141: if (loc % 2) { /* part of subdomain */
142: locals[i] = loc / 2;
143: } else {
144: PetscFindInt(NDorder[i], 2 * (subd - 1), seps + 2 * subd, &loc);
145: loc = loc < 0 ? -(loc + 1) / 2 : loc / 2;
146: locals[i] = level[loc];
147: }
148: } else locals[i] = loc / 2;
149: }
150: PetscFree3(sizes, seps, level);
151: } else {
152: if (pmetis->repartition) {
153: 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,
154: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm));
155: } else if (isImprove) {
156: 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,
157: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm));
158: } else {
159: 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,
160: (idx_t *)&pmetis->cuts, (idx_t *)locals, &comm));
161: }
162: }
163: MPI_Comm_free(&comm);
165: PetscFree(tpwgts);
166: PetscFree(ubvec);
167: if (PetscLogPrintInfo) pmetis->printout = itmp;
169: if (bs > 1) {
170: PetscInt i, j, *newlocals;
171: PetscMalloc1(bs * pmat->rmap->n, &newlocals);
172: for (i = 0; i < pmat->rmap->n; i++) {
173: for (j = 0; j < bs; j++) newlocals[bs * i + j] = locals[i];
174: }
175: PetscFree(locals);
176: ISCreateGeneral(PetscObjectComm((PetscObject)part), bs * pmat->rmap->n, newlocals, PETSC_OWN_POINTER, partitioning);
177: } else {
178: ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, locals, PETSC_OWN_POINTER, partitioning);
179: }
180: if (useND) {
181: IS ndis;
183: if (bs > 1) {
184: ISCreateBlock(PetscObjectComm((PetscObject)part), bs, pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis);
185: } else {
186: ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis);
187: }
188: ISSetPermutation(ndis);
189: PetscObjectCompose((PetscObject)(*partitioning), "_petsc_matpartitioning_ndorder", (PetscObject)ndis);
190: ISDestroy(&ndis);
191: }
192: } else {
193: ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, partitioning);
194: if (useND) {
195: IS ndis;
197: if (bs > 1) {
198: ISCreateBlock(PetscObjectComm((PetscObject)part), bs, 0, NULL, PETSC_COPY_VALUES, &ndis);
199: } else {
200: ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, &ndis);
201: }
202: ISSetPermutation(ndis);
203: PetscObjectCompose((PetscObject)(*partitioning), "_petsc_matpartitioning_ndorder", (PetscObject)ndis);
204: ISDestroy(&ndis);
205: }
206: }
207: MatDestroy(&pmat);
208: MatDestroy(&amat);
209: return 0;
210: }
212: /*
213: Uses the ParMETIS parallel matrix partitioner to compute a nested dissection ordering of the matrix in parallel
214: */
215: static PetscErrorCode MatPartitioningApplyND_Parmetis(MatPartitioning part, IS *partitioning)
216: {
217: MatPartitioningApply_Parmetis_Private(part, PETSC_TRUE, PETSC_FALSE, partitioning);
218: return 0;
219: }
221: /*
222: Uses the ParMETIS parallel matrix partitioner to partition the matrix in parallel
223: */
224: static PetscErrorCode MatPartitioningApply_Parmetis(MatPartitioning part, IS *partitioning)
225: {
226: MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_FALSE, partitioning);
227: return 0;
228: }
230: /*
231: Uses the ParMETIS to improve the quality of a partition
232: */
233: static PetscErrorCode MatPartitioningImprove_Parmetis(MatPartitioning part, IS *partitioning)
234: {
235: MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_TRUE, partitioning);
236: return 0;
237: }
239: PetscErrorCode MatPartitioningView_Parmetis(MatPartitioning part, PetscViewer viewer)
240: {
241: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
242: PetscMPIInt rank;
243: PetscBool iascii;
245: MPI_Comm_rank(PetscObjectComm((PetscObject)part), &rank);
246: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
247: if (iascii) {
248: if (pmetis->parallel == 2) {
249: PetscViewerASCIIPrintf(viewer, " Using parallel coarse grid partitioner\n");
250: } else {
251: PetscViewerASCIIPrintf(viewer, " Using sequential coarse grid partitioner\n");
252: }
253: PetscViewerASCIIPrintf(viewer, " Using %" PetscInt_FMT " fold factor\n", pmetis->foldfactor);
254: PetscViewerASCIIPushSynchronized(viewer);
255: PetscViewerASCIISynchronizedPrintf(viewer, " [%d]Number of cuts found %" PetscInt_FMT "\n", rank, pmetis->cuts);
256: PetscViewerFlush(viewer);
257: PetscViewerASCIIPopSynchronized(viewer);
258: }
259: return 0;
260: }
262: /*@
263: MatPartitioningParmetisSetCoarseSequential - Use the sequential code to
264: do the partitioning of the coarse grid.
266: Logically Collective
268: Input Parameter:
269: . part - the partitioning context
271: Level: advanced
273: .seealso: `MATPARTITIONINGPARMETIS`
274: @*/
275: PetscErrorCode MatPartitioningParmetisSetCoarseSequential(MatPartitioning part)
276: {
277: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
279: pmetis->parallel = 1;
280: return 0;
281: }
283: /*@
284: MatPartitioningParmetisSetRepartition - Repartition
285: current mesh to rebalance computation.
287: Logically Collective
289: Input Parameter:
290: . part - the partitioning context
292: Level: advanced
294: .seealso: `MATPARTITIONINGPARMETIS`
295: @*/
296: PetscErrorCode MatPartitioningParmetisSetRepartition(MatPartitioning part)
297: {
298: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
300: pmetis->repartition = PETSC_TRUE;
301: return 0;
302: }
304: /*@
305: MatPartitioningParmetisGetEdgeCut - Returns the number of edge cuts in the vertex partition.
307: Input Parameter:
308: . part - the partitioning context
310: Output Parameter:
311: . cut - the edge cut
313: Level: advanced
315: .seealso: `MATPARTITIONINGPARMETIS`
316: @*/
317: PetscErrorCode MatPartitioningParmetisGetEdgeCut(MatPartitioning part, PetscInt *cut)
318: {
319: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
321: *cut = pmetis->cuts;
322: return 0;
323: }
325: PetscErrorCode MatPartitioningSetFromOptions_Parmetis(MatPartitioning part, PetscOptionItems *PetscOptionsObject)
326: {
327: PetscBool flag = PETSC_FALSE;
329: PetscOptionsHeadBegin(PetscOptionsObject, "Set ParMeTiS partitioning options");
330: PetscOptionsBool("-mat_partitioning_parmetis_coarse_sequential", "Use sequential coarse partitioner", "MatPartitioningParmetisSetCoarseSequential", flag, &flag, NULL);
331: if (flag) MatPartitioningParmetisSetCoarseSequential(part);
332: PetscOptionsBool("-mat_partitioning_parmetis_repartition", "", "MatPartitioningParmetisSetRepartition", flag, &flag, NULL);
333: if (flag) MatPartitioningParmetisSetRepartition(part);
334: PetscOptionsHeadEnd();
335: return 0;
336: }
338: PetscErrorCode MatPartitioningDestroy_Parmetis(MatPartitioning part)
339: {
340: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
342: PetscFree(pmetis);
343: return 0;
344: }
346: /*MC
347: MATPARTITIONINGPARMETIS - Creates a partitioning context via the external package PARMETIS.
349: Collective
351: Input Parameter:
352: . part - the partitioning context
354: Options Database Keys:
355: . -mat_partitioning_parmetis_coarse_sequential - use sequential PARMETIS coarse partitioner
357: Level: beginner
359: Note:
360: See https://www-users.cs.umn.edu/~karypis/metis/
362: .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MatPartitioningParmetisSetCoarseSequential()`, `MatPartitioningParmetisSetRepartition()`,
363: `MatPartitioningParmetisGetEdgeCut()`
364: M*/
366: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Parmetis(MatPartitioning part)
367: {
368: MatPartitioning_Parmetis *pmetis;
370: PetscNew(&pmetis);
371: part->data = (void *)pmetis;
373: pmetis->cuts = 0; /* output variable */
374: pmetis->foldfactor = 150; /*folding factor */
375: pmetis->parallel = 2; /* use parallel partitioner for coarse grid */
376: pmetis->indexing = 0; /* index numbering starts from 0 */
377: pmetis->printout = 0; /* print no output while running */
378: pmetis->repartition = PETSC_FALSE;
380: part->ops->apply = MatPartitioningApply_Parmetis;
381: part->ops->applynd = MatPartitioningApplyND_Parmetis;
382: part->ops->improve = MatPartitioningImprove_Parmetis;
383: part->ops->view = MatPartitioningView_Parmetis;
384: part->ops->destroy = MatPartitioningDestroy_Parmetis;
385: part->ops->setfromoptions = MatPartitioningSetFromOptions_Parmetis;
386: return 0;
387: }
389: /*@
390: 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
391: between cells (the "dual" graph) and is suitable for partitioning with the `MatPartitioning object`. Use this to partition
392: cells of a mesh.
394: Collective
396: Input Parameters:
397: + mesh - the graph that represents the coupling of the vertices of the mesh
398: - ncommonnodes - mesh elements that share this number of common nodes are considered neighbors, use 2 for triangles and
399: quadrilaterials, 3 for tetrahedrals and 4 for hexahedrals
401: Output Parameter:
402: . dual - the dual graph
404: Notes:
405: Currently requires ParMetis to be installed and uses ParMETIS_V3_Mesh2Dual()
407: $ Each row of the mesh object represents a single cell in the mesh. For triangles it has 3 entries, quadrilaterials 4 entries,
408: $ tetrahedrals 4 entries and hexahedrals 8 entries. You can mix triangles and quadrilaterals in the same mesh, but cannot
409: $ mix tetrahedrals and hexahedrals
410: $ The columns of each row of the Mat mesh are the global vertex numbers of the vertices of that row's cell.
411: $ The number of rows in mesh is number of cells, the number of columns is the number of vertices.
413: Level: advanced
415: .seealso: `MatCreateMPIAdj()`, `MatPartitioningCreate()`
416: @*/
417: PetscErrorCode MatMeshToCellGraph(Mat mesh, PetscInt ncommonnodes, Mat *dual)
418: {
419: PetscInt *newxadj, *newadjncy;
420: PetscInt numflag = 0;
421: Mat_MPIAdj *adj = (Mat_MPIAdj *)mesh->data, *newadj;
422: PetscBool flg;
423: MPI_Comm comm;
425: PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg);
428: PetscObjectGetComm((PetscObject)mesh, &comm);
429: 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));
430: MatCreateMPIAdj(PetscObjectComm((PetscObject)mesh), mesh->rmap->n, mesh->rmap->N, newxadj, newadjncy, NULL, dual);
431: newadj = (Mat_MPIAdj *)(*dual)->data;
433: newadj->freeaijwithfree = PETSC_TRUE; /* signal the matrix should be freed with system free since space was allocated by ParMETIS */
434: return 0;
435: }