Actual source code: plexpartition.c

petsc-master 2020-04-04
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/hashseti.h>

  4: PetscClassId PETSCPARTITIONER_CLASSID = 0;

  6: PetscFunctionList PetscPartitionerList              = NULL;
  7: PetscBool         PetscPartitionerRegisterAllCalled = PETSC_FALSE;

  9: PetscBool ChacoPartitionercite = PETSC_FALSE;
 10: const char ChacoPartitionerCitation[] = "@inproceedings{Chaco95,\n"
 11:                                "  author    = {Bruce Hendrickson and Robert Leland},\n"
 12:                                "  title     = {A multilevel algorithm for partitioning graphs},\n"
 13:                                "  booktitle = {Supercomputing '95: Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM)},"
 14:                                "  isbn      = {0-89791-816-9},\n"
 15:                                "  pages     = {28},\n"
 16:                                "  doi       = {https://doi.acm.org/10.1145/224170.224228},\n"
 17:                                "  publisher = {ACM Press},\n"
 18:                                "  address   = {New York},\n"
 19:                                "  year      = {1995}\n}\n";

 21: PetscBool ParMetisPartitionercite = PETSC_FALSE;
 22: const char ParMetisPartitionerCitation[] = "@article{KarypisKumar98,\n"
 23:                                "  author  = {George Karypis and Vipin Kumar},\n"
 24:                                "  title   = {A Parallel Algorithm for Multilevel Graph Partitioning and Sparse Matrix Ordering},\n"
 25:                                "  journal = {Journal of Parallel and Distributed Computing},\n"
 26:                                "  volume  = {48},\n"
 27:                                "  pages   = {71--85},\n"
 28:                                "  year    = {1998}\n}\n";

 30: PetscBool PTScotchPartitionercite = PETSC_FALSE;
 31: const char PTScotchPartitionerCitation[] =
 32:   "@article{PTSCOTCH,\n"
 33:   "  author  = {C. Chevalier and F. Pellegrini},\n"
 34:   "  title   = {{PT-SCOTCH}: a tool for efficient parallel graph ordering},\n"
 35:   "  journal = {Parallel Computing},\n"
 36:   "  volume  = {34},\n"
 37:   "  number  = {6},\n"
 38:   "  pages   = {318--331},\n"
 39:   "  year    = {2008},\n"
 40:   "  doi     = {https://doi.org/10.1016/j.parco.2007.12.001}\n"
 41:   "}\n";


 44: PETSC_STATIC_INLINE PetscInt DMPlex_GlobalID(PetscInt point) { return point >= 0 ? point : -(point+1); }

 46: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
 47: {
 48:   PetscInt       dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
 49:   PetscInt      *adj = NULL, *vOffsets = NULL, *graph = NULL;
 50:   IS             cellNumbering;
 51:   const PetscInt *cellNum;
 52:   PetscBool      useCone, useClosure;
 53:   PetscSection   section;
 54:   PetscSegBuffer adjBuffer;
 55:   PetscSF        sfPoint;
 56:   PetscInt       *adjCells = NULL, *remoteCells = NULL;
 57:   const PetscInt *local;
 58:   PetscInt       nroots, nleaves, l;
 59:   PetscMPIInt    rank;

 63:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
 64:   DMGetDimension(dm, &dim);
 65:   DMPlexGetDepth(dm, &depth);
 66:   if (dim != depth) {
 67:     /* We do not handle the uninterpolated case here */
 68:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
 69:     /* DMPlexCreateNeighborCSR does not make a numbering */
 70:     if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
 71:     /* Different behavior for empty graphs */
 72:     if (!*numVertices) {
 73:       PetscMalloc1(1, offsets);
 74:       (*offsets)[0] = 0;
 75:     }
 76:     /* Broken in parallel */
 77:     if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
 78:     return(0);
 79:   }
 80:   DMGetPointSF(dm, &sfPoint);
 81:   DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
 82:   /* Build adjacency graph via a section/segbuffer */
 83:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &section);
 84:   PetscSectionSetChart(section, pStart, pEnd);
 85:   PetscSegBufferCreate(sizeof(PetscInt),1000,&adjBuffer);
 86:   /* Always use FVM adjacency to create partitioner graph */
 87:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
 88:   DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
 89:   DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
 90:   if (globalNumbering) {
 91:     PetscObjectReference((PetscObject)cellNumbering);
 92:     *globalNumbering = cellNumbering;
 93:   }
 94:   ISGetIndices(cellNumbering, &cellNum);
 95:   /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
 96:      Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
 97:    */
 98:   PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
 99:   if (nroots >= 0) {
100:     PetscInt fStart, fEnd, f;

102:     PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
103:     DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
104:     for (l = 0; l < nroots; ++l) adjCells[l] = -3;
105:     for (f = fStart; f < fEnd; ++f) {
106:       const PetscInt *support;
107:       PetscInt        supportSize;

109:       DMPlexGetSupport(dm, f, &support);
110:       DMPlexGetSupportSize(dm, f, &supportSize);
111:       if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
112:       else if (supportSize == 2) {
113:         PetscFindInt(support[0], nleaves, local, &p);
114:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1]]);
115:         PetscFindInt(support[1], nleaves, local, &p);
116:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
117:       }
118:       /* Handle non-conforming meshes */
119:       if (supportSize > 2) {
120:         PetscInt        numChildren, i;
121:         const PetscInt *children;

123:         DMPlexGetTreeChildren(dm, f, &numChildren, &children);
124:         for (i = 0; i < numChildren; ++i) {
125:           const PetscInt child = children[i];
126:           if (fStart <= child && child < fEnd) {
127:             DMPlexGetSupport(dm, child, &support);
128:             DMPlexGetSupportSize(dm, child, &supportSize);
129:             if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
130:             else if (supportSize == 2) {
131:               PetscFindInt(support[0], nleaves, local, &p);
132:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1]]);
133:               PetscFindInt(support[1], nleaves, local, &p);
134:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
135:             }
136:           }
137:         }
138:       }
139:     }
140:     for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
141:     PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells);
142:     PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells);
143:     PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
144:     PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
145:   }
146:   /* Combine local and global adjacencies */
147:   for (*numVertices = 0, p = pStart; p < pEnd; p++) {
148:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
149:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
150:     /* Add remote cells */
151:     if (remoteCells) {
152:       const PetscInt gp = DMPlex_GlobalID(cellNum[p]);
153:       PetscInt       coneSize, numChildren, c, i;
154:       const PetscInt *cone, *children;

156:       DMPlexGetCone(dm, p, &cone);
157:       DMPlexGetConeSize(dm, p, &coneSize);
158:       for (c = 0; c < coneSize; ++c) {
159:         const PetscInt point = cone[c];
160:         if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
161:           PetscInt *PETSC_RESTRICT pBuf;
162:           PetscSectionAddDof(section, p, 1);
163:           PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
164:           *pBuf = remoteCells[point];
165:         }
166:         /* Handle non-conforming meshes */
167:         DMPlexGetTreeChildren(dm, point, &numChildren, &children);
168:         for (i = 0; i < numChildren; ++i) {
169:           const PetscInt child = children[i];
170:           if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
171:             PetscInt *PETSC_RESTRICT pBuf;
172:             PetscSectionAddDof(section, p, 1);
173:             PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
174:             *pBuf = remoteCells[child];
175:           }
176:         }
177:       }
178:     }
179:     /* Add local cells */
180:     adjSize = PETSC_DETERMINE;
181:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
182:     for (a = 0; a < adjSize; ++a) {
183:       const PetscInt point = adj[a];
184:       if (point != p && pStart <= point && point < pEnd) {
185:         PetscInt *PETSC_RESTRICT pBuf;
186:         PetscSectionAddDof(section, p, 1);
187:         PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
188:         *pBuf = DMPlex_GlobalID(cellNum[point]);
189:       }
190:     }
191:     (*numVertices)++;
192:   }
193:   PetscFree(adj);
194:   PetscFree2(adjCells, remoteCells);
195:   DMSetBasicAdjacency(dm, useCone, useClosure);

197:   /* Derive CSR graph from section/segbuffer */
198:   PetscSectionSetUp(section);
199:   PetscSectionGetStorageSize(section, &size);
200:   PetscMalloc1(*numVertices+1, &vOffsets);
201:   for (idx = 0, p = pStart; p < pEnd; p++) {
202:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
203:     PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
204:   }
205:   vOffsets[*numVertices] = size;
206:   PetscSegBufferExtractAlloc(adjBuffer, &graph);

208:   if (nroots >= 0) {
209:     /* Filter out duplicate edges using section/segbuffer */
210:     PetscSectionReset(section);
211:     PetscSectionSetChart(section, 0, *numVertices);
212:     for (p = 0; p < *numVertices; p++) {
213:       PetscInt start = vOffsets[p], end = vOffsets[p+1];
214:       PetscInt numEdges = end-start, *PETSC_RESTRICT edges;
215:       PetscSortRemoveDupsInt(&numEdges, &graph[start]);
216:       PetscSectionSetDof(section, p, numEdges);
217:       PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
218:       PetscArraycpy(edges, &graph[start], numEdges);
219:     }
220:     PetscFree(vOffsets);
221:     PetscFree(graph);
222:     /* Derive CSR graph from section/segbuffer */
223:     PetscSectionSetUp(section);
224:     PetscSectionGetStorageSize(section, &size);
225:     PetscMalloc1(*numVertices+1, &vOffsets);
226:     for (idx = 0, p = 0; p < *numVertices; p++) {
227:       PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
228:     }
229:     vOffsets[*numVertices] = size;
230:     PetscSegBufferExtractAlloc(adjBuffer, &graph);
231:   } else {
232:     /* Sort adjacencies (not strictly necessary) */
233:     for (p = 0; p < *numVertices; p++) {
234:       PetscInt start = vOffsets[p], end = vOffsets[p+1];
235:       PetscSortInt(end-start, &graph[start]);
236:     }
237:   }

239:   if (offsets) *offsets = vOffsets;
240:   if (adjacency) *adjacency = graph;

242:   /* Cleanup */
243:   PetscSegBufferDestroy(&adjBuffer);
244:   PetscSectionDestroy(&section);
245:   ISRestoreIndices(cellNumbering, &cellNum);
246:   ISDestroy(&cellNumbering);
247:   return(0);
248: }

250: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
251: {
252:   Mat            conn, CSR;
253:   IS             fis, cis, cis_own;
254:   PetscSF        sfPoint;
255:   const PetscInt *rows, *cols, *ii, *jj;
256:   PetscInt       *idxs,*idxs2;
257:   PetscInt       dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
258:   PetscMPIInt    rank;
259:   PetscBool      flg;

263:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
264:   DMGetDimension(dm, &dim);
265:   DMPlexGetDepth(dm, &depth);
266:   if (dim != depth) {
267:     /* We do not handle the uninterpolated case here */
268:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
269:     /* DMPlexCreateNeighborCSR does not make a numbering */
270:     if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
271:     /* Different behavior for empty graphs */
272:     if (!*numVertices) {
273:       PetscMalloc1(1, offsets);
274:       (*offsets)[0] = 0;
275:     }
276:     /* Broken in parallel */
277:     if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
278:     return(0);
279:   }
280:   /* Interpolated and parallel case */

282:   /* numbering */
283:   DMGetPointSF(dm, &sfPoint);
284:   DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
285:   DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
286:   DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
287:   DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
288:   if (globalNumbering) {
289:     ISDuplicate(cis, globalNumbering);
290:   }

292:   /* get positive global ids and local sizes for facets and cells */
293:   ISGetLocalSize(fis, &m);
294:   ISGetIndices(fis, &rows);
295:   PetscMalloc1(m, &idxs);
296:   for (i = 0, floc = 0; i < m; i++) {
297:     const PetscInt p = rows[i];

299:     if (p < 0) {
300:       idxs[i] = -(p+1);
301:     } else {
302:       idxs[i] = p;
303:       floc   += 1;
304:     }
305:   }
306:   ISRestoreIndices(fis, &rows);
307:   ISDestroy(&fis);
308:   ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);

310:   ISGetLocalSize(cis, &m);
311:   ISGetIndices(cis, &cols);
312:   PetscMalloc1(m, &idxs);
313:   PetscMalloc1(m, &idxs2);
314:   for (i = 0, cloc = 0; i < m; i++) {
315:     const PetscInt p = cols[i];

317:     if (p < 0) {
318:       idxs[i] = -(p+1);
319:     } else {
320:       idxs[i]       = p;
321:       idxs2[cloc++] = p;
322:     }
323:   }
324:   ISRestoreIndices(cis, &cols);
325:   ISDestroy(&cis);
326:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
327:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);

329:   /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
330:   MatCreate(PetscObjectComm((PetscObject)dm), &conn);
331:   MatSetSizes(conn, floc, cloc, M, N);
332:   MatSetType(conn, MATMPIAIJ);
333:   DMPlexGetMaxSizes(dm, NULL, &lm);
334:   MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dm));
335:   MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);

337:   /* Assemble matrix */
338:   ISGetIndices(fis, &rows);
339:   ISGetIndices(cis, &cols);
340:   for (c = cStart; c < cEnd; c++) {
341:     const PetscInt *cone;
342:     PetscInt        coneSize, row, col, f;

344:     col  = cols[c-cStart];
345:     DMPlexGetCone(dm, c, &cone);
346:     DMPlexGetConeSize(dm, c, &coneSize);
347:     for (f = 0; f < coneSize; f++) {
348:       const PetscScalar v = 1.0;
349:       const PetscInt *children;
350:       PetscInt        numChildren, ch;

352:       row  = rows[cone[f]-fStart];
353:       MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);

355:       /* non-conforming meshes */
356:       DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
357:       for (ch = 0; ch < numChildren; ch++) {
358:         const PetscInt child = children[ch];

360:         if (child < fStart || child >= fEnd) continue;
361:         row  = rows[child-fStart];
362:         MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
363:       }
364:     }
365:   }
366:   ISRestoreIndices(fis, &rows);
367:   ISRestoreIndices(cis, &cols);
368:   ISDestroy(&fis);
369:   ISDestroy(&cis);
370:   MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
371:   MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);

373:   /* Get parallel CSR by doing conn^T * conn */
374:   MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
375:   MatDestroy(&conn);

377:   /* extract local part of the CSR */
378:   MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
379:   MatDestroy(&CSR);
380:   MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
381:   if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");

