Actual source code: plexpartition.c

petsc-master 2019-11-16
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: PETSC_STATIC_INLINE PetscInt DMPlex_GlobalID(PetscInt point) { return point >= 0 ? point : -(point+1); }

 32: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
 33: {
 34:   PetscInt       dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
 35:   PetscInt      *adj = NULL, *vOffsets = NULL, *graph = NULL;
 36:   IS             cellNumbering;
 37:   const PetscInt *cellNum;
 38:   PetscBool      useCone, useClosure;
 39:   PetscSection   section;
 40:   PetscSegBuffer adjBuffer;
 41:   PetscSF        sfPoint;
 42:   PetscInt       *adjCells = NULL, *remoteCells = NULL;
 43:   const PetscInt *local;
 44:   PetscInt       nroots, nleaves, l;
 45:   PetscMPIInt    rank;

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

 88:     PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
 89:     DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
 90:     for (l = 0; l < nroots; ++l) adjCells[l] = -3;
 91:     for (f = fStart; f < fEnd; ++f) {
 92:       const PetscInt *support;
 93:       PetscInt        supportSize;

 95:       DMPlexGetSupport(dm, f, &support);
 96:       DMPlexGetSupportSize(dm, f, &supportSize);
 97:       if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
 98:       else if (supportSize == 2) {
 99:         PetscFindInt(support[0], nleaves, local, &p);
100:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1]]);
101:         PetscFindInt(support[1], nleaves, local, &p);
102:         if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
103:       }
104:       /* Handle non-conforming meshes */
105:       if (supportSize > 2) {
106:         PetscInt        numChildren, i;
107:         const PetscInt *children;

109:         DMPlexGetTreeChildren(dm, f, &numChildren, &children);
110:         for (i = 0; i < numChildren; ++i) {
111:           const PetscInt child = children[i];
112:           if (fStart <= child && child < fEnd) {
113:             DMPlexGetSupport(dm, child, &support);
114:             DMPlexGetSupportSize(dm, child, &supportSize);
115:             if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
116:             else if (supportSize == 2) {
117:               PetscFindInt(support[0], nleaves, local, &p);
118:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1]]);
119:               PetscFindInt(support[1], nleaves, local, &p);
120:               if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
121:             }
122:           }
123:         }
124:       }
125:     }
126:     for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
127:     PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells);
128:     PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells);
129:     PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
130:     PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
131:   }
132:   /* Combine local and global adjacencies */
133:   for (*numVertices = 0, p = pStart; p < pEnd; p++) {
134:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
135:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
136:     /* Add remote cells */
137:     if (remoteCells) {
138:       const PetscInt gp = DMPlex_GlobalID(cellNum[p]);
139:       PetscInt       coneSize, numChildren, c, i;
140:       const PetscInt *cone, *children;

142:       DMPlexGetCone(dm, p, &cone);
143:       DMPlexGetConeSize(dm, p, &coneSize);
144:       for (c = 0; c < coneSize; ++c) {
145:         const PetscInt point = cone[c];
146:         if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
147:           PetscInt *PETSC_RESTRICT pBuf;
148:           PetscSectionAddDof(section, p, 1);
149:           PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
150:           *pBuf = remoteCells[point];
151:         }
152:         /* Handle non-conforming meshes */
153:         DMPlexGetTreeChildren(dm, point, &numChildren, &children);
154:         for (i = 0; i < numChildren; ++i) {
155:           const PetscInt child = children[i];
156:           if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
157:             PetscInt *PETSC_RESTRICT pBuf;
158:             PetscSectionAddDof(section, p, 1);
159:             PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
160:             *pBuf = remoteCells[child];
161:           }
162:         }
163:       }
164:     }
165:     /* Add local cells */
166:     adjSize = PETSC_DETERMINE;
167:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
168:     for (a = 0; a < adjSize; ++a) {
169:       const PetscInt point = adj[a];
170:       if (point != p && pStart <= point && point < pEnd) {
171:         PetscInt *PETSC_RESTRICT pBuf;
172:         PetscSectionAddDof(section, p, 1);
173:         PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
174:         *pBuf = DMPlex_GlobalID(cellNum[point]);
175:       }
176:     }
177:     (*numVertices)++;
178:   }
179:   PetscFree(adj);
180:   PetscFree2(adjCells, remoteCells);
181:   DMSetBasicAdjacency(dm, useCone, useClosure);

183:   /* Derive CSR graph from section/segbuffer */
184:   PetscSectionSetUp(section);
185:   PetscSectionGetStorageSize(section, &size);
186:   PetscMalloc1(*numVertices+1, &vOffsets);
187:   for (idx = 0, p = pStart; p < pEnd; p++) {
188:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
189:     PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
190:   }
191:   vOffsets[*numVertices] = size;
192:   PetscSegBufferExtractAlloc(adjBuffer, &graph);

194:   if (nroots >= 0) {
195:     /* Filter out duplicate edges using section/segbuffer */
196:     PetscSectionReset(section);
197:     PetscSectionSetChart(section, 0, *numVertices);
198:     for (p = 0; p < *numVertices; p++) {
199:       PetscInt start = vOffsets[p], end = vOffsets[p+1];
200:       PetscInt numEdges = end-start, *PETSC_RESTRICT edges;
201:       PetscSortRemoveDupsInt(&numEdges, &graph[start]);
202:       PetscSectionSetDof(section, p, numEdges);
203:       PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
204:       PetscArraycpy(edges, &graph[start], numEdges);
205:     }
206:     PetscFree(vOffsets);
207:     PetscFree(graph);
208:     /* Derive CSR graph from section/segbuffer */
209:     PetscSectionSetUp(section);
210:     PetscSectionGetStorageSize(section, &size);
211:     PetscMalloc1(*numVertices+1, &vOffsets);
212:     for (idx = 0, p = 0; p < *numVertices; p++) {
213:       PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
214:     }
215:     vOffsets[*numVertices] = size;
216:     PetscSegBufferExtractAlloc(adjBuffer, &graph);
217:   } else {
218:     /* Sort adjacencies (not strictly necessary) */
219:     for (p = 0; p < *numVertices; p++) {
220:       PetscInt start = vOffsets[p], end = vOffsets[p+1];
221:       PetscSortInt(end-start, &graph[start]);
222:     }
223:   }

225:   if (offsets) *offsets = vOffsets;
226:   if (adjacency) *adjacency = graph;

228:   /* Cleanup */
229:   PetscSegBufferDestroy(&adjBuffer);
230:   PetscSectionDestroy(&section);
231:   ISRestoreIndices(cellNumbering, &cellNum);
232:   ISDestroy(&cellNumbering);
233:   return(0);
234: }

236: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
237: {
238:   Mat            conn, CSR;
239:   IS             fis, cis, cis_own;
240:   PetscSF        sfPoint;
241:   const PetscInt *rows, *cols, *ii, *jj;
242:   PetscInt       *idxs,*idxs2;
243:   PetscInt       dim, depth, floc, cloc, i, M, N, c, m, cStart, cEnd, fStart, fEnd;
244:   PetscMPIInt    rank;
245:   PetscBool      flg;

249:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
250:   DMGetDimension(dm, &dim);
251:   DMPlexGetDepth(dm, &depth);
252:   if (dim != depth) {
253:     /* We do not handle the uninterpolated case here */
254:     DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
255:     /* DMPlexCreateNeighborCSR does not make a numbering */
256:     if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
257:     /* Different behavior for empty graphs */
258:     if (!*numVertices) {
259:       PetscMalloc1(1, offsets);
260:       (*offsets)[0] = 0;
261:     }
262:     /* Broken in parallel */
263:     if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
264:     return(0);
265:   }
266:   /* Interpolated and parallel case */

268:   /* numbering */
269:   DMGetPointSF(dm, &sfPoint);
270:   DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
271:   DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
272:   DMPlexCreateNumbering_Internal(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
273:   DMPlexCreateNumbering_Internal(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
274:   if (globalNumbering) {
275:     ISDuplicate(cis, globalNumbering);
276:   }

278:   /* get positive global ids and local sizes for facets and cells */
279:   ISGetLocalSize(fis, &m);
280:   ISGetIndices(fis, &rows);
281:   PetscMalloc1(m, &idxs);
282:   for (i = 0, floc = 0; i < m; i++) {
283:     const PetscInt p = rows[i];

285:     if (p < 0) {
286:       idxs[i] = -(p+1);
287:     } else {
288:       idxs[i] = p;
289:       floc   += 1;
290:     }
291:   }
292:   ISRestoreIndices(fis, &rows);
293:   ISDestroy(&fis);
294:   ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);

296:   ISGetLocalSize(cis, &m);
297:   ISGetIndices(cis, &cols);
298:   PetscMalloc1(m, &idxs);
299:   PetscMalloc1(m, &idxs2);
300:   for (i = 0, cloc = 0; i < m; i++) {
301:     const PetscInt p = cols[i];

303:     if (p < 0) {
304:       idxs[i] = -(p+1);
305:     } else {
306:       idxs[i]       = p;
307:       idxs2[cloc++] = p;
308:     }
309:   }
310:   ISRestoreIndices(cis, &cols);
311:   ISDestroy(&cis);
312:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
313:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);

315:   /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
316:   MatCreate(PetscObjectComm((PetscObject)dm), &conn);
317:   MatSetSizes(conn, floc, cloc, M, N);
318:   MatSetType(conn, MATMPIAIJ);
319:   DMPlexGetMaxSizes(dm, NULL, &m);
320:   MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);

322:   /* Assemble matrix */
323:   ISGetIndices(fis, &rows);
324:   ISGetIndices(cis, &cols);
325:   for (c = cStart; c < cEnd; c++) {
326:     const PetscInt *cone;
327:     PetscInt        coneSize, row, col, f;

329:     col  = cols[c-cStart];
330:     DMPlexGetCone(dm, c, &cone);
331:     DMPlexGetConeSize(dm, c, &coneSize);
332:     for (f = 0; f < coneSize; f++) {
333:       const PetscScalar v = 1.0;
334:       const PetscInt *children;
335:       PetscInt        numChildren, ch;

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

340:       /* non-conforming meshes */
341:       DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
342:       for (ch = 0; ch < numChildren; ch++) {
343:         const PetscInt child = children[ch];

345:         if (child < fStart || child >= fEnd) continue;
346:         row  = rows[child-fStart];
347:         MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
348:       }
349:     }
350:   }
351:   ISRestoreIndices(fis, &rows);
352:   ISRestoreIndices(cis, &cols);
353:   ISDestroy(&fis);
354:   ISDestroy(&cis);
355:   MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
356:   MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);

358:   /* Get parallel CSR by doing conn^T * conn */
359:   MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
360:   MatDestroy(&conn);

362:   /* extract local part of the CSR */
363:   MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
364:   MatDestroy(&CSR);
365:   MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
366:   if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");

368:   /* get back requested output */
369:   if (numVertices) *numVertices = m;
370:   if (offsets) {
371:     PetscCalloc1(m+1, &idxs);
372:     for (i = 1; i < m+1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
373:     *offsets = idxs;
374:   }
375:   if (adjacency) {
376:     PetscMalloc1(ii[m] - m, &idxs);
377:     ISGetIndices(cis_own, &rows);
378:     for (i = 0, c = 0; i < m; i++) {
379:       PetscInt j, g = rows[i];

381:       for (j = ii[i]; j < ii[i+1]; j++) {
382:         if (jj[j] == g) continue; /* again, self-connectivity */
383:         idxs[c++] = jj[j];
384:       }
385:     }
386:     if (c != ii[m] - m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %D != %D",c,ii[m]-m);
387:     ISRestoreIndices(cis_own, &rows);
388:     *adjacency = idxs;
389:   }

391:   /* cleanup */
392:   ISDestroy(&cis_own);
393:   MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
394:   if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
395:   MatDestroy(&conn);
396:   return(0);
397: }

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

402:   Input Parameters:
403: + dm      - The mesh DM dm
404: - height  - Height of the strata from which to construct the graph

406:   Output Parameter:
407: + numVertices     - Number of vertices in the graph
408: . offsets         - Point offsets in the graph
409: . adjacency       - Point connectivity in the graph
410: - 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.

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

415:   Level: developer

417: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
418: @*/
419: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
420: {
422:   PetscBool      usemat = PETSC_FALSE;

425:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_csr_via_mat", &usemat, NULL);
426:   if (usemat) {
427:     DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
428:   } else {
429:     DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
430:   }
431:   return(0);
432: }

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

437:   Collective on DM

439:   Input Arguments:
440: + dm - The DMPlex
441: - cellHeight - The height of mesh points to treat as cells (default should be 0)

443:   Output Arguments:
444: + numVertices - The number of local vertices in the graph, or cells in the mesh.
445: . offsets     - The offset to the adjacency list for each cell
446: - adjacency   - The adjacency list for all cells

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

450:   Level: advanced

452: .seealso: DMPlexCreate()
453: @*/
454: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
455: {
456:   const PetscInt maxFaceCases = 30;
457:   PetscInt       numFaceCases = 0;
458:   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
459:   PetscInt      *off, *adj;
460:   PetscInt      *neighborCells = NULL;
461:   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;

465:   /* For parallel partitioning, I think you have to communicate supports */
466:   DMGetDimension(dm, &dim);
467:   cellDim = dim - cellHeight;
468:   DMPlexGetDepth(dm, &depth);
469:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
470:   if (cEnd - cStart == 0) {
471:     if (numVertices) *numVertices = 0;
472:     if (offsets)   *offsets   = NULL;
473:     if (adjacency) *adjacency = NULL;
474:     return(0);
475:   }
476:   numCells  = cEnd - cStart;
477:   faceDepth = depth - cellHeight;
478:   if (dim == depth) {
479:     PetscInt f, fStart, fEnd;

481:     PetscCalloc1(numCells+1, &off);
482:     /* Count neighboring cells */
483:     DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
484:     for (f = fStart; f < fEnd; ++f) {
485:       const PetscInt *support;
486:       PetscInt        supportSize;
487:       DMPlexGetSupportSize(dm, f, &supportSize);
488:       DMPlexGetSupport(dm, f, &support);
489:       if (supportSize == 2) {
490:         PetscInt numChildren;

492:         DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
493:         if (!numChildren) {
494:           ++off[support[0]-cStart+1];
495:           ++off[support[1]-cStart+1];
496:         }
497:       }
498:     }
499:     /* Prefix sum */
500:     for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
501:     if (adjacency) {
502:       PetscInt *tmp;

504:       PetscMalloc1(off[numCells], &adj);
505:       PetscMalloc1(numCells+1, &tmp);
506:       PetscArraycpy(tmp, off, numCells+1);
507:       /* Get neighboring cells */
508:       for (f = fStart; f < fEnd; ++f) {
509:         const PetscInt *support;
510:         PetscInt        supportSize;
511:         DMPlexGetSupportSize(dm, f, &supportSize);
512:         DMPlexGetSupport(dm, f, &support);
513:         if (supportSize == 2) {
514:           PetscInt numChildren;

516:           DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
517:           if (!numChildren) {
518:             adj[tmp[support[0]-cStart]++] = support[1];
519:             adj[tmp[support[1]-cStart]++] = support[0];
520:           }
521:         }
522:       }
523: #if defined(PETSC_USE_DEBUG)
524:       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);
525: #endif
526:       PetscFree(tmp);
527:     }
528:     if (numVertices) *numVertices = numCells;
529:     if (offsets)   *offsets   = off;
530:     if (adjacency) *adjacency = adj;
531:     return(0);
532:   }
533:   /* Setup face recognition */
534:   if (faceDepth == 1) {
535:     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 */

537:     for (c = cStart; c < cEnd; ++c) {
538:       PetscInt corners;

540:       DMPlexGetConeSize(dm, c, &corners);
541:       if (!cornersSeen[corners]) {
542:         PetscInt nFV;

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

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

549:         numFaceVertices[numFaceCases++] = nFV;
550:       }
551:     }
552:   }
553:   PetscCalloc1(numCells+1, &off);
554:   /* Count neighboring cells */
555:   for (cell = cStart; cell < cEnd; ++cell) {
556:     PetscInt numNeighbors = PETSC_DETERMINE, n;

558:     DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
559:     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
560:     for (n = 0; n < numNeighbors; ++n) {
561:       PetscInt        cellPair[2];
562:       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
563:       PetscInt        meetSize = 0;
564:       const PetscInt *meet    = NULL;

566:       cellPair[0] = cell; cellPair[1] = neighborCells[n];
567:       if (cellPair[0] == cellPair[1]) continue;
568:       if (!found) {
569:         DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
570:         if (meetSize) {
571:           PetscInt f;

573:           for (f = 0; f < numFaceCases; ++f) {
574:             if (numFaceVertices[f] == meetSize) {
575:               found = PETSC_TRUE;
576:               break;
577:             }
578:           }
579:         }
580:         DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
581:       }
582:       if (found) ++off[cell-cStart+1];
583:     }
584:   }
585:   /* Prefix sum */
586:   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];

