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