383:   /* get back requested output */
384:   if (numVertices) *numVertices = m;
385:   if (offsets) {
386:     PetscCalloc1(m+1, &idxs);
387:     for (i = 1; i < m+1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
388:     *offsets = idxs;
389:   }
390:   if (adjacency) {
391:     PetscMalloc1(ii[m] - m, &idxs);
392:     ISGetIndices(cis_own, &rows);
393:     for (i = 0, c = 0; i < m; i++) {
394:       PetscInt j, g = rows[i];

396:       for (j = ii[i]; j < ii[i+1]; j++) {
397:         if (jj[j] == g) continue; /* again, self-connectivity */
398:         idxs[c++] = jj[j];
399:       }
400:     }
401:     if (c != ii[m] - m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %D != %D",c,ii[m]-m);
402:     ISRestoreIndices(cis_own, &rows);
403:     *adjacency = idxs;
404:   }

406:   /* cleanup */
407:   ISDestroy(&cis_own);
408:   MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
409:   if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
410:   MatDestroy(&conn);
411:   return(0);
412: }

414: /*@C
415:   DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner

417:   Input Parameters:
418: + dm      - The mesh DM dm
419: - height  - Height of the strata from which to construct the graph

421:   Output Parameter:
422: + numVertices     - Number of vertices in the graph
423: . offsets         - Point offsets in the graph
424: . adjacency       - Point connectivity in the graph
425: - globalNumbering - A map from the local cell numbering to the global numbering used in "adjacency".  Negative indicates that the cell is a duplicate from another process.

427:   The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function
428:   representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().

430:   Level: developer

432: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
433: @*/
434: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
435: {
437:   PetscBool      usemat = PETSC_FALSE;

440:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_csr_via_mat", &usemat, NULL);
441:   if (usemat) {
442:     DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
443:   } else {
444:     DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
445:   }
446:   return(0);
447: }

449: /*@C
450:   DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.

452:   Collective on DM

454:   Input Arguments:
455: + dm - The DMPlex
456: - cellHeight - The height of mesh points to treat as cells (default should be 0)

458:   Output Arguments:
459: + numVertices - The number of local vertices in the graph, or cells in the mesh.
460: . offsets     - The offset to the adjacency list for each cell
461: - adjacency   - The adjacency list for all cells

463:   Note: This is suitable for input to a mesh partitioner like ParMetis.

465:   Level: advanced

467: .seealso: DMPlexCreate()
468: @*/
469: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
470: {
471:   const PetscInt maxFaceCases = 30;
472:   PetscInt       numFaceCases = 0;
473:   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
474:   PetscInt      *off, *adj;
475:   PetscInt      *neighborCells = NULL;
476:   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;

480:   /* For parallel partitioning, I think you have to communicate supports */
481:   DMGetDimension(dm, &dim);
482:   cellDim = dim - cellHeight;
483:   DMPlexGetDepth(dm, &depth);
484:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
485:   if (cEnd - cStart == 0) {
486:     if (numVertices) *numVertices = 0;
487:     if (offsets)   *offsets   = NULL;
488:     if (adjacency) *adjacency = NULL;
489:     return(0);
490:   }
491:   numCells  = cEnd - cStart;
492:   faceDepth = depth - cellHeight;
493:   if (dim == depth) {
494:     PetscInt f, fStart, fEnd;

496:     PetscCalloc1(numCells+1, &off);
497:     /* Count neighboring cells */
498:     DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
499:     for (f = fStart; f < fEnd; ++f) {
500:       const PetscInt *support;
501:       PetscInt        supportSize;
502:       DMPlexGetSupportSize(dm, f, &supportSize);
503:       DMPlexGetSupport(dm, f, &support);
504:       if (supportSize == 2) {
505:         PetscInt numChildren;

507:         DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
508:         if (!numChildren) {
509:           ++off[support[0]-cStart+1];
510:           ++off[support[1]-cStart+1];
511:         }
512:       }
513:     }
514:     /* Prefix sum */
515:     for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
516:     if (adjacency) {
517:       PetscInt *tmp;

519:       PetscMalloc1(off[numCells], &adj);
520:       PetscMalloc1(numCells+1, &tmp);
521:       PetscArraycpy(tmp, off, numCells+1);
522:       /* Get neighboring cells */
523:       for (f = fStart; f < fEnd; ++f) {
524:         const PetscInt *support;
525:         PetscInt        supportSize;
526:         DMPlexGetSupportSize(dm, f, &supportSize);
527:         DMPlexGetSupport(dm, f, &support);
528:         if (supportSize == 2) {
529:           PetscInt numChildren;

531:           DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
532:           if (!numChildren) {
533:             adj[tmp[support[0]-cStart]++] = support[1];
534:             adj[tmp[support[1]-cStart]++] = support[0];
535:           }
536:         }
537:       }
538: #if defined(PETSC_USE_DEBUG)
539:       for (c = 0; c < cEnd-cStart; ++c) if (tmp[c] != off[c+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %d != %d for cell %d", tmp[c], off[c], c+cStart);
540: #endif
541:       PetscFree(tmp);
542:     }
543:     if (numVertices) *numVertices = numCells;
544:     if (offsets)   *offsets   = off;
545:     if (adjacency) *adjacency = adj;
546:     return(0);
547:   }
548:   /* Setup face recognition */
549:   if (faceDepth == 1) {
550:     PetscInt cornersSeen[30] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /* Could use PetscBT */

552:     for (c = cStart; c < cEnd; ++c) {
553:       PetscInt corners;

555:       DMPlexGetConeSize(dm, c, &corners);
556:       if (!cornersSeen[corners]) {
557:         PetscInt nFV;

559:         if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
560:         cornersSeen[corners] = 1;

562:         DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);

564:         numFaceVertices[numFaceCases++] = nFV;
565:       }
566:     }
567:   }
568:   PetscCalloc1(numCells+1, &off);
569:   /* Count neighboring cells */
570:   for (cell = cStart; cell < cEnd; ++cell) {
571:     PetscInt numNeighbors = PETSC_DETERMINE, n;

573:     DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
574:     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
575:     for (n = 0; n < numNeighbors; ++n) {
576:       PetscInt        cellPair[2];
577:       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
578:       PetscInt        meetSize = 0;
579:       const PetscInt *meet    = NULL;

581:       cellPair[0] = cell; cellPair[1] = neighborCells[n];
582:       if (cellPair[0] == cellPair[1]) continue;
583:       if (!found) {
584:         DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
585:         if (meetSize) {
586:           PetscInt f;

588:           for (f = 0; f < numFaceCases; ++f) {
589:             if (numFaceVertices[f] == meetSize) {
590:               found = PETSC_TRUE;
591:               break;
592:             }
593:           }
594:         }
595:         DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
596:       }
597:       if (found) ++off[cell-cStart+1];
598:     }
599:   }
600:   /* Prefix sum */
601:   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];

603:   if (adjacency) {
604:     PetscMalloc1(off[numCells], &adj);
605:     /* Get neighboring cells */
606:     for (cell = cStart; cell < cEnd; ++cell) {
607:       PetscInt numNeighbors = PETSC_DETERMINE, n;
608:       PetscInt cellOffset   = 0;

610:       DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
611:       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
612:       for (n = 0; n < numNeighbors; ++n) {
613:         PetscInt        cellPair[2];
614:         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
615:         PetscInt        meetSize = 0;
616:         const PetscInt *meet    = NULL;

618:         cellPair[0] = cell; cellPair[1] = neighborCells[n];
619:         if (cellPair[0] == cellPair[1]) continue;
620:         if (!found) {
621:           DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
622:           if (meetSize) {
623:             PetscInt f;

625:             for (f = 0; f < numFaceCases; ++f) {
626:               if (numFaceVertices[f] == meetSize) {
627:                 found = PETSC_TRUE;
628:                 break;
629:               }
630:             }
631:           }
632:           DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
633:         }
634:         if (found) {
635:           adj[off[cell-cStart]+cellOffset] = neighborCells[n];
636:           ++cellOffset;
637:         }
638:       }
639:     }
640:   }
641:   PetscFree(neighborCells);
642:   if (numVertices) *numVertices = numCells;
643:   if (offsets)   *offsets   = off;
644:   if (adjacency) *adjacency = adj;
645:   return(0);
646: }

648: /*@C
649:   PetscPartitionerRegister - Adds a new PetscPartitioner implementation

651:   Not Collective

653:   Input Parameters:
654: + name        - The name of a new user-defined creation routine
655: - create_func - The creation routine itself

657:   Notes:
658:   PetscPartitionerRegister() may be called multiple times to add several user-defined PetscPartitioners

660:   Sample usage:
661: .vb
662:     PetscPartitionerRegister("my_part", MyPetscPartitionerCreate);
663: .ve

665:   Then, your PetscPartitioner type can be chosen with the procedural interface via
666: .vb
667:     PetscPartitionerCreate(MPI_Comm, PetscPartitioner *);
668:     PetscPartitionerSetType(PetscPartitioner, "my_part");
669: .ve
670:    or at runtime via the option
671: .vb
672:     -petscpartitioner_type my_part
673: .ve

675:   Level: advanced

677: .seealso: PetscPartitionerRegisterAll(), PetscPartitionerRegisterDestroy()

679: @*/
680: PetscErrorCode PetscPartitionerRegister(const char sname[], PetscErrorCode (*function)(PetscPartitioner))
681: {

685:   PetscFunctionListAdd(&PetscPartitionerList, sname, function);
686:   return(0);
687: }

689: /*@C
690:   PetscPartitionerSetType - Builds a particular PetscPartitioner

692:   Collective on PetscPartitioner

694:   Input Parameters:
695: + part - The PetscPartitioner object
696: - name - The kind of partitioner

698:   Options Database Key:
699: . -petscpartitioner_type <type> - Sets the PetscPartitioner type; use -help for a list of available types

701:   Level: intermediate

703: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate()
704: @*/
705: PetscErrorCode PetscPartitionerSetType(PetscPartitioner part, PetscPartitionerType name)
706: {
707:   PetscErrorCode (*r)(PetscPartitioner);
708:   PetscBool      match;

713:   PetscObjectTypeCompare((PetscObject) part, name, &match);
714:   if (match) return(0);

716:   PetscFunctionListFind(PetscPartitionerList, name, &r);
717:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscPartitioner type: %s", name);

719:   if (part->ops->destroy) {
720:     (*part->ops->destroy)(part);
721:   }
722:   part->noGraph = PETSC_FALSE;
723:   PetscMemzero(part->ops, sizeof(*part->ops));
724:   PetscObjectChangeTypeName((PetscObject) part, name);
725:   (*r)(part);
726:   return(0);
727: }

729: /*@C
730:   PetscPartitionerGetType - Gets the PetscPartitioner type name (as a string) from the object.

732:   Not Collective

734:   Input Parameter:
735: . part - The PetscPartitioner

737:   Output Parameter:
738: . name - The PetscPartitioner type name

740:   Level: intermediate

742: .seealso: PetscPartitionerSetType(), PetscPartitionerCreate()
743: @*/
744: PetscErrorCode PetscPartitionerGetType(PetscPartitioner part, PetscPartitionerType *name)
745: {
749:   *name = ((PetscObject) part)->type_name;
750:   return(0);
751: }

753: /*@C
754:    PetscPartitionerViewFromOptions - View from Options

756:    Collective on PetscPartitioner

758:    Input Parameters:
759: +  A - the PetscPartitioner object
760: .  obj - Optional object
761: -  name - command line option

763:    Level: intermediate
764: .seealso:  PetscPartitionerView(), PetscObjectViewFromOptions()
765: @*/
766: PetscErrorCode PetscPartitionerViewFromOptions(PetscPartitioner A,PetscObject obj,const char name[])
767: {

772:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
773:   return(0);
774: }

776: /*@C
777:   PetscPartitionerView - Views a PetscPartitioner

779:   Collective on PetscPartitioner

781:   Input Parameter:
782: + part - the PetscPartitioner object to view
783: - v    - the viewer

785:   Level: developer

787: .seealso: PetscPartitionerDestroy()
788: @*/
789: PetscErrorCode PetscPartitionerView(PetscPartitioner part, PetscViewer v)
790: {
791:   PetscMPIInt    size;
792:   PetscBool      isascii;

797:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) part), &v);}
798:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
799:   if (isascii) {
800:     MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
801:     PetscViewerASCIIPrintf(v, "Graph Partitioner: %d MPI Process%s\n", size, size > 1 ? "es" : "");
802:     PetscViewerASCIIPrintf(v, "  type: %s\n", part->hdr.type_name);
803:     PetscViewerASCIIPrintf(v, "  edge cut: %D\n", part->edgeCut);
804:     PetscViewerASCIIPrintf(v, "  balance: %.2g\n", part->balance);
805:     PetscViewerASCIIPrintf(v, "  use vertex weights: %d\n", part->usevwgt);
806:   }
807:   if (part->ops->view) {(*part->ops->view)(part, v);}
808:   return(0);
809: }

811: static PetscErrorCode PetscPartitionerGetDefaultType(const char *currentType, const char **defaultType)
812: {
814:   if (!currentType) {
815: #if defined(PETSC_HAVE_PARMETIS)
816:     *defaultType = PETSCPARTITIONERPARMETIS;
817: #elif defined(PETSC_HAVE_PTSCOTCH)
818:     *defaultType = PETSCPARTITIONERPTSCOTCH;
819: #elif defined(PETSC_HAVE_CHACO)
820:     *defaultType = PETSCPARTITIONERCHACO;
821: #else
822:     *defaultType = PETSCPARTITIONERSIMPLE;
823: #endif
824:   } else {
825:     *defaultType = currentType;
826:   }
827:   return(0);
828: }

830: /*@
831:   PetscPartitionerSetFromOptions - sets parameters in a PetscPartitioner from the options database

833:   Collective on PetscPartitioner

835:   Input Parameter:
836: . part - the PetscPartitioner object to set options for

838:   Options Database Keys:
839: +  -petscpartitioner_type <type> - Sets the PetscPartitioner type; use -help for a list of available types
840: .  -petscpartitioner_use_vertex_weights - Uses weights associated with the graph vertices
841: -  -petscpartitioner_view_graph - View the graph each time PetscPartitionerPartition is called. Viewer can be customized, see PetscOptionsGetViewer()

843:   Level: developer

845: .seealso: PetscPartitionerView(), PetscPartitionerSetType(), PetscPartitionerPartition()
846: @*/
847: PetscErrorCode PetscPartitionerSetFromOptions(PetscPartitioner part)
848: {
849:   const char    *defaultType;
850:   char           name[256];
851:   PetscBool      flg;

856:   PetscPartitionerGetDefaultType(((PetscObject) part)->type_name,&defaultType);
857:   PetscObjectOptionsBegin((PetscObject) part);
858:   PetscOptionsFList("-petscpartitioner_type", "Graph partitioner", "PetscPartitionerSetType", PetscPartitionerList, defaultType, name, sizeof(name), &flg);
859:   if (flg) {
860:     PetscPartitionerSetType(part, name);
861:   } else if (!((PetscObject) part)->type_name) {
862:     PetscPartitionerSetType(part, defaultType);
863:   }
864:   PetscOptionsBool("-petscpartitioner_use_vertex_weights","Use vertex weights","",part->usevwgt,&part->usevwgt,NULL);
865:   if (part->ops->setfromoptions) {
866:     (*part->ops->setfromoptions)(PetscOptionsObject,part);
867:   }
868:   PetscViewerDestroy(&part->viewer);
869:   PetscViewerDestroy(&part->viewerGraph);
870:   PetscOptionsGetViewer(((PetscObject) part)->comm, ((PetscObject) part)->options, ((PetscObject) part)->prefix, "-petscpartitioner_view", &part->viewer, NULL, NULL);
871:   PetscOptionsGetViewer(((PetscObject) part)->comm, ((PetscObject) part)->options, ((PetscObject) part)->prefix, "-petscpartitioner_view_graph", &part->viewerGraph, NULL, &part->viewGraph);
872:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
873:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) part);
874:   PetscOptionsEnd();
875:   return(0);
876: }

878: /*@C
879:   PetscPartitionerSetUp - Construct data structures for the PetscPartitioner

881:   Collective on PetscPartitioner

883:   Input Parameter:
884: . part - the PetscPartitioner object to setup

886:   Level: developer

888: .seealso: PetscPartitionerView(), PetscPartitionerDestroy()
889: @*/
890: PetscErrorCode PetscPartitionerSetUp(PetscPartitioner part)
891: {

896:   if (part->ops->setup) {(*part->ops->setup)(part);}
897:   return(0);
898: }

900: /*@
901:   PetscPartitionerDestroy - Destroys a PetscPartitioner object

903:   Collective on PetscPartitioner

905:   Input Parameter:
906: . part - the PetscPartitioner object to destroy

908:   Level: developer

910: .seealso: PetscPartitionerView()
911: @*/
912: PetscErrorCode PetscPartitionerDestroy(PetscPartitioner *part)
913: {

917:   if (!*part) return(0);

920:   if (--((PetscObject)(*part))->refct > 0) {*part = 0; return(0);}
921:   ((PetscObject) (*part))->refct = 0;

923:   PetscViewerDestroy(&(*part)->viewer);
924:   PetscViewerDestroy(&(*part)->viewerGraph);
925:   if ((*part)->ops->destroy) {(*(*part)->ops->destroy)(*part);}
926:   PetscHeaderDestroy(part);
927:   return(0);
928: }

930: /*@
931:   PetscPartitionerPartition - Partition a graph

933:   Collective on PetscPartitioner

935:   Input Parameters:
936: + part    - The PetscPartitioner
937: . nparts  - Number of partitions
938: . numVertices - Number of vertices in the local part of the graph
939: . start - row pointers for the local part of the graph (CSR style)
940: . adjacency - adjacency list (CSR style)
941: . vertexSection - PetscSection describing the absolute weight of each local vertex (can be NULL)
942: - targetSection - PetscSection describing the absolute weight of each partition (can be NULL)

944:   Output Parameters:
945: + partSection     - The PetscSection giving the division of points by partition
946: - partition       - The list of points by partition

948:   Options Database:
949: + -petscpartitioner_view - View the partitioner information
950: - -petscpartitioner_view_graph - View the graph we are partitioning

952:   Notes:
953:     The chart of the vertexSection (if present) must contain [0,numVertices), with the number of dofs in the section specifying the absolute weight for each vertex.
954:     The chart of the targetSection (if present) must contain [0,nparts), with the number of dofs in the section specifying the absolute weight for each partition. This information must be the same across processes, PETSc does not check it.

956:   Level: developer

958: .seealso PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscSectionSetDof()
959: @*/
960: PetscErrorCode PetscPartitionerPartition(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertexSection, PetscSection targetSection, PetscSection partSection, IS *partition)
961: {

967:   if (nparts <= 0) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Number of parts must be positive");
968:   if (numVertices < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices must be non-negative");
969:   if (numVertices && !part->noGraph) {
973:   }
974:   if (vertexSection) {
975:     PetscInt s,e;

978:     PetscSectionGetChart(vertexSection, &s, &e);
979:     if (s > 0 || e < numVertices) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid vertexSection chart [%D,%D)",s,e);
980:   }
981:   if (targetSection) {
982:     PetscInt s,e;

985:     PetscSectionGetChart(targetSection, &s, &e);
986:     if (s > 0 || e < nparts) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid targetSection chart [%D,%D)",s,e);
987:   }

991:   PetscSectionReset(partSection);
992:   PetscSectionSetChart(partSection, 0, nparts);
993:   if (nparts == 1) { /* quick */
994:     PetscSectionSetDof(partSection, 0, numVertices);
995:     ISCreateStride(PetscObjectComm((PetscObject)part),numVertices,0,1,partition);
996:   } else {
997:     if (!part->ops->partition) SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "PetscPartitioner %s has no partitioning method", ((PetscObject)part)->type_name);
998:     (*part->ops->partition)(part, nparts, numVertices, start, adjacency, vertexSection, targetSection, partSection, partition);
999:   }
1000:   PetscSectionSetUp(partSection);
1001:   if (part->viewerGraph) {
1002:     PetscViewer viewer = part->viewerGraph;
1003:     PetscBool   isascii;
1004:     PetscInt    v, i;
1005:     PetscMPIInt rank;

1007:     MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
1008:     PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1009:     if (isascii) {
1010:       PetscViewerASCIIPushSynchronized(viewer);
1011:       PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Nv: %D\n", rank, numVertices);
1012:       for (v = 0; v < numVertices; ++v) {
1013:         const PetscInt s = start[v];
1014:         const PetscInt e = start[v+1];

1016:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]  ", rank);
1017:         for (i = s; i < e; ++i) {PetscViewerASCIISynchronizedPrintf(viewer, "%D ", adjacency[i]);}
1018:         PetscViewerASCIISynchronizedPrintf(viewer, "[%D-%D)\n", s, e);
1019:       }
1020:       PetscViewerFlush(viewer);
1021:       PetscViewerASCIIPopSynchronized(viewer);
1022:     }
1023:   }
1024:   if (part->viewer) {
1025:     PetscPartitionerView(part,part->viewer);
1026:   }
1027:   return(0);
1028: }

