Actual source code: plexpartition.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/partitionerimpl.h>
  3: #include <petsc/private/hashseti.h>

  5: const char *const DMPlexCSRAlgorithms[] = {"mat", "graph", "overlap", "DMPlexCSRAlgorithm", "DM_PLEX_CSR_", NULL};

  7: static inline PetscInt DMPlex_GlobalID(PetscInt point)
  8: {
  9:   return point >= 0 ? point : -(point + 1);
 10: }

 12: static PetscErrorCode DMPlexCreatePartitionerGraph_Overlap(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
 13: {
 14:   DM              ovdm;
 15:   PetscSF         sfPoint;
 16:   IS              cellNumbering;
 17:   const PetscInt *cellNum;
 18:   PetscInt       *adj = NULL, *vOffsets = NULL, *vAdj = NULL;
 19:   PetscBool       useCone, useClosure;
 20:   PetscInt        dim, depth, overlap, cStart, cEnd, c, v;
 21:   PetscMPIInt     rank, size;

 23:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
 24:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
 25:   DMGetDimension(dm, &dim);
 26:   DMPlexGetDepth(dm, &depth);
 27:   if (dim != depth) {
 28:     /* We do not handle the uninterpolated case here */
 29:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
 30:     /* DMPlexCreateNeighborCSR does not make a numbering */
 31:     if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
 32:     /* Different behavior for empty graphs */
 33:     if (!*numVertices) {
 34:       PetscMalloc1(1, offsets);
 35:       (*offsets)[0] = 0;
 36:     }
 37:     /* Broken in parallel */
 39:     return 0;
 40:   }
 41:   /* Always use FVM adjacency to create partitioner graph */
 42:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
 43:   DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
 44:   /* Need overlap >= 1 */
 45:   DMPlexGetOverlap(dm, &overlap);
 46:   if (size && overlap < 1) {
 47:     DMPlexDistributeOverlap(dm, 1, NULL, &ovdm);
 48:   } else {
 49:     PetscObjectReference((PetscObject)dm);
 50:     ovdm = dm;
 51:   }
 52:   DMGetPointSF(ovdm, &sfPoint);
 53:   DMPlexGetHeightStratum(ovdm, height, &cStart, &cEnd);
 54:   DMPlexCreateNumbering_Plex(ovdm, cStart, cEnd, 0, NULL, sfPoint, &cellNumbering);
 55:   if (globalNumbering) {
 56:     PetscObjectReference((PetscObject)cellNumbering);
 57:     *globalNumbering = cellNumbering;
 58:   }
 59:   ISGetIndices(cellNumbering, &cellNum);
 60:   /* Determine sizes */
 61:   for (*numVertices = 0, c = cStart; c < cEnd; ++c) {
 62:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
 63:     if (cellNum[c - cStart] < 0) continue;
 64:     (*numVertices)++;
 65:   }
 66:   PetscCalloc1(*numVertices + 1, &vOffsets);
 67:   for (c = cStart, v = 0; c < cEnd; ++c) {
 68:     PetscInt adjSize = PETSC_DETERMINE, a, vsize = 0;

 70:     if (cellNum[c - cStart] < 0) continue;
 71:     DMPlexGetAdjacency(ovdm, c, &adjSize, &adj);
 72:     for (a = 0; a < adjSize; ++a) {
 73:       const PetscInt point = adj[a];
 74:       if (point != c && cStart <= point && point < cEnd) ++vsize;
 75:     }
 76:     vOffsets[v + 1] = vOffsets[v] + vsize;
 77:     ++v;
 78:   }
 79:   /* Determine adjacency */
 80:   PetscMalloc1(vOffsets[*numVertices], &vAdj);
 81:   for (c = cStart, v = 0; c < cEnd; ++c) {
 82:     PetscInt adjSize = PETSC_DETERMINE, a, off = vOffsets[v];

 84:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
 85:     if (cellNum[c - cStart] < 0) continue;
 86:     DMPlexGetAdjacency(ovdm, c, &adjSize, &adj);
 87:     for (a = 0; a < adjSize; ++a) {
 88:       const PetscInt point = adj[a];
 89:       if (point != c && cStart <= point && point < cEnd) vAdj[off++] = DMPlex_GlobalID(cellNum[point - cStart]);
 90:     }
 92:     /* Sort adjacencies (not strictly necessary) */
 93:     PetscSortInt(off - vOffsets[v], &vAdj[vOffsets[v]]);
 94:     ++v;
 95:   }
 96:   PetscFree(adj);
 97:   ISRestoreIndices(cellNumbering, &cellNum);
 98:   ISDestroy(&cellNumbering);
 99:   DMSetBasicAdjacency(dm, useCone, useClosure);
100:   DMDestroy(&ovdm);
101:   if (offsets) {
102:     *offsets = vOffsets;
103:   } else PetscFree(vOffsets);
104:   if (adjacency) {
105:     *adjacency = vAdj;
106:   } else PetscFree(vAdj);
107:   return 0;
108: }

110: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
111: {
112:   PetscInt        dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
113:   PetscInt       *adj = NULL, *vOffsets = NULL, *graph = NULL;
114:   IS              cellNumbering;
115:   const PetscInt *cellNum;
116:   PetscBool       useCone, useClosure;
117:   PetscSection    section;
118:   PetscSegBuffer  adjBuffer;
119:   PetscSF         sfPoint;
120:   PetscInt       *adjCells = NULL, *remoteCells = NULL;
121:   const PetscInt *local;
122:   PetscInt        nroots, nleaves, l;
123:   PetscMPIInt     rank;

125:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
126:   DMGetDimension(dm, &dim);
127:   DMPlexGetDepth(dm, &depth);
128:   if (dim != depth) {
129:     /* We do not handle the uninterpolated case here */
130:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
131:     /* DMPlexCreateNeighborCSR does not make a numbering */
132:     if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
133:     /* Different behavior for empty graphs */
134:     if (!*numVertices) {
135:       PetscMalloc1(1, offsets);
136:       (*offsets)[0] = 0;
137:     }
138:     /* Broken in parallel */
140:     return 0;
141:   }
142:   DMGetPointSF(dm, &sfPoint);
143:   DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
144:   /* Build adjacency graph via a section/segbuffer */
145:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section);
146:   PetscSectionSetChart(section, pStart, pEnd);
147:   PetscSegBufferCreate(sizeof(PetscInt), 1000, &adjBuffer);
148:   /* Always use FVM adjacency to create partitioner graph */
149:   DMGetBasicAdjacency(dm, &useCone, &useClosure);
150:   DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
151:   DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
152:   if (globalNumbering) {
153:     PetscObjectReference((PetscObject)cellNumbering);
154:     *globalNumbering = cellNumbering;
155:   }
156:   ISGetIndices(cellNumbering, &cellNum);
157:   /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
158:      Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
159:    */
160:   PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
161:   if (nroots >= 0) {
162:     PetscInt fStart, fEnd, f;

164:     PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
165:     DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd);
166:     for (l = 0; l < nroots; ++l) adjCells[l] = -3;
167:     for (f = fStart; f < fEnd; ++f) {
168:       const PetscInt *support;
169:       PetscInt        supportSize;

171:       DMPlexGetSupport(dm, f, &support);
172:       DMPlexGetSupportSize(dm, f, &supportSize);
173:       if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
174:       else if (supportSize == 2) {
175:         PetscFindInt(support[0], nleaves, local, &p);
176:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
177:         PetscFindInt(support[1], nleaves, local, &p);
178:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
179:       }
180:       /* Handle non-conforming meshes */
181:       if (supportSize > 2) {
182:         PetscInt        numChildren, i;
183:         const PetscInt *children;

185:         DMPlexGetTreeChildren(dm, f, &numChildren, &children);
186:         for (i = 0; i < numChildren; ++i) {
187:           const PetscInt child = children[i];
188:           if (fStart <= child && child < fEnd) {
189:             DMPlexGetSupport(dm, child, &support);
190:             DMPlexGetSupportSize(dm, child, &supportSize);
191:             if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
192:             else if (supportSize == 2) {
193:               PetscFindInt(support[0], nleaves, local, &p);
194:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
195:               PetscFindInt(support[1], nleaves, local, &p);
196:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
197:             }
198:           }
199:         }
200:       }
201:     }
202:     for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
203:     PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE);
204:     PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE);
205:     PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
206:     PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
207:   }
208:   /* Combine local and global adjacencies */
209:   for (*numVertices = 0, p = pStart; p < pEnd; p++) {
210:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
211:     if (nroots > 0) {
212:       if (cellNum[p - pStart] < 0) continue;
213:     }
214:     /* Add remote cells */
215:     if (remoteCells) {
216:       const PetscInt  gp = DMPlex_GlobalID(cellNum[p - pStart]);
217:       PetscInt        coneSize, numChildren, c, i;
218:       const PetscInt *cone, *children;

220:       DMPlexGetCone(dm, p, &cone);
221:       DMPlexGetConeSize(dm, p, &coneSize);
222:       for (c = 0; c < coneSize; ++c) {
223:         const PetscInt point = cone[c];
224:         if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
225:           PetscInt *PETSC_RESTRICT pBuf;
226:           PetscSectionAddDof(section, p, 1);
227:           PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
228:           *pBuf = remoteCells[point];
229:         }
230:         /* Handle non-conforming meshes */
231:         DMPlexGetTreeChildren(dm, point, &numChildren, &children);
232:         for (i = 0; i < numChildren; ++i) {
233:           const PetscInt child = children[i];
234:           if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
235:             PetscInt *PETSC_RESTRICT pBuf;
236:             PetscSectionAddDof(section, p, 1);
237:             PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
238:             *pBuf = remoteCells[child];
239:           }
240:         }
241:       }
242:     }
243:     /* Add local cells */
244:     adjSize = PETSC_DETERMINE;
245:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
246:     for (a = 0; a < adjSize; ++a) {
247:       const PetscInt point = adj[a];
248:       if (point != p && pStart <= point && point < pEnd) {
249:         PetscInt *PETSC_RESTRICT pBuf;
250:         PetscSectionAddDof(section, p, 1);
251:         PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
252:         *pBuf = DMPlex_GlobalID(cellNum[point - pStart]);
253:       }
254:     }
255:     (*numVertices)++;
256:   }
257:   PetscFree(adj);
258:   PetscFree2(adjCells, remoteCells);
259:   DMSetBasicAdjacency(dm, useCone, useClosure);

