Actual source code: plexpartition.c

petsc-master 2020-01-15
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) {
972:   }
973:   if (vertexSection) {
974:     PetscInt s,e;

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

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

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

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

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

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

1032:   Collective

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

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

1040:   Level: beginner

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

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

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

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

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

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

1070:   Collective on PetscPartitioner

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

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

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

1085:   Level: developer

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

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

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

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

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

1143:       /* 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) */
1144:       /* We do this only if the local section has been set */
1145:       if (section) {
1146:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
1147:         if (!clSection) {
1148:           DMPlexCreateClosureIndex(dm,NULL);
1149:         }
1150:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
1151:         ISGetIndices(clPoints,&clIdx);
1152:       }
1153:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
1154:       PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
1155:       PetscSectionSetChart(vertSection, 0, numVertices);
1156:       if (globalNumbering) {
1157:         ISGetIndices(globalNumbering,&gid);
1158:       } else gid = NULL;
1159:       for (p = pStart, v = 0; p < pEnd; ++p) {
1160:         PetscInt dof = 1;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1335:   Level: intermediate

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

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

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

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

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

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

1361:   Collective on PetscPartitioner

1363:   Input Parameters:
1364: + part   - The PetscPartitioner
1365: . size   - The number of partitions
1366: . sizes  - array of length size (or NULL) providing the number of points in each partition
1367: - 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.)

1369:   Level: developer

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

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

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

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

1404:   Collective on PetscPartitioner

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

1410:   Level: intermediate

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

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

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

1427:   Collective on PetscPartitioner

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

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

1435:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1617:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

1690:   Level: intermediate

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

1887:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2130:   Level: intermediate

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

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

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

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

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

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

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

2164: #if defined(PETSC_HAVE_PTSCOTCH)

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

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

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

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

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

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

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

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

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

2274: #endif /* PETSC_HAVE_PTSCOTCH */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2452:   Level: intermediate

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

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

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

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

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

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

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

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

2485:   Not collective

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

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

2493:   Level: developer

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

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

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

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

2513:   logically collective on DM

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

2519:   Level: developer

2521:   Note: Any existing PetscPartitioner will be destroyed.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2703:   Level: developer

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

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

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

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

2741:   Level: developer

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

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

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

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

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

2784:   Level: developer

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

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

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

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

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

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

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

2853:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2987:   Level: developer

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

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

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

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

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

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

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

3074:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

3244: #endif

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

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

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

3258:   Level: intermediate

3260: @*/

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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