1030: /*@
1031:   PetscPartitionerCreate - Creates an empty PetscPartitioner object. The type can then be set with PetscPartitionerSetType().

1033:   Collective

1035:   Input Parameter:
1036: . comm - The communicator for the PetscPartitioner object

1038:   Output Parameter:
1039: . part - The PetscPartitioner object

1041:   Level: beginner

1043: .seealso: PetscPartitionerSetType(), PETSCPARTITIONERCHACO, PETSCPARTITIONERPARMETIS, PETSCPARTITIONERSHELL, PETSCPARTITIONERSIMPLE, PETSCPARTITIONERGATHER
1044: @*/
1045: PetscErrorCode PetscPartitionerCreate(MPI_Comm comm, PetscPartitioner *part)
1046: {
1047:   PetscPartitioner p;
1048:   const char       *partitionerType = NULL;
1049:   PetscErrorCode   ierr;

1053:   *part = NULL;
1054:   DMInitializePackage();

1056:   PetscHeaderCreate(p, PETSCPARTITIONER_CLASSID, "PetscPartitioner", "Graph Partitioner", "PetscPartitioner", comm, PetscPartitionerDestroy, PetscPartitionerView);
1057:   PetscPartitionerGetDefaultType(NULL,&partitionerType);
1058:   PetscPartitionerSetType(p,partitionerType);

1060:   p->edgeCut = 0;
1061:   p->balance = 0.0;
1062:   p->usevwgt = PETSC_TRUE;

1064:   *part = p;
1065:   return(0);
1066: }

1068: /*@
1069:   PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh

1071:   Collective on PetscPartitioner

1073:   Input Parameters:
1074: + part    - The PetscPartitioner
1075: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
1076: - dm      - The mesh DM

1078:   Output Parameters:
1079: + partSection     - The PetscSection giving the division of points by partition
1080: - partition       - The list of points by partition

1082:   Notes:
1083:     If the DM has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
1084:     by the section in the transitive closure of the point.

1086:   Level: developer

1088: .seealso DMPlexDistribute(), PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscPartitionerPartition()
1089: @*/
1090: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
1091: {
1092:   PetscMPIInt    size;
1093:   PetscBool      isplex;
1095:   PetscSection   vertSection = NULL;

1103:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isplex);
1104:   if (!isplex) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not for type %s",((PetscObject)dm)->type_name);
1105:   MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
1106:   if (size == 1) {
1107:     PetscInt *points;
1108:     PetscInt  cStart, cEnd, c;

1110:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
1111:     PetscSectionReset(partSection);
1112:     PetscSectionSetChart(partSection, 0, size);
1113:     PetscSectionSetDof(partSection, 0, cEnd-cStart);
1114:     PetscSectionSetUp(partSection);
1115:     PetscMalloc1(cEnd-cStart, &points);
1116:     for (c = cStart; c < cEnd; ++c) points[c] = c;
1117:     ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
1118:     return(0);
1119:   }
1120:   if (part->height == 0) {
1121:     PetscInt numVertices = 0;
1122:     PetscInt *start     = NULL;
1123:     PetscInt *adjacency = NULL;
1124:     IS       globalNumbering;

1126:     if (!part->noGraph || part->viewGraph) {
1127:       DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
1128:     } else { /* only compute the number of owned local vertices */
1129:       const PetscInt *idxs;
1130:       PetscInt       p, pStart, pEnd;

1132:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
1133:       DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
1134:       ISGetIndices(globalNumbering, &idxs);
1135:       for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
1136:       ISRestoreIndices(globalNumbering, &idxs);
1137:     }
1138:     if (part->usevwgt) {
1139:       PetscSection   section = dm->localSection, clSection = NULL;
1140:       IS             clPoints = NULL;
1141:       const PetscInt *gid,*clIdx;
1142:       PetscInt       v, p, pStart, pEnd;

1144:       /* dm->localSection encodes degrees of freedom per point, not per cell. We need to get the closure index to properly specify cell weights (aka dofs) */
1145:       /* We do this only if the local section has been set */
1146:       if (section) {
1147:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
1148:         if (!clSection) {
1149:           DMPlexCreateClosureIndex(dm,NULL);
1150:         }
1151:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
1152:         ISGetIndices(clPoints,&clIdx);
1153:       }
1154:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
1155:       PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
1156:       PetscSectionSetChart(vertSection, 0, numVertices);
1157:       if (globalNumbering) {
1158:         ISGetIndices(globalNumbering,&gid);
1159:       } else gid = NULL;
1160:       for (p = pStart, v = 0; p < pEnd; ++p) {
1161:         PetscInt dof = 1;

1163:         /* skip cells in the overlap */
1164:         if (gid && gid[p-pStart] < 0) continue;

1166:         if (section) {
1167:           PetscInt cl, clSize, clOff;

1169:           dof  = 0;
1170:           PetscSectionGetDof(clSection, p, &clSize);
1171:           PetscSectionGetOffset(clSection, p, &clOff);
1172:           for (cl = 0; cl < clSize; cl+=2) {
1173:             PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */

1175:             PetscSectionGetDof(section, clPoint, &clDof);
1176:             dof += clDof;
1177:           }
1178:         }
1179:         if (!dof) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of dofs for point %D in the local section should be positive",p);
1180:         PetscSectionSetDof(vertSection, v, dof);
1181:         v++;
1182:       }
1183:       if (globalNumbering) {
1184:         ISRestoreIndices(globalNumbering,&gid);
1185:       }
1186:       if (clPoints) {
1187:         ISRestoreIndices(clPoints,&clIdx);
1188:       }
1189:       PetscSectionSetUp(vertSection);
1190:     }
1191:     PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
1192:     PetscFree(start);
1193:     PetscFree(adjacency);
1194:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
1195:       const PetscInt *globalNum;
1196:       const PetscInt *partIdx;
1197:       PetscInt       *map, cStart, cEnd;
1198:       PetscInt       *adjusted, i, localSize, offset;
1199:       IS             newPartition;

1201:       ISGetLocalSize(*partition,&localSize);
1202:       PetscMalloc1(localSize,&adjusted);
1203:       ISGetIndices(globalNumbering,&globalNum);
1204:       ISGetIndices(*partition,&partIdx);
1205:       PetscMalloc1(localSize,&map);
1206:       DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
1207:       for (i = cStart, offset = 0; i < cEnd; i++) {
1208:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
1209:       }
1210:       for (i = 0; i < localSize; i++) {
1211:         adjusted[i] = map[partIdx[i]];
1212:       }
1213:       PetscFree(map);
1214:       ISRestoreIndices(*partition,&partIdx);
1215:       ISRestoreIndices(globalNumbering,&globalNum);
1216:       ISCreateGeneral(PETSC_COMM_SELF,localSize,adjusted,PETSC_OWN_POINTER,&newPartition);
1217:       ISDestroy(&globalNumbering);
1218:       ISDestroy(partition);
1219:       *partition = newPartition;
1220:     }
1221:   } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
1222:   PetscSectionDestroy(&vertSection);
1223:   return(0);
1224: }

1226: static PetscErrorCode PetscPartitionerDestroy_Shell(PetscPartitioner part)
1227: {
1228:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1229:   PetscErrorCode          ierr;

1232:   PetscSectionDestroy(&p->section);
1233:   ISDestroy(&p->partition);
1234:   PetscFree(p);
1235:   return(0);
1236: }

1238: static PetscErrorCode PetscPartitionerView_Shell_Ascii(PetscPartitioner part, PetscViewer viewer)
1239: {
1240:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1241:   PetscErrorCode          ierr;

1244:   if (p->random) {
1245:     PetscViewerASCIIPushTab(viewer);
1246:     PetscViewerASCIIPrintf(viewer, "using random partition\n");
1247:     PetscViewerASCIIPopTab(viewer);
1248:   }
1249:   return(0);
1250: }

1252: static PetscErrorCode PetscPartitionerView_Shell(PetscPartitioner part, PetscViewer viewer)
1253: {
1254:   PetscBool      iascii;

1260:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1261:   if (iascii) {PetscPartitionerView_Shell_Ascii(part, viewer);}
1262:   return(0);
1263: }

