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