588:   if (adjacency) {
589:     PetscMalloc1(off[numCells], &adj);
590:     /* Get neighboring cells */
591:     for (cell = cStart; cell < cEnd; ++cell) {
592:       PetscInt numNeighbors = PETSC_DETERMINE, n;
593:       PetscInt cellOffset   = 0;

595:       DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
596:       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
597:       for (n = 0; n < numNeighbors; ++n) {
598:         PetscInt        cellPair[2];
599:         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
600:         PetscInt        meetSize = 0;
601:         const PetscInt *meet    = NULL;

603:         cellPair[0] = cell; cellPair[1] = neighborCells[n];
604:         if (cellPair[0] == cellPair[1]) continue;
605:         if (!found) {
606:           DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
607:           if (meetSize) {
608:             PetscInt f;

610:             for (f = 0; f < numFaceCases; ++f) {
611:               if (numFaceVertices[f] == meetSize) {
612:                 found = PETSC_TRUE;
613:                 break;
614:               }
615:             }
616:           }
617:           DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
618:         }
619:         if (found) {
620:           adj[off[cell-cStart]+cellOffset] = neighborCells[n];
621:           ++cellOffset;
622:         }
623:       }
624:     }
625:   }
626:   PetscFree(neighborCells);
627:   if (numVertices) *numVertices = numCells;
628:   if (offsets)   *offsets   = off;
629:   if (adjacency) *adjacency = adj;
630:   return(0);
631: }

633: /*@C
634:   PetscPartitionerRegister - Adds a new PetscPartitioner implementation

636:   Not Collective

638:   Input Parameters:
639: + name        - The name of a new user-defined creation routine
640: - create_func - The creation routine itself

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

645:   Sample usage:
646: .vb
647:     PetscPartitionerRegister("my_part", MyPetscPartitionerCreate);
648: .ve

650:   Then, your PetscPartitioner type can be chosen with the procedural interface via
651: .vb
652:     PetscPartitionerCreate(MPI_Comm, PetscPartitioner *);
653:     PetscPartitionerSetType(PetscPartitioner, "my_part");
654: .ve
655:    or at runtime via the option
656: .vb
657:     -petscpartitioner_type my_part
658: .ve

660:   Level: advanced

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

664: @*/
665: PetscErrorCode PetscPartitionerRegister(const char sname[], PetscErrorCode (*function)(PetscPartitioner))
666: {

670:   PetscFunctionListAdd(&PetscPartitionerList, sname, function);
671:   return(0);
672: }

674: /*@C
675:   PetscPartitionerSetType - Builds a particular PetscPartitioner

677:   Collective on PetscPartitioner

679:   Input Parameters:
680: + part - The PetscPartitioner object
681: - name - The kind of partitioner

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

686:   Level: intermediate

688: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate()
689: @*/
690: PetscErrorCode PetscPartitionerSetType(PetscPartitioner part, PetscPartitionerType name)
691: {
692:   PetscErrorCode (*r)(PetscPartitioner);
693:   PetscBool      match;

698:   PetscObjectTypeCompare((PetscObject) part, name, &match);
699:   if (match) return(0);

701:   PetscPartitionerRegisterAll();
702:   PetscFunctionListFind(PetscPartitionerList, name, &r);
703:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscPartitioner type: %s", name);

705:   if (part->ops->destroy) {
706:     (*part->ops->destroy)(part);
707:   }
708:   part->noGraph = PETSC_FALSE;
709:   PetscMemzero(part->ops, sizeof(*part->ops));
710:   PetscObjectChangeTypeName((PetscObject) part, name);
711:   (*r)(part);
712:   return(0);
713: }

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

718:   Not Collective

720:   Input Parameter:
721: . part - The PetscPartitioner

723:   Output Parameter:
724: . name - The PetscPartitioner type name

726:   Level: intermediate

728: .seealso: PetscPartitionerSetType(), PetscPartitionerCreate()
729: @*/
730: PetscErrorCode PetscPartitionerGetType(PetscPartitioner part, PetscPartitionerType *name)
731: {

737:   PetscPartitionerRegisterAll();
738:   *name = ((PetscObject) part)->type_name;
739:   return(0);
740: }

742: /*@C
743:    PetscPartitionerViewFromOptions - View from Options

745:    Collective on PetscPartitioner

747:    Input Parameters:
748: +  A - the PetscPartitioner object
749: .  obj - Optional object
750: -  name - command line option

752:    Level: intermediate
753: .seealso:  PetscPartitionerView(), PetscObjectViewFromOptions()
754: @*/
755: PetscErrorCode PetscPartitionerViewFromOptions(PetscPartitioner A,PetscObject obj,const char name[])
756: {

761:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
762:   return(0);
763: }

765: /*@C
766:   PetscPartitionerView - Views a PetscPartitioner

768:   Collective on PetscPartitioner

770:   Input Parameter:
771: + part - the PetscPartitioner object to view
772: - v    - the viewer

774:   Level: developer

776: .seealso: PetscPartitionerDestroy()
777: @*/
778: PetscErrorCode PetscPartitionerView(PetscPartitioner part, PetscViewer v)
779: {
780:   PetscMPIInt    size;
781:   PetscBool      isascii;

786:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) part), &v);}
787:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
788:   if (isascii) {
789:     MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
790:     PetscViewerASCIIPrintf(v, "Graph Partitioner: %d MPI Process%s\n", size, size > 1 ? "es" : "");
791:     PetscViewerASCIIPrintf(v, "  type: %s\n", part->hdr.type_name);
792:     PetscViewerASCIIPushTab(v);
793:     PetscViewerASCIIPrintf(v, "edge cut: %D\n", part->edgeCut);
794:     PetscViewerASCIIPrintf(v, "balance:  %.2g\n", part->balance);
795:     PetscViewerASCIIPopTab(v);
796:   }
797:   if (part->ops->view) {(*part->ops->view)(part, v);}
798:   return(0);
799: }

801: static PetscErrorCode PetscPartitionerGetDefaultType(const char *currentType, const char **defaultType)
802: {
804:   if (!currentType) {
805: #if defined(PETSC_HAVE_PARMETIS)
806:     *defaultType = PETSCPARTITIONERPARMETIS;
807: #elif defined(PETSC_HAVE_PTSCOTCH)
808:     *defaultType = PETSCPARTITIONERPTSCOTCH;
809: #elif defined(PETSC_HAVE_CHACO)
810:     *defaultType = PETSCPARTITIONERCHACO;
811: #else
812:     *defaultType = PETSCPARTITIONERSIMPLE;
813: #endif
814:   } else {
815:     *defaultType = currentType;
816:   }
817:   return(0);
818: }

820: /*@
821:   PetscPartitionerSetFromOptions - sets parameters in a PetscPartitioner from the options database

823:   Collective on PetscPartitioner

825:   Input Parameter:
826: . part - the PetscPartitioner object to set options for

828:   Level: developer

830: .seealso: PetscPartitionerView()
831: @*/
832: PetscErrorCode PetscPartitionerSetFromOptions(PetscPartitioner part)
833: {
834:   const char    *defaultType;
835:   char           name[256];
836:   PetscBool      flg;

841:   PetscPartitionerRegisterAll();
842:   PetscPartitionerGetDefaultType(((PetscObject) part)->type_name,&defaultType);
843:   PetscObjectOptionsBegin((PetscObject) part);
844:   PetscOptionsFList("-petscpartitioner_type", "Graph partitioner", "PetscPartitionerSetType", PetscPartitionerList, defaultType, name, sizeof(name), &flg);
845:   if (flg) {
846:     PetscPartitionerSetType(part, name);
847:   } else if (!((PetscObject) part)->type_name) {
848:     PetscPartitionerSetType(part, defaultType);
849:   }
850:   if (part->ops->setfromoptions) {
851:     (*part->ops->setfromoptions)(PetscOptionsObject,part);
852:   }
853:   PetscViewerDestroy(&part->viewerGraph);
854:   PetscOptionsGetViewer(((PetscObject) part)->comm, ((PetscObject) part)->options, ((PetscObject) part)->prefix, "-petscpartitioner_view_graph", &part->viewerGraph, &part->formatGraph, &part->viewGraph);
855:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
856:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) part);
857:   PetscOptionsEnd();
858:   return(0);
859: }

861: /*@C
862:   PetscPartitionerSetUp - Construct data structures for the PetscPartitioner

864:   Collective on PetscPartitioner

866:   Input Parameter:
867: . part - the PetscPartitioner object to setup

869:   Level: developer

871: .seealso: PetscPartitionerView(), PetscPartitionerDestroy()
872: @*/
873: PetscErrorCode PetscPartitionerSetUp(PetscPartitioner part)
874: {

879:   if (part->ops->setup) {(*part->ops->setup)(part);}
880:   return(0);
881: }

883: /*@
884:   PetscPartitionerDestroy - Destroys a PetscPartitioner object

886:   Collective on PetscPartitioner

888:   Input Parameter:
889: . part - the PetscPartitioner object to destroy

891:   Level: developer

893: .seealso: PetscPartitionerView()
894: @*/
895: PetscErrorCode PetscPartitionerDestroy(PetscPartitioner *part)
896: {

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

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

906:   PetscViewerDestroy(&(*part)->viewerGraph);
907:   if ((*part)->ops->destroy) {(*(*part)->ops->destroy)(*part);}
908:   PetscHeaderDestroy(part);
909:   return(0);
910: }

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

915:   Collective

917:   Input Parameter:
918: . comm - The communicator for the PetscPartitioner object

920:   Output Parameter:
921: . part - The PetscPartitioner object

923:   Level: beginner

925: .seealso: PetscPartitionerSetType(), PETSCPARTITIONERCHACO, PETSCPARTITIONERPARMETIS, PETSCPARTITIONERSHELL, PETSCPARTITIONERSIMPLE, PETSCPARTITIONERGATHER
926: @*/
927: PetscErrorCode PetscPartitionerCreate(MPI_Comm comm, PetscPartitioner *part)
928: {
929:   PetscPartitioner p;
930:   const char       *partitionerType = NULL;
931:   PetscErrorCode   ierr;

935:   *part = NULL;
936:   DMInitializePackage();

938:   PetscHeaderCreate(p, PETSCPARTITIONER_CLASSID, "PetscPartitioner", "Graph Partitioner", "PetscPartitioner", comm, PetscPartitionerDestroy, PetscPartitionerView);
939:   PetscPartitionerGetDefaultType(NULL,&partitionerType);
940:   PetscPartitionerSetType(p,partitionerType);

942:   p->edgeCut = 0;
943:   p->balance = 0.0;

945:   *part = p;
946:   return(0);
947: }