1265: static PetscErrorCode PetscPartitionerSetFromOptions_Shell(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1266: {
1267:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1268:   PetscErrorCode          ierr;

1271:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner Shell Options");
1272:   PetscOptionsBool("-petscpartitioner_shell_random", "Use a random partition", "PetscPartitionerView", PETSC_FALSE, &p->random, NULL);
1273:   PetscOptionsTail();
1274:   return(0);
1275: }

1277: static PetscErrorCode PetscPartitionerPartition_Shell(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1278: {
1279:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1280:   PetscInt                np;
1281:   PetscErrorCode          ierr;

1284:   if (p->random) {
1285:     PetscRandom r;
1286:     PetscInt   *sizes, *points, v, p;
1287:     PetscMPIInt rank;

1289:     MPI_Comm_rank(PetscObjectComm((PetscObject) part), &rank);
1290:     PetscRandomCreate(PETSC_COMM_SELF, &r);
1291:     PetscRandomSetInterval(r, 0.0, (PetscScalar) nparts);
1292:     PetscRandomSetFromOptions(r);
1293:     PetscCalloc2(nparts, &sizes, numVertices, &points);
1294:     for (v = 0; v < numVertices; ++v) {points[v] = v;}
1295:     for (p = 0; p < nparts; ++p) {sizes[p] = numVertices/nparts + (PetscInt) (p < numVertices % nparts);}
1296:     for (v = numVertices-1; v > 0; --v) {
1297:       PetscReal val;
1298:       PetscInt  w, tmp;

1300:       PetscRandomSetInterval(r, 0.0, (PetscScalar) (v+1));
1301:       PetscRandomGetValueReal(r, &val);
1302:       w    = PetscFloorReal(val);
1303:       tmp       = points[v];
1304:       points[v] = points[w];
1305:       points[w] = tmp;
1306:     }
1307:     PetscRandomDestroy(&r);
1308:     PetscPartitionerShellSetPartition(part, nparts, sizes, points);
1309:     PetscFree2(sizes, points);
1310:   }
1311:   if (!p->section) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_WRONG, "Shell partitioner information not provided. Please call PetscPartitionerShellSetPartition()");
1312:   PetscSectionGetChart(p->section, NULL, &np);
1313:   if (nparts != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of requested partitions %d != configured partitions %d", nparts, np);
1314:   ISGetLocalSize(p->partition, &np);
1315:   if (numVertices != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input vertices %d != configured vertices %d", numVertices, np);
1316:   PetscSectionCopy(p->section, partSection);
1317:   *partition = p->partition;
1318:   PetscObjectReference((PetscObject) p->partition);
1319:   return(0);
1320: }

1322: static PetscErrorCode PetscPartitionerInitialize_Shell(PetscPartitioner part)
1323: {
1325:   part->noGraph             = PETSC_TRUE; /* PetscPartitionerShell cannot overload the partition call, so it is safe for now */
1326:   part->ops->view           = PetscPartitionerView_Shell;
1327:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_Shell;
1328:   part->ops->destroy        = PetscPartitionerDestroy_Shell;
1329:   part->ops->partition      = PetscPartitionerPartition_Shell;
1330:   return(0);
1331: }

1333: /*MC
1334:   PETSCPARTITIONERSHELL = "shell" - A PetscPartitioner object

1336:   Level: intermediate

1338:   Options Database Keys:
1339: .  -petscpartitioner_shell_random - Use a random partition

1341: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1342: M*/

1344: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Shell(PetscPartitioner part)
1345: {
1346:   PetscPartitioner_Shell *p;
1347:   PetscErrorCode          ierr;

1351:   PetscNewLog(part, &p);
1352:   part->data = p;

1354:   PetscPartitionerInitialize_Shell(part);
1355:   p->random = PETSC_FALSE;
1356:   return(0);
1357: }

1359: /*@C
1360:   PetscPartitionerShellSetPartition - Set an artifical partition for a mesh

1362:   Collective on PetscPartitioner

1364:   Input Parameters:
1365: + part   - The PetscPartitioner
1366: . size   - The number of partitions
1367: . sizes  - array of length size (or NULL) providing the number of points in each partition
1368: - points - array of length sum(sizes) (may be NULL iff sizes is NULL), a permutation of the points that groups those assigned to each partition in order (i.e., partition 0 first, partition 1 next, etc.)

1370:   Level: developer

1372:   Notes:
1373:     It is safe to free the sizes and points arrays after use in this routine.

1375: .seealso DMPlexDistribute(), PetscPartitionerCreate()
1376: @*/
1377: PetscErrorCode PetscPartitionerShellSetPartition(PetscPartitioner part, PetscInt size, const PetscInt sizes[], const PetscInt points[])
1378: {
1379:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1380:   PetscInt                proc, numPoints;
1381:   PetscErrorCode          ierr;

1387:   PetscSectionDestroy(&p->section);
1388:   ISDestroy(&p->partition);
1389:   PetscSectionCreate(PetscObjectComm((PetscObject) part), &p->section);
1390:   PetscSectionSetChart(p->section, 0, size);
1391:   if (sizes) {
1392:     for (proc = 0; proc < size; ++proc) {
1393:       PetscSectionSetDof(p->section, proc, sizes[proc]);
1394:     }
1395:   }
1396:   PetscSectionSetUp(p->section);
1397:   PetscSectionGetStorageSize(p->section, &numPoints);
1398:   ISCreateGeneral(PetscObjectComm((PetscObject) part), numPoints, points, PETSC_COPY_VALUES, &p->partition);
1399:   return(0);
1400: }

1402: /*@
1403:   PetscPartitionerShellSetRandom - Set the flag to use a random partition

1405:   Collective on PetscPartitioner

1407:   Input Parameters:
1408: + part   - The PetscPartitioner
1409: - random - The flag to use a random partition

1411:   Level: intermediate

1413: .seealso PetscPartitionerShellGetRandom(), PetscPartitionerCreate()
1414: @*/
1415: PetscErrorCode PetscPartitionerShellSetRandom(PetscPartitioner part, PetscBool random)
1416: {
1417:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;

1421:   p->random = random;
1422:   return(0);
1423: }

1425: /*@
1426:   PetscPartitionerShellGetRandom - get the flag to use a random partition

1428:   Collective on PetscPartitioner

1430:   Input Parameter:
1431: . part   - The PetscPartitioner

1433:   Output Parameter:
1434: . random - The flag to use a random partition

1436:   Level: intermediate

1438: .seealso PetscPartitionerShellSetRandom(), PetscPartitionerCreate()
1439: @*/
1440: PetscErrorCode PetscPartitionerShellGetRandom(PetscPartitioner part, PetscBool *random)
1441: {
1442:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;

1447:   *random = p->random;
1448:   return(0);
1449: }

1451: static PetscErrorCode PetscPartitionerDestroy_Simple(PetscPartitioner part)
1452: {
1453:   PetscPartitioner_Simple *p = (PetscPartitioner_Simple *) part->data;
1454:   PetscErrorCode          ierr;

1457:   PetscFree(p);
1458:   return(0);
1459: }

1461: static PetscErrorCode PetscPartitionerView_Simple_Ascii(PetscPartitioner part, PetscViewer viewer)
1462: {
1464:   return(0);
1465: }

1467: static PetscErrorCode PetscPartitionerView_Simple(PetscPartitioner part, PetscViewer viewer)
1468: {
1469:   PetscBool      iascii;

1475:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1476:   if (iascii) {PetscPartitionerView_Simple_Ascii(part, viewer);}
1477:   return(0);
1478: }

1480: static PetscErrorCode PetscPartitionerPartition_Simple(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1481: {
1482:   MPI_Comm       comm;
1483:   PetscInt       np, *tpwgts = NULL, sumw = 0, numVerticesGlobal  = 0;
1484:   PetscMPIInt    size;

1488:   if (vertSection) { PetscInfo(part,"PETSCPARTITIONERSIMPLE ignores vertex weights\n"); }
1489:   comm = PetscObjectComm((PetscObject)part);
1490:   MPI_Comm_size(comm,&size);
1491:   if (targetSection) {
1492:     MPIU_Allreduce(&numVertices, &numVerticesGlobal, 1, MPIU_INT, MPI_SUM, comm);
1493:     PetscCalloc1(nparts,&tpwgts);
1494:     for (np = 0; np < nparts; ++np) {
1495:       PetscSectionGetDof(targetSection,np,&tpwgts[np]);
1496:       sumw += tpwgts[np];
1497:     }
1498:     if (!sumw) {
1499:       PetscFree(tpwgts);
1500:     } else {
1501:       PetscInt m,mp;
1502:       for (np = 0; np < nparts; ++np) tpwgts[np] = (tpwgts[np]*numVerticesGlobal)/sumw;
1503:       for (np = 0, m = -1, mp = 0, sumw = 0; np < nparts; ++np) {
1504:         if (m < tpwgts[np]) { m = tpwgts[np]; mp = np; }
1505:         sumw += tpwgts[np];
1506:       }
1507:       if (sumw != numVerticesGlobal) tpwgts[mp] += numVerticesGlobal - sumw;
1508:     }
1509:   }

1511:   ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1512:   if (size == 1) {
1513:     if (tpwgts) {
1514:       for (np = 0; np < nparts; ++np) {
1515:         PetscSectionSetDof(partSection, np, tpwgts[np]);
1516:       }
1517:     } else {
1518:       for (np = 0; np < nparts; ++np) {
1519:         PetscSectionSetDof(partSection, np, numVertices/nparts + ((numVertices % nparts) > np));
1520:       }
1521:     }
1522:   } else {
1523:     if (tpwgts) {
1524:       Vec         v;
1525:       PetscScalar *array;
1526:       PetscInt    st,j;
1527:       PetscMPIInt rank;

1529:       VecCreate(comm,&v);
1530:       VecSetSizes(v,numVertices,numVerticesGlobal);
1531:       VecSetType(v,VECSTANDARD);
1532:       MPI_Comm_rank(comm,&rank);
1533:       for (np = 0,st = 0; np < nparts; ++np) {
1534:         if (rank == np || (rank == size-1 && size < nparts && np >= size)) {
1535:           for (j = 0; j < tpwgts[np]; j++) {
1536:             VecSetValue(v,st+j,np,INSERT_VALUES);
1537:           }
1538:         }
1539:         st += tpwgts[np];
1540:       }
1541:       VecAssemblyBegin(v);
1542:       VecAssemblyEnd(v);
1543:       VecGetArray(v,&array);
1544:       for (j = 0; j < numVertices; ++j) {
1545:         PetscSectionAddDof(partSection,PetscRealPart(array[j]),1);
1546:       }
1547:       VecRestoreArray(v,&array);
1548:       VecDestroy(&v);
1549:     } else {
1550:       PetscMPIInt rank;
1551:       PetscInt nvGlobal, *offsets, myFirst, myLast;

1553:       PetscMalloc1(size+1,&offsets);
1554:       offsets[0] = 0;
1555:       MPI_Allgather(&numVertices,1,MPIU_INT,&offsets[1],1,MPIU_INT,comm);
1556:       for (np = 2; np <= size; np++) {
1557:         offsets[np] += offsets[np-1];
1558:       }
1559:       nvGlobal = offsets[size];
1560:       MPI_Comm_rank(comm,&rank);
1561:       myFirst = offsets[rank];
1562:       myLast  = offsets[rank + 1] - 1;
1563:       PetscFree(offsets);
1564:       if (numVertices) {
1565:         PetscInt firstPart = 0, firstLargePart = 0;
1566:         PetscInt lastPart = 0, lastLargePart = 0;
1567:         PetscInt rem = nvGlobal % nparts;
1568:         PetscInt pSmall = nvGlobal/nparts;
1569:         PetscInt pBig = nvGlobal/nparts + 1;

1571:         if (rem) {
1572:           firstLargePart = myFirst / pBig;
1573:           lastLargePart  = myLast  / pBig;

1575:           if (firstLargePart < rem) {
1576:             firstPart = firstLargePart;
1577:           } else {
1578:             firstPart = rem + (myFirst - (rem * pBig)) / pSmall;
1579:           }
1580:           if (lastLargePart < rem) {
1581:             lastPart = lastLargePart;
1582:           } else {
1583:             lastPart = rem + (myLast - (rem * pBig)) / pSmall;
1584:           }
1585:         } else {
1586:           firstPart = myFirst / (nvGlobal/nparts);
1587:           lastPart  = myLast  / (nvGlobal/nparts);
1588:         }

1590:         for (np = firstPart; np <= lastPart; np++) {
1591:           PetscInt PartStart =  np    * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np);
1592:           PetscInt PartEnd   = (np+1) * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np+1);

1594:           PartStart = PetscMax(PartStart,myFirst);
1595:           PartEnd   = PetscMin(PartEnd,myLast+1);
1596:           PetscSectionSetDof(partSection,np,PartEnd-PartStart);
1597:         }
1598:       }
1599:     }
1600:   }
1601:   PetscFree(tpwgts);
1602:   return(0);
1603: }

1605: static PetscErrorCode PetscPartitionerInitialize_Simple(PetscPartitioner part)
1606: {
1608:   part->noGraph        = PETSC_TRUE;
1609:   part->ops->view      = PetscPartitionerView_Simple;
1610:   part->ops->destroy   = PetscPartitionerDestroy_Simple;
1611:   part->ops->partition = PetscPartitionerPartition_Simple;
1612:   return(0);
1613: }

1615: /*MC
1616:   PETSCPARTITIONERSIMPLE = "simple" - A PetscPartitioner object

1618:   Level: intermediate

1620: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1621: M*/

1623: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Simple(PetscPartitioner part)
1624: {
1625:   PetscPartitioner_Simple *p;
1626:   PetscErrorCode           ierr;

1630:   PetscNewLog(part, &p);
1631:   part->data = p;

1633:   PetscPartitionerInitialize_Simple(part);
1634:   return(0);
1635: }

1637: static PetscErrorCode PetscPartitionerDestroy_Gather(PetscPartitioner part)
1638: {
1639:   PetscPartitioner_Gather *p = (PetscPartitioner_Gather *) part->data;
1640:   PetscErrorCode          ierr;

1643:   PetscFree(p);
1644:   return(0);
1645: }

1647: static PetscErrorCode PetscPartitionerView_Gather_Ascii(PetscPartitioner part, PetscViewer viewer)
1648: {
1650:   return(0);
1651: }

1653: static PetscErrorCode PetscPartitionerView_Gather(PetscPartitioner part, PetscViewer viewer)
1654: {
1655:   PetscBool      iascii;

1661:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1662:   if (iascii) {PetscPartitionerView_Gather_Ascii(part, viewer);}
1663:   return(0);
1664: }

1666: static PetscErrorCode PetscPartitionerPartition_Gather(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1667: {
1668:   PetscInt       np;

1672:   ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1673:   PetscSectionSetDof(partSection,0,numVertices);
1674:   for (np = 1; np < nparts; ++np) {PetscSectionSetDof(partSection, np, 0);}
1675:   return(0);
1676: }

1678: static PetscErrorCode PetscPartitionerInitialize_Gather(PetscPartitioner part)
1679: {
1681:   part->noGraph        = PETSC_TRUE;
1682:   part->ops->view      = PetscPartitionerView_Gather;
1683:   part->ops->destroy   = PetscPartitionerDestroy_Gather;
1684:   part->ops->partition = PetscPartitionerPartition_Gather;
1685:   return(0);
1686: }

1688: /*MC
1689:   PETSCPARTITIONERGATHER = "gather" - A PetscPartitioner object

1691:   Level: intermediate

1693: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1694: M*/

1696: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Gather(PetscPartitioner part)
1697: {
1698:   PetscPartitioner_Gather *p;
1699:   PetscErrorCode           ierr;

1703:   PetscNewLog(part, &p);
1704:   part->data = p;

1706:   PetscPartitionerInitialize_Gather(part);
1707:   return(0);
1708: }


1711: static PetscErrorCode PetscPartitionerDestroy_Chaco(PetscPartitioner part)
1712: {
1713:   PetscPartitioner_Chaco *p = (PetscPartitioner_Chaco *) part->data;
1714:   PetscErrorCode          ierr;

1717:   PetscFree(p);
1718:   return(0);
1719: }

1721: static PetscErrorCode PetscPartitionerView_Chaco_Ascii(PetscPartitioner part, PetscViewer viewer)
1722: {
1724:   return(0);
1725: }

1727: static PetscErrorCode PetscPartitionerView_Chaco(PetscPartitioner part, PetscViewer viewer)
1728: {
1729:   PetscBool      iascii;

1735:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1736:   if (iascii) {PetscPartitionerView_Chaco_Ascii(part, viewer);}
1737:   return(0);
1738: }

1740: #if defined(PETSC_HAVE_CHACO)
1741: #if defined(PETSC_HAVE_UNISTD_H)
1742: #include <unistd.h>
1743: #endif
1744: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1745: #include <chaco.h>
1746: #else
1747: /* Older versions of Chaco do not have an include file */
1748: PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
1749:                        float *ewgts, float *x, float *y, float *z, char *outassignname,
1750:                        char *outfilename, short *assignment, int architecture, int ndims_tot,
1751:                        int mesh_dims[3], double *goal, int global_method, int local_method,
1752:                        int rqi_flag, int vmax, int ndims, double eigtol, long seed);
1753: #endif
1754: extern int FREE_GRAPH;
1755: #endif

1757: static PetscErrorCode PetscPartitionerPartition_Chaco(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1758: {
1759: #if defined(PETSC_HAVE_CHACO)
1760:   enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
1761:   MPI_Comm       comm;
1762:   int            nvtxs          = numVertices; /* number of vertices in full graph */
1763:   int           *vwgts          = NULL;   /* weights for all vertices */
1764:   float         *ewgts          = NULL;   /* weights for all edges */
1765:   float         *x              = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
1766:   char          *outassignname  = NULL;   /*  name of assignment output file */
1767:   char          *outfilename    = NULL;   /* output file name */
1768:   int            architecture   = 1;      /* 0 => hypercube, d => d-dimensional mesh */
1769:   int            ndims_tot      = 0;      /* total number of cube dimensions to divide */
1770:   int            mesh_dims[3];            /* dimensions of mesh of processors */
1771:   double        *goal          = NULL;    /* desired set sizes for each set */
1772:   int            global_method = 1;       /* global partitioning algorithm */
1773:   int            local_method  = 1;       /* local partitioning algorithm */
1774:   int            rqi_flag      = 0;       /* should I use RQI/Symmlq eigensolver? */
1775:   int            vmax          = 200;     /* how many vertices to coarsen down to? */
1776:   int            ndims         = 1;       /* number of eigenvectors (2^d sets) */
1777:   double         eigtol        = 0.001;   /* tolerance on eigenvectors */
1778:   long           seed          = 123636512; /* for random graph mutations */
1779: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1780:   int           *assignment;              /* Output partition */
1781: #else
1782:   short int     *assignment;              /* Output partition */
1783: #endif
1784:   int            fd_stdout, fd_pipe[2];
1785:   PetscInt      *points;
1786:   int            i, v, p;

1790:   PetscObjectGetComm((PetscObject)part,&comm);
1791: #if defined (PETSC_USE_DEBUG)
1792:   {
1793:     int ival,isum;
1794:     PetscBool distributed;

1796:     ival = (numVertices > 0);
1797:     MPI_Allreduce(&ival, &isum, 1, MPI_INT, MPI_SUM, comm);
1798:     distributed = (isum > 1) ? PETSC_TRUE : PETSC_FALSE;
1799:     if (distributed) SETERRQ(comm, PETSC_ERR_SUP, "Chaco cannot partition a distributed graph");
1800:   }
1801: #endif
1802:   if (!numVertices) { /* distributed case, return if not holding the graph */
1803:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
1804:     return(0);
1805:   }
1806:   FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
1807:   for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];

1809:   if (global_method == INERTIAL_METHOD) {
1810:     /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
1811:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
1812:   }
1813:   mesh_dims[0] = nparts;
1814:   mesh_dims[1] = 1;
1815:   mesh_dims[2] = 1;
1816:   PetscMalloc1(nvtxs, &assignment);
1817:   /* Chaco outputs to stdout. We redirect this to a buffer. */
1818:   /* TODO: check error codes for UNIX calls */
1819: #if defined(PETSC_HAVE_UNISTD_H)
1820:   {
1821:     int piperet;
1822:     piperet = pipe(fd_pipe);
1823:     if (piperet) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"Could not create pipe");
1824:     fd_stdout = dup(1);
1825:     close(1);
1826:     dup2(fd_pipe[1], 1);
1827:   }
1828: #endif
1829:   if (part->usevwgt) { PetscInfo(part,"PETSCPARTITIONERCHACO ignores vertex weights\n"); }
1830:   interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
1831:                    assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
1832:                    vmax, ndims, eigtol, seed);
1833: #if defined(PETSC_HAVE_UNISTD_H)
1834:   {
1835:     char msgLog[10000];
1836:     int  count;

1838:     fflush(stdout);
1839:     count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
1840:     if (count < 0) count = 0;
1841:     msgLog[count] = 0;
1842:     close(1);
1843:     dup2(fd_stdout, 1);
1844:     close(fd_stdout);
1845:     close(fd_pipe[0]);
1846:     close(fd_pipe[1]);
1847:     if (ierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
1848:   }
1849: #else
1850:   if (ierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in Chaco library: %s", "error in stdout");
1851: #endif
1852:   /* Convert to PetscSection+IS */
1853:   for (v = 0; v < nvtxs; ++v) {
1854:     PetscSectionAddDof(partSection, assignment[v], 1);
1855:   }
1856:   PetscMalloc1(nvtxs, &points);
1857:   for (p = 0, i = 0; p < nparts; ++p) {
1858:     for (v = 0; v < nvtxs; ++v) {
1859:       if (assignment[v] == p) points[i++] = v;
1860:     }
1861:   }
1862:   if (i != nvtxs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1863:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
1864:   if (global_method == INERTIAL_METHOD) {
1865:     /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
1866:   }
1867:   PetscFree(assignment);
1868:   for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
1869:   return(0);
1870: #else
1871:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco.");
1872: #endif
1873: }

1875: static PetscErrorCode PetscPartitionerInitialize_Chaco(PetscPartitioner part)
1876: {
1878:   part->noGraph        = PETSC_FALSE;
1879:   part->ops->view      = PetscPartitionerView_Chaco;
1880:   part->ops->destroy   = PetscPartitionerDestroy_Chaco;
1881:   part->ops->partition = PetscPartitionerPartition_Chaco;
1882:   return(0);
1883: }

1885: /*MC
1886:   PETSCPARTITIONERCHACO = "chaco" - A PetscPartitioner object using the Chaco library

1888:   Level: intermediate

1890: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1891: M*/

1893: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Chaco(PetscPartitioner part)
1894: {
1895:   PetscPartitioner_Chaco *p;
1896:   PetscErrorCode          ierr;

1900:   PetscNewLog(part, &p);
1901:   part->data = p;

1903:   PetscPartitionerInitialize_Chaco(part);
1904:   PetscCitationsRegister(ChacoPartitionerCitation, &ChacoPartitionercite);
1905:   return(0);
1906: }

1908: static const char *ptypes[] = {"kway", "rb"};

1910: static PetscErrorCode PetscPartitionerDestroy_ParMetis(PetscPartitioner part)
1911: {
1912:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1913:   PetscErrorCode             ierr;

1916:   MPI_Comm_free(&p->pcomm);
1917:   PetscFree(p);
1918:   return(0);
1919: }

1921: static PetscErrorCode PetscPartitionerView_ParMetis_Ascii(PetscPartitioner part, PetscViewer viewer)
1922: {
1923:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1924:   PetscErrorCode             ierr;

1927:   PetscViewerASCIIPushTab(viewer);
1928:   PetscViewerASCIIPrintf(viewer, "ParMetis type: %s\n", ptypes[p->ptype]);
1929:   PetscViewerASCIIPrintf(viewer, "load imbalance ratio %g\n", (double) p->imbalanceRatio);
1930:   PetscViewerASCIIPrintf(viewer, "debug flag %D\n", p->debugFlag);
1931:   PetscViewerASCIIPrintf(viewer, "random seed %D\n", p->randomSeed);
1932:   PetscViewerASCIIPopTab(viewer);
1933:   return(0);
1934: }

1936: static PetscErrorCode PetscPartitionerView_ParMetis(PetscPartitioner part, PetscViewer viewer)
1937: {
1938:   PetscBool      iascii;

1944:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1945:   if (iascii) {PetscPartitionerView_ParMetis_Ascii(part, viewer);}
1946:   return(0);
1947: }

1949: static PetscErrorCode PetscPartitionerSetFromOptions_ParMetis(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1950: {
1951:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1952:   PetscErrorCode            ierr;

1955:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner ParMetis Options");
1956:   PetscOptionsEList("-petscpartitioner_parmetis_type", "Partitioning method", "", ptypes, 2, ptypes[p->ptype], &p->ptype, NULL);
1957:   PetscOptionsReal("-petscpartitioner_parmetis_imbalance_ratio", "Load imbalance ratio limit", "", p->imbalanceRatio, &p->imbalanceRatio, NULL);
1958:   PetscOptionsInt("-petscpartitioner_parmetis_debug", "Debugging flag", "", p->debugFlag, &p->debugFlag, NULL);
1959:   PetscOptionsInt("-petscpartitioner_parmetis_seed", "Random seed", "", p->randomSeed, &p->randomSeed, NULL);
1960:   PetscOptionsTail();
1961:   return(0);
1962: }

1964: #if defined(PETSC_HAVE_PARMETIS)
1965: #include <parmetis.h>
1966: #endif

1968: static PetscErrorCode PetscPartitionerPartition_ParMetis(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1969: {
1970: #if defined(PETSC_HAVE_PARMETIS)
1971:   PetscPartitioner_ParMetis *pm = (PetscPartitioner_ParMetis *) part->data;
1972:   MPI_Comm       comm;
1973:   PetscInt       nvtxs       = numVertices; /* The number of vertices in full graph */
1974:   PetscInt      *vtxdist;                   /* Distribution of vertices across processes */
1975:   PetscInt      *xadj        = start;       /* Start of edge list for each vertex */
1976:   PetscInt      *adjncy      = adjacency;   /* Edge lists for all vertices */
1977:   PetscInt      *vwgt        = NULL;        /* Vertex weights */
1978:   PetscInt      *adjwgt      = NULL;        /* Edge weights */
1979:   PetscInt       wgtflag     = 0;           /* Indicates which weights are present */
1980:   PetscInt       numflag     = 0;           /* Indicates initial offset (0 or 1) */
1981:   PetscInt       ncon        = 1;           /* The number of weights per vertex */
1982:   PetscInt       metis_ptype = pm->ptype;   /* kway or recursive bisection */
1983:   real_t        *tpwgts;                    /* The fraction of vertex weights assigned to each partition */
1984:   real_t        *ubvec;                     /* The balance intolerance for vertex weights */
1985:   PetscInt       options[64];               /* Options */
1986:   PetscInt       v, i, *assignment, *points;
1987:   PetscMPIInt    p, size, rank;
1988:   PetscBool      hasempty = PETSC_FALSE;

1992:   PetscObjectGetComm((PetscObject) part, &comm);
1993:   MPI_Comm_size(comm, &size);
1994:   MPI_Comm_rank(comm, &rank);
1995:   /* Calculate vertex distribution */
1996:   PetscMalloc4(size+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment);
1997:   vtxdist[0] = 0;
1998:   MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
1999:   for (p = 2; p <= size; ++p) {
2000:     hasempty = (PetscBool)(hasempty || !vtxdist[p-1] || !vtxdist[p]);
2001:     vtxdist[p] += vtxdist[p-1];
2002:   }
2003:   /* null graph */
2004:   if (vtxdist[size] == 0) {
2005:     PetscFree4(vtxdist,tpwgts,ubvec,assignment);
2006:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
2007:     return(0);
2008:   }
2009:   /* Calculate partition weights */
2010:   if (targetSection) {
2011:     PetscInt p;
2012:     real_t   sumt = 0.0;

2014:     for (p = 0; p < nparts; ++p) {
2015:       PetscInt tpd;

2017:       PetscSectionGetDof(targetSection,p,&tpd);
2018:       sumt += tpd;
2019:       tpwgts[p] = tpd;
2020:     }
2021:     if (sumt) { /* METIS/ParMETIS do not like exactly zero weight */
2022:       for (p = 0, sumt = 0.0; p < nparts; ++p) {
2023:         tpwgts[p] = PetscMax(tpwgts[p],PETSC_SMALL);
2024:         sumt += tpwgts[p];
2025:       }
2026:       for (p = 0; p < nparts; ++p) tpwgts[p] /= sumt;
2027:       for (p = 0, sumt = 0.0; p < nparts-1; ++p) sumt += tpwgts[p];
2028:       tpwgts[nparts - 1] = 1. - sumt;
2029:     }
2030:   } else {
2031:     for (p = 0; p < nparts; ++p) tpwgts[p] = 1.0/nparts;
2032:   }
2033:   ubvec[0] = pm->imbalanceRatio;

2035:   /* Weight cells */
2036:   if (vertSection) {
2037:     PetscMalloc1(nvtxs,&vwgt);
2038:     for (v = 0; v < nvtxs; ++v) {
2039:       PetscSectionGetDof(vertSection, v, &vwgt[v]);
2040:     }
2041:     wgtflag |= 2; /* have weights on graph vertices */
2042:   }

2044:   for (p = 0; !vtxdist[p+1] && p < size; ++p);
2045:   if (vtxdist[p+1] == vtxdist[size]) {
2046:     if (rank == p) {
2047:       METIS_SetDefaultOptions(options); /* initialize all defaults */
2048:       options[METIS_OPTION_DBGLVL] = pm->debugFlag;
2049:       options[METIS_OPTION_SEED]   = pm->randomSeed;
2050:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
2051:       if (metis_ptype == 1) {
2052:         PetscStackPush("METIS_PartGraphRecursive");
2053:         METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
2054:         PetscStackPop;
2055:         if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphRecursive()");
2056:       } else {
2057:         /*
2058:          It would be nice to activate the two options below, but they would need some actual testing.
2059:          - Turning on these options may exercise path of the METIS code that have bugs and may break production runs.
2060:          - If CONTIG is set to 1, METIS will exit with error if the graph is disconnected, despite the manual saying the option is ignored in such case.
2061:         */
2062:         /* options[METIS_OPTION_CONTIG]  = 1; */ /* try to produce partitions that are contiguous */
2063:         /* options[METIS_OPTION_MINCONN] = 1; */ /* minimize the maximum degree of the subdomain graph */
2064:         PetscStackPush("METIS_PartGraphKway");
2065:         METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
2066:         PetscStackPop;
2067:         if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
2068:       }
2069:     }
2070:   } else {
2071:     MPI_Comm pcomm;

2073:     options[0] = 1; /*use options */
2074:     options[1] = pm->debugFlag;
2075:     options[2] = (pm->randomSeed == -1) ? 15 : pm->randomSeed; /* default is GLOBAL_SEED=15 from `libparmetis/defs.h` */

2077:     if (hasempty) { /* parmetis does not support empty graphs on some of the processes */
2078:       PetscInt cnt;

2080:       MPI_Comm_split(pm->pcomm,!!nvtxs,rank,&pcomm);
2081:       for (p=0,cnt=0;p<size;p++) {
2082:         if (vtxdist[p+1] != vtxdist[p]) {
2083:           vtxdist[cnt+1] = vtxdist[p+1];
2084:           cnt++;
2085:         }
2086:       }
2087:     } else pcomm = pm->pcomm;
2088:     if (nvtxs) {
2089:       PetscStackPush("ParMETIS_V3_PartKway");
2090:       ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment, &pcomm);
2091:       PetscStackPop;
2092:       if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error %d in ParMETIS_V3_PartKway()", ierr);
2093:     }
2094:     if (hasempty) {
2095:       MPI_Comm_free(&pcomm);
2096:     }
2097:   }

2099:   /* Convert to PetscSection+IS */
2100:   for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
2101:   PetscMalloc1(nvtxs, &points);
2102:   for (p = 0, i = 0; p < nparts; ++p) {
2103:     for (v = 0; v < nvtxs; ++v) {
2104:       if (assignment[v] == p) points[i++] = v;
2105:     }
2106:   }
2107:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
2108:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
2109:   PetscFree4(vtxdist,tpwgts,ubvec,assignment);
2110:   PetscFree(vwgt);
2111:   return(0);
2112: #else
2113:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2114: #endif
2115: }

2117: static PetscErrorCode PetscPartitionerInitialize_ParMetis(PetscPartitioner part)
2118: {
2120:   part->noGraph             = PETSC_FALSE;
2121:   part->ops->view           = PetscPartitionerView_ParMetis;
2122:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_ParMetis;
2123:   part->ops->destroy        = PetscPartitionerDestroy_ParMetis;
2124:   part->ops->partition      = PetscPartitionerPartition_ParMetis;
2125:   return(0);
2126: }

2128: /*MC
2129:   PETSCPARTITIONERPARMETIS = "parmetis" - A PetscPartitioner object using the ParMETIS library

2131:   Level: intermediate

2133:   Options Database Keys:
2134: +  -petscpartitioner_parmetis_type <string> - ParMETIS partitioning type. Either "kway" or "rb" (recursive bisection)
2135: .  -petscpartitioner_parmetis_imbalance_ratio <value> - Load imbalance ratio limit
2136: .  -petscpartitioner_parmetis_debug <int> - Debugging flag passed to ParMETIS/METIS routines
2137: -  -petscpartitioner_parmetis_seed <int> - Random seed

2139:   Notes: when the graph is on a single process, this partitioner actually calls METIS and not ParMETIS

2141: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
2142: M*/

2144: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMetis(PetscPartitioner part)
2145: {
2146:   PetscPartitioner_ParMetis *p;
2147:   PetscErrorCode          ierr;

2151:   PetscNewLog(part, &p);
2152:   part->data = p;

2154:   MPI_Comm_dup(PetscObjectComm((PetscObject)part),&p->pcomm);
2155:   p->ptype          = 0;
2156:   p->imbalanceRatio = 1.05;
2157:   p->debugFlag      = 0;
2158:   p->randomSeed     = -1; /* defaults to GLOBAL_SEED=15 from `libparmetis/defs.h` */

2160:   PetscPartitionerInitialize_ParMetis(part);
2161:   PetscCitationsRegister(ParMetisPartitionerCitation, &ParMetisPartitionercite);
2162:   return(0);
2163: }

2165: #if defined(PETSC_HAVE_PTSCOTCH)

2167: EXTERN_C_BEGIN
2168: #include <ptscotch.h>
2169: EXTERN_C_END

2171: #define CHKERRPTSCOTCH(ierr) do { if (ierr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error calling PT-Scotch library"); } while(0)

2173: static int PTScotch_Strategy(PetscInt strategy)
2174: {
2175:   switch (strategy) {
2176:   case  0: return SCOTCH_STRATDEFAULT;
2177:   case  1: return SCOTCH_STRATQUALITY;
2178:   case  2: return SCOTCH_STRATSPEED;
2179:   case  3: return SCOTCH_STRATBALANCE;
2180:   case  4: return SCOTCH_STRATSAFETY;
2181:   case  5: return SCOTCH_STRATSCALABILITY;
2182:   case  6: return SCOTCH_STRATRECURSIVE;
2183:   case  7: return SCOTCH_STRATREMAP;
2184:   default: return SCOTCH_STRATDEFAULT;
2185:   }
2186: }

2188: static PetscErrorCode PTScotch_PartGraph_Seq(SCOTCH_Num strategy, double imbalance, SCOTCH_Num n, SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
2189:                                              SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num tpart[], SCOTCH_Num part[])
2190: {
2191:   SCOTCH_Graph   grafdat;
2192:   SCOTCH_Strat   stradat;
2193:   SCOTCH_Num     vertnbr = n;
2194:   SCOTCH_Num     edgenbr = xadj[n];
2195:   SCOTCH_Num*    velotab = vtxwgt;
2196:   SCOTCH_Num*    edlotab = adjwgt;
2197:   SCOTCH_Num     flagval = strategy;
2198:   double         kbalval = imbalance;

2202:   {
2203:     PetscBool flg = PETSC_TRUE;
2204:     PetscOptionsDeprecatedNoObject("-petscpartititoner_ptscotch_vertex_weight",NULL,"3.13","Use -petscpartitioner_use_vertex_weights");
2205:     PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
2206:     if (!flg) velotab = NULL;
2207:   }
2208:   SCOTCH_graphInit(&grafdat);CHKERRPTSCOTCH(ierr);
2209:   SCOTCH_graphBuild(&grafdat, 0, vertnbr, xadj, xadj + 1, velotab, NULL, edgenbr, adjncy, edlotab);CHKERRPTSCOTCH(ierr);
2210:   SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
2211:   SCOTCH_stratGraphMapBuild(&stradat, flagval, nparts, kbalval);CHKERRPTSCOTCH(ierr);
2212:   if (tpart) {
2213:     SCOTCH_Arch archdat;
2214:     SCOTCH_archInit(&archdat);CHKERRPTSCOTCH(ierr);
2215:     SCOTCH_archCmpltw(&archdat, nparts, tpart);CHKERRPTSCOTCH(ierr);
2216:     SCOTCH_graphMap(&grafdat, &archdat, &stradat, part);CHKERRPTSCOTCH(ierr);
2217:     SCOTCH_archExit(&archdat);
2218:   } else {
2219:     SCOTCH_graphPart(&grafdat, nparts, &stradat, part);CHKERRPTSCOTCH(ierr);
2220:   }
2221:   SCOTCH_stratExit(&stradat);
2222:   SCOTCH_graphExit(&grafdat);
2223:   return(0);
2224: }

2226: static PetscErrorCode PTScotch_PartGraph_MPI(SCOTCH_Num strategy, double imbalance, SCOTCH_Num vtxdist[], SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
2227:                                              SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num tpart[], SCOTCH_Num part[], MPI_Comm comm)
2228: {
2229:   PetscMPIInt     procglbnbr;
2230:   PetscMPIInt     proclocnum;
2231:   SCOTCH_Arch     archdat;
2232:   SCOTCH_Dgraph   grafdat;
2233:   SCOTCH_Dmapping mappdat;
2234:   SCOTCH_Strat    stradat;
2235:   SCOTCH_Num      vertlocnbr;
2236:   SCOTCH_Num      edgelocnbr;
2237:   SCOTCH_Num*     veloloctab = vtxwgt;
2238:   SCOTCH_Num*     edloloctab = adjwgt;
2239:   SCOTCH_Num      flagval = strategy;
2240:   double          kbalval = imbalance;
2241:   PetscErrorCode  ierr;

2244:   {
2245:     PetscBool flg = PETSC_TRUE;
2246:     PetscOptionsDeprecatedNoObject("-petscpartititoner_ptscotch_vertex_weight",NULL,"3.13","Use -petscpartitioner_use_vertex_weights");
2247:     PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
2248:     if (!flg) veloloctab = NULL;
2249:   }
2250:   MPI_Comm_size(comm, &procglbnbr);
2251:   MPI_Comm_rank(comm, &proclocnum);
2252:   vertlocnbr = vtxdist[proclocnum + 1] - vtxdist[proclocnum];
2253:   edgelocnbr = xadj[vertlocnbr];

2255:   SCOTCH_dgraphInit(&grafdat, comm);CHKERRPTSCOTCH(ierr);
2256:   SCOTCH_dgraphBuild(&grafdat, 0, vertlocnbr, vertlocnbr, xadj, xadj + 1, veloloctab, NULL, edgelocnbr, edgelocnbr, adjncy, NULL, edloloctab);CHKERRPTSCOTCH(ierr);
2257:   SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
2258:   SCOTCH_stratDgraphMapBuild(&stradat, flagval, procglbnbr, nparts, kbalval);
2259:   SCOTCH_archInit(&archdat);CHKERRPTSCOTCH(ierr);
2260:   if (tpart) { /* target partition weights */
2261:     SCOTCH_archCmpltw(&archdat, nparts, tpart);CHKERRPTSCOTCH(ierr);
2262:   } else {
2263:     SCOTCH_archCmplt(&archdat, nparts);CHKERRPTSCOTCH(ierr);
2264:   }
2265:   SCOTCH_dgraphMapInit(&grafdat, &mappdat, &archdat, part);CHKERRPTSCOTCH(ierr);

2267:   SCOTCH_dgraphMapCompute(&grafdat, &mappdat, &stradat);CHKERRPTSCOTCH(ierr);
2268:   SCOTCH_dgraphMapExit(&grafdat, &mappdat);
2269:   SCOTCH_archExit(&archdat);
2270:   SCOTCH_stratExit(&stradat);
2271:   SCOTCH_dgraphExit(&grafdat);
2272:   return(0);
2273: }

2275: #endif /* PETSC_HAVE_PTSCOTCH */

2277: static PetscErrorCode PetscPartitionerDestroy_PTScotch(PetscPartitioner part)
2278: {
2279:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2280:   PetscErrorCode             ierr;

2283:   MPI_Comm_free(&p->pcomm);
2284:   PetscFree(p);
2285:   return(0);
2286: }

2288: static PetscErrorCode PetscPartitionerView_PTScotch_Ascii(PetscPartitioner part, PetscViewer viewer)
2289: {
2290:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2291:   PetscErrorCode            ierr;

2294:   PetscViewerASCIIPushTab(viewer);
2295:   PetscViewerASCIIPrintf(viewer,"using partitioning strategy %s\n",PTScotchStrategyList[p->strategy]);
2296:   PetscViewerASCIIPrintf(viewer,"using load imbalance ratio %g\n",(double)p->imbalance);
2297:   PetscViewerASCIIPopTab(viewer);
2298:   return(0);
2299: }

2301: static PetscErrorCode PetscPartitionerView_PTScotch(PetscPartitioner part, PetscViewer viewer)
2302: {
2303:   PetscBool      iascii;

2309:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2310:   if (iascii) {PetscPartitionerView_PTScotch_Ascii(part, viewer);}
2311:   return(0);
2312: }

2314: static PetscErrorCode PetscPartitionerSetFromOptions_PTScotch(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
2315: {
2316:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2317:   const char *const         *slist = PTScotchStrategyList;
2318:   PetscInt                  nlist = (PetscInt)(sizeof(PTScotchStrategyList)/sizeof(PTScotchStrategyList[0]));
2319:   PetscBool                 flag;
2320:   PetscErrorCode            ierr;

2323:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner PTScotch Options");
2324:   PetscOptionsEList("-petscpartitioner_ptscotch_strategy","Partitioning strategy","",slist,nlist,slist[p->strategy],&p->strategy,&flag);
2325:   PetscOptionsReal("-petscpartitioner_ptscotch_imbalance","Load imbalance ratio","",p->imbalance,&p->imbalance,&flag);
2326:   PetscOptionsTail();
2327:   return(0);
2328: }

2330: static PetscErrorCode PetscPartitionerPartition_PTScotch(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
2331: {
2332: #if defined(PETSC_HAVE_PTSCOTCH)
2333:   MPI_Comm       comm;
2334:   PetscInt       nvtxs = numVertices;   /* The number of vertices in full graph */
2335:   PetscInt       *vtxdist;              /* Distribution of vertices across processes */
2336:   PetscInt       *xadj   = start;       /* Start of edge list for each vertex */
2337:   PetscInt       *adjncy = adjacency;   /* Edge lists for all vertices */
2338:   PetscInt       *vwgt   = NULL;        /* Vertex weights */
2339:   PetscInt       *adjwgt = NULL;        /* Edge weights */
2340:   PetscInt       v, i, *assignment, *points;
2341:   PetscMPIInt    size, rank, p;
2342:   PetscBool      hasempty = PETSC_FALSE;
2343:   PetscInt       *tpwgts = NULL;

2347:   PetscObjectGetComm((PetscObject)part,&comm);
2348:   MPI_Comm_size(comm, &size);
2349:   MPI_Comm_rank(comm, &rank);
2350:   PetscMalloc2(size+1,&vtxdist,PetscMax(nvtxs,1),&assignment);
2351:   /* Calculate vertex distribution */
2352:   vtxdist[0] = 0;
2353:   MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
2354:   for (p = 2; p <= size; ++p) {
2355:     hasempty = (PetscBool)(hasempty || !vtxdist[p-1] || !vtxdist[p]);
2356:     vtxdist[p] += vtxdist[p-1];
2357:   }
2358:   /* null graph */
2359:   if (vtxdist[size] == 0) {
2360:     PetscFree2(vtxdist, assignment);
2361:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
2362:     return(0);
2363:   }

2365:   /* Calculate vertex weights */
2366:   if (vertSection) {
2367:     PetscMalloc1(nvtxs,&vwgt);
2368:     for (v = 0; v < nvtxs; ++v) {
2369:       PetscSectionGetDof(vertSection, v, &vwgt[v]);
2370:     }
2371:   }

2373:   /* Calculate partition weights */
2374:   if (targetSection) {
2375:     PetscInt sumw;

2377:     PetscCalloc1(nparts,&tpwgts);
2378:     for (p = 0, sumw = 0; p < nparts; ++p) {
2379:       PetscSectionGetDof(targetSection,p,&tpwgts[p]);
2380:       sumw += tpwgts[p];
2381:     }
2382:     if (!sumw) {
2383:       PetscFree(tpwgts);
2384:     }
2385:   }

2387:   {
2388:     PetscPartitioner_PTScotch *pts = (PetscPartitioner_PTScotch *) part->data;
2389:     int                       strat = PTScotch_Strategy(pts->strategy);
2390:     double                    imbal = (double)pts->imbalance;

2392:     for (p = 0; !vtxdist[p+1] && p < size; ++p);
2393:     if (vtxdist[p+1] == vtxdist[size]) {
2394:       if (rank == p) {
2395:         PTScotch_PartGraph_Seq(strat, imbal, nvtxs, xadj, adjncy, vwgt, adjwgt, nparts, tpwgts, assignment);
2396:       }
2397:     } else {
2398:       PetscInt cnt;
2399:       MPI_Comm pcomm;

2401:       if (hasempty) {
2402:         MPI_Comm_split(pts->pcomm,!!nvtxs,rank,&pcomm);
2403:         for (p=0,cnt=0;p<size;p++) {
2404:           if (vtxdist[p+1] != vtxdist[p]) {
2405:             vtxdist[cnt+1] = vtxdist[p+1];
2406:             cnt++;
2407:           }
2408:         }
2409:       } else pcomm = pts->pcomm;
2410:       if (nvtxs) {
2411:         PTScotch_PartGraph_MPI(strat, imbal, vtxdist, xadj, adjncy, vwgt, adjwgt, nparts, tpwgts, assignment, pcomm);
2412:       }
2413:       if (hasempty) {
2414:         MPI_Comm_free(&pcomm);
2415:       }
2416:     }
2417:   }
2418:   PetscFree(vwgt);
2419:   PetscFree(tpwgts);

2421:   /* Convert to PetscSection+IS */
2422:   for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
2423:   PetscMalloc1(nvtxs, &points);
2424:   for (p = 0, i = 0; p < nparts; ++p) {
2425:     for (v = 0; v < nvtxs; ++v) {
2426:       if (assignment[v] == p) points[i++] = v;
2427:     }
2428:   }
2429:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
2430:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);

2432:   PetscFree2(vtxdist,assignment);
2433:   return(0);
2434: #else
2435:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-ptscotch.");
2436: #endif
2437: }

2439: static PetscErrorCode PetscPartitionerInitialize_PTScotch(PetscPartitioner part)
2440: {
2442:   part->noGraph             = PETSC_FALSE;
2443:   part->ops->view           = PetscPartitionerView_PTScotch;
2444:   part->ops->destroy        = PetscPartitionerDestroy_PTScotch;
2445:   part->ops->partition      = PetscPartitionerPartition_PTScotch;
2446:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_PTScotch;
2447:   return(0);
2448: }

2450: /*MC
2451:   PETSCPARTITIONERPTSCOTCH = "ptscotch" - A PetscPartitioner object using the PT-Scotch library

2453:   Level: intermediate

2455:   Options Database Keys:
2456: +  -petscpartitioner_ptscotch_strategy <string> - PT-Scotch strategy. Choose one of default quality speed balance safety scalability recursive remap
2457: -  -petscpartitioner_ptscotch_imbalance <val> - Load imbalance ratio

2459:   Notes: when the graph is on a single process, this partitioner actually uses Scotch and not PT-Scotch

2461: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
2462: M*/

2464: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_PTScotch(PetscPartitioner part)
2465: {
2466:   PetscPartitioner_PTScotch *p;
2467:   PetscErrorCode          ierr;

2471:   PetscNewLog(part, &p);
2472:   part->data = p;

2474:   MPI_Comm_dup(PetscObjectComm((PetscObject)part),&p->pcomm);
2475:   p->strategy  = 0;
2476:   p->imbalance = 0.01;

2478:   PetscPartitionerInitialize_PTScotch(part);
2479:   PetscCitationsRegister(PTScotchPartitionerCitation, &PTScotchPartitionercite);
2480:   return(0);
2481: }

2483: /*@
2484:   DMPlexGetPartitioner - Get the mesh partitioner

2486:   Not collective

2488:   Input Parameter:
2489: . dm - The DM

2491:   Output Parameter:
2492: . part - The PetscPartitioner

2494:   Level: developer

2496:   Note: This gets a borrowed reference, so the user should not destroy this PetscPartitioner.

2498: .seealso DMPlexDistribute(), DMPlexSetPartitioner(), PetscPartitionerDMPlexPartition(), PetscPartitionerCreate()
2499: @*/
2500: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
2501: {
2502:   DM_Plex       *mesh = (DM_Plex *) dm->data;

2507:   *part = mesh->partitioner;
2508:   return(0);
2509: }

2511: /*@
2512:   DMPlexSetPartitioner - Set the mesh partitioner

2514:   logically collective on DM

2516:   Input Parameters:
2517: + dm - The DM
2518: - part - The partitioner

2520:   Level: developer

2522:   Note: Any existing PetscPartitioner will be destroyed.

2524: .seealso DMPlexDistribute(), DMPlexGetPartitioner(), PetscPartitionerCreate()
2525: @*/
2526: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
2527: {
2528:   DM_Plex       *mesh = (DM_Plex *) dm->data;

2534:   PetscObjectReference((PetscObject)part);
2535:   PetscPartitionerDestroy(&mesh->partitioner);
2536:   mesh->partitioner = part;
2537:   return(0);
2538: }

2540: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
2541: {
2542:   const PetscInt *cone;
2543:   PetscInt       coneSize, c;
2544:   PetscBool      missing;

2548:   PetscHSetIQueryAdd(ht, point, &missing);
2549:   if (missing) {
2550:     DMPlexGetCone(dm, point, &cone);
2551:     DMPlexGetConeSize(dm, point, &coneSize);
2552:     for (c = 0; c < coneSize; c++) {
2553:       DMPlexAddClosure_Private(dm, ht, cone[c]);
2554:     }
2555:   }
2556:   return(0);
2557: }

2559: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
2560: {

2564:   if (up) {
2565:     PetscInt parent;

2567:     DMPlexGetTreeParent(dm,point,&parent,NULL);
2568:     if (parent != point) {
2569:       PetscInt closureSize, *closure = NULL, i;

2571:       DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
2572:       for (i = 0; i < closureSize; i++) {
2573:         PetscInt cpoint = closure[2*i];

2575:         PetscHSetIAdd(ht, cpoint);
2576:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
2577:       }
2578:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
2579:     }
2580:   }
2581:   if (down) {
2582:     PetscInt numChildren;
2583:     const PetscInt *children;

2585:     DMPlexGetTreeChildren(dm,point,&numChildren,&children);
2586:     if (numChildren) {
2587:       PetscInt i;

2589:       for (i = 0; i < numChildren; i++) {
2590:         PetscInt cpoint = children[i];

2592:         PetscHSetIAdd(ht, cpoint);
2593:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
2594:       }
2595:     }
2596:   }
2597:   return(0);
2598: }

2600: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
2601: {
2602:   PetscInt       parent;

2606:   DMPlexGetTreeParent(dm, point, &parent,NULL);
2607:   if (point != parent) {
2608:     const PetscInt *cone;
2609:     PetscInt       coneSize, c;

2611:     DMPlexAddClosureTree_Up_Private(dm, ht, parent);
2612:     DMPlexAddClosure_Private(dm, ht, parent);
2613:     DMPlexGetCone(dm, parent, &cone);
2614:     DMPlexGetConeSize(dm, parent, &coneSize);
2615:     for (c = 0; c < coneSize; c++) {
2616:       const PetscInt cp = cone[c];

2618:       DMPlexAddClosureTree_Up_Private(dm, ht, cp);
2619:     }
2620:   }
2621:   return(0);
2622: }

2624: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
2625: {
2626:   PetscInt       i, numChildren;
2627:   const PetscInt *children;

2631:   DMPlexGetTreeChildren(dm, point, &numChildren, &children);
2632:   for (i = 0; i < numChildren; i++) {
2633:     PetscHSetIAdd(ht, children[i]);
2634:   }
2635:   return(0);
2636: }

2638: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
2639: {
2640:   const PetscInt *cone;
2641:   PetscInt       coneSize, c;

2645:   PetscHSetIAdd(ht, point);
2646:   DMPlexAddClosureTree_Up_Private(dm, ht, point);
2647:   DMPlexAddClosureTree_Down_Private(dm, ht, point);
2648:   DMPlexGetCone(dm, point, &cone);
2649:   DMPlexGetConeSize(dm, point, &coneSize);
2650:   for (c = 0; c < coneSize; c++) {
2651:     DMPlexAddClosureTree_Private(dm, ht, cone[c]);
2652:   }
2653:   return(0);
2654: }

2656: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
2657: {
2658:   DM_Plex         *mesh = (DM_Plex *)dm->data;
2659:   const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
2660:   PetscInt        nelems, *elems, off = 0, p;
2661:   PetscHSetI      ht;
2662:   PetscErrorCode  ierr;

2665:   PetscHSetICreate(&ht);
2666:   PetscHSetIResize(ht, numPoints*16);
2667:   if (!hasTree) {
2668:     for (p = 0; p < numPoints; ++p) {
2669:       DMPlexAddClosure_Private(dm, ht, points[p]);
2670:     }
2671:   } else {
2672: #if 1
2673:     for (p = 0; p < numPoints; ++p) {
2674:       DMPlexAddClosureTree_Private(dm, ht, points[p]);
2675:     }
2676: #else
2677:     PetscInt  *closure = NULL, closureSize, c;
2678:     for (p = 0; p < numPoints; ++p) {
2679:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2680:       for (c = 0; c < closureSize*2; c += 2) {
2681:         PetscHSetIAdd(ht, closure[c]);
2682:         if (hasTree) {DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);}
2683:       }
2684:     }
2685:     if (closure) {DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);}
2686: #endif
2687:   }
2688:   PetscHSetIGetSize(ht, &nelems);
2689:   PetscMalloc1(nelems, &elems);
2690:   PetscHSetIGetElems(ht, &off, elems);
2691:   PetscHSetIDestroy(&ht);
2692:   PetscSortInt(nelems, elems);
2693:   ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
2694:   return(0);
2695: }

2697: /*@
2698:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

2700:   Input Parameters:
2701: + dm     - The DM
2702: - label  - DMLabel assinging ranks to remote roots

2704:   Level: developer

2706: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
2707: @*/
2708: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
2709: {
2710:   IS              rankIS,   pointIS, closureIS;
2711:   const PetscInt *ranks,   *points;
2712:   PetscInt        numRanks, numPoints, r;
2713:   PetscErrorCode  ierr;

2716:   DMLabelGetValueIS(label, &rankIS);
2717:   ISGetLocalSize(rankIS, &numRanks);
2718:   ISGetIndices(rankIS, &ranks);
2719:   for (r = 0; r < numRanks; ++r) {
2720:     const PetscInt rank = ranks[r];
2721:     DMLabelGetStratumIS(label, rank, &pointIS);
2722:     ISGetLocalSize(pointIS, &numPoints);
2723:     ISGetIndices(pointIS, &points);
2724:     DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
2725:     ISRestoreIndices(pointIS, &points);
2726:     ISDestroy(&pointIS);
2727:     DMLabelSetStratumIS(label, rank, closureIS);
2728:     ISDestroy(&closureIS);
2729:   }
2730:   ISRestoreIndices(rankIS, &ranks);
2731:   ISDestroy(&rankIS);
2732:   return(0);
2733: }

2735: /*@
2736:   DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label

2738:   Input Parameters:
2739: + dm     - The DM
2740: - label  - DMLabel assinging ranks to remote roots

2742:   Level: developer

2744: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2745: @*/
2746: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
2747: {
2748:   IS              rankIS,   pointIS;
2749:   const PetscInt *ranks,   *points;
2750:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
2751:   PetscInt       *adj = NULL;
2752:   PetscErrorCode  ierr;

2755:   DMLabelGetValueIS(label, &rankIS);
2756:   ISGetLocalSize(rankIS, &numRanks);
2757:   ISGetIndices(rankIS, &ranks);
2758:   for (r = 0; r < numRanks; ++r) {
2759:     const PetscInt rank = ranks[r];

2761:     DMLabelGetStratumIS(label, rank, &pointIS);
2762:     ISGetLocalSize(pointIS, &numPoints);
2763:     ISGetIndices(pointIS, &points);
2764:     for (p = 0; p < numPoints; ++p) {
2765:       adjSize = PETSC_DETERMINE;
2766:       DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
2767:       for (a = 0; a < adjSize; ++a) {DMLabelSetValue(label, adj[a], rank);}
2768:     }
2769:     ISRestoreIndices(pointIS, &points);
2770:     ISDestroy(&pointIS);
2771:   }
2772:   ISRestoreIndices(rankIS, &ranks);
2773:   ISDestroy(&rankIS);
2774:   PetscFree(adj);
2775:   return(0);
2776: }

2778: /*@
2779:   DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point SF

2781:   Input Parameters:
2782: + dm     - The DM
2783: - label  - DMLabel assinging ranks to remote roots

2785:   Level: developer

2787:   Note: This is required when generating multi-level overlaps to capture
2788:   overlap points from non-neighbouring partitions.

2790: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2791: @*/
2792: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
2793: {
2794:   MPI_Comm        comm;
2795:   PetscMPIInt     rank;
2796:   PetscSF         sfPoint;
2797:   DMLabel         lblRoots, lblLeaves;
2798:   IS              rankIS, pointIS;
2799:   const PetscInt *ranks;
2800:   PetscInt        numRanks, r;
2801:   PetscErrorCode  ierr;

2804:   PetscObjectGetComm((PetscObject) dm, &comm);
2805:   MPI_Comm_rank(comm, &rank);
2806:   DMGetPointSF(dm, &sfPoint);
2807:   /* Pull point contributions from remote leaves into local roots */
2808:   DMLabelGather(label, sfPoint, &lblLeaves);
2809:   DMLabelGetValueIS(lblLeaves, &rankIS);
2810:   ISGetLocalSize(rankIS, &numRanks);
2811:   ISGetIndices(rankIS, &ranks);
2812:   for (r = 0; r < numRanks; ++r) {
2813:     const PetscInt remoteRank = ranks[r];
2814:     if (remoteRank == rank) continue;
2815:     DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
2816:     DMLabelInsertIS(label, pointIS, remoteRank);
2817:     ISDestroy(&pointIS);
2818:   }
2819:   ISRestoreIndices(rankIS, &ranks);
2820:   ISDestroy(&rankIS);
2821:   DMLabelDestroy(&lblLeaves);
2822:   /* Push point contributions from roots into remote leaves */
2823:   DMLabelDistribute(label, sfPoint, &lblRoots);
2824:   DMLabelGetValueIS(lblRoots, &rankIS);
2825:   ISGetLocalSize(rankIS, &numRanks);
2826:   ISGetIndices(rankIS, &ranks);
2827:   for (r = 0; r < numRanks; ++r) {
2828:     const PetscInt remoteRank = ranks[r];
2829:     if (remoteRank == rank) continue;
2830:     DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
2831:     DMLabelInsertIS(label, pointIS, remoteRank);
2832:     ISDestroy(&pointIS);
2833:   }
2834:   ISRestoreIndices(rankIS, &ranks);
2835:   ISDestroy(&rankIS);
2836:   DMLabelDestroy(&lblRoots);
2837:   return(0);
2838: }

2840: /*@
2841:   DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label

2843:   Input Parameters:
2844: + dm        - The DM
2845: . rootLabel - DMLabel assinging ranks to local roots
2846: - processSF - A star forest mapping into the local index on each remote rank

2848:   Output Parameter:
2849: . leafLabel - DMLabel assinging ranks to remote roots

2851:   Note: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
2852:   resulting leafLabel is a receiver mapping of remote roots to their parent rank.

2854:   Level: developer

2856: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2857: @*/
2858: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
2859: {
2860:   MPI_Comm           comm;
2861:   PetscMPIInt        rank, size, r;
2862:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
2863:   PetscSF            sfPoint;
2864:   PetscSection       rootSection;
2865:   PetscSFNode       *rootPoints, *leafPoints;
2866:   const PetscSFNode *remote;
2867:   const PetscInt    *local, *neighbors;
2868:   IS                 valueIS;
2869:   PetscBool          mpiOverflow = PETSC_FALSE;
2870:   PetscErrorCode     ierr;

2873:   PetscLogEventBegin(DMPLEX_PartLabelInvert,dm,0,0,0);
2874:   PetscObjectGetComm((PetscObject) dm, &comm);
2875:   MPI_Comm_rank(comm, &rank);
2876:   MPI_Comm_size(comm, &size);
2877:   DMGetPointSF(dm, &sfPoint);

2879:   /* Convert to (point, rank) and use actual owners */
2880:   PetscSectionCreate(comm, &rootSection);
2881:   PetscSectionSetChart(rootSection, 0, size);
2882:   DMLabelGetValueIS(rootLabel, &valueIS);
2883:   ISGetLocalSize(valueIS, &numNeighbors);
2884:   ISGetIndices(valueIS, &neighbors);
2885:   for (n = 0; n < numNeighbors; ++n) {
2886:     DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
2887:     PetscSectionAddDof(rootSection, neighbors[n], numPoints);
2888:   }
2889:   PetscSectionSetUp(rootSection);
2890:   PetscSectionGetStorageSize(rootSection, &rootSize);
2891:   PetscMalloc1(rootSize, &rootPoints);
2892:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
2893:   for (n = 0; n < numNeighbors; ++n) {
2894:     IS              pointIS;
2895:     const PetscInt *points;

2897:     PetscSectionGetOffset(rootSection, neighbors[n], &off);
2898:     DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
2899:     ISGetLocalSize(pointIS, &numPoints);
2900:     ISGetIndices(pointIS, &points);
2901:     for (p = 0; p < numPoints; ++p) {
2902:       if (local) {PetscFindInt(points[p], nleaves, local, &l);}
2903:       else       {l = -1;}
2904:       if (l >= 0) {rootPoints[off+p] = remote[l];}
2905:       else        {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
2906:     }
2907:     ISRestoreIndices(pointIS, &points);
2908:     ISDestroy(&pointIS);
2909:   }

2911:   /* Try to communicate overlap using All-to-All */
2912:   if (!processSF) {
2913:     PetscInt64  counter = 0;
2914:     PetscBool   locOverflow = PETSC_FALSE;
2915:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

2917:     PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
2918:     for (n = 0; n < numNeighbors; ++n) {
2919:       PetscSectionGetDof(rootSection, neighbors[n], &dof);
2920:       PetscSectionGetOffset(rootSection, neighbors[n], &off);
2921: #if defined(PETSC_USE_64BIT_INDICES)
2922:       if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2923:       if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2924: #endif
2925:       scounts[neighbors[n]] = (PetscMPIInt) dof;
2926:       sdispls[neighbors[n]] = (PetscMPIInt) off;
2927:     }
2928:     MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
2929:     for (r = 0; r < size; ++r) { rdispls[r] = (int)counter; counter += rcounts[r]; }
2930:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
2931:     MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
2932:     if (!mpiOverflow) {
2933:       PetscInfo(dm,"Using Alltoallv for mesh distribution\n");
2934:       leafSize = (PetscInt) counter;
2935:       PetscMalloc1(leafSize, &leafPoints);
2936:       MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
2937:     }
2938:     PetscFree4(scounts, sdispls, rcounts, rdispls);
2939:   }

2941:   /* Communicate overlap using process star forest */
2942:   if (processSF || mpiOverflow) {
2943:     PetscSF      procSF;
2944:     PetscSection leafSection;

2946:     if (processSF) {
2947:       PetscInfo(dm,"Using processSF for mesh distribution\n");
2948:       PetscObjectReference((PetscObject)processSF);
2949:       procSF = processSF;
2950:     } else {
2951:       PetscInfo(dm,"Using processSF for mesh distribution (MPI overflow)\n");
2952:       PetscSFCreate(comm,&procSF);
2953:       PetscSFSetGraphWithPattern(procSF,NULL,PETSCSF_PATTERN_ALLTOALL);
2954:     }

2956:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
2957:     DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
2958:     PetscSectionGetStorageSize(leafSection, &leafSize);
2959:     PetscSectionDestroy(&leafSection);
2960:     PetscSFDestroy(&procSF);
2961:   }

2963:   for (p = 0; p < leafSize; p++) {
2964:     DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
2965:   }

2967:   ISRestoreIndices(valueIS, &neighbors);
2968:   ISDestroy(&valueIS);
2969:   PetscSectionDestroy(&rootSection);
2970:   PetscFree(rootPoints);
2971:   PetscFree(leafPoints);
2972:   PetscLogEventEnd(DMPLEX_PartLabelInvert,dm,0,0,0);
2973:   return(0);
2974: }

2976: /*@
2977:   DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points

2979:   Input Parameters:
2980: + dm    - The DM
2981: - label - DMLabel assinging ranks to remote roots

2983:   Output Parameter:
2984: . sf    - The star forest communication context encapsulating the defined mapping

2986:   Note: The incoming label is a receiver mapping of remote points to their parent rank.

2988:   Level: developer

2990: .seealso: DMPlexDistribute(), DMPlexCreateOverlap()
2991: @*/
2992: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
2993: {
2994:   PetscMPIInt     rank;
2995:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
2996:   PetscSFNode    *remotePoints;
2997:   IS              remoteRootIS, neighborsIS;
2998:   const PetscInt *remoteRoots, *neighbors;

3002:   PetscLogEventBegin(DMPLEX_PartLabelCreateSF,dm,0,0,0);
3003:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);

3005:   DMLabelGetValueIS(label, &neighborsIS);
3006: #if 0
3007:   {
3008:     IS is;
3009:     ISDuplicate(neighborsIS, &is);
3010:     ISSort(is);
3011:     ISDestroy(&neighborsIS);
3012:     neighborsIS = is;
3013:   }
3014: #endif
3015:   ISGetLocalSize(neighborsIS, &nNeighbors);
3016:   ISGetIndices(neighborsIS, &neighbors);
3017:   for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
3018:     DMLabelGetStratumSize(label, neighbors[n], &numPoints);
3019:     numRemote += numPoints;
3020:   }
3021:   PetscMalloc1(numRemote, &remotePoints);
3022:   /* Put owned points first */
3023:   DMLabelGetStratumSize(label, rank, &numPoints);
3024:   if (numPoints > 0) {
3025:     DMLabelGetStratumIS(label, rank, &remoteRootIS);
3026:     ISGetIndices(remoteRootIS, &remoteRoots);
3027:     for (p = 0; p < numPoints; p++) {
3028:       remotePoints[idx].index = remoteRoots[p];
3029:       remotePoints[idx].rank = rank;
3030:       idx++;
3031:     }
3032:     ISRestoreIndices(remoteRootIS, &remoteRoots);
3033:     ISDestroy(&remoteRootIS);
3034:   }
3035:   /* Now add remote points */
3036:   for (n = 0; n < nNeighbors; ++n) {
3037:     const PetscInt nn = neighbors[n];

3039:     DMLabelGetStratumSize(label, nn, &numPoints);
3040:     if (nn == rank || numPoints <= 0) continue;
3041:     DMLabelGetStratumIS(label, nn, &remoteRootIS);
3042:     ISGetIndices(remoteRootIS, &remoteRoots);
3043:     for (p = 0; p < numPoints; p++) {
3044:       remotePoints[idx].index = remoteRoots[p];
3045:       remotePoints[idx].rank = nn;
3046:       idx++;
3047:     }
3048:     ISRestoreIndices(remoteRootIS, &remoteRoots);
3049:     ISDestroy(&remoteRootIS);
3050:   }
3051:   PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
3052:   DMPlexGetChart(dm, &pStart, &pEnd);
3053:   PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
3054:   PetscSFSetUp(*sf);
3055:   ISDestroy(&neighborsIS);
3056:   PetscLogEventEnd(DMPLEX_PartLabelCreateSF,dm,0,0,0);
3057:   return(0);
3058: }

3060: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
3061:  * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
3062:  * them out in that case. */
3063: #if defined(PETSC_HAVE_PARMETIS)
3064: /*@C

3066:   DMPlexRewriteSF - Rewrites the ownership of the SF of a DM (in place).

3068:   Input parameters:
3069: + dm                - The DMPlex object.
3070: . n                 - The number of points.
3071: . pointsToRewrite   - The points in the SF whose ownership will change.
3072: . targetOwners      - New owner for each element in pointsToRewrite.
3073: - degrees           - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.

3075:   Level: developer

3077: @*/
3078: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
3079: {
3080:   PetscInt      ierr, pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
3081:   PetscInt     *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
3082:   PetscSFNode  *leafLocationsNew;
3083:   const         PetscSFNode *iremote;
3084:   const         PetscInt *ilocal;
3085:   PetscBool    *isLeaf;
3086:   PetscSF       sf;
3087:   MPI_Comm      comm;
3088:   PetscMPIInt   rank, size;

3091:   PetscObjectGetComm((PetscObject) dm, &comm);
3092:   MPI_Comm_rank(comm, &rank);
3093:   MPI_Comm_size(comm, &size);
3094:   DMPlexGetChart(dm, &pStart, &pEnd);

3096:   DMGetPointSF(dm, &sf);
3097:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
3098:   PetscMalloc1(pEnd-pStart, &isLeaf);
3099:   for (i=0; i<pEnd-pStart; i++) {
3100:     isLeaf[i] = PETSC_FALSE;
3101:   }
3102:   for (i=0; i<nleafs; i++) {
3103:     isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
3104:   }

3106:   PetscMalloc1(pEnd-pStart+1, &cumSumDegrees);
3107:   cumSumDegrees[0] = 0;
3108:   for (i=1; i<=pEnd-pStart; i++) {
3109:     cumSumDegrees[i] = cumSumDegrees[i-1] + degrees[i-1];
3110:   }
3111:   sumDegrees = cumSumDegrees[pEnd-pStart];
3112:   /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */

3114:   PetscMalloc1(sumDegrees, &locationsOfLeafs);
3115:   PetscMalloc1(pEnd-pStart, &rankOnLeafs);
3116:   for (i=0; i<pEnd-pStart; i++) {
3117:     rankOnLeafs[i] = rank;
3118:   }
3119:   PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
3120:   PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
3121:   PetscFree(rankOnLeafs);

3123:   /* get the remote local points of my leaves */
3124:   PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
3125:   PetscMalloc1(pEnd-pStart, &points);
3126:   for (i=0; i<pEnd-pStart; i++) {
3127:     points[i] = pStart+i;
3128:   }
3129:   PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
3130:   PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
3131:   PetscFree(points);
3132:   /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
3133:   PetscMalloc1(pEnd-pStart, &newOwners);
3134:   PetscMalloc1(pEnd-pStart, &newNumbers);
3135:   for (i=0; i<pEnd-pStart; i++) {
3136:     newOwners[i] = -1;
3137:     newNumbers[i] = -1;
3138:   }
3139:   {
3140:     PetscInt oldNumber, newNumber, oldOwner, newOwner;
3141:     for (i=0; i<n; i++) {
3142:       oldNumber = pointsToRewrite[i];
3143:       newNumber = -1;
3144:       oldOwner = rank;
3145:       newOwner = targetOwners[i];
3146:       if (oldOwner == newOwner) {
3147:         newNumber = oldNumber;
3148:       } else {
3149:         for (j=0; j<degrees[oldNumber]; j++) {
3150:           if (locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) {
3151:             newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j];
3152:             break;
3153:           }
3154:         }
3155:       }
3156:       if (newNumber == -1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");

3158:       newOwners[oldNumber] = newOwner;
3159:       newNumbers[oldNumber] = newNumber;
3160:     }
3161:   }
3162:   PetscFree(cumSumDegrees);
3163:   PetscFree(locationsOfLeafs);
3164:   PetscFree(remoteLocalPointOfLeafs);

3166:   PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners);
3167:   PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners);
3168:   PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers);
3169:   PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers);