261:   /* Derive CSR graph from section/segbuffer */
262:   PetscSectionSetUp(section);
263:   PetscSectionGetStorageSize(section, &size);
264:   PetscMalloc1(*numVertices + 1, &vOffsets);
265:   for (idx = 0, p = pStart; p < pEnd; p++) {
266:     if (nroots > 0) {
267:       if (cellNum[p - pStart] < 0) continue;
268:     }
269:     PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
270:   }
271:   vOffsets[*numVertices] = size;
272:   PetscSegBufferExtractAlloc(adjBuffer, &graph);

274:   if (nroots >= 0) {
275:     /* Filter out duplicate edges using section/segbuffer */
276:     PetscSectionReset(section);
277:     PetscSectionSetChart(section, 0, *numVertices);
278:     for (p = 0; p < *numVertices; p++) {
279:       PetscInt start = vOffsets[p], end = vOffsets[p + 1];
280:       PetscInt numEdges = end - start, *PETSC_RESTRICT edges;
281:       PetscSortRemoveDupsInt(&numEdges, &graph[start]);
282:       PetscSectionSetDof(section, p, numEdges);
283:       PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
284:       PetscArraycpy(edges, &graph[start], numEdges);
285:     }
286:     PetscFree(vOffsets);
287:     PetscFree(graph);
288:     /* Derive CSR graph from section/segbuffer */
289:     PetscSectionSetUp(section);
290:     PetscSectionGetStorageSize(section, &size);
291:     PetscMalloc1(*numVertices + 1, &vOffsets);
292:     for (idx = 0, p = 0; p < *numVertices; p++) PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
293:     vOffsets[*numVertices] = size;
294:     PetscSegBufferExtractAlloc(adjBuffer, &graph);
295:   } else {
296:     /* Sort adjacencies (not strictly necessary) */
297:     for (p = 0; p < *numVertices; p++) {
298:       PetscInt start = vOffsets[p], end = vOffsets[p + 1];
299:       PetscSortInt(end - start, &graph[start]);
300:     }
301:   }

303:   if (offsets) *offsets = vOffsets;
304:   if (adjacency) *adjacency = graph;

306:   /* Cleanup */
307:   PetscSegBufferDestroy(&adjBuffer);
308:   PetscSectionDestroy(&section);
309:   ISRestoreIndices(cellNumbering, &cellNum);
310:   ISDestroy(&cellNumbering);
311:   return 0;
312: }

314: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
315: {
316:   Mat             conn, CSR;
317:   IS              fis, cis, cis_own;
318:   PetscSF         sfPoint;
319:   const PetscInt *rows, *cols, *ii, *jj;
320:   PetscInt       *idxs, *idxs2;
321:   PetscInt        dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
322:   PetscMPIInt     rank;
323:   PetscBool       flg;

325:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
326:   DMGetDimension(dm, &dim);
327:   DMPlexGetDepth(dm, &depth);
328:   if (dim != depth) {
329:     /* We do not handle the uninterpolated case here */
330:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
331:     /* DMPlexCreateNeighborCSR does not make a numbering */
332:     if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
333:     /* Different behavior for empty graphs */
334:     if (!*numVertices) {
335:       PetscMalloc1(1, offsets);
336:       (*offsets)[0] = 0;
337:     }
338:     /* Broken in parallel */
340:     return 0;
341:   }
342:   /* Interpolated and parallel case */

344:   /* numbering */
345:   DMGetPointSF(dm, &sfPoint);
346:   DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
347:   DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd);
348:   DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
349:   DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
350:   if (globalNumbering) ISDuplicate(cis, globalNumbering);

352:   /* get positive global ids and local sizes for facets and cells */
353:   ISGetLocalSize(fis, &m);
354:   ISGetIndices(fis, &rows);
355:   PetscMalloc1(m, &idxs);
356:   for (i = 0, floc = 0; i < m; i++) {
357:     const PetscInt p = rows[i];

359:     if (p < 0) {
360:       idxs[i] = -(p + 1);
361:     } else {
362:       idxs[i] = p;
363:       floc += 1;
364:     }
365:   }
366:   ISRestoreIndices(fis, &rows);
367:   ISDestroy(&fis);
368:   ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);

370:   ISGetLocalSize(cis, &m);
371:   ISGetIndices(cis, &cols);
372:   PetscMalloc1(m, &idxs);
373:   PetscMalloc1(m, &idxs2);
374:   for (i = 0, cloc = 0; i < m; i++) {
375:     const PetscInt p = cols[i];

377:     if (p < 0) {
378:       idxs[i] = -(p + 1);
379:     } else {
380:       idxs[i]       = p;
381:       idxs2[cloc++] = p;
382:     }
383:   }
384:   ISRestoreIndices(cis, &cols);
385:   ISDestroy(&cis);
386:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
387:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);

389:   /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
390:   MatCreate(PetscObjectComm((PetscObject)dm), &conn);
391:   MatSetSizes(conn, floc, cloc, M, N);
392:   MatSetType(conn, MATMPIAIJ);
393:   DMPlexGetMaxSizes(dm, NULL, &lm);
394:   MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
395:   MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);

397:   /* Assemble matrix */
398:   ISGetIndices(fis, &rows);
399:   ISGetIndices(cis, &cols);
400:   for (c = cStart; c < cEnd; c++) {
401:     const PetscInt *cone;
402:     PetscInt        coneSize, row, col, f;

404:     col = cols[c - cStart];
405:     DMPlexGetCone(dm, c, &cone);
406:     DMPlexGetConeSize(dm, c, &coneSize);
407:     for (f = 0; f < coneSize; f++) {
408:       const PetscScalar v = 1.0;
409:       const PetscInt   *children;
410:       PetscInt          numChildren, ch;

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

415:       /* non-conforming meshes */
416:       DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
417:       for (ch = 0; ch < numChildren; ch++) {
418:         const PetscInt child = children[ch];

420:         if (child < fStart || child >= fEnd) continue;
421:         row = rows[child - fStart];
422:         MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
423:       }
424:     }
425:   }
426:   ISRestoreIndices(fis, &rows);
427:   ISRestoreIndices(cis, &cols);
428:   ISDestroy(&fis);
429:   ISDestroy(&cis);
430:   MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
431:   MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);

433:   /* Get parallel CSR by doing conn^T * conn */
434:   MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
435:   MatDestroy(&conn);

437:   /* extract local part of the CSR */
438:   MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
439:   MatDestroy(&CSR);
440:   MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);

443:   /* get back requested output */
444:   if (numVertices) *numVertices = m;
445:   if (offsets) {
446:     PetscCalloc1(m + 1, &idxs);
447:     for (i = 1; i < m + 1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
448:     *offsets = idxs;
449:   }
450:   if (adjacency) {
451:     PetscMalloc1(ii[m] - m, &idxs);
452:     ISGetIndices(cis_own, &rows);
453:     for (i = 0, c = 0; i < m; i++) {
454:       PetscInt j, g = rows[i];

456:       for (j = ii[i]; j < ii[i + 1]; j++) {
457:         if (jj[j] == g) continue; /* again, self-connectivity */
458:         idxs[c++] = jj[j];
459:       }
460:     }
462:     ISRestoreIndices(cis_own, &rows);
463:     *adjacency = idxs;
464:   }

466:   /* cleanup */
467:   ISDestroy(&cis_own);
468:   MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
470:   MatDestroy(&conn);
471:   return 0;
472: }

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

477:   Input Parameters:
478: + dm      - The mesh DM dm
479: - height  - Height of the strata from which to construct the graph

481:   Output Parameters:
482: + numVertices     - Number of vertices in the graph
483: . offsets         - Point offsets in the graph
484: . adjacency       - Point connectivity in the graph
485: - 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.

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

490:   Options Database Keys:
491: . -dm_plex_csr_alg <mat,graph,overlap> - Choose the algorithm for computing the CSR graph

493:   Level: developer

495: .seealso: `PetscPartitionerGetType()`, `PetscPartitionerCreate()`, `DMSetAdjacency()`
496: @*/
497: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
498: {
499:   DMPlexCSRAlgorithm alg = DM_PLEX_CSR_GRAPH;

501:   PetscOptionsGetEnum(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *)&alg, NULL);
502:   switch (alg) {
503:   case DM_PLEX_CSR_MAT:
504:     DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
505:     break;
506:   case DM_PLEX_CSR_GRAPH:
507:     DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
508:     break;
509:   case DM_PLEX_CSR_OVERLAP:
510:     DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering);
511:     break;
512:   }
513:   return 0;
514: }

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