949: /*@
950:   PetscPartitionerPartition - Create a non-overlapping partition of the cells in the mesh

952:   Collective on PetscPartitioner

954:   Input Parameters:
955: + part    - The PetscPartitioner
956: - dm      - The mesh DM

958:   Output Parameters:
959: + partSection     - The PetscSection giving the division of points by partition
960: - partition       - The list of points by partition

962:   Options Database:
963: . -petscpartitioner_view_graph - View the graph we are partitioning

965:   Note: Instead of cells, points at a given height can be partitioned by calling PetscPartitionerSetPointHeight()

967:   Level: developer

969: .seealso DMPlexDistribute(), PetscPartitionerSetPointHeight(), PetscPartitionerCreate()
970: @*/
971: PetscErrorCode PetscPartitionerPartition(PetscPartitioner part, DM dm, PetscSection partSection, IS *partition)
972: {
973:   PetscMPIInt    size;

981:   MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
982:   if (size == 1) {
983:     PetscInt *points;
984:     PetscInt  cStart, cEnd, c;

986:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
987:     PetscSectionSetChart(partSection, 0, size);
988:     PetscSectionSetDof(partSection, 0, cEnd-cStart);
989:     PetscSectionSetUp(partSection);
990:     PetscMalloc1(cEnd-cStart, &points);
991:     for (c = cStart; c < cEnd; ++c) points[c] = c;
992:     ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
993:     return(0);
994:   }
995:   if (part->height == 0) {
996:     PetscInt numVertices = 0;
997:     PetscInt *start     = NULL;
998:     PetscInt *adjacency = NULL;
999:     IS       globalNumbering;

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

1007:       DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
1008:       DMPlexCreateNumbering_Internal(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
1009:       ISGetIndices(globalNumbering, &idxs);
1010:       for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
1011:       ISRestoreIndices(globalNumbering, &idxs);
1012:     }
1013:     if (part->viewGraph) {
1014:       PetscViewer viewer = part->viewerGraph;
1015:       PetscBool   isascii;
1016:       PetscInt    v, i;
1017:       PetscMPIInt rank;

1019:       MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
1020:       PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1021:       if (isascii) {
1022:         PetscViewerASCIIPushSynchronized(viewer);
1023:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Nv: %D\n", rank, numVertices);
1024:         for (v = 0; v < numVertices; ++v) {
1025:           const PetscInt s = start[v];
1026:           const PetscInt e = start[v+1];

1028:           PetscViewerASCIISynchronizedPrintf(viewer, "[%d]  ", rank);
1029:           for (i = s; i < e; ++i) {PetscViewerASCIISynchronizedPrintf(viewer, "%D ", adjacency[i]);}
1030:           PetscViewerASCIISynchronizedPrintf(viewer, "[%D-%D)\n", s, e);
1031:         }
1032:         PetscViewerFlush(viewer);
1033:         PetscViewerASCIIPopSynchronized(viewer);
1034:       }
1035:     }
1036:     if (!part->ops->partition) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_WRONGSTATE, "PetscPartitioner has no partitioning method");
1037:     (*part->ops->partition)(part, dm, size, numVertices, start, adjacency, partSection, partition);
1038:     PetscFree(start);
1039:     PetscFree(adjacency);
1040:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
1041:       const PetscInt *globalNum;
1042:       const PetscInt *partIdx;
1043:       PetscInt       *map, cStart, cEnd;
1044:       PetscInt       *adjusted, i, localSize, offset;
1045:       IS             newPartition;

1047:       ISGetLocalSize(*partition,&localSize);
1048:       PetscMalloc1(localSize,&adjusted);
1049:       ISGetIndices(globalNumbering,&globalNum);
1050:       ISGetIndices(*partition,&partIdx);
1051:       PetscMalloc1(localSize,&map);
1052:       DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
1053:       for (i = cStart, offset = 0; i < cEnd; i++) {
1054:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
1055:       }
1056:       for (i = 0; i < localSize; i++) {
1057:         adjusted[i] = map[partIdx[i]];
1058:       }
1059:       PetscFree(map);
1060:       ISRestoreIndices(*partition,&partIdx);
1061:       ISRestoreIndices(globalNumbering,&globalNum);
1062:       ISCreateGeneral(PETSC_COMM_SELF,localSize,adjusted,PETSC_OWN_POINTER,&newPartition);
1063:       ISDestroy(&globalNumbering);
1064:       ISDestroy(partition);
1065:       *partition = newPartition;
1066:     }
1067:   } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
1068:   PetscPartitionerViewFromOptions(part, NULL, "-petscpartitioner_view");
1069:   return(0);
1070: }

1072: static PetscErrorCode PetscPartitionerDestroy_Shell(PetscPartitioner part)
1073: {
1074:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1075:   PetscErrorCode          ierr;

1078:   PetscSectionDestroy(&p->section);
1079:   ISDestroy(&p->partition);
1080:   PetscFree(p);
1081:   return(0);
1082: }

1084: static PetscErrorCode PetscPartitionerView_Shell_Ascii(PetscPartitioner part, PetscViewer viewer)
1085: {
1086:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1087:   PetscErrorCode          ierr;

1090:   if (p->random) {
1091:     PetscViewerASCIIPushTab(viewer);
1092:     PetscViewerASCIIPrintf(viewer, "using random partition\n");
1093:     PetscViewerASCIIPopTab(viewer);
1094:   }
1095:   return(0);
1096: }

1098: static PetscErrorCode PetscPartitionerView_Shell(PetscPartitioner part, PetscViewer viewer)
1099: {
1100:   PetscBool      iascii;

1106:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1107:   if (iascii) {PetscPartitionerView_Shell_Ascii(part, viewer);}
1108:   return(0);
1109: }