3171:   /* Now count how many leafs we have on each processor. */
3172:   leafCounter=0;
3173:   for (i=0; i<pEnd-pStart; i++) {
3174:     if (newOwners[i] >= 0) {
3175:       if (newOwners[i] != rank) {
3176:         leafCounter++;
3177:       }
3178:     } else {
3179:       if (isLeaf[i]) {
3180:         leafCounter++;
3181:       }
3182:     }
3183:   }

3185:   /* Now set up the new sf by creating the leaf arrays */
3186:   PetscMalloc1(leafCounter, &leafsNew);
3187:   PetscMalloc1(leafCounter, &leafLocationsNew);

3189:   leafCounter = 0;
3190:   counter = 0;
3191:   for (i=0; i<pEnd-pStart; i++) {
3192:     if (newOwners[i] >= 0) {
3193:       if (newOwners[i] != rank) {
3194:         leafsNew[leafCounter] = i;
3195:         leafLocationsNew[leafCounter].rank = newOwners[i];
3196:         leafLocationsNew[leafCounter].index = newNumbers[i];
3197:         leafCounter++;
3198:       }
3199:     } else {
3200:       if (isLeaf[i]) {
3201:         leafsNew[leafCounter] = i;
3202:         leafLocationsNew[leafCounter].rank = iremote[counter].rank;
3203:         leafLocationsNew[leafCounter].index = iremote[counter].index;
3204:         leafCounter++;
3205:       }
3206:     }
3207:     if (isLeaf[i]) {
3208:       counter++;
3209:     }
3210:   }