519:   Collective on DM

521:   Input Parameters:
522: + dm - The DMPlex
523: - cellHeight - The height of mesh points to treat as cells (default should be 0)

525:   Output Parameters:
526: + numVertices - The number of local vertices in the graph, or cells in the mesh.
527: . offsets     - The offset to the adjacency list for each cell
528: - adjacency   - The adjacency list for all cells

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

532:   Level: advanced

534: .seealso: `DMPlexCreate()`
535: @*/
536: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
537: {
538:   const PetscInt maxFaceCases = 30;
539:   PetscInt       numFaceCases = 0;
540:   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
541:   PetscInt      *off, *adj;
542:   PetscInt      *neighborCells = NULL;
543:   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;

545:   /* For parallel partitioning, I think you have to communicate supports */
546:   DMGetDimension(dm, &dim);
547:   cellDim = dim - cellHeight;
548:   DMPlexGetDepth(dm, &depth);
549:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
550:   if (cEnd - cStart == 0) {
551:     if (numVertices) *numVertices = 0;
552:     if (offsets) *offsets = NULL;
553:     if (adjacency) *adjacency = NULL;
554:     return 0;
555:   }
556:   numCells  = cEnd - cStart;
557:   faceDepth = depth - cellHeight;
558:   if (dim == depth) {
559:     PetscInt f, fStart, fEnd;

561:     PetscCalloc1(numCells + 1, &off);
562:     /* Count neighboring cells */
563:     DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd);
564:     for (f = fStart; f < fEnd; ++f) {
565:       const PetscInt *support;
566:       PetscInt        supportSize;
567:       DMPlexGetSupportSize(dm, f, &supportSize);
568:       DMPlexGetSupport(dm, f, &support);
569:       if (supportSize == 2) {
570:         PetscInt numChildren;

572:         DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
573:         if (!numChildren) {
574:           ++off[support[0] - cStart + 1];
575:           ++off[support[1] - cStart + 1];
576:         }
577:       }
578:     }
579:     /* Prefix sum */
580:     for (c = 1; c <= numCells; ++c) off[c] += off[c - 1];
581:     if (adjacency) {
582:       PetscInt *tmp;

584:       PetscMalloc1(off[numCells], &adj);
585:       PetscMalloc1(numCells + 1, &tmp);
586:       PetscArraycpy(tmp, off, numCells + 1);
587:       /* Get neighboring cells */
588:       for (f = fStart; f < fEnd; ++f) {
589:         const PetscInt *support;
590:         PetscInt        supportSize;
591:         DMPlexGetSupportSize(dm, f, &supportSize);
592:         DMPlexGetSupport(dm, f, &support);
593:         if (supportSize == 2) {
594:           PetscInt numChildren;

596:           DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
597:           if (!numChildren) {
598:             adj[tmp[support[0] - cStart]++] = support[1];
599:             adj[tmp[support[1] - cStart]++] = support[0];
600:           }
601:         }
602:       }
603:       for (c = 0; c < cEnd - cStart; ++c) PetscAssert(tmp[c] == off[c + 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %" PetscInt_FMT " != %" PetscInt_FMT " for cell %" PetscInt_FMT, tmp[c], off[c], c + cStart);
604:       PetscFree(tmp);
605:     }
606:     if (numVertices) *numVertices = numCells;
607:     if (offsets) *offsets = off;
608:     if (adjacency) *adjacency = adj;
609:     return 0;
610:   }
611:   /* Setup face recognition */
612:   if (faceDepth == 1) {
613:     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 */

615:     for (c = cStart; c < cEnd; ++c) {
616:       PetscInt corners;

618:       DMPlexGetConeSize(dm, c, &corners);
619:       if (!cornersSeen[corners]) {
620:         PetscInt nFV;

623:         cornersSeen[corners] = 1;

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

627:         numFaceVertices[numFaceCases++] = nFV;
628:       }
629:     }
630:   }
631:   PetscCalloc1(numCells + 1, &off);
632:   /* Count neighboring cells */
633:   for (cell = cStart; cell < cEnd; ++cell) {
634:     PetscInt numNeighbors = PETSC_DETERMINE, n;

636:     DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
637:     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
638:     for (n = 0; n < numNeighbors; ++n) {
639:       PetscInt        cellPair[2];
640:       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
641:       PetscInt        meetSize = 0;
642:       const PetscInt *meet     = NULL;

644:       cellPair[0] = cell;
645:       cellPair[1] = neighborCells[n];
646:       if (cellPair[0] == cellPair[1]) continue;
647:       if (!found) {
648:         DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
649:         if (meetSize) {
650:           PetscInt f;

652:           for (f = 0; f < numFaceCases; ++f) {
653:             if (numFaceVertices[f] == meetSize) {
654:               found = PETSC_TRUE;
655:               break;
656:             }
657:           }
658:         }
659:         DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
660:       }
661:       if (found) ++off[cell - cStart + 1];
662:     }
663:   }
664:   /* Prefix sum */
665:   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell - 1];

667:   if (adjacency) {
668:     PetscMalloc1(off[numCells], &adj);
669:     /* Get neighboring cells */
670:     for (cell = cStart; cell < cEnd; ++cell) {
671:       PetscInt numNeighbors = PETSC_DETERMINE, n;
672:       PetscInt cellOffset   = 0;

674:       DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
675:       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
676:       for (n = 0; n < numNeighbors; ++n) {
677:         PetscInt        cellPair[2];
678:         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
679:         PetscInt        meetSize = 0;
680:         const PetscInt *meet     = NULL;

682:         cellPair[0] = cell;
683:         cellPair[1] = neighborCells[n];
684:         if (cellPair[0] == cellPair[1]) continue;
685:         if (!found) {
686:           DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
687:           if (meetSize) {
688:             PetscInt f;

690:             for (f = 0; f < numFaceCases; ++f) {
691:               if (numFaceVertices[f] == meetSize) {
692:                 found = PETSC_TRUE;
693:                 break;
694:               }
695:             }
696:           }
697:           DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
698:         }
699:         if (found) {
700:           adj[off[cell - cStart] + cellOffset] = neighborCells[n];
701:           ++cellOffset;
702:         }
703:       }
704:     }
705:   }
706:   PetscFree(neighborCells);
707:   if (numVertices) *numVertices = numCells;
708:   if (offsets) *offsets = off;
709:   if (adjacency) *adjacency = adj;
710:   return 0;
711: }

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

716:   Collective on PetscPartitioner

718:   Input Parameters:
719: + part    - The PetscPartitioner
720: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
721: - dm      - The mesh DM

723:   Output Parameters:
724: + partSection     - The PetscSection giving the division of points by partition
725: - partition       - The list of points by partition

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

731:   Level: developer

733: .seealso `DMPlexDistribute()`, `PetscPartitionerCreate()`, `PetscSectionCreate()`, `PetscSectionSetChart()`, `PetscPartitionerPartition()`
734: @*/
735: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
736: {
737:   PetscMPIInt  size;
738:   PetscBool    isplex;
739:   PetscSection vertSection = NULL;

746:   PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex);
748:   MPI_Comm_size(PetscObjectComm((PetscObject)part), &size);
749:   if (size == 1) {
750:     PetscInt *points;
751:     PetscInt  cStart, cEnd, c;

753:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
754:     PetscSectionReset(partSection);
755:     PetscSectionSetChart(partSection, 0, size);
756:     PetscSectionSetDof(partSection, 0, cEnd - cStart);
757:     PetscSectionSetUp(partSection);
758:     PetscMalloc1(cEnd - cStart, &points);
759:     for (c = cStart; c < cEnd; ++c) points[c] = c;
760:     ISCreateGeneral(PetscObjectComm((PetscObject)part), cEnd - cStart, points, PETSC_OWN_POINTER, partition);
761:     return 0;
762:   }
763:   if (part->height == 0) {
764:     PetscInt  numVertices = 0;
765:     PetscInt *start       = NULL;
766:     PetscInt *adjacency   = NULL;
767:     IS        globalNumbering;

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

775:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
776:       DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
777:       ISGetIndices(globalNumbering, &idxs);
778:       for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
779:       ISRestoreIndices(globalNumbering, &idxs);
780:     }
781:     if (part->usevwgt) {
782:       PetscSection    section = dm->localSection, clSection = NULL;
783:       IS              clPoints = NULL;
784:       const PetscInt *gid, *clIdx;
785:       PetscInt        v, p, pStart, pEnd;

787:       /* 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) */
788:       /* We do this only if the local section has been set */
789:       if (section) {
790:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
791:         if (!clSection) DMPlexCreateClosureIndex(dm, NULL);
792:         PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
793:         ISGetIndices(clPoints, &clIdx);
794:       }
795:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
796:       PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
797:       PetscSectionSetChart(vertSection, 0, numVertices);
798:       if (globalNumbering) {
799:         ISGetIndices(globalNumbering, &gid);
800:       } else gid = NULL;
801:       for (p = pStart, v = 0; p < pEnd; ++p) {
802:         PetscInt dof = 1;

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

807:         if (section) {
808:           PetscInt cl, clSize, clOff;

810:           dof = 0;
811:           PetscSectionGetDof(clSection, p, &clSize);
812:           PetscSectionGetOffset(clSection, p, &clOff);
813:           for (cl = 0; cl < clSize; cl += 2) {
814:             PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */

816:             PetscSectionGetDof(section, clPoint, &clDof);
817:             dof += clDof;
818:           }
819:         }
821:         PetscSectionSetDof(vertSection, v, dof);
822:         v++;
823:       }
824:       if (globalNumbering) ISRestoreIndices(globalNumbering, &gid);
825:       if (clPoints) ISRestoreIndices(clPoints, &clIdx);
826:       PetscSectionSetUp(vertSection);
827:     }
828:     PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
829:     PetscFree(start);
830:     PetscFree(adjacency);
831:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
832:       const PetscInt *globalNum;
833:       const PetscInt *partIdx;
834:       PetscInt       *map, cStart, cEnd;
835:       PetscInt       *adjusted, i, localSize, offset;
836:       IS              newPartition;