1111: static PetscErrorCode PetscPartitionerSetFromOptions_Shell(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1112: {
1113:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1114:   PetscErrorCode          ierr;

1117:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner Shell Options");
1118:   PetscOptionsBool("-petscpartitioner_shell_random", "Use a random partition", "PetscPartitionerView", PETSC_FALSE, &p->random, NULL);
1119:   PetscOptionsTail();
1120:   return(0);
1121: }

1123: static PetscErrorCode PetscPartitionerPartition_Shell(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1124: {
1125:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1126:   PetscInt                np;
1127:   PetscErrorCode          ierr;

1130:   if (p->random) {
1131:     PetscRandom r;
1132:     PetscInt   *sizes, *points, v, p;
1133:     PetscMPIInt rank;

1135:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1136:     PetscRandomCreate(PETSC_COMM_SELF, &r);
1137:     PetscRandomSetInterval(r, 0.0, (PetscScalar) nparts);
1138:     PetscRandomSetFromOptions(r);
1139:     PetscCalloc2(nparts, &sizes, numVertices, &points);
1140:     for (v = 0; v < numVertices; ++v) {points[v] = v;}
1141:     for (p = 0; p < nparts; ++p) {sizes[p] = numVertices/nparts + (PetscInt) (p < numVertices % nparts);}
1142:     for (v = numVertices-1; v > 0; --v) {
1143:       PetscReal val;
1144:       PetscInt  w, tmp;

1146:       PetscRandomSetInterval(r, 0.0, (PetscScalar) (v+1));
1147:       PetscRandomGetValueReal(r, &val);
1148:       w    = PetscFloorReal(val);
1149:       tmp       = points[v];
1150:       points[v] = points[w];
1151:       points[w] = tmp;
1152:     }
1153:     PetscRandomDestroy(&r);
1154:     PetscPartitionerShellSetPartition(part, nparts, sizes, points);
1155:     PetscFree2(sizes, points);
1156:   }
1157:   if (!p->section) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Shell partitioner information not provided. Please call PetscPartitionerShellSetPartition()");
1158:   PetscSectionGetChart(p->section, NULL, &np);
1159:   if (nparts != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of requested partitions %d != configured partitions %d", nparts, np);
1160:   ISGetLocalSize(p->partition, &np);
1161:   if (numVertices != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input vertices %d != configured vertices %d", numVertices, np);
1162:   PetscSectionCopy(p->section, partSection);
1163:   *partition = p->partition;
1164:   PetscObjectReference((PetscObject) p->partition);
1165:   return(0);
1166: }

1168: static PetscErrorCode PetscPartitionerInitialize_Shell(PetscPartitioner part)
1169: {
1171:   part->noGraph             = PETSC_TRUE; /* PetscPartitionerShell cannot overload the partition call, so it is safe for now */
1172:   part->ops->view           = PetscPartitionerView_Shell;
1173:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_Shell;
1174:   part->ops->destroy        = PetscPartitionerDestroy_Shell;
1175:   part->ops->partition      = PetscPartitionerPartition_Shell;
1176:   return(0);
1177: }

1179: /*MC
1180:   PETSCPARTITIONERSHELL = "shell" - A PetscPartitioner object

1182:   Level: intermediate

1184: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1185: M*/

1187: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Shell(PetscPartitioner part)
1188: {
1189:   PetscPartitioner_Shell *p;
1190:   PetscErrorCode          ierr;

1194:   PetscNewLog(part, &p);
1195:   part->data = p;

1197:   PetscPartitionerInitialize_Shell(part);
1198:   p->random = PETSC_FALSE;
1199:   return(0);
1200: }

1202: /*@C
1203:   PetscPartitionerShellSetPartition - Set an artifical partition for a mesh

1205:   Collective on PetscPartitioner

1207:   Input Parameters:
1208: + part   - The PetscPartitioner
1209: . size   - The number of partitions
1210: . sizes  - array of size size (or NULL) providing the number of points in each partition
1211: - points - array of size 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.)

1213:   Level: developer

1215:   Notes:

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

1219: .seealso DMPlexDistribute(), PetscPartitionerCreate()
1220: @*/
1221: PetscErrorCode PetscPartitionerShellSetPartition(PetscPartitioner part, PetscInt size, const PetscInt sizes[], const PetscInt points[])
1222: {
1223:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1224:   PetscInt                proc, numPoints;
1225:   PetscErrorCode          ierr;

1231:   PetscSectionDestroy(&p->section);
1232:   ISDestroy(&p->partition);
1233:   PetscSectionCreate(PetscObjectComm((PetscObject) part), &p->section);
1234:   PetscSectionSetChart(p->section, 0, size);
1235:   if (sizes) {
1236:     for (proc = 0; proc < size; ++proc) {
1237:       PetscSectionSetDof(p->section, proc, sizes[proc]);
1238:     }
1239:   }
1240:   PetscSectionSetUp(p->section);
1241:   PetscSectionGetStorageSize(p->section, &numPoints);
1242:   ISCreateGeneral(PetscObjectComm((PetscObject) part), numPoints, points, PETSC_COPY_VALUES, &p->partition);
1243:   return(0);
1244: }

1246: /*@
1247:   PetscPartitionerShellSetRandom - Set the flag to use a random partition

1249:   Collective on PetscPartitioner

1251:   Input Parameters:
1252: + part   - The PetscPartitioner
1253: - random - The flag to use a random partition

1255:   Level: intermediate

1257: .seealso PetscPartitionerShellGetRandom(), PetscPartitionerCreate()
1258: @*/
1259: PetscErrorCode PetscPartitionerShellSetRandom(PetscPartitioner part, PetscBool random)
1260: {
1261:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;

1265:   p->random = random;
1266:   return(0);
1267: }

1269: /*@
1270:   PetscPartitionerShellGetRandom - get the flag to use a random partition

1272:   Collective on PetscPartitioner

1274:   Input Parameter:
1275: . part   - The PetscPartitioner

1277:   Output Parameter
1278: . random - The flag to use a random partition

1280:   Level: intermediate

1282: .seealso PetscPartitionerShellSetRandom(), PetscPartitionerCreate()
1283: @*/
1284: PetscErrorCode PetscPartitionerShellGetRandom(PetscPartitioner part, PetscBool *random)
1285: {
1286:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;

1291:   *random = p->random;
1292:   return(0);
1293: }

1295: static PetscErrorCode PetscPartitionerDestroy_Simple(PetscPartitioner part)
1296: {
1297:   PetscPartitioner_Simple *p = (PetscPartitioner_Simple *) part->data;
1298:   PetscErrorCode          ierr;

1301:   PetscFree(p);
1302:   return(0);
1303: }

1305: static PetscErrorCode PetscPartitionerView_Simple_Ascii(PetscPartitioner part, PetscViewer viewer)
1306: {
1308:   return(0);
1309: }

1311: static PetscErrorCode PetscPartitionerView_Simple(PetscPartitioner part, PetscViewer viewer)
1312: {
1313:   PetscBool      iascii;

1319:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1320:   if (iascii) {PetscPartitionerView_Simple_Ascii(part, viewer);}
1321:   return(0);
1322: }

1324: static PetscErrorCode PetscPartitionerPartition_Simple(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1325: {
1326:   MPI_Comm       comm;
1327:   PetscInt       np;
1328:   PetscMPIInt    size;

1332:   comm = PetscObjectComm((PetscObject)part);
1333:   MPI_Comm_size(comm,&size);
1334:   PetscSectionSetChart(partSection, 0, nparts);
1335:   ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1336:   if (size == 1) {
1337:     for (np = 0; np < nparts; ++np) {PetscSectionSetDof(partSection, np, numVertices/nparts + ((numVertices % nparts) > np));}
1338:   } else {
1339:     PetscMPIInt rank;
1340:     PetscInt nvGlobal, *offsets, myFirst, myLast;

1342:     PetscMalloc1(size+1,&offsets);
1343:     offsets[0] = 0;
1344:     MPI_Allgather(&numVertices,1,MPIU_INT,&offsets[1],1,MPIU_INT,comm);
1345:     for (np = 2; np <= size; np++) {
1346:       offsets[np] += offsets[np-1];
1347:     }
1348:     nvGlobal = offsets[size];
1349:     MPI_Comm_rank(comm,&rank);
1350:     myFirst = offsets[rank];
1351:     myLast  = offsets[rank + 1] - 1;
1352:     PetscFree(offsets);
1353:     if (numVertices) {
1354:       PetscInt firstPart = 0, firstLargePart = 0;
1355:       PetscInt lastPart = 0, lastLargePart = 0;
1356:       PetscInt rem = nvGlobal % nparts;
1357:       PetscInt pSmall = nvGlobal/nparts;
1358:       PetscInt pBig = nvGlobal/nparts + 1;


1361:       if (rem) {
1362:         firstLargePart = myFirst / pBig;
1363:         lastLargePart  = myLast  / pBig;

1365:         if (firstLargePart < rem) {
1366:           firstPart = firstLargePart;
1367:         } else {
1368:           firstPart = rem + (myFirst - (rem * pBig)) / pSmall;
1369:         }
1370:         if (lastLargePart < rem) {
1371:           lastPart = lastLargePart;
1372:         } else {
1373:           lastPart = rem + (myLast - (rem * pBig)) / pSmall;
1374:         }
1375:       } else {
1376:         firstPart = myFirst / (nvGlobal/nparts);
1377:         lastPart  = myLast  / (nvGlobal/nparts);
1378:       }

1380:       for (np = firstPart; np <= lastPart; np++) {
1381:         PetscInt PartStart =  np    * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np);
1382:         PetscInt PartEnd   = (np+1) * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np+1);

1384:         PartStart = PetscMax(PartStart,myFirst);
1385:         PartEnd   = PetscMin(PartEnd,myLast+1);
1386:         PetscSectionSetDof(partSection,np,PartEnd-PartStart);
1387:       }
1388:     }
1389:   }
1390:   PetscSectionSetUp(partSection);
1391:   return(0);
1392: }

1394: static PetscErrorCode PetscPartitionerInitialize_Simple(PetscPartitioner part)
1395: {
1397:   part->noGraph        = PETSC_TRUE;
1398:   part->ops->view      = PetscPartitionerView_Simple;
1399:   part->ops->destroy   = PetscPartitionerDestroy_Simple;
1400:   part->ops->partition = PetscPartitionerPartition_Simple;
1401:   return(0);
1402: }

1404: /*MC
1405:   PETSCPARTITIONERSIMPLE = "simple" - A PetscPartitioner object

1407:   Level: intermediate

1409: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1410: M*/

1412: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Simple(PetscPartitioner part)
1413: {
1414:   PetscPartitioner_Simple *p;
1415:   PetscErrorCode           ierr;

1419:   PetscNewLog(part, &p);
1420:   part->data = p;

1422:   PetscPartitionerInitialize_Simple(part);
1423:   return(0);
1424: }

1426: static PetscErrorCode PetscPartitionerDestroy_Gather(PetscPartitioner part)
1427: {
1428:   PetscPartitioner_Gather *p = (PetscPartitioner_Gather *) part->data;
1429:   PetscErrorCode          ierr;

1432:   PetscFree(p);
1433:   return(0);
1434: }

1436: static PetscErrorCode PetscPartitionerView_Gather_Ascii(PetscPartitioner part, PetscViewer viewer)
1437: {
1439:   return(0);
1440: }

1442: static PetscErrorCode PetscPartitionerView_Gather(PetscPartitioner part, PetscViewer viewer)
1443: {
1444:   PetscBool      iascii;

1450:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1451:   if (iascii) {PetscPartitionerView_Gather_Ascii(part, viewer);}
1452:   return(0);
1453: }

1455: static PetscErrorCode PetscPartitionerPartition_Gather(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1456: {
1457:   PetscInt       np;

1461:   PetscSectionSetChart(partSection, 0, nparts);
1462:   ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1463:   PetscSectionSetDof(partSection,0,numVertices);
1464:   for (np = 1; np < nparts; ++np) {PetscSectionSetDof(partSection, np, 0);}
1465:   PetscSectionSetUp(partSection);
1466:   return(0);
1467: }

1469: static PetscErrorCode PetscPartitionerInitialize_Gather(PetscPartitioner part)
1470: {
1472:   part->noGraph        = PETSC_TRUE;
1473:   part->ops->view      = PetscPartitionerView_Gather;
1474:   part->ops->destroy   = PetscPartitionerDestroy_Gather;
1475:   part->ops->partition = PetscPartitionerPartition_Gather;
1476:   return(0);
1477: }

1479: /*MC
1480:   PETSCPARTITIONERGATHER = "gather" - A PetscPartitioner object

1482:   Level: intermediate

1484: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1485: M*/

1487: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Gather(PetscPartitioner part)
1488: {
1489:   PetscPartitioner_Gather *p;
1490:   PetscErrorCode           ierr;

1494:   PetscNewLog(part, &p);
1495:   part->data = p;

1497:   PetscPartitionerInitialize_Gather(part);
1498:   return(0);
1499: }


1502: static PetscErrorCode PetscPartitionerDestroy_Chaco(PetscPartitioner part)
1503: {
1504:   PetscPartitioner_Chaco *p = (PetscPartitioner_Chaco *) part->data;
1505:   PetscErrorCode          ierr;

1508:   PetscFree(p);
1509:   return(0);
1510: }

1512: static PetscErrorCode PetscPartitionerView_Chaco_Ascii(PetscPartitioner part, PetscViewer viewer)
1513: {
1515:   return(0);
1516: }

1518: static PetscErrorCode PetscPartitionerView_Chaco(PetscPartitioner part, PetscViewer viewer)
1519: {
1520:   PetscBool      iascii;

1526:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1527:   if (iascii) {PetscPartitionerView_Chaco_Ascii(part, viewer);}
1528:   return(0);
1529: }

1531: #if defined(PETSC_HAVE_CHACO)
1532: #if defined(PETSC_HAVE_UNISTD_H)
1533: #include <unistd.h>
1534: #endif
1535: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1536: #include <chaco.h>
1537: #else
1538: /* Older versions of Chaco do not have an include file */
1539: PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
1540:                        float *ewgts, float *x, float *y, float *z, char *outassignname,
1541:                        char *outfilename, short *assignment, int architecture, int ndims_tot,
1542:                        int mesh_dims[3], double *goal, int global_method, int local_method,
1543:                        int rqi_flag, int vmax, int ndims, double eigtol, long seed);
1544: #endif
1545: extern int FREE_GRAPH;
1546: #endif

1548: static PetscErrorCode PetscPartitionerPartition_Chaco(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1549: {
1550: #if defined(PETSC_HAVE_CHACO)
1551:   enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
1552:   MPI_Comm       comm;
1553:   int            nvtxs          = numVertices; /* number of vertices in full graph */
1554:   int           *vwgts          = NULL;   /* weights for all vertices */
1555:   float         *ewgts          = NULL;   /* weights for all edges */
1556:   float         *x              = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
1557:   char          *outassignname  = NULL;   /*  name of assignment output file */
1558:   char          *outfilename    = NULL;   /* output file name */
1559:   int            architecture   = 1;      /* 0 => hypercube, d => d-dimensional mesh */
1560:   int            ndims_tot      = 0;      /* total number of cube dimensions to divide */
1561:   int            mesh_dims[3];            /* dimensions of mesh of processors */
1562:   double        *goal          = NULL;    /* desired set sizes for each set */
1563:   int            global_method = 1;       /* global partitioning algorithm */
1564:   int            local_method  = 1;       /* local partitioning algorithm */
1565:   int            rqi_flag      = 0;       /* should I use RQI/Symmlq eigensolver? */
1566:   int            vmax          = 200;     /* how many vertices to coarsen down to? */
1567:   int            ndims         = 1;       /* number of eigenvectors (2^d sets) */
1568:   double         eigtol        = 0.001;   /* tolerance on eigenvectors */
1569:   long           seed          = 123636512; /* for random graph mutations */
1570: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1571:   int           *assignment;              /* Output partition */
1572: #else
1573:   short int     *assignment;              /* Output partition */
1574: #endif
1575:   int            fd_stdout, fd_pipe[2];
1576:   PetscInt      *points;
1577:   int            i, v, p;

1581:   PetscObjectGetComm((PetscObject)dm,&comm);
1582: #if defined (PETSC_USE_DEBUG)
1583:   {
1584:     int ival,isum;
1585:     PetscBool distributed;

1587:     ival = (numVertices > 0);
1588:     MPI_Allreduce(&ival, &isum, 1, MPI_INT, MPI_SUM, comm);
1589:     distributed = (isum > 1) ? PETSC_TRUE : PETSC_FALSE;
1590:     if (distributed) SETERRQ(comm, PETSC_ERR_SUP, "Chaco cannot partition a distributed graph");
1591:   }
1592: #endif
1593:   if (!numVertices) {
1594:     PetscSectionSetChart(partSection, 0, nparts);
1595:     PetscSectionSetUp(partSection);
1596:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
1597:     return(0);
1598:   }
1599:   FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
1600:   for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];

1602:   if (global_method == INERTIAL_METHOD) {
1603:     /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
1604:     SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
1605:   }
1606:   mesh_dims[0] = nparts;
1607:   mesh_dims[1] = 1;
1608:   mesh_dims[2] = 1;
1609:   PetscMalloc1(nvtxs, &assignment);
1610:   /* Chaco outputs to stdout. We redirect this to a buffer. */
1611:   /* TODO: check error codes for UNIX calls */
1612: #if defined(PETSC_HAVE_UNISTD_H)
1613:   {
1614:     int piperet;
1615:     piperet = pipe(fd_pipe);
1616:     if (piperet) SETERRQ(comm,PETSC_ERR_SYS,"Could not create pipe");
1617:     fd_stdout = dup(1);
1618:     close(1);
1619:     dup2(fd_pipe[1], 1);
1620:   }
1621: #endif
1622:   interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
1623:                    assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
1624:                    vmax, ndims, eigtol, seed);
1625: #if defined(PETSC_HAVE_UNISTD_H)
1626:   {
1627:     char msgLog[10000];
1628:     int  count;

1630:     fflush(stdout);
1631:     count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
1632:     if (count < 0) count = 0;
1633:     msgLog[count] = 0;
1634:     close(1);
1635:     dup2(fd_stdout, 1);
1636:     close(fd_stdout);
1637:     close(fd_pipe[0]);
1638:     close(fd_pipe[1]);
1639:     if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
1640:   }
1641: #else
1642:   if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", "error in stdout");
1643: #endif
1644:   /* Convert to PetscSection+IS */
1645:   PetscSectionSetChart(partSection, 0, nparts);
1646:   for (v = 0; v < nvtxs; ++v) {
1647:     PetscSectionAddDof(partSection, assignment[v], 1);
1648:   }
1649:   PetscSectionSetUp(partSection);
1650:   PetscMalloc1(nvtxs, &points);
1651:   for (p = 0, i = 0; p < nparts; ++p) {
1652:     for (v = 0; v < nvtxs; ++v) {
1653:       if (assignment[v] == p) points[i++] = v;
1654:     }
1655:   }
1656:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1657:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
1658:   if (global_method == INERTIAL_METHOD) {
1659:     /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
1660:   }
1661:   PetscFree(assignment);
1662:   for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
1663:   return(0);
1664: #else
1665:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco.");
1666: #endif
1667: }

1669: static PetscErrorCode PetscPartitionerInitialize_Chaco(PetscPartitioner part)
1670: {
1672:   part->noGraph        = PETSC_FALSE;
1673:   part->ops->view      = PetscPartitionerView_Chaco;
1674:   part->ops->destroy   = PetscPartitionerDestroy_Chaco;
1675:   part->ops->partition = PetscPartitionerPartition_Chaco;
1676:   return(0);
1677: }

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

1682:   Level: intermediate

1684: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1685: M*/

1687: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Chaco(PetscPartitioner part)
1688: {
1689:   PetscPartitioner_Chaco *p;
1690:   PetscErrorCode          ierr;

1694:   PetscNewLog(part, &p);
1695:   part->data = p;

1697:   PetscPartitionerInitialize_Chaco(part);
1698:   PetscCitationsRegister(ChacoPartitionerCitation, &ChacoPartitionercite);
1699:   return(0);
1700: }

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

1704: static PetscErrorCode PetscPartitionerDestroy_ParMetis(PetscPartitioner part)
1705: {
1706:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1707:   PetscErrorCode             ierr;

1710:   PetscFree(p);
1711:   return(0);
1712: }

1714: static PetscErrorCode PetscPartitionerView_ParMetis_Ascii(PetscPartitioner part, PetscViewer viewer)
1715: {
1716:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1717:   PetscErrorCode             ierr;

1720:   PetscViewerASCIIPushTab(viewer);
1721:   PetscViewerASCIIPrintf(viewer, "ParMetis type: %s\n", ptypes[p->ptype]);
1722:   PetscViewerASCIIPrintf(viewer, "load imbalance ratio %g\n", (double) p->imbalanceRatio);
1723:   PetscViewerASCIIPrintf(viewer, "debug flag %D\n", p->debugFlag);
1724:   PetscViewerASCIIPrintf(viewer, "random seed %D\n", p->randomSeed);
1725:   PetscViewerASCIIPopTab(viewer);
1726:   return(0);
1727: }

1729: static PetscErrorCode PetscPartitionerView_ParMetis(PetscPartitioner part, PetscViewer viewer)
1730: {
1731:   PetscBool      iascii;

1737:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1738:   if (iascii) {PetscPartitionerView_ParMetis_Ascii(part, viewer);}
1739:   return(0);
1740: }

1742: static PetscErrorCode PetscPartitionerSetFromOptions_ParMetis(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1743: {
1744:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1745:   PetscErrorCode            ierr;

1748:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner ParMetis Options");
1749:   PetscOptionsEList("-petscpartitioner_parmetis_type", "Partitioning method", "", ptypes, 2, ptypes[p->ptype], &p->ptype, NULL);
1750:   PetscOptionsReal("-petscpartitioner_parmetis_imbalance_ratio", "Load imbalance ratio limit", "", p->imbalanceRatio, &p->imbalanceRatio, NULL);
1751:   PetscOptionsInt("-petscpartitioner_parmetis_debug", "Debugging flag", "", p->debugFlag, &p->debugFlag, NULL);
1752:   PetscOptionsInt("-petscpartitioner_parmetis_seed", "Random seed", "", p->randomSeed, &p->randomSeed, NULL);
1753:   PetscOptionsTail();
1754:   return(0);
1755: }

1757: #if defined(PETSC_HAVE_PARMETIS)
1758: #include <parmetis.h>
1759: #endif

1761: static PetscErrorCode PetscPartitionerPartition_ParMetis(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1762: {
1763: #if defined(PETSC_HAVE_PARMETIS)
1764:   PetscPartitioner_ParMetis *pm = (PetscPartitioner_ParMetis *) part->data;
1765:   MPI_Comm       comm;
1766:   PetscSection   section;
1767:   PetscInt       nvtxs       = numVertices; /* The number of vertices in full graph */
1768:   PetscInt      *vtxdist;                   /* Distribution of vertices across processes */
1769:   PetscInt      *xadj        = start;       /* Start of edge list for each vertex */
1770:   PetscInt      *adjncy      = adjacency;   /* Edge lists for all vertices */
1771:   PetscInt      *vwgt        = NULL;        /* Vertex weights */
1772:   PetscInt      *adjwgt      = NULL;        /* Edge weights */
1773:   PetscInt       wgtflag     = 0;           /* Indicates which weights are present */
1774:   PetscInt       numflag     = 0;           /* Indicates initial offset (0 or 1) */
1775:   PetscInt       ncon        = 1;           /* The number of weights per vertex */
1776:   PetscInt       metis_ptype = pm->ptype;   /* kway or recursive bisection */
1777:   real_t        *tpwgts;                    /* The fraction of vertex weights assigned to each partition */
1778:   real_t        *ubvec;                     /* The balance intolerance for vertex weights */
1779:   PetscInt       options[64];               /* Options */
1780:   /* Outputs */
1781:   PetscInt       v, i, *assignment, *points;
1782:   PetscMPIInt    size, rank, p;

1786:   PetscObjectGetComm((PetscObject) part, &comm);
1787:   MPI_Comm_size(comm, &size);
1788:   MPI_Comm_rank(comm, &rank);
1789:   /* Calculate vertex distribution */
1790:   PetscMalloc5(size+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment,nvtxs,&vwgt);
1791:   vtxdist[0] = 0;
1792:   MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
1793:   for (p = 2; p <= size; ++p) {
1794:     vtxdist[p] += vtxdist[p-1];
1795:   }
1796:   /* Calculate partition weights */
1797:   for (p = 0; p < nparts; ++p) {
1798:     tpwgts[p] = 1.0/nparts;
1799:   }
1800:   ubvec[0] = pm->imbalanceRatio;
1801:   /* Weight cells by dofs on cell by default */
1802:   DMGetLocalSection(dm, &section);
1803:   for (v = 0; v < nvtxs; ++v) vwgt[v] = 1;
1804:   if (section) {
1805:     PetscInt cStart, cEnd, dof;

1807:     /* WARNING: Assumes that meshes with overlap have the overlapped cells at the end of the stratum. */
1808:     /* To do this properly, we should use the cell numbering created in DMPlexCreatePartitionerGraph. */
1809:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
1810:     for (v = cStart; v < cStart + numVertices; ++v) {
1811:       PetscSectionGetDof(section, v, &dof);
1812:       vwgt[v-cStart] = PetscMax(dof, 1);
1813:     }
1814:   }
1815:   wgtflag |= 2; /* have weights on graph vertices */

1817:   if (nparts == 1) {
1818:     PetscArrayzero(assignment, nvtxs);
1819:   } else {
1820:     for (p = 0; !vtxdist[p+1] && p < size; ++p);
1821:     if (vtxdist[p+1] == vtxdist[size]) {
1822:       if (rank == p) {
1823:         METIS_SetDefaultOptions(options); /* initialize all defaults */
1824:         options[METIS_OPTION_DBGLVL] = pm->debugFlag;
1825:         options[METIS_OPTION_SEED]   = pm->randomSeed;
1826:         if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
1827:         if (metis_ptype == 1) {
1828:           PetscStackPush("METIS_PartGraphRecursive");
1829:           METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
1830:           PetscStackPop;
1831:           if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphRecursive()");
1832:         } else {
1833:           /*
1834:            It would be nice to activate the two options below, but they would need some actual testing.
1835:            - Turning on these options may exercise path of the METIS code that have bugs and may break production runs.
1836:            - 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.
1837:           */
1838:           /* options[METIS_OPTION_CONTIG]  = 1; */ /* try to produce partitions that are contiguous */
1839:           /* options[METIS_OPTION_MINCONN] = 1; */ /* minimize the maximum degree of the subdomain graph */
1840:           PetscStackPush("METIS_PartGraphKway");
1841:           METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
1842:           PetscStackPop;
1843:           if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
1844:         }
1845:       }
1846:     } else {
1847:       options[0] = 1; /*use options */
1848:       options[1] = pm->debugFlag;
1849:       options[2] = (pm->randomSeed == -1) ? 15 : pm->randomSeed; /* default is GLOBAL_SEED=15 from `libparmetis/defs.h` */
1850:       PetscStackPush("ParMETIS_V3_PartKway");
1851:       ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment, &comm);
1852:       PetscStackPop;
1853:       if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error %d in ParMETIS_V3_PartKway()", ierr);
1854:     }
1855:   }
1856:   /* Convert to PetscSection+IS */
1857:   PetscSectionSetChart(partSection, 0, nparts);
1858:   for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
1859:   PetscSectionSetUp(partSection);
1860:   PetscMalloc1(nvtxs, &points);
1861:   for (p = 0, i = 0; p < nparts; ++p) {
1862:     for (v = 0; v < nvtxs; ++v) {
1863:       if (assignment[v] == p) points[i++] = v;
1864:     }
1865:   }
1866:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1867:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
1868:   PetscFree5(vtxdist,tpwgts,ubvec,assignment,vwgt);
1869:   return(0);
1870: #else
1871:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
1872: #endif
1873: }

1875: static PetscErrorCode PetscPartitionerInitialize_ParMetis(PetscPartitioner part)
1876: {
1878:   part->noGraph             = PETSC_FALSE;
1879:   part->ops->view           = PetscPartitionerView_ParMetis;
1880:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_ParMetis;
1881:   part->ops->destroy        = PetscPartitionerDestroy_ParMetis;
1882:   part->ops->partition      = PetscPartitionerPartition_ParMetis;
1883:   return(0);
1884: }

1886: /*MC
1887:   PETSCPARTITIONERPARMETIS = "parmetis" - A PetscPartitioner object using the ParMetis library

1889:   Level: intermediate

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

1894: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMetis(PetscPartitioner part)
1895: {
1896:   PetscPartitioner_ParMetis *p;
1897:   PetscErrorCode          ierr;

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

1904:   p->ptype          = 0;
1905:   p->imbalanceRatio = 1.05;
1906:   p->debugFlag      = 0;
1907:   p->randomSeed     = -1; /* defaults to GLOBAL_SEED=15 from `libparmetis/defs.h` */

1909:   PetscPartitionerInitialize_ParMetis(part);
1910:   PetscCitationsRegister(ParMetisPartitionerCitation, &ParMetisPartitionercite);
1911:   return(0);
1912: }

1914: PetscBool PTScotchPartitionercite = PETSC_FALSE;
1915: const char PTScotchPartitionerCitation[] =
1916:   "@article{PTSCOTCH,\n"
1917:   "  author  = {C. Chevalier and F. Pellegrini},\n"
1918:   "  title   = {{PT-SCOTCH}: a tool for efficient parallel graph ordering},\n"
1919:   "  journal = {Parallel Computing},\n"
1920:   "  volume  = {34},\n"
1921:   "  number  = {6},\n"
1922:   "  pages   = {318--331},\n"
1923:   "  year    = {2008},\n"
1924:   "  doi     = {https://doi.org/10.1016/j.parco.2007.12.001}\n"
1925:   "}\n";

1927: #if defined(PETSC_HAVE_PTSCOTCH)

1929: EXTERN_C_BEGIN
1930: #include <ptscotch.h>
1931: EXTERN_C_END

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

1935: static int PTScotch_Strategy(PetscInt strategy)
1936: {
1937:   switch (strategy) {
1938:   case  0: return SCOTCH_STRATDEFAULT;
1939:   case  1: return SCOTCH_STRATQUALITY;
1940:   case  2: return SCOTCH_STRATSPEED;
1941:   case  3: return SCOTCH_STRATBALANCE;
1942:   case  4: return SCOTCH_STRATSAFETY;
1943:   case  5: return SCOTCH_STRATSCALABILITY;
1944:   case  6: return SCOTCH_STRATRECURSIVE;
1945:   case  7: return SCOTCH_STRATREMAP;
1946:   default: return SCOTCH_STRATDEFAULT;
1947:   }
1948: }


1951: static PetscErrorCode PTScotch_PartGraph_Seq(SCOTCH_Num strategy, double imbalance, SCOTCH_Num n, SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
1952:                                              SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num part[])
1953: {
1954:   SCOTCH_Graph   grafdat;
1955:   SCOTCH_Strat   stradat;
1956:   SCOTCH_Num     vertnbr = n;
1957:   SCOTCH_Num     edgenbr = xadj[n];
1958:   SCOTCH_Num*    velotab = vtxwgt;
1959:   SCOTCH_Num*    edlotab = adjwgt;
1960:   SCOTCH_Num     flagval = strategy;
1961:   double         kbalval = imbalance;

1965:   {
1966:     PetscBool flg = PETSC_TRUE;
1967:     PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
1968:     if (!flg) velotab = NULL;
1969:   }
1970:   SCOTCH_graphInit(&grafdat);CHKERRPTSCOTCH(ierr);
1971:   SCOTCH_graphBuild(&grafdat, 0, vertnbr, xadj, xadj + 1, velotab, NULL, edgenbr, adjncy, edlotab);CHKERRPTSCOTCH(ierr);
1972:   SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
1973:   SCOTCH_stratGraphMapBuild(&stradat, flagval, nparts, kbalval);CHKERRPTSCOTCH(ierr);
1974: #if defined(PETSC_USE_DEBUG)
1975:   SCOTCH_graphCheck(&grafdat);CHKERRPTSCOTCH(ierr);
1976: #endif
1977:   SCOTCH_graphPart(&grafdat, nparts, &stradat, part);CHKERRPTSCOTCH(ierr);
1978:   SCOTCH_stratExit(&stradat);
1979:   SCOTCH_graphExit(&grafdat);
1980:   return(0);
1981: }

1983: static PetscErrorCode PTScotch_PartGraph_MPI(SCOTCH_Num strategy, double imbalance, SCOTCH_Num vtxdist[], SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
1984:                                              SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num part[], MPI_Comm comm)
1985: {
1986:   PetscMPIInt     procglbnbr;
1987:   PetscMPIInt     proclocnum;
1988:   SCOTCH_Arch     archdat;
1989:   SCOTCH_Dgraph   grafdat;
1990:   SCOTCH_Dmapping mappdat;
1991:   SCOTCH_Strat    stradat;
1992:   SCOTCH_Num      vertlocnbr;
1993:   SCOTCH_Num      edgelocnbr;
1994:   SCOTCH_Num*     veloloctab = vtxwgt;
1995:   SCOTCH_Num*     edloloctab = adjwgt;
1996:   SCOTCH_Num      flagval = strategy;
1997:   double          kbalval = imbalance;
1998:   PetscErrorCode  ierr;

2001:   {
2002:     PetscBool flg = PETSC_TRUE;
2003:     PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
2004:     if (!flg) veloloctab = NULL;
2005:   }
2006:   MPI_Comm_size(comm, &procglbnbr);
2007:   MPI_Comm_rank(comm, &proclocnum);
2008:   vertlocnbr = vtxdist[proclocnum + 1] - vtxdist[proclocnum];
2009:   edgelocnbr = xadj[vertlocnbr];

2011:   SCOTCH_dgraphInit(&grafdat, comm);CHKERRPTSCOTCH(ierr);
2012:   SCOTCH_dgraphBuild(&grafdat, 0, vertlocnbr, vertlocnbr, xadj, xadj + 1, veloloctab, NULL, edgelocnbr, edgelocnbr, adjncy, NULL, edloloctab);CHKERRPTSCOTCH(ierr);
2013: #if defined(PETSC_USE_DEBUG)
2014:   SCOTCH_dgraphCheck(&grafdat);CHKERRPTSCOTCH(ierr);
2015: #endif
2016:   SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
2017:   SCOTCH_stratDgraphMapBuild(&stradat, flagval, procglbnbr, nparts, kbalval);
2018:   SCOTCH_archInit(&archdat);CHKERRPTSCOTCH(ierr);
2019:   SCOTCH_archCmplt(&archdat, nparts);CHKERRPTSCOTCH(ierr);
2020:   SCOTCH_dgraphMapInit(&grafdat, &mappdat, &archdat, part);CHKERRPTSCOTCH(ierr);

2022:   SCOTCH_dgraphMapCompute(&grafdat, &mappdat, &stradat);CHKERRPTSCOTCH(ierr);
2023:   SCOTCH_dgraphMapExit(&grafdat, &mappdat);
2024:   SCOTCH_archExit(&archdat);
2025:   SCOTCH_stratExit(&stradat);
2026:   SCOTCH_dgraphExit(&grafdat);
2027:   return(0);
2028: }

2030: #endif /* PETSC_HAVE_PTSCOTCH */

2032: static PetscErrorCode PetscPartitionerDestroy_PTScotch(PetscPartitioner part)
2033: {
2034:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2035:   PetscErrorCode             ierr;

2038:   PetscFree(p);
2039:   return(0);
2040: }

2042: static PetscErrorCode PetscPartitionerView_PTScotch_Ascii(PetscPartitioner part, PetscViewer viewer)
2043: {
2044:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2045:   PetscErrorCode            ierr;

2048:   PetscViewerASCIIPushTab(viewer);
2049:   PetscViewerASCIIPrintf(viewer, "using partitioning strategy %s\n",PTScotchStrategyList[p->strategy]);
2050:   PetscViewerASCIIPrintf(viewer, "using load imbalance ratio %g\n",(double)p->imbalance);
2051:   PetscViewerASCIIPopTab(viewer);
2052:   return(0);
2053: }

2055: static PetscErrorCode PetscPartitionerView_PTScotch(PetscPartitioner part, PetscViewer viewer)
2056: {
2057:   PetscBool      iascii;

2063:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2064:   if (iascii) {PetscPartitionerView_PTScotch_Ascii(part, viewer);}
2065:   return(0);
2066: }

2068: static PetscErrorCode PetscPartitionerSetFromOptions_PTScotch(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
2069: {
2070:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2071:   const char *const         *slist = PTScotchStrategyList;
2072:   PetscInt                  nlist = (PetscInt)(sizeof(PTScotchStrategyList)/sizeof(PTScotchStrategyList[0]));
2073:   PetscBool                 flag;
2074:   PetscErrorCode            ierr;

2077:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner PTScotch Options");
2078:   PetscOptionsEList("-petscpartitioner_ptscotch_strategy","Partitioning strategy","",slist,nlist,slist[p->strategy],&p->strategy,&flag);
2079:   PetscOptionsReal("-petscpartitioner_ptscotch_imbalance","Load imbalance ratio","",p->imbalance,&p->imbalance,&flag);
2080:   PetscOptionsTail();
2081:   return(0);
2082: }

2084: static PetscErrorCode PetscPartitionerPartition_PTScotch(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
2085: {
2086: #if defined(PETSC_HAVE_PTSCOTCH)
2087:   MPI_Comm       comm       = PetscObjectComm((PetscObject)part);
2088:   PetscInt       nvtxs      = numVertices; /* The number of vertices in full graph */
2089:   PetscInt      *vtxdist;                  /* Distribution of vertices across processes */
2090:   PetscInt      *xadj       = start;       /* Start of edge list for each vertex */
2091:   PetscInt      *adjncy     = adjacency;   /* Edge lists for all vertices */
2092:   PetscInt      *vwgt       = NULL;        /* Vertex weights */
2093:   PetscInt      *adjwgt     = NULL;        /* Edge weights */
2094:   PetscInt       v, i, *assignment, *points;
2095:   PetscMPIInt    size, rank, p;

2099:   MPI_Comm_size(comm, &size);
2100:   MPI_Comm_rank(comm, &rank);
2101:   PetscMalloc2(size+1,&vtxdist,PetscMax(nvtxs,1),&assignment);

2103:   /* Calculate vertex distribution */
2104:   vtxdist[0] = 0;
2105:   MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
2106:   for (p = 2; p <= size; ++p) {
2107:     vtxdist[p] += vtxdist[p-1];
2108:   }

2110:   if (nparts == 1) {
2111:     PetscArrayzero(assignment, nvtxs);
2112:   } else { /* Weight cells by dofs on cell by default */
2113:     PetscSection section;

2115:     /* WARNING: Assumes that meshes with overlap have the overlapped cells at the end of the stratum. */
2116:     /* To do this properly, we should use the cell numbering created in DMPlexCreatePartitionerGraph. */
2117:     PetscMalloc1(PetscMax(nvtxs,1),&vwgt);
2118:     for (v = 0; v < PetscMax(nvtxs,1); ++v) vwgt[v] = 1;
2119:     DMGetLocalSection(dm, &section);
2120:     if (section) {
2121:       PetscInt vStart, vEnd, dof;
2122:       DMPlexGetHeightStratum(dm, part->height, &vStart, &vEnd);
2123:       for (v = vStart; v < vStart + numVertices; ++v) {
2124:         PetscSectionGetDof(section, v, &dof);
2125:         vwgt[v-vStart] = PetscMax(dof, 1);
2126:       }
2127:     }
2128:     {
2129:       PetscPartitioner_PTScotch *pts = (PetscPartitioner_PTScotch *) part->data;
2130:       int                       strat = PTScotch_Strategy(pts->strategy);
2131:       double                    imbal = (double)pts->imbalance;

2133:       for (p = 0; !vtxdist[p+1] && p < size; ++p);
2134:       if (vtxdist[p+1] == vtxdist[size]) {
2135:         if (rank == p) {
2136:           PTScotch_PartGraph_Seq(strat, imbal, nvtxs, xadj, adjncy, vwgt, adjwgt, nparts, assignment);
2137:         }
2138:       } else {
2139:         PTScotch_PartGraph_MPI(strat, imbal, vtxdist, xadj, adjncy, vwgt, adjwgt, nparts, assignment, comm);
2140:       }
2141:     }
2142:     PetscFree(vwgt);
2143:   }

2145:   /* Convert to PetscSection+IS */
2146:   PetscSectionSetChart(partSection, 0, nparts);
2147:   for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
2148:   PetscSectionSetUp(partSection);
2149:   PetscMalloc1(nvtxs, &points);
2150:   for (p = 0, i = 0; p < nparts; ++p) {
2151:     for (v = 0; v < nvtxs; ++v) {
2152:       if (assignment[v] == p) points[i++] = v;
2153:     }
2154:   }
2155:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
2156:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);

2158:   PetscFree2(vtxdist,assignment);
2159:   return(0);
2160: #else
2161:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-ptscotch.");
2162: #endif
2163: }

2165: static PetscErrorCode PetscPartitionerInitialize_PTScotch(PetscPartitioner part)
2166: {
2168:   part->noGraph             = PETSC_FALSE;
2169:   part->ops->view           = PetscPartitionerView_PTScotch;
2170:   part->ops->destroy        = PetscPartitionerDestroy_PTScotch;
2171:   part->ops->partition      = PetscPartitionerPartition_PTScotch;
2172:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_PTScotch;
2173:   return(0);
2174: }

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

2179:   Level: intermediate

2181: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
2182: M*/

2184: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_PTScotch(PetscPartitioner part)
2185: {
2186:   PetscPartitioner_PTScotch *p;
2187:   PetscErrorCode          ierr;

2191:   PetscNewLog(part, &p);
2192:   part->data = p;

2194:   p->strategy  = 0;
2195:   p->imbalance = 0.01;

2197:   PetscPartitionerInitialize_PTScotch(part);
2198:   PetscCitationsRegister(PTScotchPartitionerCitation, &PTScotchPartitionercite);
2199:   return(0);
2200: }


2203: /*@
2204:   DMPlexGetPartitioner - Get the mesh partitioner

2206:   Not collective

2208:   Input Parameter:
2209: . dm - The DM

2211:   Output Parameter:
2212: . part - The PetscPartitioner

2214:   Level: developer

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

2218: .seealso DMPlexDistribute(), DMPlexSetPartitioner(), PetscPartitionerCreate()
2219: @*/
2220: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
2221: {
2222:   DM_Plex       *mesh = (DM_Plex *) dm->data;

2227:   *part = mesh->partitioner;
2228:   return(0);
2229: }

2231: /*@
2232:   DMPlexSetPartitioner - Set the mesh partitioner

2234:   logically collective on DM

2236:   Input Parameters:
2237: + dm - The DM
2238: - part - The partitioner

2240:   Level: developer

2242:   Note: Any existing PetscPartitioner will be destroyed.

2244: .seealso DMPlexDistribute(), DMPlexGetPartitioner(), PetscPartitionerCreate()
2245: @*/
2246: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
2247: {
2248:   DM_Plex       *mesh = (DM_Plex *) dm->data;

2254:   PetscObjectReference((PetscObject)part);
2255:   PetscPartitionerDestroy(&mesh->partitioner);
2256:   mesh->partitioner = part;
2257:   return(0);
2258: }

2260: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
2261: {
2262:   const PetscInt *cone;
2263:   PetscInt       coneSize, c;
2264:   PetscBool      missing;

2268:   PetscHSetIQueryAdd(ht, point, &missing);
2269:   if (missing) {
2270:     DMPlexGetCone(dm, point, &cone);
2271:     DMPlexGetConeSize(dm, point, &coneSize);
2272:     for (c = 0; c < coneSize; c++) {
2273:       DMPlexAddClosure_Private(dm, ht, cone[c]);
2274:     }
2275:   }
2276:   return(0);
2277: }

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

2284:   if (up) {
2285:     PetscInt parent;

2287:     DMPlexGetTreeParent(dm,point,&parent,NULL);
2288:     if (parent != point) {
2289:       PetscInt closureSize, *closure = NULL, i;

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

2295:         PetscHSetIAdd(ht, cpoint);
2296:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
2297:       }
2298:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
2299:     }
2300:   }
2301:   if (down) {
2302:     PetscInt numChildren;
2303:     const PetscInt *children;

2305:     DMPlexGetTreeChildren(dm,point,&numChildren,&children);
2306:     if (numChildren) {
2307:       PetscInt i;

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

2312:         PetscHSetIAdd(ht, cpoint);
2313:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
2314:       }
2315:     }
2316:   }
2317:   return(0);
2318: }