3212:   PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
3213:   PetscFree(newOwners);
3214:   PetscFree(newNumbers);
3215:   PetscFree(isLeaf);
3216:   return(0);
3217: }

3219: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
3220: {
3221:   PetscInt *distribution, min, max, sum, i, ierr;
3222:   PetscMPIInt rank, size;
3224:   MPI_Comm_size(comm, &size);
3225:   MPI_Comm_rank(comm, &rank);
3226:   PetscCalloc1(size, &distribution);
3227:   for (i=0; i<n; i++) {
3228:     if (part) distribution[part[i]] += vtxwgt[skip*i];
3229:     else distribution[rank] += vtxwgt[skip*i];
3230:   }
3231:   MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
3232:   min = distribution[0];
3233:   max = distribution[0];
3234:   sum = distribution[0];
3235:   for (i=1; i<size; i++) {
3236:     if (distribution[i]<min) min=distribution[i];
3237:     if (distribution[i]>max) max=distribution[i];
3238:     sum += distribution[i];
3239:   }
3240:   PetscViewerASCIIPrintf(viewer, "Min: %D, Avg: %D, Max: %D, Balance: %f\n", min, sum/size, max, (max*1.*size)/sum);
3241:   PetscFree(distribution);
3242:   return(0);
3243: }

3245: #endif