838:       ISGetLocalSize(*partition, &localSize);
839:       PetscMalloc1(localSize, &adjusted);
840:       ISGetIndices(globalNumbering, &globalNum);
841:       ISGetIndices(*partition, &partIdx);
842:       PetscMalloc1(localSize, &map);
843:       DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
844:       for (i = cStart, offset = 0; i < cEnd; i++) {
845:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
846:       }
847:       for (i = 0; i < localSize; i++) adjusted[i] = map[partIdx[i]];
848:       PetscFree(map);
849:       ISRestoreIndices(*partition, &partIdx);
850:       ISRestoreIndices(globalNumbering, &globalNum);
851:       ISCreateGeneral(PETSC_COMM_SELF, localSize, adjusted, PETSC_OWN_POINTER, &newPartition);
852:       ISDestroy(&globalNumbering);
853:       ISDestroy(partition);
854:       *partition = newPartition;
855:     }
856:   } else SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %" PetscInt_FMT " for points to partition", part->height);
857:   PetscSectionDestroy(&vertSection);
858:   return 0;
859: }

861: /*@
862:   DMPlexGetPartitioner - Get the mesh partitioner

864:   Not collective

866:   Input Parameter:
867: . dm - The DM

869:   Output Parameter:
870: . part - The PetscPartitioner

872:   Level: developer

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

876: .seealso `DMPlexDistribute()`, `DMPlexSetPartitioner()`, `PetscPartitionerDMPlexPartition()`, `PetscPartitionerCreate()`
877: @*/
878: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
879: {
880:   DM_Plex *mesh = (DM_Plex *)dm->data;

884:   *part = mesh->partitioner;
885:   return 0;
886: }

888: /*@
889:   DMPlexSetPartitioner - Set the mesh partitioner

891:   logically collective on DM

893:   Input Parameters:
894: + dm - The DM
895: - part - The partitioner

897:   Level: developer

899:   Note: Any existing PetscPartitioner will be destroyed.

901: .seealso `DMPlexDistribute()`, `DMPlexGetPartitioner()`, `PetscPartitionerCreate()`
902: @*/
903: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
904: {
905:   DM_Plex *mesh = (DM_Plex *)dm->data;

909:   PetscObjectReference((PetscObject)part);
910:   PetscPartitionerDestroy(&mesh->partitioner);
911:   mesh->partitioner = part;
912:   return 0;
913: }

915: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
916: {
917:   const PetscInt *cone;
918:   PetscInt        coneSize, c;
919:   PetscBool       missing;

922:   PetscHSetIQueryAdd(ht, point, &missing);
923:   if (missing) {
924:     DMPlexGetCone(dm, point, &cone);
925:     DMPlexGetConeSize(dm, point, &coneSize);
926:     for (c = 0; c < coneSize; c++) DMPlexAddClosure_Private(dm, ht, cone[c]);
927:   }
928:   return 0;
929: }

931: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
932: {
933:   if (up) {
934:     PetscInt parent;

936:     DMPlexGetTreeParent(dm, point, &parent, NULL);
937:     if (parent != point) {
938:       PetscInt closureSize, *closure = NULL, i;

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

944:         PetscHSetIAdd(ht, cpoint);
945:         DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_TRUE, PETSC_FALSE);
946:       }
947:       DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
948:     }
949:   }
950:   if (down) {
951:     PetscInt        numChildren;
952:     const PetscInt *children;

954:     DMPlexGetTreeChildren(dm, point, &numChildren, &children);
955:     if (numChildren) {
956:       PetscInt i;

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

961:         PetscHSetIAdd(ht, cpoint);
962:         DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_FALSE, PETSC_TRUE);
963:       }
964:     }
965:   }
966:   return 0;
967: }

969: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
970: {
971:   PetscInt parent;

974:   DMPlexGetTreeParent(dm, point, &parent, NULL);
975:   if (point != parent) {
976:     const PetscInt *cone;
977:     PetscInt        coneSize, c;

979:     DMPlexAddClosureTree_Up_Private(dm, ht, parent);
980:     DMPlexAddClosure_Private(dm, ht, parent);
981:     DMPlexGetCone(dm, parent, &cone);
982:     DMPlexGetConeSize(dm, parent, &coneSize);
983:     for (c = 0; c < coneSize; c++) {
984:       const PetscInt cp = cone[c];

986:       DMPlexAddClosureTree_Up_Private(dm, ht, cp);
987:     }
988:   }
989:   return 0;
990: }

992: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
993: {
994:   PetscInt        i, numChildren;
995:   const PetscInt *children;

998:   DMPlexGetTreeChildren(dm, point, &numChildren, &children);
999:   for (i = 0; i < numChildren; i++) PetscHSetIAdd(ht, children[i]);
1000:   return 0;
1001: }

1003: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
1004: {
1005:   const PetscInt *cone;
1006:   PetscInt        coneSize, c;

1009:   PetscHSetIAdd(ht, point);
1010:   DMPlexAddClosureTree_Up_Private(dm, ht, point);
1011:   DMPlexAddClosureTree_Down_Private(dm, ht, point);
1012:   DMPlexGetCone(dm, point, &cone);
1013:   DMPlexGetConeSize(dm, point, &coneSize);
1014:   for (c = 0; c < coneSize; c++) DMPlexAddClosureTree_Private(dm, ht, cone[c]);
1015:   return 0;
1016: }

1018: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
1019: {
1020:   DM_Plex        *mesh    = (DM_Plex *)dm->data;
1021:   const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
1022:   PetscInt        nelems, *elems, off = 0, p;
1023:   PetscHSetI      ht = NULL;

1025:   PetscHSetICreate(&ht);
1026:   PetscHSetIResize(ht, numPoints * 16);
1027:   if (!hasTree) {
1028:     for (p = 0; p < numPoints; ++p) DMPlexAddClosure_Private(dm, ht, points[p]);
1029:   } else {
1030: #if 1
1031:     for (p = 0; p < numPoints; ++p) DMPlexAddClosureTree_Private(dm, ht, points[p]);
1032: #else
1033:     PetscInt *closure = NULL, closureSize, c;
1034:     for (p = 0; p < numPoints; ++p) {
1035:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
1036:       for (c = 0; c < closureSize * 2; c += 2) {
1037:         PetscHSetIAdd(ht, closure[c]);
1038:         if (hasTree) DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);
1039:       }
1040:     }
1041:     if (closure) DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);
1042: #endif
1043:   }
1044:   PetscHSetIGetSize(ht, &nelems);
1045:   PetscMalloc1(nelems, &elems);
1046:   PetscHSetIGetElems(ht, &off, elems);
1047:   PetscHSetIDestroy(&ht);
1048:   PetscSortInt(nelems, elems);
1049:   ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
1050:   return 0;
1051: }

1053: /*@
1054:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

1056:   Input Parameters:
1057: + dm     - The DM
1058: - label  - DMLabel assigning ranks to remote roots

1060:   Level: developer

1062: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1063: @*/
1064: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
1065: {
1066:   IS              rankIS, pointIS, closureIS;
1067:   const PetscInt *ranks, *points;
1068:   PetscInt        numRanks, numPoints, r;

1070:   DMLabelGetValueIS(label, &rankIS);
1071:   ISGetLocalSize(rankIS, &numRanks);
1072:   ISGetIndices(rankIS, &ranks);
1073:   for (r = 0; r < numRanks; ++r) {
1074:     const PetscInt rank = ranks[r];
1075:     DMLabelGetStratumIS(label, rank, &pointIS);
1076:     ISGetLocalSize(pointIS, &numPoints);
1077:     ISGetIndices(pointIS, &points);
1078:     DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
1079:     ISRestoreIndices(pointIS, &points);
1080:     ISDestroy(&pointIS);
1081:     DMLabelSetStratumIS(label, rank, closureIS);
1082:     ISDestroy(&closureIS);
1083:   }
1084:   ISRestoreIndices(rankIS, &ranks);
1085:   ISDestroy(&rankIS);
1086:   return 0;
1087: }

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

1092:   Input Parameters:
1093: + dm     - The DM
1094: - label  - DMLabel assigning ranks to remote roots

1096:   Level: developer