2320: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
2321: {
2322:   PetscInt       parent;

2326:   DMPlexGetTreeParent(dm, point, &parent,NULL);
2327:   if (point != parent) {
2328:     const PetscInt *cone;
2329:     PetscInt       coneSize, c;

2331:     DMPlexAddClosureTree_Up_Private(dm, ht, parent);
2332:     DMPlexAddClosure_Private(dm, ht, parent);
2333:     DMPlexGetCone(dm, parent, &cone);
2334:     DMPlexGetConeSize(dm, parent, &coneSize);
2335:     for (c = 0; c < coneSize; c++) {
2336:       const PetscInt cp = cone[c];

2338:       DMPlexAddClosureTree_Up_Private(dm, ht, cp);
2339:     }
2340:   }
2341:   return(0);
2342: }

2344: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
2345: {
2346:   PetscInt       i, numChildren;
2347:   const PetscInt *children;

2351:   DMPlexGetTreeChildren(dm, point, &numChildren, &children);
2352:   for (i = 0; i < numChildren; i++) {
2353:     PetscHSetIAdd(ht, children[i]);
2354:   }
2355:   return(0);
2356: }

2358: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
2359: {
2360:   const PetscInt *cone;
2361:   PetscInt       coneSize, c;

2365:   PetscHSetIAdd(ht, point);
2366:   DMPlexAddClosureTree_Up_Private(dm, ht, point);
2367:   DMPlexAddClosureTree_Down_Private(dm, ht, point);
2368:   DMPlexGetCone(dm, point, &cone);
2369:   DMPlexGetConeSize(dm, point, &coneSize);
2370:   for (c = 0; c < coneSize; c++) {
2371:     DMPlexAddClosureTree_Private(dm, ht, cone[c]);
2372:   }
2373:   return(0);
2374: }