3247: /*@
3248:   DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the PointSF of the DM inplace.

3250:   Input parameters:
3251: + dm               - The DMPlex object.
3252: . entityDepth      - depth of the entity to balance (0 -> balance vertices).
3253: . useInitialGuess  - whether to use the current distribution as initial guess (only used by ParMETIS).
3254: - parallel         - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.

3256:   Output parameters:
3257: . success          - whether the graph partitioning was successful or not. If not, try useInitialGuess=True and parallel=True.

3259:   Level: intermediate

3261: @*/

3263: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
3264: {
3265: #if defined(PETSC_HAVE_PARMETIS)
3266:   PetscSF     sf;
3267:   PetscInt    ierr, i, j, idx, jdx;
3268:   PetscInt    eBegin, eEnd, nroots, nleafs, pStart, pEnd;
3269:   const       PetscInt *degrees, *ilocal;
3270:   const       PetscSFNode *iremote;
3271:   PetscBool   *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
3272:   PetscInt    numExclusivelyOwned, numNonExclusivelyOwned;
3273:   PetscMPIInt rank, size;
3274:   PetscInt    *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
3275:   const       PetscInt *cumSumVertices;
3276:   PetscInt    offset, counter;
3277:   PetscInt    lenadjncy;
3278:   PetscInt    *xadj, *adjncy, *vtxwgt;
3279:   PetscInt    lenxadj;
3280:   PetscInt    *adjwgt = NULL;
3281:   PetscInt    *part, *options;
3282:   PetscInt    nparts, wgtflag, numflag, ncon, edgecut;
3283:   real_t      *ubvec;
3284:   PetscInt    *firstVertices, *renumbering;
3285:   PetscInt    failed, failedGlobal;
3286:   MPI_Comm    comm;
3287:   Mat         A, Apre;
3288:   const char *prefix = NULL;
3289:   PetscViewer       viewer;
3290:   PetscViewerFormat format;
3291:   PetscLayout layout;

3294:   if (success) *success = PETSC_FALSE;
3295:   PetscObjectGetComm((PetscObject) dm, &comm);
3296:   MPI_Comm_rank(comm, &rank);
3297:   MPI_Comm_size(comm, &size);
3298:   if (size==1) return(0);

3300:   PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);

