Actual source code: plexpartition.c

petsc-master 2020-06-03
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 (PetscDefined(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:       }
541:       PetscFree(tmp);
542:     }
543:     if (numVertices) *numVertices = numCells;
544:     if (offsets)   *offsets   = off;
545:     if (adjacency) *adjacency = adj;
546:     return(0);
547:   }
548:   /* Setup face recognition */
549:   if (faceDepth == 1) {
550:     PetscInt cornersSeen[30] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /* Could use PetscBT */

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

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

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

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

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

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

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

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

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

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

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

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

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

651:   Not Collective

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

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

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

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

675:   Level: advanced

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

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

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

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

692:   Collective on PetscPartitioner

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

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

701:   Level: intermediate

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

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

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

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

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

732:   Not Collective

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

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

740:   Level: intermediate

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

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

756:    Collective on PetscPartitioner

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

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

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

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

779:   Collective on PetscPartitioner

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

785:   Level: developer

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

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

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

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

833:   Collective on PetscPartitioner

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

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

843:   Level: developer

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

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

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

881:   Collective on PetscPartitioner

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

886:   Level: developer

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

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

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

903:   Collective on PetscPartitioner

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

908:   Level: developer

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

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

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

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

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

933:   Collective on PetscPartitioner

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

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

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

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

956:   Level: developer

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

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

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

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

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

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

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

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

1033:   Collective

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

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

1041:   Level: beginner

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

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

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

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

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

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

1071:   Collective on PetscPartitioner

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

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

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

1086:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1336:   Level: intermediate

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

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

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

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

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

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

1362:   Collective on PetscPartitioner

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

1370:   Level: developer

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

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

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

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

1405:   Collective on PetscPartitioner

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

1411:   Level: intermediate

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

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

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

1428:   Collective on PetscPartitioner

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

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

1436:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1618:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

1691:   Level: intermediate

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

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

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

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


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

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

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

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

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

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

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

1790:   PetscObjectGetComm((PetscObject)part,&comm);
1791:   if (PetscDefined (USE_DEBUG)) {
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:   if (!numVertices) { /* distributed case, return if not holding the graph */
1801:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
1802:     return(0);
1803:   }
1804:   FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
1805:   for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];

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

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

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

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

1886:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2129:   Level: intermediate

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

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

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

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

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

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

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

2163: #if defined(PETSC_HAVE_PTSCOTCH)

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

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

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

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

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

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

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

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

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

2273: #endif /* PETSC_HAVE_PTSCOTCH */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2451:   Level: intermediate

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

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

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

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

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

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

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

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

2484:   Not collective

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

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

2492:   Level: developer

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

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

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

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

2512:   logically collective on DM

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

2518:   Level: developer

2520:   Note: Any existing PetscPartitioner will be destroyed.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2702:   Level: developer

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

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

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

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

2740:   Level: developer

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

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

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

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

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

2783:   Level: developer

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

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

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

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

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

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

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

2852:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2986:   Level: developer

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

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

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

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

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

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

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

3073:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

3243: #endif

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

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

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

3257:   Level: intermediate

3259: @*/

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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