2376: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
2377: {
2378:   DM_Plex         *mesh = (DM_Plex *)dm->data;
2379:   const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
2380:   PetscInt        nelems, *elems, off = 0, p;
2381:   PetscHSetI      ht;
2382:   PetscErrorCode  ierr;

2385:   PetscHSetICreate(&ht);
2386:   PetscHSetIResize(ht, numPoints*16);
2387:   if (!hasTree) {
2388:     for (p = 0; p < numPoints; ++p) {
2389:       DMPlexAddClosure_Private(dm, ht, points[p]);
2390:     }
2391:   } else {
2392: #if 1
2393:     for (p = 0; p < numPoints; ++p) {
2394:       DMPlexAddClosureTree_Private(dm, ht, points[p]);
2395:     }
2396: #else
2397:     PetscInt  *closure = NULL, closureSize, c;
2398:     for (p = 0; p < numPoints; ++p) {
2399:       DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2400:       for (c = 0; c < closureSize*2; c += 2) {
2401:         PetscHSetIAdd(ht, closure[c]);
2402:         if (hasTree) {DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);}
2403:       }
2404:     }
2405:     if (closure) {DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);}
2406: #endif
2407:   }
2408:   PetscHSetIGetSize(ht, &nelems);
2409:   PetscMalloc1(nelems, &elems);
2410:   PetscHSetIGetElems(ht, &off, elems);
2411:   PetscHSetIDestroy(&ht);
2412:   PetscSortInt(nelems, elems);
2413:   ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
2414:   return(0);
2415: }

2417: /*@
2418:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

2420:   Input Parameters:
2421: + dm     - The DM
2422: - label  - DMLabel assinging ranks to remote roots

2424:   Level: developer

2426: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
2427: @*/
2428: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
2429: {
2430:   IS              rankIS,   pointIS, closureIS;
2431:   const PetscInt *ranks,   *points;
2432:   PetscInt        numRanks, numPoints, r;
2433:   PetscErrorCode  ierr;

2436:   DMLabelGetValueIS(label, &rankIS);
2437:   ISGetLocalSize(rankIS, &numRanks);
2438:   ISGetIndices(rankIS, &ranks);
2439:   for (r = 0; r < numRanks; ++r) {
2440:     const PetscInt rank = ranks[r];
2441:     DMLabelGetStratumIS(label, rank, &pointIS);
2442:     ISGetLocalSize(pointIS, &numPoints);
2443:     ISGetIndices(pointIS, &points);
2444:     DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
2445:     ISRestoreIndices(pointIS, &points);
2446:     ISDestroy(&pointIS);
2447:     DMLabelSetStratumIS(label, rank, closureIS);
2448:     ISDestroy(&closureIS);
2449:   }
2450:   ISRestoreIndices(rankIS, &ranks);
2451:   ISDestroy(&rankIS);
2452:   return(0);
2453: }

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

2458:   Input Parameters:
2459: + dm     - The DM
2460: - label  - DMLabel assinging ranks to remote roots

2462:   Level: developer

2464: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2465: @*/
2466: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
2467: {
2468:   IS              rankIS,   pointIS;
2469:   const PetscInt *ranks,   *points;
2470:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
2471:   PetscInt       *adj = NULL;
2472:   PetscErrorCode  ierr;

2475:   DMLabelGetValueIS(label, &rankIS);
2476:   ISGetLocalSize(rankIS, &numRanks);
2477:   ISGetIndices(rankIS, &ranks);
2478:   for (r = 0; r < numRanks; ++r) {
2479:     const PetscInt rank = ranks[r];

2481:     DMLabelGetStratumIS(label, rank, &pointIS);
2482:     ISGetLocalSize(pointIS, &numPoints);
2483:     ISGetIndices(pointIS, &points);
2484:     for (p = 0; p < numPoints; ++p) {
2485:       adjSize = PETSC_DETERMINE;
2486:       DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
2487:       for (a = 0; a < adjSize; ++a) {DMLabelSetValue(label, adj[a], rank);}
2488:     }
2489:     ISRestoreIndices(pointIS, &points);
2490:     ISDestroy(&pointIS);
2491:   }
2492:   ISRestoreIndices(rankIS, &ranks);
2493:   ISDestroy(&rankIS);
2494:   PetscFree(adj);
2495:   return(0);
2496: }

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

2501:   Input Parameters:
2502: + dm     - The DM
2503: - label  - DMLabel assinging ranks to remote roots

2505:   Level: developer

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

2510: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2511: @*/
2512: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
2513: {
2514:   MPI_Comm        comm;
2515:   PetscMPIInt     rank;
2516:   PetscSF         sfPoint;
2517:   DMLabel         lblRoots, lblLeaves;
2518:   IS              rankIS, pointIS;
2519:   const PetscInt *ranks;
2520:   PetscInt        numRanks, r;
2521:   PetscErrorCode  ierr;

2524:   PetscObjectGetComm((PetscObject) dm, &comm);
2525:   MPI_Comm_rank(comm, &rank);
2526:   DMGetPointSF(dm, &sfPoint);
2527:   /* Pull point contributions from remote leaves into local roots */
2528:   DMLabelGather(label, sfPoint, &lblLeaves);
2529:   DMLabelGetValueIS(lblLeaves, &rankIS);
2530:   ISGetLocalSize(rankIS, &numRanks);
2531:   ISGetIndices(rankIS, &ranks);
2532:   for (r = 0; r < numRanks; ++r) {
2533:     const PetscInt remoteRank = ranks[r];
2534:     if (remoteRank == rank) continue;
2535:     DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
2536:     DMLabelInsertIS(label, pointIS, remoteRank);
2537:     ISDestroy(&pointIS);
2538:   }
2539:   ISRestoreIndices(rankIS, &ranks);
2540:   ISDestroy(&rankIS);
2541:   DMLabelDestroy(&lblLeaves);
2542:   /* Push point contributions from roots into remote leaves */
2543:   DMLabelDistribute(label, sfPoint, &lblRoots);
2544:   DMLabelGetValueIS(lblRoots, &rankIS);
2545:   ISGetLocalSize(rankIS, &numRanks);
2546:   ISGetIndices(rankIS, &ranks);
2547:   for (r = 0; r < numRanks; ++r) {
2548:     const PetscInt remoteRank = ranks[r];
2549:     if (remoteRank == rank) continue;
2550:     DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
2551:     DMLabelInsertIS(label, pointIS, remoteRank);
2552:     ISDestroy(&pointIS);
2553:   }
2554:   ISRestoreIndices(rankIS, &ranks);
2555:   ISDestroy(&rankIS);
2556:   DMLabelDestroy(&lblRoots);
2557:   return(0);
2558: }

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

2563:   Input Parameters:
2564: + dm        - The DM
2565: . rootLabel - DMLabel assinging ranks to local roots
2566: - processSF - A star forest mapping into the local index on each remote rank

2568:   Output Parameter:
2569: . leafLabel - DMLabel assinging ranks to remote roots

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

2574:   Level: developer

