Actual source code: pmetis.c

  1: #include <../src/mat/impls/adj/mpi/mpiadj.h>

  3: /*
  4:    Currently using ParMetis-4.0.2
  5: */

  7: #include <parmetis.h>

  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;

 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)

 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)

 36: #define PetscCallParmetis(func, args) PetscCallParmetis_(PetscStringize(func), func, args)

 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;

 46:   PetscFunctionBegin;
 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)));

 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:     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:     }

 82:     PetscCall(PetscMalloc1(pmat->rmap->n, &locals));

 84:     if (isImprove) {
 85:       PetscInt        i;
 86:       const PetscInt *part_indices;
 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:     }

 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:     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] = part->part_weights[i * nparts + j];
107:         } else {
108:           tpwgts[i * nparts + j] = 1. / nparts;
109:         }
110:       }
111:     }
112:     PetscCall(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:     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 = 1.05;

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;

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));

165:     PetscCall(PetscFree(tpwgts));
166:     PetscCall(PetscFree(ubvec));
167:     if (PetscLogPrintInfo) pmetis->printout = itmp;

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;

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;

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: }

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:   PetscFunctionBegin;
218:   PetscCall(MatPartitioningApply_Parmetis_Private(part, PETSC_TRUE, PETSC_FALSE, partitioning));
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: /*
223:    Uses the ParMETIS parallel matrix partitioner to partition the matrix in parallel
224: */
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: }

232: /*
233:    Uses the ParMETIS to improve the quality  of a partition
234: */
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: }

242: static PetscErrorCode MatPartitioningView_Parmetis(MatPartitioning part, PetscViewer viewer)
243: {
244:   MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;
245:   PetscMPIInt               rank;
246:   PetscBool                 iascii;

248:   PetscFunctionBegin;
249:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)part), &rank));
250:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
251:   if (iascii) {
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: }

266: /*@
267:   MatPartitioningParmetisSetCoarseSequential - Use the sequential code to
268:   do the partitioning of the coarse grid.

270:   Logically Collective

272:   Input Parameter:
273: . part - the partitioning context

275:   Level: advanced

277: .seealso: `MATPARTITIONINGPARMETIS`
278: @*/
279: PetscErrorCode MatPartitioningParmetisSetCoarseSequential(MatPartitioning part)
280: {
281:   MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;

283:   PetscFunctionBegin;
284:   pmetis->parallel = 1;
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: /*@
289:   MatPartitioningParmetisSetRepartition - Repartition
290:   current mesh to rebalance computation.

292:   Logically Collective

294:   Input Parameter:
295: . part - the partitioning context

297:   Level: advanced

299: .seealso: `MATPARTITIONINGPARMETIS`
300: @*/
301: PetscErrorCode MatPartitioningParmetisSetRepartition(MatPartitioning part)
302: {
303:   MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;

305:   PetscFunctionBegin;
306:   pmetis->repartition = PETSC_TRUE;
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

310: /*@
311:   MatPartitioningParmetisGetEdgeCut - Returns the number of edge cuts in the vertex partition.

313:   Input Parameter:
314: . part - the partitioning context

316:   Output Parameter:
317: . cut - the edge cut

319:   Level: advanced

321: .seealso: `MATPARTITIONINGPARMETIS`
322: @*/
323: PetscErrorCode MatPartitioningParmetisGetEdgeCut(MatPartitioning part, PetscInt *cut)
324: {
325:   MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;

327:   PetscFunctionBegin;
328:   *cut = pmetis->cuts;
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: static PetscErrorCode MatPartitioningSetFromOptions_Parmetis(MatPartitioning part, PetscOptionItems *PetscOptionsObject)
333: {
334:   PetscBool flag = PETSC_FALSE;

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: }

346: static PetscErrorCode MatPartitioningDestroy_Parmetis(MatPartitioning part)
347: {
348:   MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis *)part->data;

350:   PetscFunctionBegin;
351:   PetscCall(PetscFree(pmetis));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: /*MC
356:    MATPARTITIONINGPARMETIS - Creates a partitioning context via the external package PARMETIS.

358:    Collective

360:    Input Parameter:
361: .  part - the partitioning context

363:    Options Database Key:
364: .  -mat_partitioning_parmetis_coarse_sequential - use sequential PARMETIS coarse partitioner

366:    Level: beginner

368:    Note:
369:     See https://www-users.cs.umn.edu/~karypis/metis/

371: .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MatPartitioningParmetisSetCoarseSequential()`, `MatPartitioningParmetisSetRepartition()`,
372:           `MatPartitioningParmetisGetEdgeCut()`
373: M*/

375: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Parmetis(MatPartitioning part)
376: {
377:   MatPartitioning_Parmetis *pmetis;

379:   PetscFunctionBegin;
380:   PetscCall(PetscNew(&pmetis));
381:   part->data = (void *)pmetis;

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;

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: }

399: /*@
400:   MatMeshToCellGraph - Convert a mesh to a cell graph.

402:   Collective

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

409:   Output Parameter:
410: . dual - the dual graph

412:   Level: advanced

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.

420:   Currently requires ParMetis to be installed and uses ParMETIS_V3_Mesh2Dual()

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.

428: .seealso: `MatCreateMPIAdj()`, `MatPartitioningCreate()`
429: @*/
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;

438:   PetscFunctionBegin;
439:   PetscCall(PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg));
440:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Must use MPIAdj matrix type");

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;

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: }