1098: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1099: @*/
1100: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1101: {
1102:   IS              rankIS, pointIS;
1103:   const PetscInt *ranks, *points;
1104:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
1105:   PetscInt       *adj = NULL;

1107:   DMLabelGetValueIS(label, &rankIS);
1108:   ISGetLocalSize(rankIS, &numRanks);
1109:   ISGetIndices(rankIS, &ranks);
1110:   for (r = 0; r < numRanks; ++r) {
1111:     const PetscInt rank = ranks[r];

1113:     DMLabelGetStratumIS(label, rank, &pointIS);
1114:     ISGetLocalSize(pointIS, &numPoints);
1115:     ISGetIndices(pointIS, &points);
1116:     for (p = 0; p < numPoints; ++p) {
1117:       adjSize = PETSC_DETERMINE;
1118:       DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
1119:       for (a = 0; a < adjSize; ++a) DMLabelSetValue(label, adj[a], rank);
1120:     }
1121:     ISRestoreIndices(pointIS, &points);
1122:     ISDestroy(&pointIS);
1123:   }
1124:   ISRestoreIndices(rankIS, &ranks);
1125:   ISDestroy(&rankIS);
1126:   PetscFree(adj);
1127:   return 0;
1128: }

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

1133:   Input Parameters:
1134: + dm     - The DM
1135: - label  - DMLabel assigning ranks to remote roots

1137:   Level: developer

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

1142: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1143: @*/
1144: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1145: {
1146:   MPI_Comm        comm;
1147:   PetscMPIInt     rank;
1148:   PetscSF         sfPoint;
1149:   DMLabel         lblRoots, lblLeaves;
1150:   IS              rankIS, pointIS;
1151:   const PetscInt *ranks;
1152:   PetscInt        numRanks, r;

1154:   PetscObjectGetComm((PetscObject)dm, &comm);
1155:   MPI_Comm_rank(comm, &rank);
1156:   DMGetPointSF(dm, &sfPoint);
1157:   /* Pull point contributions from remote leaves into local roots */
1158:   DMLabelGather(label, sfPoint, &lblLeaves);
1159:   DMLabelGetValueIS(lblLeaves, &rankIS);
1160:   ISGetLocalSize(rankIS, &numRanks);
1161:   ISGetIndices(rankIS, &ranks);
1162:   for (r = 0; r < numRanks; ++r) {
1163:     const PetscInt remoteRank = ranks[r];
1164:     if (remoteRank == rank) continue;
1165:     DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
1166:     DMLabelInsertIS(label, pointIS, remoteRank);
1167:     ISDestroy(&pointIS);
1168:   }
1169:   ISRestoreIndices(rankIS, &ranks);
1170:   ISDestroy(&rankIS);
1171:   DMLabelDestroy(&lblLeaves);
1172:   /* Push point contributions from roots into remote leaves */
1173:   DMLabelDistribute(label, sfPoint, &lblRoots);
1174:   DMLabelGetValueIS(lblRoots, &rankIS);
1175:   ISGetLocalSize(rankIS, &numRanks);
1176:   ISGetIndices(rankIS, &ranks);
1177:   for (r = 0; r < numRanks; ++r) {
1178:     const PetscInt remoteRank = ranks[r];
1179:     if (remoteRank == rank) continue;
1180:     DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
1181:     DMLabelInsertIS(label, pointIS, remoteRank);
1182:     ISDestroy(&pointIS);
1183:   }
1184:   ISRestoreIndices(rankIS, &ranks);
1185:   ISDestroy(&rankIS);
1186:   DMLabelDestroy(&lblRoots);
1187:   return 0;
1188: }

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

1193:   Input Parameters:
1194: + dm        - The DM
1195: . rootLabel - DMLabel assigning ranks to local roots
1196: - processSF - A star forest mapping into the local index on each remote rank

1198:   Output Parameter:
1199: . leafLabel - DMLabel assigning ranks to remote roots

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

1204:   Level: developer

1206: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1207: @*/
1208: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1209: {
1210:   MPI_Comm           comm;
1211:   PetscMPIInt        rank, size, r;
1212:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1213:   PetscSF            sfPoint;
1214:   PetscSection       rootSection;
1215:   PetscSFNode       *rootPoints, *leafPoints;
1216:   const PetscSFNode *remote;
1217:   const PetscInt    *local, *neighbors;
1218:   IS                 valueIS;
1219:   PetscBool          mpiOverflow = PETSC_FALSE;

1221:   PetscLogEventBegin(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1222:   PetscObjectGetComm((PetscObject)dm, &comm);
1223:   MPI_Comm_rank(comm, &rank);
1224:   MPI_Comm_size(comm, &size);
1225:   DMGetPointSF(dm, &sfPoint);

1227:   /* Convert to (point, rank) and use actual owners */
1228:   PetscSectionCreate(comm, &rootSection);
1229:   PetscSectionSetChart(rootSection, 0, size);
1230:   DMLabelGetValueIS(rootLabel, &valueIS);
1231:   ISGetLocalSize(valueIS, &numNeighbors);
1232:   ISGetIndices(valueIS, &neighbors);
1233:   for (n = 0; n < numNeighbors; ++n) {
1234:     DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
1235:     PetscSectionAddDof(rootSection, neighbors[n], numPoints);
1236:   }
1237:   PetscSectionSetUp(rootSection);
1238:   PetscSectionGetStorageSize(rootSection, &rootSize);
1239:   PetscMalloc1(rootSize, &rootPoints);
1240:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
1241:   for (n = 0; n < numNeighbors; ++n) {
1242:     IS              pointIS;
1243:     const PetscInt *points;

1245:     PetscSectionGetOffset(rootSection, neighbors[n], &off);
1246:     DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
1247:     ISGetLocalSize(pointIS, &numPoints);
1248:     ISGetIndices(pointIS, &points);
1249:     for (p = 0; p < numPoints; ++p) {
1250:       if (local) PetscFindInt(points[p], nleaves, local, &l);
1251:       else l = -1;
1252:       if (l >= 0) {
1253:         rootPoints[off + p] = remote[l];
1254:       } else {
1255:         rootPoints[off + p].index = points[p];
1256:         rootPoints[off + p].rank  = rank;
1257:       }
1258:     }
1259:     ISRestoreIndices(pointIS, &points);
1260:     ISDestroy(&pointIS);
1261:   }

1263:   /* Try to communicate overlap using All-to-All */
1264:   if (!processSF) {
1265:     PetscInt64   counter     = 0;
1266:     PetscBool    locOverflow = PETSC_FALSE;
1267:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

1269:     PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
1270:     for (n = 0; n < numNeighbors; ++n) {
1271:       PetscSectionGetDof(rootSection, neighbors[n], &dof);
1272:       PetscSectionGetOffset(rootSection, neighbors[n], &off);
1273: #if defined(PETSC_USE_64BIT_INDICES)
1274:       if (dof > PETSC_MPI_INT_MAX) {
1275:         locOverflow = PETSC_TRUE;
1276:         break;
1277:       }
1278:       if (off > PETSC_MPI_INT_MAX) {
1279:         locOverflow = PETSC_TRUE;
1280:         break;
1281:       }
1282: #endif
1283:       scounts[neighbors[n]] = (PetscMPIInt)dof;
1284:       sdispls[neighbors[n]] = (PetscMPIInt)off;
1285:     }
1286:     MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
1287:     for (r = 0; r < size; ++r) {
1288:       rdispls[r] = (int)counter;
1289:       counter += rcounts[r];
1290:     }
1291:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
1292:     MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
1293:     if (!mpiOverflow) {
1294:       PetscInfo(dm, "Using Alltoallv for mesh distribution\n");
1295:       leafSize = (PetscInt)counter;
1296:       PetscMalloc1(leafSize, &leafPoints);
1297:       MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
1298:     }
1299:     PetscFree4(scounts, sdispls, rcounts, rdispls);
1300:   }

1302:   /* Communicate overlap using process star forest */
1303:   if (processSF || mpiOverflow) {
1304:     PetscSF      procSF;
1305:     PetscSection leafSection;

1307:     if (processSF) {
1308:       PetscInfo(dm, "Using processSF for mesh distribution\n");
1309:       PetscObjectReference((PetscObject)processSF);
1310:       procSF = processSF;
1311:     } else {
1312:       PetscInfo(dm, "Using processSF for mesh distribution (MPI overflow)\n");
1313:       PetscSFCreate(comm, &procSF);
1314:       PetscSFSetGraphWithPattern(procSF, NULL, PETSCSF_PATTERN_ALLTOALL);
1315:     }

1317:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
1318:     DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void **)&leafPoints);
1319:     PetscSectionGetStorageSize(leafSection, &leafSize);
1320:     PetscSectionDestroy(&leafSection);
1321:     PetscSFDestroy(&procSF);
1322:   }

1324:   for (p = 0; p < leafSize; p++) DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);

1326:   ISRestoreIndices(valueIS, &neighbors);
1327:   ISDestroy(&valueIS);
1328:   PetscSectionDestroy(&rootSection);
1329:   PetscFree(rootPoints);
1330:   PetscFree(leafPoints);
1331:   PetscLogEventEnd(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1332:   return 0;
1333: }

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

1338:   Input Parameters:
1339: + dm    - The DM
1340: - label - DMLabel assigning ranks to remote roots

1342:   Output Parameter:
1343: . sf    - The star forest communication context encapsulating the defined mapping

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

1347:   Level: developer