2576: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2577: @*/
2578: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
2579: {
2580:   MPI_Comm           comm;
2581:   PetscMPIInt        rank, size, r;
2582:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
2583:   PetscSF            sfPoint;
2584:   PetscSection       rootSection;
2585:   PetscSFNode       *rootPoints, *leafPoints;
2586:   const PetscSFNode *remote;
2587:   const PetscInt    *local, *neighbors;
2588:   IS                 valueIS;
2589:   PetscBool          mpiOverflow = PETSC_FALSE;
2590:   PetscErrorCode     ierr;

2593:   PetscLogEventBegin(DMPLEX_PartLabelInvert,dm,0,0,0);
2594:   PetscObjectGetComm((PetscObject) dm, &comm);
2595:   MPI_Comm_rank(comm, &rank);
2596:   MPI_Comm_size(comm, &size);
2597:   DMGetPointSF(dm, &sfPoint);

2599:   /* Convert to (point, rank) and use actual owners */
2600:   PetscSectionCreate(comm, &rootSection);
2601:   PetscSectionSetChart(rootSection, 0, size);
2602:   DMLabelGetValueIS(rootLabel, &valueIS);
2603:   ISGetLocalSize(valueIS, &numNeighbors);
2604:   ISGetIndices(valueIS, &neighbors);
2605:   for (n = 0; n < numNeighbors; ++n) {
2606:     DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
2607:     PetscSectionAddDof(rootSection, neighbors[n], numPoints);
2608:   }
2609:   PetscSectionSetUp(rootSection);
2610:   PetscSectionGetStorageSize(rootSection, &rootSize);
2611:   PetscMalloc1(rootSize, &rootPoints);
2612:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
2613:   for (n = 0; n < numNeighbors; ++n) {
2614:     IS              pointIS;
2615:     const PetscInt *points;

2617:     PetscSectionGetOffset(rootSection, neighbors[n], &off);
2618:     DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
2619:     ISGetLocalSize(pointIS, &numPoints);
2620:     ISGetIndices(pointIS, &points);
2621:     for (p = 0; p < numPoints; ++p) {
2622:       if (local) {PetscFindInt(points[p], nleaves, local, &l);}
2623:       else       {l = -1;}
2624:       if (l >= 0) {rootPoints[off+p] = remote[l];}
2625:       else        {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
2626:     }
2627:     ISRestoreIndices(pointIS, &points);
2628:     ISDestroy(&pointIS);
2629:   }

2631:   /* Try to communicate overlap using All-to-All */
2632:   if (!processSF) {
2633:     PetscInt64  counter = 0;
2634:     PetscBool   locOverflow = PETSC_FALSE;
2635:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

2637:     PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
2638:     for (n = 0; n < numNeighbors; ++n) {
2639:       PetscSectionGetDof(rootSection, neighbors[n], &dof);
2640:       PetscSectionGetOffset(rootSection, neighbors[n], &off);
2641: #if defined(PETSC_USE_64BIT_INDICES)
2642:       if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2643:       if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2644: #endif
2645:       scounts[neighbors[n]] = (PetscMPIInt) dof;
2646:       sdispls[neighbors[n]] = (PetscMPIInt) off;
2647:     }
2648:     MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
2649:     for (r = 0; r < size; ++r) { rdispls[r] = (int)counter; counter += rcounts[r]; }
2650:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
2651:     MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
2652:     if (!mpiOverflow) {
2653:       PetscInfo(dm,"Using Alltoallv for mesh distribution\n");
2654:       leafSize = (PetscInt) counter;
2655:       PetscMalloc1(leafSize, &leafPoints);
2656:       MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
2657:     }
2658:     PetscFree4(scounts, sdispls, rcounts, rdispls);
2659:   }

2661:   /* Communicate overlap using process star forest */
2662:   if (processSF || mpiOverflow) {
2663:     PetscSF      procSF;
2664:     PetscSection leafSection;

2666:     if (processSF) {
2667:       PetscInfo(dm,"Using processSF for mesh distribution\n");
2668:       PetscObjectReference((PetscObject)processSF);
2669:       procSF = processSF;
2670:     } else {
2671:       PetscInfo(dm,"Using processSF for mesh distribution (MPI overflow)\n");
2672:       PetscSFCreate(comm,&procSF);
2673:       PetscSFSetGraphWithPattern(procSF,NULL,PETSCSF_PATTERN_ALLTOALL);
2674:     }

2676:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
2677:     DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
2678:     PetscSectionGetStorageSize(leafSection, &leafSize);
2679:     PetscSectionDestroy(&leafSection);
2680:     PetscSFDestroy(&procSF);
2681:   }

2683:   for (p = 0; p < leafSize; p++) {
2684:     DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
2685:   }

2687:   ISRestoreIndices(valueIS, &neighbors);
2688:   ISDestroy(&valueIS);
2689:   PetscSectionDestroy(&rootSection);
2690:   PetscFree(rootPoints);
2691:   PetscFree(leafPoints);
2692:   PetscLogEventEnd(DMPLEX_PartLabelInvert,dm,0,0,0);
2693:   return(0);
2694: }

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

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

2703:   Output Parameter:
2704: . sf    - The star forest communication context encapsulating the defined mapping

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

2708:   Level: developer

2710: .seealso: DMPlexDistribute(), DMPlexCreateOverlap()
2711: @*/
2712: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
2713: {
2714:   PetscMPIInt     rank;
2715:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
2716:   PetscSFNode    *remotePoints;
2717:   IS              remoteRootIS, neighborsIS;
2718:   const PetscInt *remoteRoots, *neighbors;

2722:   PetscLogEventBegin(DMPLEX_PartLabelCreateSF,dm,0,0,0);
2723:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);

2725:   DMLabelGetValueIS(label, &neighborsIS);
2726: #if 0
2727:   {
2728:     IS is;
2729:     ISDuplicate(neighborsIS, &is);
2730:     ISSort(is);
2731:     ISDestroy(&neighborsIS);
2732:     neighborsIS = is;
2733:   }
2734: #endif
2735:   ISGetLocalSize(neighborsIS, &nNeighbors);
2736:   ISGetIndices(neighborsIS, &neighbors);
2737:   for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
2738:     DMLabelGetStratumSize(label, neighbors[n], &numPoints);
2739:     numRemote += numPoints;
2740:   }
2741:   PetscMalloc1(numRemote, &remotePoints);
2742:   /* Put owned points first */
2743:   DMLabelGetStratumSize(label, rank, &numPoints);
2744:   if (numPoints > 0) {
2745:     DMLabelGetStratumIS(label, rank, &remoteRootIS);
2746:     ISGetIndices(remoteRootIS, &remoteRoots);
2747:     for (p = 0; p < numPoints; p++) {
2748:       remotePoints[idx].index = remoteRoots[p];
2749:       remotePoints[idx].rank = rank;
2750:       idx++;
2751:     }
2752:     ISRestoreIndices(remoteRootIS, &remoteRoots);
2753:     ISDestroy(&remoteRootIS);
2754:   }
2755:   /* Now add remote points */
2756:   for (n = 0; n < nNeighbors; ++n) {
2757:     const PetscInt nn = neighbors[n];

2759:     DMLabelGetStratumSize(label, nn, &numPoints);
2760:     if (nn == rank || numPoints <= 0) continue;
2761:     DMLabelGetStratumIS(label, nn, &remoteRootIS);
2762:     ISGetIndices(remoteRootIS, &remoteRoots);
2763:     for (p = 0; p < numPoints; p++) {
2764:       remotePoints[idx].index = remoteRoots[p];
2765:       remotePoints[idx].rank = nn;
2766:       idx++;
2767:     }
2768:     ISRestoreIndices(remoteRootIS, &remoteRoots);
2769:     ISDestroy(&remoteRootIS);
2770:   }
2771:   PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
2772:   DMPlexGetChart(dm, &pStart, &pEnd);
2773:   PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
2774:   PetscSFSetUp(*sf);
2775:   ISDestroy(&neighborsIS);
2776:   PetscLogEventEnd(DMPLEX_PartLabelCreateSF,dm,0,0,0);
2777:   return(0);
2778: }

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

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

2788:   Input parameters:
2789: + dm                - The DMPlex object.
2790: . n                 - The number of points.
2791: . pointsToRewrite   - The points in the SF whose ownership will change.
2792: . targetOwners      - New owner for each element in pointsToRewrite.
2793: - degrees           - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.

2795:   Level: developer

2797: @*/
2798: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
2799: {
2800:   PetscInt      ierr, pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
2801:   PetscInt     *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
2802:   PetscSFNode  *leafLocationsNew;
2803:   const         PetscSFNode *iremote;
2804:   const         PetscInt *ilocal;
2805:   PetscBool    *isLeaf;
2806:   PetscSF       sf;
2807:   MPI_Comm      comm;
2808:   PetscMPIInt   rank, size;

2811:   PetscObjectGetComm((PetscObject) dm, &comm);
2812:   MPI_Comm_rank(comm, &rank);
2813:   MPI_Comm_size(comm, &size);
2814:   DMPlexGetChart(dm, &pStart, &pEnd);

2816:   DMGetPointSF(dm, &sf);
2817:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
2818:   PetscMalloc1(pEnd-pStart, &isLeaf);
2819:   for (i=0; i<pEnd-pStart; i++) {
2820:     isLeaf[i] = PETSC_FALSE;
2821:   }
2822:   for (i=0; i<nleafs; i++) {
2823:     isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
2824:   }

2826:   PetscMalloc1(pEnd-pStart+1, &cumSumDegrees);
2827:   cumSumDegrees[0] = 0;
2828:   for (i=1; i<=pEnd-pStart; i++) {
2829:     cumSumDegrees[i] = cumSumDegrees[i-1] + degrees[i-1];
2830:   }
2831:   sumDegrees = cumSumDegrees[pEnd-pStart];
2832:   /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */

2834:   PetscMalloc1(sumDegrees, &locationsOfLeafs);
2835:   PetscMalloc1(pEnd-pStart, &rankOnLeafs);
2836:   for (i=0; i<pEnd-pStart; i++) {
2837:     rankOnLeafs[i] = rank;
2838:   }
2839:   PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
2840:   PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
2841:   PetscFree(rankOnLeafs);

2843:   /* get the remote local points of my leaves */
2844:   PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
2845:   PetscMalloc1(pEnd-pStart, &points);
2846:   for (i=0; i<pEnd-pStart; i++) {
2847:     points[i] = pStart+i;
2848:   }
2849:   PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
2850:   PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
2851:   PetscFree(points);
2852:   /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
2853:   PetscMalloc1(pEnd-pStart, &newOwners);
2854:   PetscMalloc1(pEnd-pStart, &newNumbers);
2855:   for (i=0; i<pEnd-pStart; i++) {
2856:     newOwners[i] = -1;
2857:     newNumbers[i] = -1;
2858:   }
2859:   {
2860:     PetscInt oldNumber, newNumber, oldOwner, newOwner;
2861:     for (i=0; i<n; i++) {
2862:       oldNumber = pointsToRewrite[i];
2863:       newNumber = -1;
2864:       oldOwner = rank;
2865:       newOwner = targetOwners[i];
2866:       if (oldOwner == newOwner) {
2867:         newNumber = oldNumber;
2868:       } else {
2869:         for (j=0; j<degrees[oldNumber]; j++) {
2870:           if (locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) {
2871:             newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j];
2872:             break;
2873:           }
2874:         }
2875:       }
2876:       if (newNumber == -1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");

2878:       newOwners[oldNumber] = newOwner;
2879:       newNumbers[oldNumber] = newNumber;
2880:     }
2881:   }
2882:   PetscFree(cumSumDegrees);
2883:   PetscFree(locationsOfLeafs);
2884:   PetscFree(remoteLocalPointOfLeafs);

2886:   PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners);
2887:   PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners);
2888:   PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers);
2889:   PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers);

2891:   /* Now count how many leafs we have on each processor. */
2892:   leafCounter=0;
2893:   for (i=0; i<pEnd-pStart; i++) {
2894:     if (newOwners[i] >= 0) {
2895:       if (newOwners[i] != rank) {
2896:         leafCounter++;
2897:       }
2898:     } else {
2899:       if (isLeaf[i]) {
2900:         leafCounter++;
2901:       }
2902:     }
2903:   }

2905:   /* Now set up the new sf by creating the leaf arrays */
2906:   PetscMalloc1(leafCounter, &leafsNew);
2907:   PetscMalloc1(leafCounter, &leafLocationsNew);

2909:   leafCounter = 0;
2910:   counter = 0;
2911:   for (i=0; i<pEnd-pStart; i++) {
2912:     if (newOwners[i] >= 0) {
2913:       if (newOwners[i] != rank) {
2914:         leafsNew[leafCounter] = i;
2915:         leafLocationsNew[leafCounter].rank = newOwners[i];
2916:         leafLocationsNew[leafCounter].index = newNumbers[i];
2917:         leafCounter++;
2918:       }
2919:     } else {
2920:       if (isLeaf[i]) {
2921:         leafsNew[leafCounter] = i;
2922:         leafLocationsNew[leafCounter].rank = iremote[counter].rank;
2923:         leafLocationsNew[leafCounter].index = iremote[counter].index;
2924:         leafCounter++;
2925:       }
2926:     }
2927:     if (isLeaf[i]) {
2928:       counter++;
2929:     }
2930:   }

2932:   PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
2933:   PetscFree(newOwners);
2934:   PetscFree(newNumbers);
2935:   PetscFree(isLeaf);
2936:   return(0);
2937: }

2939: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
2940: {
2941:   PetscInt *distribution, min, max, sum, i, ierr;
2942:   PetscMPIInt rank, size;
2944:   MPI_Comm_size(comm, &size);
2945:   MPI_Comm_rank(comm, &rank);
2946:   PetscCalloc1(size, &distribution);
2947:   for (i=0; i<n; i++) {
2948:     if (part) distribution[part[i]] += vtxwgt[skip*i];
2949:     else distribution[rank] += vtxwgt[skip*i];
2950:   }
2951:   MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
2952:   min = distribution[0];
2953:   max = distribution[0];
2954:   sum = distribution[0];
2955:   for (i=1; i<size; i++) {
2956:     if (distribution[i]<min) min=distribution[i];
2957:     if (distribution[i]>max) max=distribution[i];
2958:     sum += distribution[i];
2959:   }
2960:   PetscViewerASCIIPrintf(viewer, "Min: %D, Avg: %D, Max: %D, Balance: %f\n", min, sum/size, max, (max*1.*size)/sum);
2961:   PetscFree(distribution);
2962:   return(0);
2963: }

2965: #endif

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

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

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

2979:   Level: intermediate

2981: @*/

