1 #include <../src/mat/impls/adj/mpi/mpiadj.h> /*I "petscmat.h" I*/
2
3 /*
4 Currently using ParMetis-4.0.2
5 */
6
7 #include <parmetis.h>
8
9 /*
10 The first 5 elements of this structure are the input control array to Metis
11 */
12 typedef struct {
13 PetscInt cuts; /* number of cuts made (output) */
14 PetscInt foldfactor;
15 PetscInt parallel; /* use parallel partitioner for coarse problem */
16 PetscInt indexing; /* 0 indicates C indexing, 1 Fortran */
17 PetscInt printout; /* indicates if one wishes Metis to print info */
18 PetscBool repartition;
19 } MatPartitioning_Parmetis;
20
21 #define PetscCallPARMETIS(n, func) \
22 do { \
23 PetscCheck(n != METIS_ERROR_INPUT, PETSC_COMM_SELF, PETSC_ERR_LIB, "ParMETIS error due to wrong inputs and/or options for %s", func); \
24 PetscCheck(n != METIS_ERROR_MEMORY, PETSC_COMM_SELF, PETSC_ERR_LIB, "ParMETIS error due to insufficient memory in %s", func); \
25 PetscCheck(n != METIS_ERROR, PETSC_COMM_SELF, PETSC_ERR_LIB, "ParMETIS general error in %s", func); \
26 } while (0)
27
28 #define PetscCallParmetis_(name, func, args) \
29 do { \
30 PetscStackPushExternal(name); \
31 int status = func args; \
32 PetscStackPop; \
33 PetscCallPARMETIS(status, name); \
34 } while (0)
35
36 #define PetscCallParmetis(func, args) PetscCallParmetis_(PetscStringize(func), func, args)
37
MatPartitioningApply_Parmetis_Private(MatPartitioning part,PetscBool useND,PetscBool isImprove,IS * partitioning)38 static PetscErrorCode MatPartitioningApply_Parmetis_Private(MatPartitioning part, PetscBool useND, PetscBool isImprove, IS *partitioning)
39 {
40 MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
41 PetscInt *locals = NULL;
42 Mat mat = part->adj, amat, pmat;
43 PetscBool flg;
44 PetscInt bs = 1;
45
46 PetscFunctionBegin;
47 PetscValidHeaderSpecific(part, MAT_PARTITIONING_CLASSID, 1);
48 PetscAssertPointer(partitioning, 4);
49 PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIADJ, &flg));
50 if (flg) {
51 amat = mat;
52 PetscCall(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 PetscCall(MatConvert(mat, MATMPIADJ, MAT_INITIAL_MATRIX, &amat));
57 if (amat->rmap->n > 0) bs = mat->rmap->n / amat->rmap->n;
58 }
59 PetscCall(MatMPIAdjCreateNonemptySubcommMat(amat, &pmat));
60 PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)part)));
61
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 = (real_t)0.1;
71
72 PetscCall(PetscObjectGetComm((PetscObject)pmat, &pcomm));
73 if (PetscDefined(USE_DEBUG)) {
74 /* check that matrix has no diagonal entries */
75 PetscInt rstart;
76 PetscCall(MatGetOwnershipRange(pmat, &rstart, NULL));
77 for (i = 0; i < pmat->rmap->n; i++) {
78 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);
79 }
80 }
81
82 PetscCall(PetscMalloc1(pmat->rmap->n, &locals));
83
84 if (isImprove) {
85 PetscInt i;
86 const PetscInt *part_indices;
87 PetscValidHeaderSpecific(*partitioning, IS_CLASSID, 4);
88 PetscCall(ISGetIndices(*partitioning, &part_indices));
89 for (i = 0; i < pmat->rmap->n; i++) locals[i] = part_indices[i * bs];
90 PetscCall(ISRestoreIndices(*partitioning, &part_indices));
91 PetscCall(ISDestroy(partitioning));
92 }
93
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;
97
98 if (PetscLogPrintInfo) {
99 itmp = pmetis->printout;
100 pmetis->printout = 127;
101 }
102 PetscCall(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] = (real_t)part->part_weights[i * nparts + j];
107 } else {
108 tpwgts[i * nparts + j] = (real_t)1. / nparts;
109 }
110 }
111 }
112 PetscCall(PetscMalloc1(ncon, &ubvec));
113 for (i = 0; i < ncon; i++) ubvec[i] = (real_t)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 PetscCallMPI(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 = (real_t)1.05;
127
128 PetscCallMPI(MPI_Comm_size(comm, &size));
129 PetscCall(PetscMalloc1(pmat->rmap->n, &NDorder));
130 PetscCall(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 PetscCall(MatPartitioningSizesToSep_Private(subd, sizes, seps, level));
135 for (i = 0; i < pmat->rmap->n; i++) {
136 PetscInt loc;
137
138 PetscCall(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 PetscCall(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 PetscCall(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 PetscCallMPI(MPI_Comm_free(&comm));
164
165 PetscCall(PetscFree(tpwgts));
166 PetscCall(PetscFree(ubvec));
167 if (PetscLogPrintInfo) pmetis->printout = itmp;
168
169 if (bs > 1) {
170 PetscInt i, j, *newlocals;
171 PetscCall(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 PetscCall(PetscFree(locals));
176 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), bs * pmat->rmap->n, newlocals, PETSC_OWN_POINTER, partitioning));
177 } else {
178 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, locals, PETSC_OWN_POINTER, partitioning));
179 }
180 if (useND) {
181 IS ndis;
182
183 if (bs > 1) {
184 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)part), bs, pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis));
185 } else {
186 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), pmat->rmap->n, NDorder, PETSC_OWN_POINTER, &ndis));
187 }
188 PetscCall(ISSetPermutation(ndis));
189 PetscCall(PetscObjectCompose((PetscObject)*partitioning, "_petsc_matpartitioning_ndorder", (PetscObject)ndis));
190 PetscCall(ISDestroy(&ndis));
191 }
192 } else {
193 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, partitioning));
194 if (useND) {
195 IS ndis;
196
197 if (bs > 1) {
198 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)part), bs, 0, NULL, PETSC_COPY_VALUES, &ndis));
199 } else {
200 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)part), 0, NULL, PETSC_COPY_VALUES, &ndis));
201 }
202 PetscCall(ISSetPermutation(ndis));
203 PetscCall(PetscObjectCompose((PetscObject)*partitioning, "_petsc_matpartitioning_ndorder", (PetscObject)ndis));
204 PetscCall(ISDestroy(&ndis));
205 }
206 }
207 PetscCall(MatDestroy(&pmat));
208 PetscCall(MatDestroy(&amat));
209 PetscFunctionReturn(PETSC_SUCCESS);
210 }
211
212 /*
213 Uses the ParMETIS parallel matrix partitioner to compute a nested dissection ordering of the matrix in parallel
214 */
MatPartitioningApplyND_Parmetis(MatPartitioning part,IS * partitioning)215 static PetscErrorCode MatPartitioningApplyND_Parmetis(MatPartitioning part, IS *partitioning)
216 {
217 PetscFunctionBegin;
218 PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_TRUE, PETSC_FALSE, partitioning));
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
222 /*
223 Uses the ParMETIS parallel matrix partitioner to partition the matrix in parallel
224 */
MatPartitioningApply_Parmetis(MatPartitioning part,IS * partitioning)225 static PetscErrorCode MatPartitioningApply_Parmetis(MatPartitioning part, IS *partitioning)
226 {
227 PetscFunctionBegin;
228 PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_FALSE, partitioning));
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
232 /*
233 Uses the ParMETIS to improve the quality of a partition
234 */
MatPartitioningImprove_Parmetis(MatPartitioning part,IS * partitioning)235 static PetscErrorCode MatPartitioningImprove_Parmetis(MatPartitioning part, IS *partitioning)
236 {
237 PetscFunctionBegin;
238 PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_TRUE, partitioning));
239 PetscFunctionReturn(PETSC_SUCCESS);
240 }
241
MatPartitioningView_Parmetis(MatPartitioning part,PetscViewer viewer)242 static PetscErrorCode MatPartitioningView_Parmetis(MatPartitioning part, PetscViewer viewer)
243 {
244 MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
245 PetscMPIInt rank;
246 PetscBool isascii;
247
248 PetscFunctionBegin;
249 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)part), &rank));
250 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
251 if (isascii) {
252 if (pmetis->parallel == 2) {
253 PetscCall(PetscViewerASCIIPrintf(viewer, " Using parallel coarse grid partitioner\n"));
254 } else {
255 PetscCall(PetscViewerASCIIPrintf(viewer, " Using sequential coarse grid partitioner\n"));
256 }
257 PetscCall(PetscViewerASCIIPrintf(viewer, " Using %" PetscInt_FMT " fold factor\n", pmetis->foldfactor));
258 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
259 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d]Number of cuts found %" PetscInt_FMT "\n", rank, pmetis->cuts));
260 PetscCall(PetscViewerFlush(viewer));
261 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
262 }
263 PetscFunctionReturn(PETSC_SUCCESS);
264 }
265
266 /*@
267 MatPartitioningParmetisSetCoarseSequential - Use the sequential code to
268 do the partitioning of the coarse grid.
269
270 Logically Collective
271
272 Input Parameter:
273 . part - the partitioning context
274
275 Level: advanced
276
277 .seealso: `MATPARTITIONINGPARMETIS`
278 @*/
MatPartitioningParmetisSetCoarseSequential(MatPartitioning part)279 PetscErrorCode MatPartitioningParmetisSetCoarseSequential(MatPartitioning part)
280 {
281 MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
282
283 PetscFunctionBegin;
284 pmetis->parallel = 1;
285 PetscFunctionReturn(PETSC_SUCCESS);
286 }
287
288 /*@
289 MatPartitioningParmetisSetRepartition - Repartition
290 current mesh to rebalance computation.
291
292 Logically Collective
293
294 Input Parameter:
295 . part - the partitioning context
296
297 Level: advanced
298
299 .seealso: `MATPARTITIONINGPARMETIS`
300 @*/
MatPartitioningParmetisSetRepartition(MatPartitioning part)301 PetscErrorCode MatPartitioningParmetisSetRepartition(MatPartitioning part)
302 {
303 MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
304
305 PetscFunctionBegin;
306 pmetis->repartition = PETSC_TRUE;
307 PetscFunctionReturn(PETSC_SUCCESS);
308 }
309
310 /*@
311 MatPartitioningParmetisGetEdgeCut - Returns the number of edge cuts in the vertex partition.
312
313 Input Parameter:
314 . part - the partitioning context
315
316 Output Parameter:
317 . cut - the edge cut
318
319 Level: advanced
320
321 .seealso: `MATPARTITIONINGPARMETIS`
322 @*/
MatPartitioningParmetisGetEdgeCut(MatPartitioning part,PetscInt * cut)323 PetscErrorCode MatPartitioningParmetisGetEdgeCut(MatPartitioning part, PetscInt *cut)
324 {
325 MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
326
327 PetscFunctionBegin;
328 *cut = pmetis->cuts;
329 PetscFunctionReturn(PETSC_SUCCESS);
330 }
331
MatPartitioningSetFromOptions_Parmetis(MatPartitioning part,PetscOptionItems PetscOptionsObject)332 static PetscErrorCode MatPartitioningSetFromOptions_Parmetis(MatPartitioning part, PetscOptionItems PetscOptionsObject)
333 {
334 PetscBool flag = PETSC_FALSE;
335
336 PetscFunctionBegin;
337 PetscOptionsHeadBegin(PetscOptionsObject, "Set ParMeTiS partitioning options");
338 PetscCall(PetscOptionsBool("-mat_partitioning_parmetis_coarse_sequential", "Use sequential coarse partitioner", "MatPartitioningParmetisSetCoarseSequential", flag, &flag, NULL));
339 if (flag) PetscCall(MatPartitioningParmetisSetCoarseSequential(part));
340 PetscCall(PetscOptionsBool("-mat_partitioning_parmetis_repartition", "", "MatPartitioningParmetisSetRepartition", flag, &flag, NULL));
341 if (flag) PetscCall(MatPartitioningParmetisSetRepartition(part));
342 PetscOptionsHeadEnd();
343 PetscFunctionReturn(PETSC_SUCCESS);
344 }
345
MatPartitioningDestroy_Parmetis(MatPartitioning part)346 static PetscErrorCode MatPartitioningDestroy_Parmetis(MatPartitioning part)
347 {
348 MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
349
350 PetscFunctionBegin;
351 PetscCall(PetscFree(pmetis));
352 PetscFunctionReturn(PETSC_SUCCESS);
353 }
354
355 /*MC
356 MATPARTITIONINGPARMETIS - Creates a partitioning context via the external package PARMETIS.
357
358 Collective
359
360 Input Parameter:
361 . part - the partitioning context
362
363 Options Database Key:
364 . -mat_partitioning_parmetis_coarse_sequential - use sequential PARMETIS coarse partitioner
365
366 Level: beginner
367
368 Note:
369 See https://www-users.cs.umn.edu/~karypis/metis/
370
371 .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MatPartitioningParmetisSetCoarseSequential()`, `MatPartitioningParmetisSetRepartition()`,
372 `MatPartitioningParmetisGetEdgeCut()`
373 M*/
374
MatPartitioningCreate_Parmetis(MatPartitioning part)375 PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Parmetis(MatPartitioning part)
376 {
377 MatPartitioning_Parmetis *pmetis;
378
379 PetscFunctionBegin;
380 PetscCall(PetscNew(&pmetis));
381 part->data = (void *)pmetis;
382
383 pmetis->cuts = 0; /* output variable */
384 pmetis->foldfactor = 150; /*folding factor */
385 pmetis->parallel = 2; /* use parallel partitioner for coarse grid */
386 pmetis->indexing = 0; /* index numbering starts from 0 */
387 pmetis->printout = 0; /* print no output while running */
388 pmetis->repartition = PETSC_FALSE;
389
390 part->ops->apply = MatPartitioningApply_Parmetis;
391 part->ops->applynd = MatPartitioningApplyND_Parmetis;
392 part->ops->improve = MatPartitioningImprove_Parmetis;
393 part->ops->view = MatPartitioningView_Parmetis;
394 part->ops->destroy = MatPartitioningDestroy_Parmetis;
395 part->ops->setfromoptions = MatPartitioningSetFromOptions_Parmetis;
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
399 /*@
400 MatMeshToCellGraph - Convert a mesh to a cell graph.
401
402 Collective
403
404 Input Parameters:
405 + mesh - the graph that represents the coupling of the vertices of the mesh
406 - ncommonnodes - mesh elements that share this number of common nodes are considered neighbors, use 2 for triangles and
407 quadrilaterials, 3 for tetrahedrals and 4 for hexahedrals
408
409 Output Parameter:
410 . dual - the dual graph
411
412 Level: advanced
413
414 Notes:
415 Uses the ParMETIS package to convert a `Mat` that represents coupling of vertices of a mesh
416 to a `Mat` the represents the graph of the coupling between cells (the "dual" graph) and is
417 suitable for partitioning with the `MatPartitioning` object. Use this to partition cells of a
418 mesh.
419
420 Currently requires ParMetis to be installed and uses ParMETIS_V3_Mesh2Dual()
421
422 Each row of the mesh object represents a single cell in the mesh. For triangles it has 3 entries, quadrilaterials 4 entries,
423 tetrahedrals 4 entries and hexahedrals 8 entries. You can mix triangles and quadrilaterals in the same mesh, but cannot
424 mix tetrahedrals and hexahedrals
425 The columns of each row of the `Mat` mesh are the global vertex numbers of the vertices of that row's cell.
426 The number of rows in mesh is number of cells, the number of columns is the number of vertices.
427
428 .seealso: `MatCreateMPIAdj()`, `MatPartitioningCreate()`
429 @*/
MatMeshToCellGraph(Mat mesh,PetscInt ncommonnodes,Mat * dual)430 PetscErrorCode MatMeshToCellGraph(Mat mesh, PetscInt ncommonnodes, Mat *dual)
431 {
432 PetscInt *newxadj, *newadjncy;
433 PetscInt numflag = 0;
434 Mat_MPIAdj *adj = (Mat_MPIAdj *)mesh->data, *newadj;
435 PetscBool flg;
436 MPI_Comm comm;
437
438 PetscFunctionBegin;
439 PetscCall(PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg));
440 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Must use MPIAdj matrix type");
441
442 PetscCall(PetscObjectGetComm((PetscObject)mesh, &comm));
443 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));
444 PetscCall(MatCreateMPIAdj(PetscObjectComm((PetscObject)mesh), mesh->rmap->n, mesh->rmap->N, newxadj, newadjncy, NULL, dual));
445 newadj = (Mat_MPIAdj *)(*dual)->data;
446
447 newadj->freeaijwithfree = PETSC_TRUE; /* signal the matrix should be freed with system free since space was allocated by ParMETIS */
448 PetscFunctionReturn(PETSC_SUCCESS);
449 }
450