1349: .seealso: `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1350: @*/
1351: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
1352: {
1353:   PetscMPIInt     rank;
1354:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1355:   PetscSFNode    *remotePoints;
1356:   IS              remoteRootIS, neighborsIS;
1357:   const PetscInt *remoteRoots, *neighbors;

1359:   PetscLogEventBegin(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1360:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);

1362:   DMLabelGetValueIS(label, &neighborsIS);
1363: #if 0
1364:   {
1365:     IS is;
1366:     ISDuplicate(neighborsIS, &is);
1367:     ISSort(is);
1368:     ISDestroy(&neighborsIS);
1369:     neighborsIS = is;
1370:   }
1371: #endif
1372:   ISGetLocalSize(neighborsIS, &nNeighbors);
1373:   ISGetIndices(neighborsIS, &neighbors);
1374:   for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1375:     DMLabelGetStratumSize(label, neighbors[n], &numPoints);
1376:     numRemote += numPoints;
1377:   }
1378:   PetscMalloc1(numRemote, &remotePoints);
1379:   /* Put owned points first */
1380:   DMLabelGetStratumSize(label, rank, &numPoints);
1381:   if (numPoints > 0) {
1382:     DMLabelGetStratumIS(label, rank, &remoteRootIS);
1383:     ISGetIndices(remoteRootIS, &remoteRoots);
1384:     for (p = 0; p < numPoints; p++) {
1385:       remotePoints[idx].index = remoteRoots[p];
1386:       remotePoints[idx].rank  = rank;
1387:       idx++;
1388:     }
1389:     ISRestoreIndices(remoteRootIS, &remoteRoots);
1390:     ISDestroy(&remoteRootIS);
1391:   }
1392:   /* Now add remote points */
1393:   for (n = 0; n < nNeighbors; ++n) {
1394:     const PetscInt nn = neighbors[n];

1396:     DMLabelGetStratumSize(label, nn, &numPoints);
1397:     if (nn == rank || numPoints <= 0) continue;
1398:     DMLabelGetStratumIS(label, nn, &remoteRootIS);
1399:     ISGetIndices(remoteRootIS, &remoteRoots);
1400:     for (p = 0; p < numPoints; p++) {
1401:       remotePoints[idx].index = remoteRoots[p];
1402:       remotePoints[idx].rank  = nn;
1403:       idx++;
1404:     }
1405:     ISRestoreIndices(remoteRootIS, &remoteRoots);
1406:     ISDestroy(&remoteRootIS);
1407:   }
1408:   PetscSFCreate(PetscObjectComm((PetscObject)dm), sf);
1409:   DMPlexGetChart(dm, &pStart, &pEnd);
1410:   PetscSFSetGraph(*sf, pEnd - pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1411:   PetscSFSetUp(*sf);
1412:   ISDestroy(&neighborsIS);
1413:   PetscLogEventEnd(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1414:   return 0;
1415: }

1417: #if defined(PETSC_HAVE_PARMETIS)
1418:   #include <parmetis.h>
1419: #endif

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

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

1429:   Input parameters:
1430: + dm                - The DMPlex object.
1431: . n                 - The number of points.
1432: . pointsToRewrite   - The points in the SF whose ownership will change.
1433: . targetOwners      - New owner for each element in pointsToRewrite.
1434: - degrees           - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.

1436:   Level: developer

1438: @*/
1439: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1440: {
1441:   PetscInt           pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1442:   PetscInt          *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1443:   PetscSFNode       *leafLocationsNew;
1444:   const PetscSFNode *iremote;
1445:   const PetscInt    *ilocal;
1446:   PetscBool         *isLeaf;
1447:   PetscSF            sf;
1448:   MPI_Comm           comm;
1449:   PetscMPIInt        rank, size;

1451:   PetscObjectGetComm((PetscObject)dm, &comm);
1452:   MPI_Comm_rank(comm, &rank);
1453:   MPI_Comm_size(comm, &size);
1454:   DMPlexGetChart(dm, &pStart, &pEnd);

1456:   DMGetPointSF(dm, &sf);
1457:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1458:   PetscMalloc1(pEnd - pStart, &isLeaf);
1459:   for (i = 0; i < pEnd - pStart; i++) isLeaf[i] = PETSC_FALSE;
1460:   for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;

1462:   PetscMalloc1(pEnd - pStart + 1, &cumSumDegrees);
1463:   cumSumDegrees[0] = 0;
1464:   for (i = 1; i <= pEnd - pStart; i++) cumSumDegrees[i] = cumSumDegrees[i - 1] + degrees[i - 1];
1465:   sumDegrees = cumSumDegrees[pEnd - pStart];
1466:   /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */

1468:   PetscMalloc1(sumDegrees, &locationsOfLeafs);
1469:   PetscMalloc1(pEnd - pStart, &rankOnLeafs);
1470:   for (i = 0; i < pEnd - pStart; i++) rankOnLeafs[i] = rank;
1471:   PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1472:   PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1473:   PetscFree(rankOnLeafs);

1475:   /* get the remote local points of my leaves */
1476:   PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
1477:   PetscMalloc1(pEnd - pStart, &points);
1478:   for (i = 0; i < pEnd - pStart; i++) points[i] = pStart + i;
1479:   PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1480:   PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1481:   PetscFree(points);
1482:   /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1483:   PetscMalloc1(pEnd - pStart, &newOwners);
1484:   PetscMalloc1(pEnd - pStart, &newNumbers);
1485:   for (i = 0; i < pEnd - pStart; i++) {
1486:     newOwners[i]  = -1;
1487:     newNumbers[i] = -1;
1488:   }
1489:   {
1490:     PetscInt oldNumber, newNumber, oldOwner, newOwner;
1491:     for (i = 0; i < n; i++) {
1492:       oldNumber = pointsToRewrite[i];
1493:       newNumber = -1;
1494:       oldOwner  = rank;
1495:       newOwner  = targetOwners[i];
1496:       if (oldOwner == newOwner) {
1497:         newNumber = oldNumber;
1498:       } else {
1499:         for (j = 0; j < degrees[oldNumber]; j++) {
1500:           if (locationsOfLeafs[cumSumDegrees[oldNumber] + j] == newOwner) {
1501:             newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber] + j];
1502:             break;
1503:           }
1504:         }
1505:       }

1508:       newOwners[oldNumber]  = newOwner;
1509:       newNumbers[oldNumber] = newNumber;
1510:     }
1511:   }
1512:   PetscFree(cumSumDegrees);
1513:   PetscFree(locationsOfLeafs);
1514:   PetscFree(remoteLocalPointOfLeafs);

1516:   PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1517:   PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1518:   PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);
1519:   PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);

1521:   /* Now count how many leafs we have on each processor. */
1522:   leafCounter = 0;
1523:   for (i = 0; i < pEnd - pStart; i++) {
1524:     if (newOwners[i] >= 0) {
1525:       if (newOwners[i] != rank) leafCounter++;
1526:     } else {
1527:       if (isLeaf[i]) leafCounter++;
1528:     }
1529:   }

1531:   /* Now set up the new sf by creating the leaf arrays */
1532:   PetscMalloc1(leafCounter, &leafsNew);
1533:   PetscMalloc1(leafCounter, &leafLocationsNew);

1535:   leafCounter = 0;
1536:   counter     = 0;
1537:   for (i = 0; i < pEnd - pStart; i++) {
1538:     if (newOwners[i] >= 0) {
1539:       if (newOwners[i] != rank) {
1540:         leafsNew[leafCounter]               = i;
1541:         leafLocationsNew[leafCounter].rank  = newOwners[i];
1542:         leafLocationsNew[leafCounter].index = newNumbers[i];
1543:         leafCounter++;
1544:       }
1545:     } else {
1546:       if (isLeaf[i]) {
1547:         leafsNew[leafCounter]               = i;
1548:         leafLocationsNew[leafCounter].rank  = iremote[counter].rank;
1549:         leafLocationsNew[leafCounter].index = iremote[counter].index;
1550:         leafCounter++;
1551:       }
1552:     }
1553:     if (isLeaf[i]) counter++;
1554:   }

1556:   PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
1557:   PetscFree(newOwners);
1558:   PetscFree(newNumbers);
1559:   PetscFree(isLeaf);
1560:   return 0;
1561: }

1563: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1564: {
1565:   PetscInt   *distribution, min, max, sum;
1566:   PetscMPIInt rank, size;

1568:   MPI_Comm_size(comm, &size);
1569:   MPI_Comm_rank(comm, &rank);
1570:   PetscCalloc1(size, &distribution);
1571:   for (PetscInt i = 0; i < n; i++) {
1572:     if (part) distribution[part[i]] += vtxwgt[skip * i];
1573:     else distribution[rank] += vtxwgt[skip * i];
1574:   }
1575:   MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
1576:   min = distribution[0];
1577:   max = distribution[0];
1578:   sum = distribution[0];
1579:   for (PetscInt i = 1; i < size; i++) {
1580:     if (distribution[i] < min) min = distribution[i];
1581:     if (distribution[i] > max) max = distribution[i];
1582:     sum += distribution[i];
1583:   }
1584:   PetscViewerASCIIPrintf(viewer, "Min: %" PetscInt_FMT ", Avg: %" PetscInt_FMT ", Max: %" PetscInt_FMT ", Balance: %f\n", min, sum / size, max, (max * 1. * size) / sum);
1585:   PetscFree(distribution);
1586:   return 0;
1587: }

1589: #endif

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

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

1600:   Output parameters:
1601: . success          - whether the graph partitioning was successful or not, optional. Unsuccessful simply means no change to the partitioning

1603:   Options Database:
1604: +  -dm_plex_rebalance_shared_points_parmetis - Use ParMetis instead of Metis for the partitioner
1605: .  -dm_plex_rebalance_shared_points_use_initial_guess - Use current partition to bootstrap ParMetis partition
1606: .  -dm_plex_rebalance_shared_points_use_mat_partitioning - Use the MatPartitioning object to perform the partition, the prefix for those operations is -dm_plex_rebalance_shared_points_
1607: -  -dm_plex_rebalance_shared_points_monitor - Monitor the shared points rebalance process

1609:   Developer Notes:
1610:   This should use MatPartitioning to allow the use of any partitioner and not be hardwired to use ParMetis

1612:   Level: intermediate
1613: @*/
1614: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1615: {
1616: #if defined(PETSC_HAVE_PARMETIS)
1617:   PetscSF            sf;
1618:   PetscInt           ierr, i, j, idx, jdx;
1619:   PetscInt           eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1620:   const PetscInt    *degrees, *ilocal;
1621:   const PetscSFNode *iremote;
1622:   PetscBool         *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1623:   PetscInt           numExclusivelyOwned, numNonExclusivelyOwned;
1624:   PetscMPIInt        rank, size;
1625:   PetscInt          *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1626:   const PetscInt    *cumSumVertices;
1627:   PetscInt           offset, counter;
1628:   PetscInt          *vtxwgt;
1629:   const PetscInt    *xadj, *adjncy;
1630:   PetscInt          *part, *options;
1631:   PetscInt           nparts, wgtflag, numflag, ncon, edgecut;
1632:   real_t            *ubvec;
1633:   PetscInt          *firstVertices, *renumbering;
1634:   PetscInt           failed, failedGlobal;
1635:   MPI_Comm           comm;
1636:   Mat                A;
1637:   PetscViewer        viewer;
1638:   PetscViewerFormat  format;
1639:   PetscLayout        layout;
1640:   real_t            *tpwgts;
1641:   PetscMPIInt       *counts, *mpiCumSumVertices;
1642:   PetscInt          *pointsToRewrite;
1643:   PetscInt           numRows;
1644:   PetscBool          done, usematpartitioning = PETSC_FALSE;
1645:   IS                 ispart = NULL;
1646:   MatPartitioning    mp;
1647:   const char        *prefix;

1649:   PetscObjectGetComm((PetscObject)dm, &comm);
1650:   MPI_Comm_size(comm, &size);
1651:   if (size == 1) {
1652:     if (success) *success = PETSC_TRUE;
1653:     return 0;
1654:   }
1655:   if (success) *success = PETSC_FALSE;
1656:   MPI_Comm_rank(comm, &rank);

1658:   parallel        = PETSC_FALSE;
1659:   useInitialGuess = PETSC_FALSE;
1660:   PetscObjectOptionsBegin((PetscObject)dm);
1661:   PetscOptionsName("-dm_plex_rebalance_shared_points_parmetis", "Use ParMetis instead of Metis for the partitioner", "DMPlexRebalanceSharedPoints", &parallel);
1662:   PetscOptionsBool("-dm_plex_rebalance_shared_points_use_initial_guess", "Use current partition to bootstrap ParMetis partition", "DMPlexRebalanceSharedPoints", useInitialGuess, &useInitialGuess, NULL);
1663:   PetscOptionsBool("-dm_plex_rebalance_shared_points_use_mat_partitioning", "Use the MatPartitioning object to partition", "DMPlexRebalanceSharedPoints", usematpartitioning, &usematpartitioning, NULL);
1664:   PetscOptionsViewer("-dm_plex_rebalance_shared_points_monitor", "Monitor the shared points rebalance process", "DMPlexRebalanceSharedPoints", &viewer, &format, NULL);
1665:   PetscOptionsEnd();
1666:   if (viewer) PetscViewerPushFormat(viewer, format);

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

1670:   DMGetOptionsPrefix(dm, &prefix);
1671:   PetscOptionsGetViewer(comm, ((PetscObject)dm)->options, prefix, "-dm_rebalance_partition_view", &viewer, &format, NULL);
1672:   if (viewer) PetscViewerPushFormat(viewer, format);

1674:   /* Figure out all points in the plex that we are interested in balancing. */
1675:   DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
1676:   DMPlexGetChart(dm, &pStart, &pEnd);
1677:   PetscMalloc1(pEnd - pStart, &toBalance);
1678:   for (i = 0; i < pEnd - pStart; i++) toBalance[i] = (PetscBool)(i >= eBegin && i < eEnd);

1680:   /* There are three types of points:
1681:    * exclusivelyOwned: points that are owned by this process and only seen by this process
1682:    * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1683:    * leaf: a point that is seen by this process but owned by a different process
1684:    */
1685:   DMGetPointSF(dm, &sf);
1686:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1687:   PetscMalloc1(pEnd - pStart, &isLeaf);
1688:   PetscMalloc1(pEnd - pStart, &isNonExclusivelyOwned);
1689:   PetscMalloc1(pEnd - pStart, &isExclusivelyOwned);
1690:   for (i = 0; i < pEnd - pStart; i++) {
1691:     isNonExclusivelyOwned[i] = PETSC_FALSE;
1692:     isExclusivelyOwned[i]    = PETSC_FALSE;
1693:     isLeaf[i]                = PETSC_FALSE;
1694:   }

1696:   /* mark all the leafs */
1697:   for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;

1699:   /* for an owned point, we can figure out whether another processor sees it or
1700:    * not by calculating its degree */
1701:   PetscSFComputeDegreeBegin(sf, &degrees);
1702:   PetscSFComputeDegreeEnd(sf, &degrees);
1703:   numExclusivelyOwned    = 0;
1704:   numNonExclusivelyOwned = 0;
1705:   for (i = 0; i < pEnd - pStart; i++) {
1706:     if (toBalance[i]) {
1707:       if (degrees[i] > 0) {
1708:         isNonExclusivelyOwned[i] = PETSC_TRUE;
1709:         numNonExclusivelyOwned += 1;
1710:       } else {
1711:         if (!isLeaf[i]) {
1712:           isExclusivelyOwned[i] = PETSC_TRUE;
1713:           numExclusivelyOwned += 1;
1714:         }
1715:       }
1716:     }
1717:   }

1719:   /* Build a graph with one vertex per core representing the
1720:    * exclusively owned points and then one vertex per nonExclusively owned
1721:    * point. */
1722:   PetscLayoutCreate(comm, &layout);
1723:   PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
1724:   PetscLayoutSetUp(layout);
1725:   PetscLayoutGetRanges(layout, &cumSumVertices);
1726:   PetscMalloc1(pEnd - pStart, &globalNumbersOfLocalOwnedVertices);
1727:   for (i = 0; i < pEnd - pStart; i++) globalNumbersOfLocalOwnedVertices[i] = pStart - 1;
1728:   offset  = cumSumVertices[rank];
1729:   counter = 0;
1730:   for (i = 0; i < pEnd - pStart; i++) {
1731:     if (toBalance[i]) {
1732:       if (degrees[i] > 0) {
1733:         globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1734:         counter++;
1735:       }
1736:     }
1737:   }

1739:   /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
1740:   PetscMalloc1(pEnd - pStart, &leafGlobalNumbers);
1741:   PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE);
1742:   PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE);

1744:   /* Build the graph for partitioning */
1745:   numRows = 1 + numNonExclusivelyOwned;
1746:   PetscLogEventBegin(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);
1747:   MatCreate(comm, &A);
1748:   MatSetType(A, MATMPIADJ);
1749:   MatSetSizes(A, numRows, numRows, cumSumVertices[size], cumSumVertices[size]);
1750:   idx = cumSumVertices[rank];
1751:   for (i = 0; i < pEnd - pStart; i++) {
1752:     if (toBalance[i]) {
1753:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1754:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1755:       else continue;
1756:       MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
1757:       MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
1758:     }
1759:   }
1760:   PetscFree(globalNumbersOfLocalOwnedVertices);
1761:   PetscFree(leafGlobalNumbers);
1762:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1763:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1764:   PetscLogEventEnd(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);

1766:   nparts = size;
1767:   ncon   = 1;
1768:   PetscMalloc1(ncon * nparts, &tpwgts);
1769:   for (i = 0; i < ncon * nparts; i++) tpwgts[i] = 1. / (nparts);
1770:   PetscMalloc1(ncon, &ubvec);
1771:   ubvec[0] = 1.05;
1772:   ubvec[1] = 1.05;

1774:   PetscMalloc1(ncon * (1 + numNonExclusivelyOwned), &vtxwgt);
1775:   if (ncon == 2) {
1776:     vtxwgt[0] = numExclusivelyOwned;
1777:     vtxwgt[1] = 1;
1778:     for (i = 0; i < numNonExclusivelyOwned; i++) {
1779:       vtxwgt[ncon * (i + 1)]     = 1;
1780:       vtxwgt[ncon * (i + 1) + 1] = 0;
1781:     }
1782:   } else {
1783:     PetscInt base, ms;
1784:     MPI_Allreduce(&numExclusivelyOwned, &base, 1, MPIU_INT, MPIU_MAX, PetscObjectComm((PetscObject)dm));
1785:     MatGetSize(A, &ms, NULL);
1786:     ms -= size;
1787:     base      = PetscMax(base, ms);
1788:     vtxwgt[0] = base + numExclusivelyOwned;
1789:     for (i = 0; i < numNonExclusivelyOwned; i++) vtxwgt[i + 1] = 1;
1790:   }

1792:   if (viewer) {
1793:     PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %" PetscInt_FMT " on interface of mesh distribution.\n", entityDepth);
1794:     PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %" PetscInt_FMT "\n", cumSumVertices[size]);
1795:   }
1796:   /* TODO: Drop the parallel/sequential choice here and just use MatPartioner for much more flexibility */
1797:   if (usematpartitioning) {
1798:     const char *prefix;

1800:     MatPartitioningCreate(PetscObjectComm((PetscObject)dm), &mp);
1801:     PetscObjectSetOptionsPrefix((PetscObject)mp, "dm_plex_rebalance_shared_points_");
1802:     PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
1803:     PetscObjectPrependOptionsPrefix((PetscObject)mp, prefix);
1804:     MatPartitioningSetAdjacency(mp, A);
1805:     MatPartitioningSetNumberVertexWeights(mp, ncon);
1806:     MatPartitioningSetVertexWeights(mp, vtxwgt);
1807:     MatPartitioningSetFromOptions(mp);
1808:     MatPartitioningApply(mp, &ispart);
1809:     ISGetIndices(ispart, (const PetscInt **)&part);
1810:   } else if (parallel) {
1811:     if (viewer) PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n");
1812:     PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1813:     MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1815:     PetscMalloc1(4, &options);
1816:     options[0] = 1;
1817:     options[1] = 0; /* Verbosity */
1818:     if (viewer) options[1] = 1;
1819:     options[2] = 0;                    /* Seed */
1820:     options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1821:     wgtflag    = 2;
1822:     numflag    = 0;
1823:     if (useInitialGuess) {
1824:       PetscViewerASCIIPrintf(viewer, "THIS DOES NOT WORK! I don't know why. Using current distribution of points as initial guess.\n");
1825:       for (i = 0; i < numRows; i++) part[i] = rank;
1826:       if (viewer) PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n");
1827:       PetscStackPushExternal("ParMETIS_V3_RefineKway");
1828:       PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1829:       ParMETIS_V3_RefineKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1830:       PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1831:       PetscStackPop;
1833:     } else {
1834:       PetscStackPushExternal("ParMETIS_V3_PartKway");
1835:       PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1836:       ParMETIS_V3_PartKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1837:       PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1838:       PetscStackPop;
1840:     }
1841:     MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1842:     PetscFree(options);
1843:   } else {
1844:     if (viewer) PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n");
1845:     Mat       As;
1846:     PetscInt *partGlobal;
1847:     PetscInt *numExclusivelyOwnedAll;

1849:     PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1850:     MatGetSize(A, &numRows, NULL);
1851:     PetscLogEventBegin(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);
1852:     MatMPIAdjToSeqRankZero(A, &As);
1853:     PetscLogEventEnd(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);

1855:     PetscMalloc1(size, &numExclusivelyOwnedAll);
1856:     numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1857:     MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, numExclusivelyOwnedAll, 1, MPIU_INT, comm);

1859:     PetscMalloc1(numRows, &partGlobal);
1860:     PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1861:     if (rank == 0) {
1862:       PetscInt *vtxwgt_g, numRows_g;

1864:       MatGetRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1865:       PetscMalloc1(2 * numRows_g, &vtxwgt_g);
1866:       for (i = 0; i < size; i++) {
1867:         vtxwgt_g[ncon * cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1868:         if (ncon > 1) vtxwgt_g[ncon * cumSumVertices[i] + 1] = 1;
1869:         for (j = cumSumVertices[i] + 1; j < cumSumVertices[i + 1]; j++) {
1870:           vtxwgt_g[ncon * j] = 1;
1871:           if (ncon > 1) vtxwgt_g[2 * j + 1] = 0;
1872:         }
1873:       }

1875:       PetscMalloc1(64, &options);
1876:       METIS_SetDefaultOptions(options); /* initialize all defaults */
1878:       options[METIS_OPTION_CONTIG] = 1;
1879:       PetscStackPushExternal("METIS_PartGraphKway");
1880:       METIS_PartGraphKway(&numRows_g, &ncon, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1881:       PetscStackPop;
1883:       PetscFree(options);
1884:       PetscFree(vtxwgt_g);
1885:       MatRestoreRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1886:       MatDestroy(&As);
1887:     }
1888:     PetscBarrier((PetscObject)dm);
1889:     PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1890:     PetscFree(numExclusivelyOwnedAll);

1892:     /* scatter the partitioning information to ranks */
1893:     PetscLogEventBegin(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1894:     PetscMalloc1(size, &counts);
1895:     PetscMalloc1(size + 1, &mpiCumSumVertices);
1896:     for (i = 0; i < size; i++) PetscMPIIntCast(cumSumVertices[i + 1] - cumSumVertices[i], &(counts[i]));
1897:     for (i = 0; i <= size; i++) PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
1898:     MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
1899:     PetscFree(counts);
1900:     PetscFree(mpiCumSumVertices);
1901:     PetscFree(partGlobal);
1902:     PetscLogEventEnd(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1903:   }
1904:   PetscFree(ubvec);
1905:   PetscFree(tpwgts);

1907:   /* Rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1908:   PetscMalloc2(size, &firstVertices, size, &renumbering);
1909:   firstVertices[rank] = part[0];
1910:   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, firstVertices, 1, MPIU_INT, comm);
1911:   for (i = 0; i < size; i++) renumbering[firstVertices[i]] = i;
1912:   for (i = 0; i < cumSumVertices[rank + 1] - cumSumVertices[rank]; i++) part[i] = renumbering[part[i]];
1913:   PetscFree2(firstVertices, renumbering);

1915:   /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1916:   failed = (PetscInt)(part[0] != rank);
1917:   MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
1918:   if (failedGlobal > 0) {
1920:     PetscFree(vtxwgt);
1921:     PetscFree(toBalance);
1922:     PetscFree(isLeaf);
1923:     PetscFree(isNonExclusivelyOwned);
1924:     PetscFree(isExclusivelyOwned);
1925:     if (usematpartitioning) {
1926:       ISRestoreIndices(ispart, (const PetscInt **)&part);
1927:       ISDestroy(&ispart);
1928:     } else PetscFree(part);
1929:     if (viewer) {
1930:       PetscViewerPopFormat(viewer);
1931:       PetscViewerDestroy(&viewer);
1932:     }
1933:     PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1934:     return 0;
1935:   }

1937:   /* Check how well we did distributing points*/
1938:   if (viewer) {
1939:     PetscViewerASCIIPrintf(viewer, "Number of owned entities of depth %" PetscInt_FMT ".\n", entityDepth);
1940:     PetscViewerASCIIPrintf(viewer, "Initial      ");
1941:     DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
1942:     PetscViewerASCIIPrintf(viewer, "Rebalanced   ");
1943:     DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1944:   }

1946:   /* Check that every vertex is owned by a process that it is actually connected to. */
1947:   MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1948:   for (i = 1; i <= numNonExclusivelyOwned; i++) {
1949:     PetscInt loc = 0;
1950:     PetscFindInt(cumSumVertices[part[i]], xadj[i + 1] - xadj[i], &adjncy[xadj[i]], &loc);
1951:     /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1952:     if (loc < 0) part[i] = rank;
1953:   }
1954:   MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1955:   MatDestroy(&A);

1957:   /* See how significant the influences of the previous fixing up step was.*/
1958:   if (viewer) {
1959:     PetscViewerASCIIPrintf(viewer, "After fix    ");
1960:     DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1961:   }
1962:   if (!usematpartitioning) PetscFree(vtxwgt);
1963:   else MatPartitioningDestroy(&mp);

1965:   PetscLayoutDestroy(&layout);

1967:   PetscLogEventBegin(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);
1968:   /* Rewrite the SF to reflect the new ownership. */
1969:   PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
1970:   counter = 0;
1971:   for (i = 0; i < pEnd - pStart; i++) {
1972:     if (toBalance[i]) {
1973:       if (isNonExclusivelyOwned[i]) {
1974:         pointsToRewrite[counter] = i + pStart;
1975:         counter++;
1976:       }
1977:     }
1978:   }
1979:   DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part + 1, degrees);
1980:   PetscFree(pointsToRewrite);
1981:   PetscLogEventEnd(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);

1983:   PetscFree(toBalance);
1984:   PetscFree(isLeaf);
1985:   PetscFree(isNonExclusivelyOwned);
1986:   PetscFree(isExclusivelyOwned);
1987:   if (usematpartitioning) {
1988:     ISRestoreIndices(ispart, (const PetscInt **)&part);
1989:     ISDestroy(&ispart);
1990:   } else PetscFree(part);
1991:   if (viewer) {
1992:     PetscViewerPopFormat(viewer);
1993:     PetscViewerDestroy(&viewer);
1994:   }
1995:   if (success) *success = PETSC_TRUE;
1996:   PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1997:   return 0;
1998: #else
1999:   SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2000: #endif
2001: }