2983: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
2984: {
2985: #if defined(PETSC_HAVE_PARMETIS)
2986:   PetscSF     sf;
2987:   PetscInt    ierr, i, j, idx, jdx;
2988:   PetscInt    eBegin, eEnd, nroots, nleafs, pStart, pEnd;
2989:   const       PetscInt *degrees, *ilocal;
2990:   const       PetscSFNode *iremote;
2991:   PetscBool   *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
2992:   PetscInt    numExclusivelyOwned, numNonExclusivelyOwned;
2993:   PetscMPIInt rank, size;
2994:   PetscInt    *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
2995:   const       PetscInt *cumSumVertices;
2996:   PetscInt    offset, counter;
2997:   PetscInt    lenadjncy;
2998:   PetscInt    *xadj, *adjncy, *vtxwgt;
2999:   PetscInt    lenxadj;
3000:   PetscInt    *adjwgt = NULL;
3001:   PetscInt    *part, *options;
3002:   PetscInt    nparts, wgtflag, numflag, ncon, edgecut;
3003:   real_t      *ubvec;
3004:   PetscInt    *firstVertices, *renumbering;
3005:   PetscInt    failed, failedGlobal;
3006:   MPI_Comm    comm;
3007:   Mat         A, Apre;
3008:   const char *prefix = NULL;
3009:   PetscViewer       viewer;
3010:   PetscViewerFormat format;
3011:   PetscLayout layout;

3014:   if (success) *success = PETSC_FALSE;
3015:   PetscObjectGetComm((PetscObject) dm, &comm);
3016:   MPI_Comm_rank(comm, &rank);
3017:   MPI_Comm_size(comm, &size);
3018:   if (size==1) return(0);

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

3022:   PetscOptionsGetViewer(comm,((PetscObject)dm)->options, prefix,"-dm_rebalance_partition_view",&viewer,&format,NULL);
3023:   if (viewer) {
3024:     PetscViewerPushFormat(viewer,format);
3025:   }

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

3032:   for (i=0; i<pEnd-pStart; i++) {
3033:     toBalance[i] = (PetscBool)(i-pStart>=eBegin && i-pStart<eEnd);
3034:   }

3036:   /* There are three types of points:
3037:    * exclusivelyOwned: points that are owned by this process and only seen by this process
3038:    * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
3039:    * leaf: a point that is seen by this process but owned by a different process
3040:    */
3041:   DMGetPointSF(dm, &sf);
3042:   PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
3043:   PetscMalloc1(pEnd-pStart, &isLeaf);
3044:   PetscMalloc1(pEnd-pStart, &isNonExclusivelyOwned);
3045:   PetscMalloc1(pEnd-pStart, &isExclusivelyOwned);
3046:   for (i=0; i<pEnd-pStart; i++) {
3047:     isNonExclusivelyOwned[i] = PETSC_FALSE;
3048:     isExclusivelyOwned[i] = PETSC_FALSE;
3049:     isLeaf[i] = PETSC_FALSE;
3050:   }

3052:   /* start by marking all the leafs */
3053:   for (i=0; i<nleafs; i++) {
3054:     isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
3055:   }

3057:   /* for an owned point, we can figure out whether another processor sees it or
3058:    * not by calculating its degree */
3059:   PetscSFComputeDegreeBegin(sf, &degrees);
3060:   PetscSFComputeDegreeEnd(sf, &degrees);

3062:   numExclusivelyOwned = 0;
3063:   numNonExclusivelyOwned = 0;
3064:   for (i=0; i<pEnd-pStart; i++) {
3065:     if (toBalance[i]) {
3066:       if (degrees[i] > 0) {
3067:         isNonExclusivelyOwned[i] = PETSC_TRUE;
3068:         numNonExclusivelyOwned += 1;
3069:       } else {
3070:         if (!isLeaf[i]) {
3071:           isExclusivelyOwned[i] = PETSC_TRUE;
3072:           numExclusivelyOwned += 1;
3073:         }
3074:       }
3075:     }
3076:   }

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

3082:   PetscLayoutCreate(comm, &layout);
3083:   PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
3084:   PetscLayoutSetUp(layout);
3085:   PetscLayoutGetRanges(layout, &cumSumVertices);

3087:   PetscMalloc1(pEnd-pStart, &globalNumbersOfLocalOwnedVertices);
3088:   for (i=0; i<pEnd-pStart; i++) {globalNumbersOfLocalOwnedVertices[i] = pStart - 1;}
3089:   offset = cumSumVertices[rank];
3090:   counter = 0;
3091:   for (i=0; i<pEnd-pStart; i++) {
3092:     if (toBalance[i]) {
3093:       if (degrees[i] > 0) {
3094:         globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
3095:         counter++;
3096:       }
3097:     }
3098:   }

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

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

3107:   MatCreate(comm, &Apre);
3108:   MatSetType(Apre, MATPREALLOCATOR);
3109:   MatSetSizes(Apre, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
3110:   MatSetUp(Apre);

3112:   for (i=0; i<pEnd-pStart; i++) {
3113:     if (toBalance[i]) {
3114:       idx = cumSumVertices[rank];
3115:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
3116:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
3117:       else continue;
3118:       MatSetValue(Apre, idx, jdx, 1, INSERT_VALUES);
3119:       MatSetValue(Apre, jdx, idx, 1, INSERT_VALUES);
3120:     }
3121:   }

3123:   MatAssemblyBegin(Apre, MAT_FINAL_ASSEMBLY);
3124:   MatAssemblyEnd(Apre, MAT_FINAL_ASSEMBLY);

3126:   MatCreate(comm, &A);
3127:   MatSetType(A, MATMPIAIJ);
3128:   MatSetSizes(A, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
3129:   MatPreallocatorPreallocate(Apre, PETSC_FALSE, A);
3130:   MatDestroy(&Apre);

3132:   for (i=0; i<pEnd-pStart; i++) {
3133:     if (toBalance[i]) {
3134:       idx = cumSumVertices[rank];
3135:       if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
3136:       else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
3137:       else continue;
3138:       MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
3139:       MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
3140:     }
3141:   }

3143:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
3144:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
3145:   PetscFree(leafGlobalNumbers);
3146:   PetscFree(globalNumbersOfLocalOwnedVertices);

3148:   nparts = size;
3149:   wgtflag = 2;
3150:   numflag = 0;
3151:   ncon = 2;
3152:   real_t *tpwgts;
3153:   PetscMalloc1(ncon * nparts, &tpwgts);
3154:   for (i=0; i<ncon*nparts; i++) {
3155:     tpwgts[i] = 1./(nparts);
3156:   }

3158:   PetscMalloc1(ncon, &ubvec);
3159:   ubvec[0] = 1.01;
3160:   ubvec[1] = 1.01;
3161:   lenadjncy = 0;
3162:   for (i=0; i<1+numNonExclusivelyOwned; i++) {
3163:     PetscInt temp=0;
3164:     MatGetRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
3165:     lenadjncy += temp;
3166:     MatRestoreRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
3167:   }
3168:   PetscMalloc1(lenadjncy, &adjncy);
3169:   lenxadj = 2 + numNonExclusivelyOwned;
3170:   PetscMalloc1(lenxadj, &xadj);
3171:   xadj[0] = 0;
3172:   counter = 0;
3173:   for (i=0; i<1+numNonExclusivelyOwned; i++) {
3174:     PetscInt        temp=0;
3175:     const PetscInt *cols;
3176:     MatGetRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
3177:     PetscArraycpy(&adjncy[counter], cols, temp);
3178:     counter += temp;
3179:     xadj[i+1] = counter;
3180:     MatRestoreRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
3181:   }

3183:   PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);
3184:   PetscMalloc1(ncon*(1 + numNonExclusivelyOwned), &vtxwgt);
3185:   vtxwgt[0] = numExclusivelyOwned;
3186:   if (ncon>1) vtxwgt[1] = 1;
3187:   for (i=0; i<numNonExclusivelyOwned; i++) {
3188:     vtxwgt[ncon*(i+1)] = 1;
3189:     if (ncon>1) vtxwgt[ncon*(i+1)+1] = 0;
3190:   }

3192:   if (viewer) {
3193:     PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %D on interface of mesh distribution.\n", entityDepth);
3194:     PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %D\n", cumSumVertices[size]);
3195:   }
3196:   if (parallel) {
3197:     PetscMalloc1(4, &options);
3198:     options[0] = 1;
3199:     options[1] = 0; /* Verbosity */
3200:     options[2] = 0; /* Seed */
3201:     options[3] = PARMETIS_PSR_COUPLED; /* Seed */
3202:     if (viewer) { PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"); }
3203:     if (useInitialGuess) {
3204:       if (viewer) { PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"); }
3205:       PetscStackPush("ParMETIS_V3_RefineKway");
3206:       ParMETIS_V3_RefineKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
3207:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_RefineKway()");
3208:       PetscStackPop;
3209:     } else {
3210:       PetscStackPush("ParMETIS_V3_PartKway");
3211:       ParMETIS_V3_PartKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
3212:       PetscStackPop;
3213:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
3214:     }
3215:     PetscFree(options);
3216:   } else {
3217:     if (viewer) { PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"); }
3218:     Mat As;
3219:     PetscInt numRows;
3220:     PetscInt *partGlobal;
3221:     MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &As);

3223:     PetscInt *numExclusivelyOwnedAll;
3224:     PetscMalloc1(size, &numExclusivelyOwnedAll);
3225:     numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
3226:     MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,numExclusivelyOwnedAll,1,MPIU_INT,comm);

3228:     MatGetSize(As, &numRows, NULL);
3229:     PetscMalloc1(numRows, &partGlobal);
3230:     if (!rank) {
3231:       PetscInt *adjncy_g, *xadj_g, *vtxwgt_g;
3232:       lenadjncy = 0;

3234:       for (i=0; i<numRows; i++) {
3235:         PetscInt temp=0;
3236:         MatGetRow(As, i, &temp, NULL, NULL);
3237:         lenadjncy += temp;
3238:         MatRestoreRow(As, i, &temp, NULL, NULL);
3239:       }
3240:       PetscMalloc1(lenadjncy, &adjncy_g);
3241:       lenxadj = 1 + numRows;
3242:       PetscMalloc1(lenxadj, &xadj_g);
3243:       xadj_g[0] = 0;
3244:       counter = 0;
3245:       for (i=0; i<numRows; i++) {
3246:         PetscInt        temp=0;
3247:         const PetscInt *cols;
3248:         MatGetRow(As, i, &temp, &cols, NULL);
3249:         PetscArraycpy(&adjncy_g[counter], cols, temp);
3250:         counter += temp;
3251:         xadj_g[i+1] = counter;
3252:         MatRestoreRow(As, i, &temp, &cols, NULL);
3253:       }
3254:       PetscMalloc1(2*numRows, &vtxwgt_g);
3255:       for (i=0; i<size; i++){
3256:         vtxwgt_g[ncon*cumSumVertices[i]] = numExclusivelyOwnedAll[i];
3257:         if (ncon>1) vtxwgt_g[ncon*cumSumVertices[i]+1] = 1;
3258:         for (j=cumSumVertices[i]+1; j<cumSumVertices[i+1]; j++) {
3259:           vtxwgt_g[ncon*j] = 1;
3260:           if (ncon>1) vtxwgt_g[2*j+1] = 0;
3261:         }
3262:       }
3263:       PetscMalloc1(64, &options);
3264:       METIS_SetDefaultOptions(options); /* initialize all defaults */
3265:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
3266:       options[METIS_OPTION_CONTIG] = 1;
3267:       PetscStackPush("METIS_PartGraphKway");
3268:       METIS_PartGraphKway(&numRows, &ncon, xadj_g, adjncy_g, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
3269:       PetscStackPop;
3270:       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
3271:       PetscFree(options);
3272:       PetscFree(xadj_g);
3273:       PetscFree(adjncy_g);
3274:       PetscFree(vtxwgt_g);
3275:     }
3276:     PetscFree(numExclusivelyOwnedAll);

3278:     /* Now scatter the parts array. */
3279:     {
3280:       PetscMPIInt *counts, *mpiCumSumVertices;
3281:       PetscMalloc1(size, &counts);
3282:       PetscMalloc1(size+1, &mpiCumSumVertices);
3283:       for(i=0; i<size; i++) {
3284:         PetscMPIIntCast(cumSumVertices[i+1] - cumSumVertices[i], &(counts[i]));
3285:       }
3286:       for(i=0; i<=size; i++) {
3287:         PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
3288:       }
3289:       MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
3290:       PetscFree(counts);
3291:       PetscFree(mpiCumSumVertices);
3292:     }

3294:     PetscFree(partGlobal);
3295:     MatDestroy(&As);
3296:   }

3298:   MatDestroy(&A);
3299:   PetscFree(ubvec);
3300:   PetscFree(tpwgts);

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

3304:   PetscMalloc1(size, &firstVertices);
3305:   PetscMalloc1(size, &renumbering);
3306:   firstVertices[rank] = part[0];
3307:   MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);
3308:   for (i=0; i<size; i++) {
3309:     renumbering[firstVertices[i]] = i;
3310:   }
3311:   for (i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) {
3312:     part[i] = renumbering[part[i]];
3313:   }
3314:   /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
3315:   failed = (PetscInt)(part[0] != rank);
3316:   MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);

3318:   PetscFree(firstVertices);
3319:   PetscFree(renumbering);

3321:   if (failedGlobal > 0) {
3322:     PetscLayoutDestroy(&layout);
3323:     PetscFree(xadj);
3324:     PetscFree(adjncy);
3325:     PetscFree(vtxwgt);
3326:     PetscFree(toBalance);
3327:     PetscFree(isLeaf);
3328:     PetscFree(isNonExclusivelyOwned);
3329:     PetscFree(isExclusivelyOwned);
3330:     PetscFree(part);
3331:     if (viewer) {
3332:       PetscViewerPopFormat(viewer);
3333:       PetscViewerDestroy(&viewer);
3334:     }
3335:     PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3336:     return(0);
3337:   }

3339:   /*Let's check how well we did distributing points*/
3340:   if (viewer) {
3341:     PetscViewerASCIIPrintf(viewer, "Comparing number of owned entities of depth %D on each process before rebalancing, after rebalancing, and after consistency checks.\n", entityDepth);
3342:     PetscViewerASCIIPrintf(viewer, "Initial.     ");
3343:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
3344:     PetscViewerASCIIPrintf(viewer, "Rebalanced.  ");
3345:     DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
3346:   }

3348:   /* Now check that every vertex is owned by a process that it is actually connected to. */
3349:   for (i=1; i<=numNonExclusivelyOwned; i++) {
3350:     PetscInt loc = 0;
3351:     PetscFindInt(cumSumVertices[part[i]], xadj[i+1]-xadj[i], &adjncy[xadj[i]], &loc);
3352:     /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
3353:     if (loc<0) {
3354:       part[i] = rank;
3355:     }
3356:   }

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

3364:   PetscLayoutDestroy(&layout);
3365:   PetscFree(xadj);
3366:   PetscFree(adjncy);
3367:   PetscFree(vtxwgt);

3369:   /* Almost done, now rewrite the SF to reflect the new ownership. */
3370:   {
3371:     PetscInt *pointsToRewrite;
3372:     PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
3373:     counter = 0;
3374:     for(i=0; i<pEnd-pStart; i++) {
3375:       if (toBalance[i]) {
3376:         if (isNonExclusivelyOwned[i]) {
3377:           pointsToRewrite[counter] = i + pStart;
3378:           counter++;
3379:         }
3380:       }
3381:     }
3382:     DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part+1, degrees);
3383:     PetscFree(pointsToRewrite);
3384:   }

3386:   PetscFree(toBalance);
3387:   PetscFree(isLeaf);
3388:   PetscFree(isNonExclusivelyOwned);
3389:   PetscFree(isExclusivelyOwned);
3390:   PetscFree(part);
3391:   if (viewer) {
3392:     PetscViewerPopFormat(viewer);
3393:     PetscViewerDestroy(&viewer);
3394:   }
3395:   if (success) *success = PETSC_TRUE;
3396:   PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3397: #else
3398:   SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
3399: #endif
3400:   return(0);
3401: }