3302:   PetscOptionsGetViewer(comm,((PetscObject)dm)->options, prefix,"-dm_rebalance_partition_view",&viewer,&format,NULL);
3303:   if (viewer) {
3304:     PetscViewerPushFormat(viewer,format);
3305:   }

3307:   /* Figure out all points in the plex that we are interested in balancing. */
3308:   DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
3309:   DMPlexGetChart(dm, &pStart, &pEnd);
3310:   PetscMalloc1(pEnd-pStart, &toBalance);

3312:   for (i=0; i<pEnd-pStart; i++) {
3313:     toBalance[i] = (PetscBool)(i-pStart>=eBegin && i-pStart<eEnd);
3314:   }

3316:   /* There are three types of points:
3317:    * exclusivelyOwned: points that are owned by this process and only seen by this process
3318:    * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
3319:    * leaf: a point that is seen by this process but owned by a different process
3320:    */
3321:   DMGetPointSF(dm, &sf);
3322:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
3323:   PetscMalloc1(pEnd-pStart, &isLeaf);
3324:   PetscMalloc1(pEnd-pStart, &isNonExclusivelyOwned);
3325:   PetscMalloc1(pEnd-pStart, &isExclusivelyOwned);
3326:   for (i=0; i<pEnd-pStart; i++) {
3327:     isNonExclusivelyOwned[i] = PETSC_FALSE;
3328:     isExclusivelyOwned[i] = PETSC_FALSE;
3329:     isLeaf[i] = PETSC_FALSE;
3330:   }

3332:   /* start by marking all the leafs */
3333:   for (i=0; i<nleafs; i++) {
3334:     isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
3335:   }

3337:   /* for an owned point, we can figure out whether another processor sees it or
3338:    * not by calculating its degree */
3339:   PetscSFComputeDegreeBegin(sf, &degrees);
3340:   PetscSFComputeDegreeEnd(sf, &degrees);

3342:   numExclusivelyOwned = 0;
3343:   numNonExclusivelyOwned = 0;
3344:   for (i=0; i<pEnd-pStart; i++) {
3345:     if (toBalance[i]) {
3346:       if (degrees[i] > 0) {
3347:         isNonExclusivelyOwned[i] = PETSC_TRUE;
3348:         numNonExclusivelyOwned += 1;
3349:       } else {
3350:         if (!isLeaf[i]) {
3351:           isExclusivelyOwned[i] = PETSC_TRUE;
3352:           numExclusivelyOwned += 1;
3353:         }
3354:       }
3355:     }
3356:   }

3358:   /* We are going to build a graph with one vertex per core representing the
3359:    * exclusively owned points and then one vertex per nonExclusively owned
3360:    * point. */

3362:   PetscLayoutCreate(comm, &layout);
3363:   PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
3364:   PetscLayoutSetUp(layout);
3365:   PetscLayoutGetRanges(layout, &cumSumVertices);

3367:   PetscMalloc1(pEnd-pStart, &globalNumbersOfLocalOwnedVertices);
3368:   for (i=0; i<pEnd-pStart; i++) {globalNumbersOfLocalOwnedVertices[i] = pStart - 1;}
3369:   offset = cumSumVertices[rank];
3370:   counter = 0;
3371:   for (i=0; i<pEnd-pStart; i++) {
3372:     if (toBalance[i]) {
3373:       if (degrees[i] > 0) {
3374:         globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
3375:         counter++;
3376:       }
3377:     }
3378:   }

3380:   /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
3381:   PetscMalloc1(pEnd-pStart, &leafGlobalNumbers);
3382:   PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);
3383:   PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);

3385:   /* Now start building the data structure for ParMETIS */

3387:   MatCreate(comm, &Apre);
3388:   MatSetType(Apre, MATPREALLOCATOR);
3389:   MatSetSizes(Apre, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
3390:   MatSetUp(Apre);

3392:   for (i=0; i<pEnd-pStart; i++) {
3393:     if (toBalance[i]) {
3394:       idx = cumSumVertices[rank];
3395:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
3396:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
3397:       else continue;
3398:       MatSetValue(Apre, idx, jdx, 1, INSERT_VALUES);
3399:       MatSetValue(Apre, jdx, idx, 1, INSERT_VALUES);
3400:     }
3401:   }

3403:   MatAssemblyBegin(Apre, MAT_FINAL_ASSEMBLY);
3404:   MatAssemblyEnd(Apre, MAT_FINAL_ASSEMBLY);

3406:   MatCreate(comm, &A);
3407:   MatSetType(A, MATMPIAIJ);
3408:   MatSetSizes(A, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
3409:   MatPreallocatorPreallocate(Apre, PETSC_FALSE, A);
3410:   MatDestroy(&Apre);

3412:   for (i=0; i<pEnd-pStart; i++) {
3413:     if (toBalance[i]) {
3414:       idx = cumSumVertices[rank];
3415:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
3416:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
3417:       else continue;
3418:       MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
3419:       MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
3420:     }
3421:   }

3423:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
3424:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
3425:   PetscFree(leafGlobalNumbers);
3426:   PetscFree(globalNumbersOfLocalOwnedVertices);

3428:   nparts = size;
3429:   wgtflag = 2;
3430:   numflag = 0;
3431:   ncon = 2;
3432:   real_t *tpwgts;
3433:   PetscMalloc1(ncon * nparts, &tpwgts);
3434:   for (i=0; i<ncon*nparts; i++) {
3435:     tpwgts[i] = 1./(nparts);
3436:   }

3438:   PetscMalloc1(ncon, &ubvec);
3439:   ubvec[0] = 1.01;
3440:   ubvec[1] = 1.01;
3441:   lenadjncy = 0;
3442:   for (i=0; i<1+numNonExclusivelyOwned; i++) {
3443:     PetscInt temp=0;
3444:     MatGetRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
3445:     lenadjncy += temp;
3446:     MatRestoreRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
3447:   }
3448:   PetscMalloc1(lenadjncy, &adjncy);
3449:   lenxadj = 2 + numNonExclusivelyOwned;
3450:   PetscMalloc1(lenxadj, &xadj);
3451:   xadj[0] = 0;
3452:   counter = 0;
3453:   for (i=0; i<1+numNonExclusivelyOwned; i++) {
3454:     PetscInt        temp=0;
3455:     const PetscInt *cols;
3456:     MatGetRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
3457:     PetscArraycpy(&adjncy[counter], cols, temp);
3458:     counter += temp;
3459:     xadj[i+1] = counter;
3460:     MatRestoreRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
3461:   }

3463:   PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);
3464:   PetscMalloc1(ncon*(1 + numNonExclusivelyOwned), &vtxwgt);
3465:   vtxwgt[0] = numExclusivelyOwned;
3466:   if (ncon>1) vtxwgt[1] = 1;
3467:   for (i=0; i<numNonExclusivelyOwned; i++) {
3468:     vtxwgt[ncon*(i+1)] = 1;
3469:     if (ncon>1) vtxwgt[ncon*(i+1)+1] = 0;
3470:   }

3472:   if (viewer) {
3473:     PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %D on interface of mesh distribution.\n", entityDepth);
3474:     PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %D\n", cumSumVertices[size]);
3475:   }
3476:   if (parallel) {
3477:     PetscMalloc1(4, &options);
3478:     options[0] = 1;
3479:     options[1] = 0; /* Verbosity */
3480:     options[2] = 0; /* Seed */
3481:     options[3] = PARMETIS_PSR_COUPLED; /* Seed */
3482:     if (viewer) { PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"); }
3483:     if (useInitialGuess) {
3484:       if (viewer) { PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"); }
3485:       PetscStackPush("ParMETIS_V3_RefineKway");
3486:       ParMETIS_V3_RefineKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
3487:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_RefineKway()");
3488:       PetscStackPop;
3489:     } else {
3490:       PetscStackPush("ParMETIS_V3_PartKway");
3491:       ParMETIS_V3_PartKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
3492:       PetscStackPop;
3493:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
3494:     }
3495:     PetscFree(options);
3496:   } else {
3497:     if (viewer) { PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"); }
3498:     Mat As;
3499:     PetscInt numRows;
3500:     PetscInt *partGlobal;
3501:     MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &As);

3503:     PetscInt *numExclusivelyOwnedAll;
3504:     PetscMalloc1(size, &numExclusivelyOwnedAll);
3505:     numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
3506:     MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,numExclusivelyOwnedAll,1,MPIU_INT,comm);

3508:     MatGetSize(As, &numRows, NULL);
3509:     PetscMalloc1(numRows, &partGlobal);
3510:     if (!rank) {
3511:       PetscInt *adjncy_g, *xadj_g, *vtxwgt_g;
3512:       lenadjncy = 0;

3514:       for (i=0; i<numRows; i++) {
3515:         PetscInt temp=0;
3516:         MatGetRow(As, i, &temp, NULL, NULL);
3517:         lenadjncy += temp;
3518:         MatRestoreRow(As, i, &temp, NULL, NULL);
3519:       }
3520:       PetscMalloc1(lenadjncy, &adjncy_g);
3521:       lenxadj = 1 + numRows;
3522:       PetscMalloc1(lenxadj, &xadj_g);
3523:       xadj_g[0] = 0;
3524:       counter = 0;
3525:       for (i=0; i<numRows; i++) {
3526:         PetscInt        temp=0;
3527:         const PetscInt *cols;
3528:         MatGetRow(As, i, &temp, &cols, NULL);
3529:         PetscArraycpy(&adjncy_g[counter], cols, temp);
3530:         counter += temp;
3531:         xadj_g[i+1] = counter;
3532:         MatRestoreRow(As, i, &temp, &cols, NULL);
3533:       }
3534:       PetscMalloc1(2*numRows, &vtxwgt_g);
3535:       for (i=0; i<size; i++){
3536:         vtxwgt_g[ncon*cumSumVertices[i]] = numExclusivelyOwnedAll[i];
3537:         if (ncon>1) vtxwgt_g[ncon*cumSumVertices[i]+1] = 1;
3538:         for (j=cumSumVertices[i]+1; j<cumSumVertices[i+1]; j++) {
3539:           vtxwgt_g[ncon*j] = 1;
3540:           if (ncon>1) vtxwgt_g[2*j+1] = 0;
3541:         }
3542:       }
3543:       PetscMalloc1(64, &options);
3544:       METIS_SetDefaultOptions(options); /* initialize all defaults */
3545:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
3546:       options[METIS_OPTION_CONTIG] = 1;
3547:       PetscStackPush("METIS_PartGraphKway");
3548:       METIS_PartGraphKway(&numRows, &ncon, xadj_g, adjncy_g, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
3549:       PetscStackPop;
3550:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
3551:       PetscFree(options);
3552:       PetscFree(xadj_g);
3553:       PetscFree(adjncy_g);
3554:       PetscFree(vtxwgt_g);
3555:     }
3556:     PetscFree(numExclusivelyOwnedAll);

3558:     /* Now scatter the parts array. */
3559:     {
3560:       PetscMPIInt *counts, *mpiCumSumVertices;
3561:       PetscMalloc1(size, &counts);
3562:       PetscMalloc1(size+1, &mpiCumSumVertices);
3563:       for(i=0; i<size; i++) {
3564:         PetscMPIIntCast(cumSumVertices[i+1] - cumSumVertices[i], &(counts[i]));
3565:       }
3566:       for(i=0; i<=size; i++) {
3567:         PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
3568:       }
3569:       MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
3570:       PetscFree(counts);
3571:       PetscFree(mpiCumSumVertices);
3572:     }

3574:     PetscFree(partGlobal);
3575:     MatDestroy(&As);
3576:   }

3578:   MatDestroy(&A);
3579:   PetscFree(ubvec);
3580:   PetscFree(tpwgts);

3582:   /* Now rename the result so that the vertex resembling the exclusively owned points stays on the same rank */

3584:   PetscMalloc1(size, &firstVertices);
3585:   PetscMalloc1(size, &renumbering);
3586:   firstVertices[rank] = part[0];
3587:   MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);
3588:   for (i=0; i<size; i++) {
3589:     renumbering[firstVertices[i]] = i;
3590:   }
3591:   for (i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) {
3592:     part[i] = renumbering[part[i]];
3593:   }
3594:   /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
3595:   failed = (PetscInt)(part[0] != rank);
3596:   MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);

3598:   PetscFree(firstVertices);
3599:   PetscFree(renumbering);

3601:   if (failedGlobal > 0) {
3602:     PetscLayoutDestroy(&layout);
3603:     PetscFree(xadj);
3604:     PetscFree(adjncy);
3605:     PetscFree(vtxwgt);
3606:     PetscFree(toBalance);
3607:     PetscFree(isLeaf);
3608:     PetscFree(isNonExclusivelyOwned);
3609:     PetscFree(isExclusivelyOwned);
3610:     PetscFree(part);
3611:     if (viewer) {
3612:       PetscViewerPopFormat(viewer);
3613:       PetscViewerDestroy(&viewer);
3614:     }
3615:     PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3616:     return(0);
3617:   }

3619:   /*Let's check how well we did distributing points*/
3620:   if (viewer) {
3621:     PetscViewerASCIIPrintf(viewer, "Comparing number of owned entities of depth %D on each process before rebalancing, after rebalancing, and after consistency checks.\n", entityDepth);
3622:     PetscViewerASCIIPrintf(viewer, "Initial.     ");
3623:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
3624:     PetscViewerASCIIPrintf(viewer, "Rebalanced.  ");
3625:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
3626:   }

3628:   /* Now check that every vertex is owned by a process that it is actually connected to. */
3629:   for (i=1; i<=numNonExclusivelyOwned; i++) {
3630:     PetscInt loc = 0;
3631:     PetscFindInt(cumSumVertices[part[i]], xadj[i+1]-xadj[i], &adjncy[xadj[i]], &loc);
3632:     /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
3633:     if (loc<0) {
3634:       part[i] = rank;
3635:     }
3636:   }

3638:   /* Let's see how significant the influences of the previous fixing up step was.*/
3639:   if (viewer) {
3640:     PetscViewerASCIIPrintf(viewer, "After.       ");
3641:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
3642:   }

3644:   PetscLayoutDestroy(&layout);
3645:   PetscFree(xadj);
3646:   PetscFree(adjncy);
3647:   PetscFree(vtxwgt);

3649:   /* Almost done, now rewrite the SF to reflect the new ownership. */
3650:   {
3651:     PetscInt *pointsToRewrite;
3652:     PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
3653:     counter = 0;
3654:     for(i=0; i<pEnd-pStart; i++) {
3655:       if (toBalance[i]) {
3656:         if (isNonExclusivelyOwned[i]) {
3657:           pointsToRewrite[counter] = i + pStart;
3658:           counter++;
3659:         }
3660:       }
3661:     }
3662:     DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part+1, degrees);
3663:     PetscFree(pointsToRewrite);
3664:   }

3666:   PetscFree(toBalance);
3667:   PetscFree(isLeaf);
3668:   PetscFree(isNonExclusivelyOwned);
3669:   PetscFree(isExclusivelyOwned);
3670:   PetscFree(part);
3671:   if (viewer) {
3672:     PetscViewerPopFormat(viewer);
3673:     PetscViewerDestroy(&viewer);
3674:   }
3675:   if (success) *success = PETSC_TRUE;
3676:   PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3677: #else
3678:   SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
3679: #endif
3680:   return(0);
3681: }