Actual source code: plexpartition.c

petsc-3.11.1 2019-04-17
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       = {http://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: /*@C
 31:   DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner

 33:   Input Parameters:
 34: + dm      - The mesh DM dm
 35: - height  - Height of the strata from which to construct the graph

 37:   Output Parameter:
 38: + numVertices     - Number of vertices in the graph
 39: . offsets         - Point offsets in the graph
 40: . adjacency       - Point connectivity in the graph
 41: - 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.

 43:   The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function
 44:   representation on the mesh.

 46:   Level: developer

 48: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
 49: @*/
 50: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
 51: {
 52:   PetscInt       dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
 53:   PetscInt      *adj = NULL, *vOffsets = NULL, *graph = NULL;
 54:   IS             cellNumbering;
 55:   const PetscInt *cellNum;
 56:   PetscBool      useCone, useClosure;
 57:   PetscSection   section;
 58:   PetscSegBuffer adjBuffer;
 59:   PetscSF        sfPoint;
 60:   PetscInt       *adjCells = NULL, *remoteCells = NULL;
 61:   const PetscInt *local;
 62:   PetscInt       nroots, nleaves, l;
 63:   PetscMPIInt    rank;

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

107:     PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
108:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
109:     for (l = 0; l < nroots; ++l) adjCells[l] = -3;
110:     for (f = fStart; f < fEnd; ++f) {
111:       const PetscInt *support;
112:       PetscInt        supportSize;

114:       DMPlexGetSupport(dm, f, &support);
115:       DMPlexGetSupportSize(dm, f, &supportSize);
116:       if (supportSize == 1) adjCells[f] = cellNum[support[0]];
117:       else if (supportSize == 2) {
118:         PetscFindInt(support[0], nleaves, local, &p);
119:         if (p >= 0) adjCells[f] = cellNum[support[1]];
120:         PetscFindInt(support[1], nleaves, local, &p);
121:         if (p >= 0) adjCells[f] = cellNum[support[0]];
122:       }
123:     }
124:     for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
125:     PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells);
126:     PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells);
127:     PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
128:     PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
129:   }
130:   /* Combine local and global adjacencies */
131:   for (*numVertices = 0, p = pStart; p < pEnd; p++) {
132:     const PetscInt *cone;
133:     PetscInt        coneSize, c;

135:     /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
136:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
137:     /* Add remote cells */
138:     if (remoteCells) {
139:       DMPlexGetCone(dm, p, &cone);
140:       DMPlexGetConeSize(dm, p, &coneSize);
141:       for (c = 0; c < coneSize; ++c) {
142:         if (remoteCells[cone[c]] != -1) {
143:           PetscInt *PETSC_RESTRICT pBuf;

145:           PetscSectionAddDof(section, p, 1);
146:           PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
147:           *pBuf = remoteCells[cone[c]];
148:         }
149:       }
150:     }
151:     /* Add local cells */
152:     adjSize = PETSC_DETERMINE;
153:     DMPlexGetAdjacency(dm, p, &adjSize, &adj);
154:     for (a = 0; a < adjSize; ++a) {
155:       const PetscInt point = adj[a];
156:       if (point != p && pStart <= point && point < pEnd) {
157:         PetscInt *PETSC_RESTRICT pBuf;
158:         PetscSectionAddDof(section, p, 1);
159:         PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
160:         *pBuf = cellNum[point];
161:       }
162:     }
163:     (*numVertices)++;
164:   }
165:   PetscFree2(adjCells, remoteCells);
166:   DMSetBasicAdjacency(dm, useCone, useClosure);
167:   /* Derive CSR graph from section/segbuffer */
168:   PetscSectionSetUp(section);
169:   PetscSectionGetStorageSize(section, &size);
170:   PetscMalloc1(*numVertices+1, &vOffsets);
171:   for (idx = 0, p = pStart; p < pEnd; p++) {
172:     if (nroots > 0) {if (cellNum[p] < 0) continue;}
173:     PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
174:   }
175:   vOffsets[*numVertices] = size;
176:   if (offsets) *offsets = vOffsets;
177:   PetscSegBufferExtractAlloc(adjBuffer, &graph);
178:   ISRestoreIndices(cellNumbering, &cellNum);
179:   ISDestroy(&cellNumbering);
180:   if (adjacency) *adjacency = graph;
181:   /* Clean up */
182:   PetscSegBufferDestroy(&adjBuffer);
183:   PetscSectionDestroy(&section);
184:   PetscFree(adj);
185:   return(0);
186: }

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

191:   Collective

193:   Input Arguments:
194: + dm - The DMPlex
195: - cellHeight - The height of mesh points to treat as cells (default should be 0)

197:   Output Arguments:
198: + numVertices - The number of local vertices in the graph, or cells in the mesh.
199: . offsets     - The offset to the adjacency list for each cell
200: - adjacency   - The adjacency list for all cells

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

204:   Level: advanced

206: .seealso: DMPlexCreate()
207: @*/
208: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
209: {
210:   const PetscInt maxFaceCases = 30;
211:   PetscInt       numFaceCases = 0;
212:   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
213:   PetscInt      *off, *adj;
214:   PetscInt      *neighborCells = NULL;
215:   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;

219:   /* For parallel partitioning, I think you have to communicate supports */
220:   DMGetDimension(dm, &dim);
221:   cellDim = dim - cellHeight;
222:   DMPlexGetDepth(dm, &depth);
223:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
224:   if (cEnd - cStart == 0) {
225:     if (numVertices) *numVertices = 0;
226:     if (offsets)   *offsets   = NULL;
227:     if (adjacency) *adjacency = NULL;
228:     return(0);
229:   }
230:   numCells  = cEnd - cStart;
231:   faceDepth = depth - cellHeight;
232:   if (dim == depth) {
233:     PetscInt f, fStart, fEnd;

235:     PetscCalloc1(numCells+1, &off);
236:     /* Count neighboring cells */
237:     DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
238:     for (f = fStart; f < fEnd; ++f) {
239:       const PetscInt *support;
240:       PetscInt        supportSize;
241:       DMPlexGetSupportSize(dm, f, &supportSize);
242:       DMPlexGetSupport(dm, f, &support);
243:       if (supportSize == 2) {
244:         PetscInt numChildren;

246:         DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
247:         if (!numChildren) {
248:           ++off[support[0]-cStart+1];
249:           ++off[support[1]-cStart+1];
250:         }
251:       }
252:     }
253:     /* Prefix sum */
254:     for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
255:     if (adjacency) {
256:       PetscInt *tmp;

258:       PetscMalloc1(off[numCells], &adj);
259:       PetscMalloc1(numCells+1, &tmp);
260:       PetscMemcpy(tmp, off, (numCells+1) * sizeof(PetscInt));
261:       /* Get neighboring cells */
262:       for (f = fStart; f < fEnd; ++f) {
263:         const PetscInt *support;
264:         PetscInt        supportSize;
265:         DMPlexGetSupportSize(dm, f, &supportSize);
266:         DMPlexGetSupport(dm, f, &support);
267:         if (supportSize == 2) {
268:           PetscInt numChildren;

270:           DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
271:           if (!numChildren) {
272:             adj[tmp[support[0]-cStart]++] = support[1];
273:             adj[tmp[support[1]-cStart]++] = support[0];
274:           }
275:         }
276:       }
277: #if defined(PETSC_USE_DEBUG)
278:       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);
279: #endif
280:       PetscFree(tmp);
281:     }
282:     if (numVertices) *numVertices = numCells;
283:     if (offsets)   *offsets   = off;
284:     if (adjacency) *adjacency = adj;
285:     return(0);
286:   }
287:   /* Setup face recognition */
288:   if (faceDepth == 1) {
289:     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 */

291:     for (c = cStart; c < cEnd; ++c) {
292:       PetscInt corners;

294:       DMPlexGetConeSize(dm, c, &corners);
295:       if (!cornersSeen[corners]) {
296:         PetscInt nFV;

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

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

303:         numFaceVertices[numFaceCases++] = nFV;
304:       }
305:     }
306:   }
307:   PetscCalloc1(numCells+1, &off);
308:   /* Count neighboring cells */
309:   for (cell = cStart; cell < cEnd; ++cell) {
310:     PetscInt numNeighbors = PETSC_DETERMINE, n;

312:     DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
313:     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
314:     for (n = 0; n < numNeighbors; ++n) {
315:       PetscInt        cellPair[2];
316:       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
317:       PetscInt        meetSize = 0;
318:       const PetscInt *meet    = NULL;

320:       cellPair[0] = cell; cellPair[1] = neighborCells[n];
321:       if (cellPair[0] == cellPair[1]) continue;
322:       if (!found) {
323:         DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
324:         if (meetSize) {
325:           PetscInt f;

327:           for (f = 0; f < numFaceCases; ++f) {
328:             if (numFaceVertices[f] == meetSize) {
329:               found = PETSC_TRUE;
330:               break;
331:             }
332:           }
333:         }
334:         DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
335:       }
336:       if (found) ++off[cell-cStart+1];
337:     }
338:   }
339:   /* Prefix sum */
340:   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];

342:   if (adjacency) {
343:     PetscMalloc1(off[numCells], &adj);
344:     /* Get neighboring cells */
345:     for (cell = cStart; cell < cEnd; ++cell) {
346:       PetscInt numNeighbors = PETSC_DETERMINE, n;
347:       PetscInt cellOffset   = 0;

349:       DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
350:       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
351:       for (n = 0; n < numNeighbors; ++n) {
352:         PetscInt        cellPair[2];
353:         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
354:         PetscInt        meetSize = 0;
355:         const PetscInt *meet    = NULL;

357:         cellPair[0] = cell; cellPair[1] = neighborCells[n];
358:         if (cellPair[0] == cellPair[1]) continue;
359:         if (!found) {
360:           DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
361:           if (meetSize) {
362:             PetscInt f;

364:             for (f = 0; f < numFaceCases; ++f) {
365:               if (numFaceVertices[f] == meetSize) {
366:                 found = PETSC_TRUE;
367:                 break;
368:               }
369:             }
370:           }
371:           DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
372:         }
373:         if (found) {
374:           adj[off[cell-cStart]+cellOffset] = neighborCells[n];
375:           ++cellOffset;
376:         }
377:       }
378:     }
379:   }
380:   PetscFree(neighborCells);
381:   if (numVertices) *numVertices = numCells;
382:   if (offsets)   *offsets   = off;
383:   if (adjacency) *adjacency = adj;
384:   return(0);
385: }

387: /*@C
388:   PetscPartitionerRegister - Adds a new PetscPartitioner implementation

390:   Not Collective

392:   Input Parameters:
393: + name        - The name of a new user-defined creation routine
394: - create_func - The creation routine itself

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

399:   Sample usage:
400: .vb
401:     PetscPartitionerRegister("my_part", MyPetscPartitionerCreate);
402: .ve

404:   Then, your PetscPartitioner type can be chosen with the procedural interface via
405: .vb
406:     PetscPartitionerCreate(MPI_Comm, PetscPartitioner *);
407:     PetscPartitionerSetType(PetscPartitioner, "my_part");
408: .ve
409:    or at runtime via the option
410: .vb
411:     -petscpartitioner_type my_part
412: .ve

414:   Level: advanced

416: .keywords: PetscPartitioner, register
417: .seealso: PetscPartitionerRegisterAll(), PetscPartitionerRegisterDestroy()

419: @*/
420: PetscErrorCode PetscPartitionerRegister(const char sname[], PetscErrorCode (*function)(PetscPartitioner))
421: {

425:   PetscFunctionListAdd(&PetscPartitionerList, sname, function);
426:   return(0);
427: }

429: /*@C
430:   PetscPartitionerSetType - Builds a particular PetscPartitioner

432:   Collective on PetscPartitioner

434:   Input Parameters:
435: + part - The PetscPartitioner object
436: - name - The kind of partitioner

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

441:   Level: intermediate

443: .keywords: PetscPartitioner, set, type
444: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate()
445: @*/
446: PetscErrorCode PetscPartitionerSetType(PetscPartitioner part, PetscPartitionerType name)
447: {
448:   PetscErrorCode (*r)(PetscPartitioner);
449:   PetscBool      match;

454:   PetscObjectTypeCompare((PetscObject) part, name, &match);
455:   if (match) return(0);

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

461:   if (part->ops->destroy) {
462:     (*part->ops->destroy)(part);
463:   }
464:   PetscMemzero(part->ops, sizeof(struct _PetscPartitionerOps));
465:   (*r)(part);
466:   PetscObjectChangeTypeName((PetscObject) part, name);
467:   return(0);
468: }

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

473:   Not Collective

475:   Input Parameter:
476: . part - The PetscPartitioner

478:   Output Parameter:
479: . name - The PetscPartitioner type name

481:   Level: intermediate

483: .keywords: PetscPartitioner, get, type, name
484: .seealso: PetscPartitionerSetType(), PetscPartitionerCreate()
485: @*/
486: PetscErrorCode PetscPartitionerGetType(PetscPartitioner part, PetscPartitionerType *name)
487: {

493:   PetscPartitionerRegisterAll();
494:   *name = ((PetscObject) part)->type_name;
495:   return(0);
496: }

498: /*@C
499:   PetscPartitionerView - Views a PetscPartitioner

501:   Collective on PetscPartitioner

503:   Input Parameter:
504: + part - the PetscPartitioner object to view
505: - v    - the viewer

507:   Level: developer

509: .seealso: PetscPartitionerDestroy()
510: @*/
511: PetscErrorCode PetscPartitionerView(PetscPartitioner part, PetscViewer v)
512: {
513:   PetscMPIInt    size;
514:   PetscBool      isascii;

519:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) part), &v);}
520:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
521:   if (isascii) {
522:     MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
523:     PetscViewerASCIIPrintf(v, "Graph Partitioner: %d MPI Process%s\n", size, size > 1 ? "es" : "");
524:     PetscViewerASCIIPrintf(v, "  type: %s\n", part->hdr.type_name);
525:     PetscViewerASCIIPushTab(v);
526:     PetscViewerASCIIPrintf(v, "edge cut: %D\n", part->edgeCut);
527:     PetscViewerASCIIPrintf(v, "balance:  %.2g\n", part->balance);
528:     PetscViewerASCIIPopTab(v);
529:   }
530:   if (part->ops->view) {(*part->ops->view)(part, v);}
531:   return(0);
532: }

534: static PetscErrorCode PetscPartitionerGetDefaultType(const char *currentType, const char **defaultType)
535: {
537:   if (!currentType) {
538: #if defined(PETSC_HAVE_CHACO)
539:     *defaultType = PETSCPARTITIONERCHACO;
540: #elif defined(PETSC_HAVE_PARMETIS)
541:     *defaultType = PETSCPARTITIONERPARMETIS;
542: #elif defined(PETSC_HAVE_PTSCOTCH)
543:     *defaultType = PETSCPARTITIONERPTSCOTCH;
544: #else
545:     *defaultType = PETSCPARTITIONERSIMPLE;
546: #endif
547:   } else {
548:     *defaultType = currentType;
549:   }
550:   return(0);
551: }

553: /*@
554:   PetscPartitionerSetFromOptions - sets parameters in a PetscPartitioner from the options database

556:   Collective on PetscPartitioner

558:   Input Parameter:
559: . part - the PetscPartitioner object to set options for

561:   Level: developer

563: .seealso: PetscPartitionerView()
564: @*/
565: PetscErrorCode PetscPartitionerSetFromOptions(PetscPartitioner part)
566: {
567:   const char    *defaultType;
568:   char           name[256];
569:   PetscBool      flg;

574:   PetscPartitionerRegisterAll();
575:   PetscPartitionerGetDefaultType(((PetscObject) part)->type_name,&defaultType);
576:   PetscObjectOptionsBegin((PetscObject) part);
577:   PetscOptionsFList("-petscpartitioner_type", "Graph partitioner", "PetscPartitionerSetType", PetscPartitionerList, defaultType, name, sizeof(name), &flg);
578:   if (flg) {
579:     PetscPartitionerSetType(part, name);
580:   } else if (!((PetscObject) part)->type_name) {
581:     PetscPartitionerSetType(part, defaultType);
582:   }
583:   if (part->ops->setfromoptions) {
584:     (*part->ops->setfromoptions)(PetscOptionsObject,part);
585:   }
586:   PetscOptionsGetViewer(((PetscObject) part)->comm, ((PetscObject) part)->options, ((PetscObject) part)->prefix, "-petscpartitioner_view_graph", &part->viewerGraph, &part->formatGraph, &part->viewGraph);
587:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
588:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) part);
589:   PetscOptionsEnd();
590:   return(0);
591: }

593: /*@C
594:   PetscPartitionerSetUp - Construct data structures for the PetscPartitioner

596:   Collective on PetscPartitioner

598:   Input Parameter:
599: . part - the PetscPartitioner object to setup

601:   Level: developer

603: .seealso: PetscPartitionerView(), PetscPartitionerDestroy()
604: @*/
605: PetscErrorCode PetscPartitionerSetUp(PetscPartitioner part)
606: {

611:   if (part->ops->setup) {(*part->ops->setup)(part);}
612:   return(0);
613: }

615: /*@
616:   PetscPartitionerDestroy - Destroys a PetscPartitioner object

618:   Collective on PetscPartitioner

620:   Input Parameter:
621: . part - the PetscPartitioner object to destroy

623:   Level: developer

625: .seealso: PetscPartitionerView()
626: @*/
627: PetscErrorCode PetscPartitionerDestroy(PetscPartitioner *part)
628: {

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

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

638:   PetscViewerDestroy(&(*part)->viewerGraph);
639:   if ((*part)->ops->destroy) {(*(*part)->ops->destroy)(*part);}
640:   PetscHeaderDestroy(part);
641:   return(0);
642: }

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

647:   Collective on MPI_Comm

649:   Input Parameter:
650: . comm - The communicator for the PetscPartitioner object

652:   Output Parameter:
653: . part - The PetscPartitioner object

655:   Level: beginner

657: .seealso: PetscPartitionerSetType(), PETSCPARTITIONERCHACO, PETSCPARTITIONERPARMETIS, PETSCPARTITIONERSHELL, PETSCPARTITIONERSIMPLE, PETSCPARTITIONERGATHER
658: @*/
659: PetscErrorCode PetscPartitionerCreate(MPI_Comm comm, PetscPartitioner *part)
660: {
661:   PetscPartitioner p;
662:   const char       *partitionerType = NULL;
663:   PetscErrorCode   ierr;

667:   *part = NULL;
668:   DMInitializePackage();

670:   PetscHeaderCreate(p, PETSCPARTITIONER_CLASSID, "PetscPartitioner", "Graph Partitioner", "PetscPartitioner", comm, PetscPartitionerDestroy, PetscPartitionerView);
671:   PetscPartitionerGetDefaultType(NULL,&partitionerType);
672:   PetscPartitionerSetType(p,partitionerType);

674:   p->edgeCut = 0;
675:   p->balance = 0.0;

677:   *part = p;
678:   return(0);
679: }

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

684:   Collective on DM

686:   Input Parameters:
687: + part    - The PetscPartitioner
688: - dm      - The mesh DM

690:   Output Parameters:
691: + partSection     - The PetscSection giving the division of points by partition
692: - partition       - The list of points by partition

694:   Options Database:
695: . -petscpartitioner_view_graph - View the graph we are partitioning

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

699:   Level: developer

701: .seealso DMPlexDistribute(), PetscPartitionerSetPointHeight(), PetscPartitionerCreate()
702: @*/
703: PetscErrorCode PetscPartitionerPartition(PetscPartitioner part, DM dm, PetscSection partSection, IS *partition)
704: {
705:   PetscMPIInt    size;

713:   MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
714:   if (size == 1) {
715:     PetscInt *points;
716:     PetscInt  cStart, cEnd, c;

718:     DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
719:     PetscSectionSetChart(partSection, 0, size);
720:     PetscSectionSetDof(partSection, 0, cEnd-cStart);
721:     PetscSectionSetUp(partSection);
722:     PetscMalloc1(cEnd-cStart, &points);
723:     for (c = cStart; c < cEnd; ++c) points[c] = c;
724:     ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
725:     return(0);
726:   }
727:   if (part->height == 0) {
728:     PetscInt numVertices;
729:     PetscInt *start     = NULL;
730:     PetscInt *adjacency = NULL;
731:     IS       globalNumbering;

733:     DMPlexCreatePartitionerGraph(dm, 0, &numVertices, &start, &adjacency, &globalNumbering);
734:     if (part->viewGraph) {
735:       PetscViewer viewer = part->viewerGraph;
736:       PetscBool   isascii;
737:       PetscInt    v, i;
738:       PetscMPIInt rank;

740:       MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
741:       PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
742:       if (isascii) {
743:         PetscViewerASCIIPushSynchronized(viewer);
744:         PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Nv: %D\n", rank, numVertices);
745:         for (v = 0; v < numVertices; ++v) {
746:           const PetscInt s = start[v];
747:           const PetscInt e = start[v+1];

749:           PetscViewerASCIISynchronizedPrintf(viewer, "[%d]  ", rank);
750:           for (i = s; i < e; ++i) {PetscViewerASCIISynchronizedPrintf(viewer, "%D ", adjacency[i]);}
751:           PetscViewerASCIISynchronizedPrintf(viewer, "[%D-%D)\n", s, e);
752:         }
753:         PetscViewerFlush(viewer);
754:         PetscViewerASCIIPopSynchronized(viewer);
755:       }
756:     }
757:     if (!part->ops->partition) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_WRONGSTATE, "PetscPartitioner has no type");
758:     (*part->ops->partition)(part, dm, size, numVertices, start, adjacency, partSection, partition);
759:     PetscFree(start);
760:     PetscFree(adjacency);
761:     if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
762:       const PetscInt *globalNum;
763:       const PetscInt *partIdx;
764:       PetscInt *map, cStart, cEnd;
765:       PetscInt *adjusted, i, localSize, offset;
766:       IS    newPartition;

768:       ISGetLocalSize(*partition,&localSize);
769:       PetscMalloc1(localSize,&adjusted);
770:       ISGetIndices(globalNumbering,&globalNum);
771:       ISGetIndices(*partition,&partIdx);
772:       PetscMalloc1(localSize,&map);
773:       DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
774:       for (i = cStart, offset = 0; i < cEnd; i++) {
775:         if (globalNum[i - cStart] >= 0) map[offset++] = i;
776:       }
777:       for (i = 0; i < localSize; i++) {
778:         adjusted[i] = map[partIdx[i]];
779:       }
780:       PetscFree(map);
781:       ISRestoreIndices(*partition,&partIdx);
782:       ISRestoreIndices(globalNumbering,&globalNum);
783:       ISCreateGeneral(PETSC_COMM_SELF,localSize,adjusted,PETSC_OWN_POINTER,&newPartition);
784:       ISDestroy(&globalNumbering);
785:       ISDestroy(partition);
786:       *partition = newPartition;
787:     }
788:   } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
789:   PetscPartitionerViewFromOptions(part, NULL, "-petscpartitioner_view");
790:   return(0);
791: }

793: static PetscErrorCode PetscPartitionerDestroy_Shell(PetscPartitioner part)
794: {
795:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
796:   PetscErrorCode          ierr;

799:   PetscSectionDestroy(&p->section);
800:   ISDestroy(&p->partition);
801:   PetscFree(p);
802:   return(0);
803: }

805: static PetscErrorCode PetscPartitionerView_Shell_Ascii(PetscPartitioner part, PetscViewer viewer)
806: {
807:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
808:   PetscErrorCode          ierr;

811:   if (p->random) {
812:     PetscViewerASCIIPushTab(viewer);
813:     PetscViewerASCIIPrintf(viewer, "using random partition\n");
814:     PetscViewerASCIIPopTab(viewer);
815:   }
816:   return(0);
817: }

819: static PetscErrorCode PetscPartitionerView_Shell(PetscPartitioner part, PetscViewer viewer)
820: {
821:   PetscBool      iascii;

827:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
828:   if (iascii) {PetscPartitionerView_Shell_Ascii(part, viewer);}
829:   return(0);
830: }

832: static PetscErrorCode PetscPartitionerSetFromOptions_Shell(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
833: {
834:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
835:   PetscErrorCode          ierr;

838:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner Shell Options");
839:   PetscOptionsBool("-petscpartitioner_shell_random", "Use a random partition", "PetscPartitionerView", PETSC_FALSE, &p->random, NULL);
840:   PetscOptionsTail();
841:   return(0);
842: }

844: static PetscErrorCode PetscPartitionerPartition_Shell(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
845: {
846:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
847:   PetscInt                np;
848:   PetscErrorCode          ierr;

851:   if (p->random) {
852:     PetscRandom r;
853:     PetscInt   *sizes, *points, v, p;
854:     PetscMPIInt rank;

856:     MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
857:     PetscRandomCreate(PETSC_COMM_SELF, &r);
858:     PetscRandomSetInterval(r, 0.0, (PetscScalar) nparts);
859:     PetscRandomSetFromOptions(r);
860:     PetscCalloc2(nparts, &sizes, numVertices, &points);
861:     for (v = 0; v < numVertices; ++v) {points[v] = v;}
862:     for (p = 0; p < nparts; ++p) {sizes[p] = numVertices/nparts + (PetscInt) (p < numVertices % nparts);}
863:     for (v = numVertices-1; v > 0; --v) {
864:       PetscReal val;
865:       PetscInt  w, tmp;

867:       PetscRandomSetInterval(r, 0.0, (PetscScalar) (v+1));
868:       PetscRandomGetValueReal(r, &val);
869:       w    = PetscFloorReal(val);
870:       tmp       = points[v];
871:       points[v] = points[w];
872:       points[w] = tmp;
873:     }
874:     PetscRandomDestroy(&r);
875:     PetscPartitionerShellSetPartition(part, nparts, sizes, points);
876:     PetscFree2(sizes, points);
877:   }
878:   if (!p->section) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Shell partitioner information not provided. Please call PetscPartitionerShellSetPartition()");
879:   PetscSectionGetChart(p->section, NULL, &np);
880:   if (nparts != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of requested partitions %d != configured partitions %d", nparts, np);
881:   ISGetLocalSize(p->partition, &np);
882:   if (numVertices != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input vertices %d != configured vertices %d", numVertices, np);
883:   PetscSectionCopy(p->section, partSection);
884:   *partition = p->partition;
885:   PetscObjectReference((PetscObject) p->partition);
886:   return(0);
887: }

889: static PetscErrorCode PetscPartitionerInitialize_Shell(PetscPartitioner part)
890: {
892:   part->ops->view           = PetscPartitionerView_Shell;
893:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_Shell;
894:   part->ops->destroy        = PetscPartitionerDestroy_Shell;
895:   part->ops->partition      = PetscPartitionerPartition_Shell;
896:   return(0);
897: }

899: /*MC
900:   PETSCPARTITIONERSHELL = "shell" - A PetscPartitioner object

902:   Level: intermediate

904: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
905: M*/

907: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Shell(PetscPartitioner part)
908: {
909:   PetscPartitioner_Shell *p;
910:   PetscErrorCode          ierr;

914:   PetscNewLog(part, &p);
915:   part->data = p;

917:   PetscPartitionerInitialize_Shell(part);
918:   p->random = PETSC_FALSE;
919:   return(0);
920: }

922: /*@C
923:   PetscPartitionerShellSetPartition - Set an artifical partition for a mesh

925:   Collective on Part

927:   Input Parameters:
928: + part     - The PetscPartitioner
929: . size - The number of partitions
930: . sizes    - array of size size (or NULL) providing the number of points in each partition
931: - 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.)

933:   Level: developer

935:   Notes:

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

939: .seealso DMPlexDistribute(), PetscPartitionerCreate()
940: @*/
941: PetscErrorCode PetscPartitionerShellSetPartition(PetscPartitioner part, PetscInt size, const PetscInt sizes[], const PetscInt points[])
942: {
943:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
944:   PetscInt                proc, numPoints;
945:   PetscErrorCode          ierr;

951:   PetscSectionDestroy(&p->section);
952:   ISDestroy(&p->partition);
953:   PetscSectionCreate(PetscObjectComm((PetscObject) part), &p->section);
954:   PetscSectionSetChart(p->section, 0, size);
955:   if (sizes) {
956:     for (proc = 0; proc < size; ++proc) {
957:       PetscSectionSetDof(p->section, proc, sizes[proc]);
958:     }
959:   }
960:   PetscSectionSetUp(p->section);
961:   PetscSectionGetStorageSize(p->section, &numPoints);
962:   ISCreateGeneral(PetscObjectComm((PetscObject) part), numPoints, points, PETSC_COPY_VALUES, &p->partition);
963:   return(0);
964: }

966: /*@
967:   PetscPartitionerShellSetRandom - Set the flag to use a random partition

969:   Collective on Part

971:   Input Parameters:
972: + part   - The PetscPartitioner
973: - random - The flag to use a random partition

975:   Level: intermediate

977: .seealso PetscPartitionerShellGetRandom(), PetscPartitionerCreate()
978: @*/
979: PetscErrorCode PetscPartitionerShellSetRandom(PetscPartitioner part, PetscBool random)
980: {
981:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;

985:   p->random = random;
986:   return(0);
987: }

989: /*@
990:   PetscPartitionerShellGetRandom - get the flag to use a random partition

992:   Collective on Part

994:   Input Parameter:
995: . part   - The PetscPartitioner

997:   Output Parameter
998: . random - The flag to use a random partition

1000:   Level: intermediate

1002: .seealso PetscPartitionerShellSetRandom(), PetscPartitionerCreate()
1003: @*/
1004: PetscErrorCode PetscPartitionerShellGetRandom(PetscPartitioner part, PetscBool *random)
1005: {
1006:   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;

1011:   *random = p->random;
1012:   return(0);
1013: }

1015: static PetscErrorCode PetscPartitionerDestroy_Simple(PetscPartitioner part)
1016: {
1017:   PetscPartitioner_Simple *p = (PetscPartitioner_Simple *) part->data;
1018:   PetscErrorCode          ierr;

1021:   PetscFree(p);
1022:   return(0);
1023: }

1025: static PetscErrorCode PetscPartitionerView_Simple_Ascii(PetscPartitioner part, PetscViewer viewer)
1026: {
1028:   return(0);
1029: }

1031: static PetscErrorCode PetscPartitionerView_Simple(PetscPartitioner part, PetscViewer viewer)
1032: {
1033:   PetscBool      iascii;

1039:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1040:   if (iascii) {PetscPartitionerView_Simple_Ascii(part, viewer);}
1041:   return(0);
1042: }

1044: static PetscErrorCode PetscPartitionerPartition_Simple(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1045: {
1046:   MPI_Comm       comm;
1047:   PetscInt       np;
1048:   PetscMPIInt    size;

1052:   comm = PetscObjectComm((PetscObject)dm);
1053:   MPI_Comm_size(comm,&size);
1054:   PetscSectionSetChart(partSection, 0, nparts);
1055:   ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1056:   if (size == 1) {
1057:     for (np = 0; np < nparts; ++np) {PetscSectionSetDof(partSection, np, numVertices/nparts + ((numVertices % nparts) > np));}
1058:   }
1059:   else {
1060:     PetscMPIInt rank;
1061:     PetscInt nvGlobal, *offsets, myFirst, myLast;

1063:     PetscMalloc1(size+1,&offsets);
1064:     offsets[0] = 0;
1065:     MPI_Allgather(&numVertices,1,MPIU_INT,&offsets[1],1,MPIU_INT,comm);
1066:     for (np = 2; np <= size; np++) {
1067:       offsets[np] += offsets[np-1];
1068:     }
1069:     nvGlobal = offsets[size];
1070:     MPI_Comm_rank(comm,&rank);
1071:     myFirst = offsets[rank];
1072:     myLast  = offsets[rank + 1] - 1;
1073:     PetscFree(offsets);
1074:     if (numVertices) {
1075:       PetscInt firstPart = 0, firstLargePart = 0;
1076:       PetscInt lastPart = 0, lastLargePart = 0;
1077:       PetscInt rem = nvGlobal % nparts;
1078:       PetscInt pSmall = nvGlobal/nparts;
1079:       PetscInt pBig = nvGlobal/nparts + 1;


1082:       if (rem) {
1083:         firstLargePart = myFirst / pBig;
1084:         lastLargePart  = myLast  / pBig;

1086:         if (firstLargePart < rem) {
1087:           firstPart = firstLargePart;
1088:         }
1089:         else {
1090:           firstPart = rem + (myFirst - (rem * pBig)) / pSmall;
1091:         }
1092:         if (lastLargePart < rem) {
1093:           lastPart = lastLargePart;
1094:         }
1095:         else {
1096:           lastPart = rem + (myLast - (rem * pBig)) / pSmall;
1097:         }
1098:       }
1099:       else {
1100:         firstPart = myFirst / (nvGlobal/nparts);
1101:         lastPart  = myLast  / (nvGlobal/nparts);
1102:       }

1104:       for (np = firstPart; np <= lastPart; np++) {
1105:         PetscInt PartStart =  np    * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np);
1106:         PetscInt PartEnd   = (np+1) * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np+1);

1108:         PartStart = PetscMax(PartStart,myFirst);
1109:         PartEnd   = PetscMin(PartEnd,myLast+1);
1110:         PetscSectionSetDof(partSection,np,PartEnd-PartStart);
1111:       }
1112:     }
1113:   }
1114:   PetscSectionSetUp(partSection);
1115:   return(0);
1116: }

1118: static PetscErrorCode PetscPartitionerInitialize_Simple(PetscPartitioner part)
1119: {
1121:   part->ops->view      = PetscPartitionerView_Simple;
1122:   part->ops->destroy   = PetscPartitionerDestroy_Simple;
1123:   part->ops->partition = PetscPartitionerPartition_Simple;
1124:   return(0);
1125: }

1127: /*MC
1128:   PETSCPARTITIONERSIMPLE = "simple" - A PetscPartitioner object

1130:   Level: intermediate

1132: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1133: M*/

1135: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Simple(PetscPartitioner part)
1136: {
1137:   PetscPartitioner_Simple *p;
1138:   PetscErrorCode           ierr;

1142:   PetscNewLog(part, &p);
1143:   part->data = p;

1145:   PetscPartitionerInitialize_Simple(part);
1146:   return(0);
1147: }

1149: static PetscErrorCode PetscPartitionerDestroy_Gather(PetscPartitioner part)
1150: {
1151:   PetscPartitioner_Gather *p = (PetscPartitioner_Gather *) part->data;
1152:   PetscErrorCode          ierr;

1155:   PetscFree(p);
1156:   return(0);
1157: }

1159: static PetscErrorCode PetscPartitionerView_Gather_Ascii(PetscPartitioner part, PetscViewer viewer)
1160: {
1162:   return(0);
1163: }

1165: static PetscErrorCode PetscPartitionerView_Gather(PetscPartitioner part, PetscViewer viewer)
1166: {
1167:   PetscBool      iascii;

1173:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1174:   if (iascii) {PetscPartitionerView_Gather_Ascii(part, viewer);}
1175:   return(0);
1176: }

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

1184:   PetscSectionSetChart(partSection, 0, nparts);
1185:   ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1186:   PetscSectionSetDof(partSection,0,numVertices);
1187:   for (np = 1; np < nparts; ++np) {PetscSectionSetDof(partSection, np, 0);}
1188:   PetscSectionSetUp(partSection);
1189:   return(0);
1190: }

1192: static PetscErrorCode PetscPartitionerInitialize_Gather(PetscPartitioner part)
1193: {
1195:   part->ops->view      = PetscPartitionerView_Gather;
1196:   part->ops->destroy   = PetscPartitionerDestroy_Gather;
1197:   part->ops->partition = PetscPartitionerPartition_Gather;
1198:   return(0);
1199: }

1201: /*MC
1202:   PETSCPARTITIONERGATHER = "gather" - A PetscPartitioner object

1204:   Level: intermediate

1206: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1207: M*/

1209: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Gather(PetscPartitioner part)
1210: {
1211:   PetscPartitioner_Gather *p;
1212:   PetscErrorCode           ierr;

1216:   PetscNewLog(part, &p);
1217:   part->data = p;

1219:   PetscPartitionerInitialize_Gather(part);
1220:   return(0);
1221: }


1224: static PetscErrorCode PetscPartitionerDestroy_Chaco(PetscPartitioner part)
1225: {
1226:   PetscPartitioner_Chaco *p = (PetscPartitioner_Chaco *) part->data;
1227:   PetscErrorCode          ierr;

1230:   PetscFree(p);
1231:   return(0);
1232: }

1234: static PetscErrorCode PetscPartitionerView_Chaco_Ascii(PetscPartitioner part, PetscViewer viewer)
1235: {
1237:   return(0);
1238: }

1240: static PetscErrorCode PetscPartitionerView_Chaco(PetscPartitioner part, PetscViewer viewer)
1241: {
1242:   PetscBool      iascii;

1248:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1249:   if (iascii) {PetscPartitionerView_Chaco_Ascii(part, viewer);}
1250:   return(0);
1251: }

1253: #if defined(PETSC_HAVE_CHACO)
1254: #if defined(PETSC_HAVE_UNISTD_H)
1255: #include <unistd.h>
1256: #endif
1257: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1258: #include <chaco.h>
1259: #else
1260: /* Older versions of Chaco do not have an include file */
1261: PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
1262:                        float *ewgts, float *x, float *y, float *z, char *outassignname,
1263:                        char *outfilename, short *assignment, int architecture, int ndims_tot,
1264:                        int mesh_dims[3], double *goal, int global_method, int local_method,
1265:                        int rqi_flag, int vmax, int ndims, double eigtol, long seed);
1266: #endif
1267: extern int FREE_GRAPH;
1268: #endif

1270: static PetscErrorCode PetscPartitionerPartition_Chaco(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1271: {
1272: #if defined(PETSC_HAVE_CHACO)
1273:   enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
1274:   MPI_Comm       comm;
1275:   int            nvtxs          = numVertices; /* number of vertices in full graph */
1276:   int           *vwgts          = NULL;   /* weights for all vertices */
1277:   float         *ewgts          = NULL;   /* weights for all edges */
1278:   float         *x              = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
1279:   char          *outassignname  = NULL;   /*  name of assignment output file */
1280:   char          *outfilename    = NULL;   /* output file name */
1281:   int            architecture   = 1;      /* 0 => hypercube, d => d-dimensional mesh */
1282:   int            ndims_tot      = 0;      /* total number of cube dimensions to divide */
1283:   int            mesh_dims[3];            /* dimensions of mesh of processors */
1284:   double        *goal          = NULL;    /* desired set sizes for each set */
1285:   int            global_method = 1;       /* global partitioning algorithm */
1286:   int            local_method  = 1;       /* local partitioning algorithm */
1287:   int            rqi_flag      = 0;       /* should I use RQI/Symmlq eigensolver? */
1288:   int            vmax          = 200;     /* how many vertices to coarsen down to? */
1289:   int            ndims         = 1;       /* number of eigenvectors (2^d sets) */
1290:   double         eigtol        = 0.001;   /* tolerance on eigenvectors */
1291:   long           seed          = 123636512; /* for random graph mutations */
1292: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1293:   int           *assignment;              /* Output partition */
1294: #else
1295:   short int     *assignment;              /* Output partition */
1296: #endif
1297:   int            fd_stdout, fd_pipe[2];
1298:   PetscInt      *points;
1299:   int            i, v, p;

1303:   PetscObjectGetComm((PetscObject)dm,&comm);
1304: #if defined (PETSC_USE_DEBUG)
1305:   {
1306:     int ival,isum;
1307:     PetscBool distributed;

1309:     ival = (numVertices > 0);
1310:     MPI_Allreduce(&ival, &isum, 1, MPI_INT, MPI_SUM, comm);
1311:     distributed = (isum > 1) ? PETSC_TRUE : PETSC_FALSE;
1312:     if (distributed) SETERRQ(comm, PETSC_ERR_SUP, "Chaco cannot partition a distributed graph");
1313:   }
1314: #endif
1315:   if (!numVertices) {
1316:     PetscSectionSetChart(partSection, 0, nparts);
1317:     PetscSectionSetUp(partSection);
1318:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
1319:     return(0);
1320:   }
1321:   FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
1322:   for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];

1324:   if (global_method == INERTIAL_METHOD) {
1325:     /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
1326:     SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
1327:   }
1328:   mesh_dims[0] = nparts;
1329:   mesh_dims[1] = 1;
1330:   mesh_dims[2] = 1;
1331:   PetscMalloc1(nvtxs, &assignment);
1332:   /* Chaco outputs to stdout. We redirect this to a buffer. */
1333:   /* TODO: check error codes for UNIX calls */
1334: #if defined(PETSC_HAVE_UNISTD_H)
1335:   {
1336:     int piperet;
1337:     piperet = pipe(fd_pipe);
1338:     if (piperet) SETERRQ(comm,PETSC_ERR_SYS,"Could not create pipe");
1339:     fd_stdout = dup(1);
1340:     close(1);
1341:     dup2(fd_pipe[1], 1);
1342:   }
1343: #endif
1344:   interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
1345:                    assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
1346:                    vmax, ndims, eigtol, seed);
1347: #if defined(PETSC_HAVE_UNISTD_H)
1348:   {
1349:     char msgLog[10000];
1350:     int  count;

1352:     fflush(stdout);
1353:     count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
1354:     if (count < 0) count = 0;
1355:     msgLog[count] = 0;
1356:     close(1);
1357:     dup2(fd_stdout, 1);
1358:     close(fd_stdout);
1359:     close(fd_pipe[0]);
1360:     close(fd_pipe[1]);
1361:     if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
1362:   }
1363: #else
1364:   if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", "error in stdout");
1365: #endif
1366:   /* Convert to PetscSection+IS */
1367:   PetscSectionSetChart(partSection, 0, nparts);
1368:   for (v = 0; v < nvtxs; ++v) {
1369:     PetscSectionAddDof(partSection, assignment[v], 1);
1370:   }
1371:   PetscSectionSetUp(partSection);
1372:   PetscMalloc1(nvtxs, &points);
1373:   for (p = 0, i = 0; p < nparts; ++p) {
1374:     for (v = 0; v < nvtxs; ++v) {
1375:       if (assignment[v] == p) points[i++] = v;
1376:     }
1377:   }
1378:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1379:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
1380:   if (global_method == INERTIAL_METHOD) {
1381:     /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
1382:   }
1383:   PetscFree(assignment);
1384:   for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
1385:   return(0);
1386: #else
1387:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco.");
1388: #endif
1389: }

1391: static PetscErrorCode PetscPartitionerInitialize_Chaco(PetscPartitioner part)
1392: {
1394:   part->ops->view      = PetscPartitionerView_Chaco;
1395:   part->ops->destroy   = PetscPartitionerDestroy_Chaco;
1396:   part->ops->partition = PetscPartitionerPartition_Chaco;
1397:   return(0);
1398: }

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

1403:   Level: intermediate

1405: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1406: M*/

1408: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Chaco(PetscPartitioner part)
1409: {
1410:   PetscPartitioner_Chaco *p;
1411:   PetscErrorCode          ierr;

1415:   PetscNewLog(part, &p);
1416:   part->data = p;

1418:   PetscPartitionerInitialize_Chaco(part);
1419:   PetscCitationsRegister(ChacoPartitionerCitation, &ChacoPartitionercite);
1420:   return(0);
1421: }

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

1425: static PetscErrorCode PetscPartitionerDestroy_ParMetis(PetscPartitioner part)
1426: {
1427:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1428:   PetscErrorCode             ierr;

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

1435: static PetscErrorCode PetscPartitionerView_ParMetis_Ascii(PetscPartitioner part, PetscViewer viewer)
1436: {
1437:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1438:   PetscErrorCode             ierr;

1441:   PetscViewerASCIIPushTab(viewer);
1442:   PetscViewerASCIIPrintf(viewer, "ParMetis type: %s\n", ptypes[p->ptype]);
1443:   PetscViewerASCIIPrintf(viewer, "load imbalance ratio %g\n", (double) p->imbalanceRatio);
1444:   PetscViewerASCIIPrintf(viewer, "debug flag %D\n", p->debugFlag);
1445:   PetscViewerASCIIPopTab(viewer);
1446:   return(0);
1447: }

1449: static PetscErrorCode PetscPartitionerView_ParMetis(PetscPartitioner part, PetscViewer viewer)
1450: {
1451:   PetscBool      iascii;

1457:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1458:   if (iascii) {PetscPartitionerView_ParMetis_Ascii(part, viewer);}
1459:   return(0);
1460: }

1462: static PetscErrorCode PetscPartitionerSetFromOptions_ParMetis(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1463: {
1464:   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1465:   PetscErrorCode            ierr;

1468:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner ParMetis Options");
1469:   PetscOptionsEList("-petscpartitioner_parmetis_type", "Partitioning method", "", ptypes, 2, ptypes[p->ptype], &p->ptype, NULL);
1470:   PetscOptionsReal("-petscpartitioner_parmetis_imbalance_ratio", "Load imbalance ratio limit", "", p->imbalanceRatio, &p->imbalanceRatio, NULL);
1471:   PetscOptionsInt("-petscpartitioner_parmetis_debug", "Debugging flag", "", p->debugFlag, &p->debugFlag, NULL);
1472:   PetscOptionsTail();
1473:   return(0);
1474: }

1476: #if defined(PETSC_HAVE_PARMETIS)
1477: #include <parmetis.h>
1478: #endif

1480: static PetscErrorCode PetscPartitionerPartition_ParMetis(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1481: {
1482: #if defined(PETSC_HAVE_PARMETIS)
1483:   PetscPartitioner_ParMetis *pm = (PetscPartitioner_ParMetis *) part->data;
1484:   MPI_Comm       comm;
1485:   PetscSection   section;
1486:   PetscInt       nvtxs       = numVertices; /* The number of vertices in full graph */
1487:   PetscInt      *vtxdist;                   /* Distribution of vertices across processes */
1488:   PetscInt      *xadj        = start;       /* Start of edge list for each vertex */
1489:   PetscInt      *adjncy      = adjacency;   /* Edge lists for all vertices */
1490:   PetscInt      *vwgt        = NULL;        /* Vertex weights */
1491:   PetscInt      *adjwgt      = NULL;        /* Edge weights */
1492:   PetscInt       wgtflag     = 0;           /* Indicates which weights are present */
1493:   PetscInt       numflag     = 0;           /* Indicates initial offset (0 or 1) */
1494:   PetscInt       ncon        = 1;           /* The number of weights per vertex */
1495:   PetscInt       metis_ptype = pm->ptype;   /* kway or recursive bisection */
1496:   real_t        *tpwgts;                    /* The fraction of vertex weights assigned to each partition */
1497:   real_t        *ubvec;                     /* The balance intolerance for vertex weights */
1498:   PetscInt       options[64];               /* Options */
1499:   /* Outputs */
1500:   PetscInt       v, i, *assignment, *points;
1501:   PetscMPIInt    size, rank, p;

1505:   PetscObjectGetComm((PetscObject) part, &comm);
1506:   MPI_Comm_size(comm, &size);
1507:   MPI_Comm_rank(comm, &rank);
1508:   /* Calculate vertex distribution */
1509:   PetscMalloc5(size+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment,nvtxs,&vwgt);
1510:   vtxdist[0] = 0;
1511:   MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
1512:   for (p = 2; p <= size; ++p) {
1513:     vtxdist[p] += vtxdist[p-1];
1514:   }
1515:   /* Calculate partition weights */
1516:   for (p = 0; p < nparts; ++p) {
1517:     tpwgts[p] = 1.0/nparts;
1518:   }
1519:   ubvec[0] = pm->imbalanceRatio;
1520:   /* Weight cells by dofs on cell by default */
1521:   DMGetSection(dm, &section);
1522:   if (section) {
1523:     PetscInt cStart, cEnd, dof;

1525:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1526:     for (v = cStart; v < cEnd; ++v) {
1527:       PetscSectionGetDof(section, v, &dof);
1528:       /* WARNING: Assumes that meshes with overlap have the overlapped cells at the end of the stratum. */
1529:       /* To do this properly, we should use the cell numbering created in DMPlexCreatePartitionerGraph. */
1530:       if (v-cStart < numVertices) vwgt[v-cStart] = PetscMax(dof, 1);
1531:     }
1532:   } else {
1533:     for (v = 0; v < nvtxs; ++v) vwgt[v] = 1;
1534:   }
1535:   wgtflag |= 2; /* have weights on graph vertices */

1537:   if (nparts == 1) {
1538:     PetscMemzero(assignment, nvtxs * sizeof(PetscInt));
1539:   } else {
1540:     for (p = 0; !vtxdist[p+1] && p < size; ++p);
1541:     if (vtxdist[p+1] == vtxdist[size]) {
1542:       if (rank == p) {
1543:         METIS_SetDefaultOptions(options); /* initialize all defaults */
1544:         if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
1545:         if (metis_ptype == 1) {
1546:           PetscStackPush("METIS_PartGraphRecursive");
1547:           METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
1548:           PetscStackPop;
1549:           if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphRecursive()");
1550:         } else {
1551:           /*
1552:            It would be nice to activate the two options below, but they would need some actual testing.
1553:            - Turning on these options may exercise path of the METIS code that have bugs and may break production runs.
1554:            - 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.
1555:           */
1556:           /* options[METIS_OPTION_CONTIG]  = 1; */ /* try to produce partitions that are contiguous */
1557:           /* options[METIS_OPTION_MINCONN] = 1; */ /* minimize the maximum degree of the subdomain graph */
1558:           PetscStackPush("METIS_PartGraphKway");
1559:           METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
1560:           PetscStackPop;
1561:           if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
1562:         }
1563:       }
1564:     } else {
1565:       options[0] = 1;
1566:       options[1] = pm->debugFlag;
1567:       PetscStackPush("ParMETIS_V3_PartKway");
1568:       ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment, &comm);
1569:       PetscStackPop;
1570:       if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error %d in ParMETIS_V3_PartKway()", ierr);
1571:     }
1572:   }
1573:   /* Convert to PetscSection+IS */
1574:   PetscSectionSetChart(partSection, 0, nparts);
1575:   for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
1576:   PetscSectionSetUp(partSection);
1577:   PetscMalloc1(nvtxs, &points);
1578:   for (p = 0, i = 0; p < nparts; ++p) {
1579:     for (v = 0; v < nvtxs; ++v) {
1580:       if (assignment[v] == p) points[i++] = v;
1581:     }
1582:   }
1583:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1584:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
1585:   PetscFree5(vtxdist,tpwgts,ubvec,assignment,vwgt);
1586:   return(0);
1587: #else
1588:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
1589: #endif
1590: }

1592: static PetscErrorCode PetscPartitionerInitialize_ParMetis(PetscPartitioner part)
1593: {
1595:   part->ops->view           = PetscPartitionerView_ParMetis;
1596:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_ParMetis;
1597:   part->ops->destroy        = PetscPartitionerDestroy_ParMetis;
1598:   part->ops->partition      = PetscPartitionerPartition_ParMetis;
1599:   return(0);
1600: }

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

1605:   Level: intermediate

1607: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1608: M*/

1610: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMetis(PetscPartitioner part)
1611: {
1612:   PetscPartitioner_ParMetis *p;
1613:   PetscErrorCode          ierr;

1617:   PetscNewLog(part, &p);
1618:   part->data = p;

1620:   p->ptype          = 0;
1621:   p->imbalanceRatio = 1.05;
1622:   p->debugFlag      = 0;

1624:   PetscPartitionerInitialize_ParMetis(part);
1625:   PetscCitationsRegister(ParMetisPartitionerCitation, &ParMetisPartitionercite);
1626:   return(0);
1627: }


1630: PetscBool PTScotchPartitionercite = PETSC_FALSE;
1631: const char PTScotchPartitionerCitation[] =
1632:   "@article{PTSCOTCH,\n"
1633:   "  author  = {C. Chevalier and F. Pellegrini},\n"
1634:   "  title   = {{PT-SCOTCH}: a tool for efficient parallel graph ordering},\n"
1635:   "  journal = {Parallel Computing},\n"
1636:   "  volume  = {34},\n"
1637:   "  number  = {6},\n"
1638:   "  pages   = {318--331},\n"
1639:   "  year    = {2008},\n"
1640:   "  doi     = {https://doi.org/10.1016/j.parco.2007.12.001}\n"
1641:   "}\n";

1643: typedef struct {
1644:   PetscInt  strategy;
1645:   PetscReal imbalance;
1646: } PetscPartitioner_PTScotch;

1648: static const char *const
1649: PTScotchStrategyList[] = {
1650:   "DEFAULT",
1651:   "QUALITY",
1652:   "SPEED",
1653:   "BALANCE",
1654:   "SAFETY",
1655:   "SCALABILITY",
1656:   "RECURSIVE",
1657:   "REMAP"
1658: };

1660: #if defined(PETSC_HAVE_PTSCOTCH)

1662: EXTERN_C_BEGIN
1663: #include <ptscotch.h>
1664: EXTERN_C_END

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

1668: static int PTScotch_Strategy(PetscInt strategy)
1669: {
1670:   switch (strategy) {
1671:   case  0: return SCOTCH_STRATDEFAULT;
1672:   case  1: return SCOTCH_STRATQUALITY;
1673:   case  2: return SCOTCH_STRATSPEED;
1674:   case  3: return SCOTCH_STRATBALANCE;
1675:   case  4: return SCOTCH_STRATSAFETY;
1676:   case  5: return SCOTCH_STRATSCALABILITY;
1677:   case  6: return SCOTCH_STRATRECURSIVE;
1678:   case  7: return SCOTCH_STRATREMAP;
1679:   default: return SCOTCH_STRATDEFAULT;
1680:   }
1681: }


1684: static PetscErrorCode PTScotch_PartGraph_Seq(SCOTCH_Num strategy, double imbalance, SCOTCH_Num n, SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
1685:                                              SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num part[])
1686: {
1687:   SCOTCH_Graph   grafdat;
1688:   SCOTCH_Strat   stradat;
1689:   SCOTCH_Num     vertnbr = n;
1690:   SCOTCH_Num     edgenbr = xadj[n];
1691:   SCOTCH_Num*    velotab = vtxwgt;
1692:   SCOTCH_Num*    edlotab = adjwgt;
1693:   SCOTCH_Num     flagval = strategy;
1694:   double         kbalval = imbalance;

1698:   {
1699:     PetscBool flg = PETSC_TRUE;
1700:     PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
1701:     if (!flg) velotab = NULL;
1702:   }
1703:   SCOTCH_graphInit(&grafdat);CHKERRPTSCOTCH(ierr);
1704:   SCOTCH_graphBuild(&grafdat, 0, vertnbr, xadj, xadj + 1, velotab, NULL, edgenbr, adjncy, edlotab);CHKERRPTSCOTCH(ierr);
1705:   SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
1706:   SCOTCH_stratGraphMapBuild(&stradat, flagval, nparts, kbalval);CHKERRPTSCOTCH(ierr);
1707: #if defined(PETSC_USE_DEBUG)
1708:   SCOTCH_graphCheck(&grafdat);CHKERRPTSCOTCH(ierr);
1709: #endif
1710:   SCOTCH_graphPart(&grafdat, nparts, &stradat, part);CHKERRPTSCOTCH(ierr);
1711:   SCOTCH_stratExit(&stradat);
1712:   SCOTCH_graphExit(&grafdat);
1713:   return(0);
1714: }

1716: static PetscErrorCode PTScotch_PartGraph_MPI(SCOTCH_Num strategy, double imbalance, SCOTCH_Num vtxdist[], SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
1717:                                              SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num part[], MPI_Comm comm)
1718: {
1719:   PetscMPIInt     procglbnbr;
1720:   PetscMPIInt     proclocnum;
1721:   SCOTCH_Arch     archdat;
1722:   SCOTCH_Dgraph   grafdat;
1723:   SCOTCH_Dmapping mappdat;
1724:   SCOTCH_Strat    stradat;
1725:   SCOTCH_Num      vertlocnbr;
1726:   SCOTCH_Num      edgelocnbr;
1727:   SCOTCH_Num*     veloloctab = vtxwgt;
1728:   SCOTCH_Num*     edloloctab = adjwgt;
1729:   SCOTCH_Num      flagval = strategy;
1730:   double          kbalval = imbalance;
1731:   PetscErrorCode  ierr;

1734:   {
1735:     PetscBool flg = PETSC_TRUE;
1736:     PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
1737:     if (!flg) veloloctab = NULL;
1738:   }
1739:   MPI_Comm_size(comm, &procglbnbr);
1740:   MPI_Comm_rank(comm, &proclocnum);
1741:   vertlocnbr = vtxdist[proclocnum + 1] - vtxdist[proclocnum];
1742:   edgelocnbr = xadj[vertlocnbr];

1744:   SCOTCH_dgraphInit(&grafdat, comm);CHKERRPTSCOTCH(ierr);
1745:   SCOTCH_dgraphBuild(&grafdat, 0, vertlocnbr, vertlocnbr, xadj, xadj + 1, veloloctab, NULL, edgelocnbr, edgelocnbr, adjncy, NULL, edloloctab);CHKERRPTSCOTCH(ierr);
1746: #if defined(PETSC_USE_DEBUG)
1747:   SCOTCH_dgraphCheck(&grafdat);CHKERRPTSCOTCH(ierr);
1748: #endif
1749:   SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
1750:   SCOTCH_stratDgraphMapBuild(&stradat, flagval, procglbnbr, nparts, kbalval);
1751:   SCOTCH_archInit(&archdat);CHKERRPTSCOTCH(ierr);
1752:   SCOTCH_archCmplt(&archdat, nparts);CHKERRPTSCOTCH(ierr);
1753:   SCOTCH_dgraphMapInit(&grafdat, &mappdat, &archdat, part);CHKERRPTSCOTCH(ierr);
1754:   SCOTCH_dgraphMapCompute(&grafdat, &mappdat, &stradat);CHKERRPTSCOTCH(ierr);
1755:   SCOTCH_dgraphMapExit(&grafdat, &mappdat);
1756:   SCOTCH_archExit(&archdat);
1757:   SCOTCH_stratExit(&stradat);
1758:   SCOTCH_dgraphExit(&grafdat);
1759:   return(0);
1760: }

1762: #endif /* PETSC_HAVE_PTSCOTCH */

1764: static PetscErrorCode PetscPartitionerDestroy_PTScotch(PetscPartitioner part)
1765: {
1766:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
1767:   PetscErrorCode             ierr;

1770:   PetscFree(p);
1771:   return(0);
1772: }

1774: static PetscErrorCode PetscPartitionerView_PTScotch_Ascii(PetscPartitioner part, PetscViewer viewer)
1775: {
1776:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
1777:   PetscErrorCode            ierr;

1780:   PetscViewerASCIIPushTab(viewer);
1781:   PetscViewerASCIIPrintf(viewer, "using partitioning strategy %s\n",PTScotchStrategyList[p->strategy]);
1782:   PetscViewerASCIIPrintf(viewer, "using load imbalance ratio %g\n",(double)p->imbalance);
1783:   PetscViewerASCIIPopTab(viewer);
1784:   return(0);
1785: }

1787: static PetscErrorCode PetscPartitionerView_PTScotch(PetscPartitioner part, PetscViewer viewer)
1788: {
1789:   PetscBool      iascii;

1795:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1796:   if (iascii) {PetscPartitionerView_PTScotch_Ascii(part, viewer);}
1797:   return(0);
1798: }

1800: static PetscErrorCode PetscPartitionerSetFromOptions_PTScotch(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1801: {
1802:   PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
1803:   const char *const         *slist = PTScotchStrategyList;
1804:   PetscInt                  nlist = (PetscInt)(sizeof(PTScotchStrategyList)/sizeof(PTScotchStrategyList[0]));
1805:   PetscBool                 flag;
1806:   PetscErrorCode            ierr;

1809:   PetscOptionsHead(PetscOptionsObject, "PetscPartitioner PTScotch Options");
1810:   PetscOptionsEList("-petscpartitioner_ptscotch_strategy","Partitioning strategy","",slist,nlist,slist[p->strategy],&p->strategy,&flag);
1811:   PetscOptionsReal("-petscpartitioner_ptscotch_imbalance","Load imbalance ratio","",p->imbalance,&p->imbalance,&flag);
1812:   PetscOptionsTail();
1813:   return(0);
1814: }

1816: static PetscErrorCode PetscPartitionerPartition_PTScotch(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1817: {
1818: #if defined(PETSC_HAVE_PTSCOTCH)
1819:   MPI_Comm       comm       = PetscObjectComm((PetscObject)part);
1820:   PetscInt       nvtxs      = numVertices; /* The number of vertices in full graph */
1821:   PetscInt      *vtxdist;                  /* Distribution of vertices across processes */
1822:   PetscInt      *xadj       = start;       /* Start of edge list for each vertex */
1823:   PetscInt      *adjncy     = adjacency;   /* Edge lists for all vertices */
1824:   PetscInt      *vwgt       = NULL;        /* Vertex weights */
1825:   PetscInt      *adjwgt     = NULL;        /* Edge weights */
1826:   PetscInt       v, i, *assignment, *points;
1827:   PetscMPIInt    size, rank, p;

1831:   MPI_Comm_size(comm, &size);
1832:   MPI_Comm_rank(comm, &rank);
1833:   PetscMalloc2(nparts+1,&vtxdist,PetscMax(nvtxs,1),&assignment);

1835:   /* Calculate vertex distribution */
1836:   vtxdist[0] = 0;
1837:   MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
1838:   for (p = 2; p <= size; ++p) {
1839:     vtxdist[p] += vtxdist[p-1];
1840:   }

1842:   if (nparts == 1) {
1843:     PetscMemzero(assignment, nvtxs * sizeof(PetscInt));
1844:   } else {
1845:     PetscSection section;
1846:     /* Weight cells by dofs on cell by default */
1847:     PetscMalloc1(PetscMax(nvtxs,1),&vwgt);
1848:     DMGetSection(dm, &section);
1849:     if (section) {
1850:       PetscInt vStart, vEnd, dof;
1851:       DMPlexGetHeightStratum(dm, 0, &vStart, &vEnd);
1852:       for (v = vStart; v < vEnd; ++v) {
1853:         PetscSectionGetDof(section, v, &dof);
1854:         /* WARNING: Assumes that meshes with overlap have the overlapped cells at the end of the stratum. */
1855:         /* To do this properly, we should use the cell numbering created in DMPlexCreatePartitionerGraph. */
1856:         if (v-vStart < numVertices) vwgt[v-vStart] = PetscMax(dof, 1);
1857:       }
1858:     } else {
1859:       for (v = 0; v < nvtxs; ++v) vwgt[v] = 1;
1860:     }
1861:     {
1862:       PetscPartitioner_PTScotch *pts = (PetscPartitioner_PTScotch *) part->data;
1863:       int                       strat = PTScotch_Strategy(pts->strategy);
1864:       double                    imbal = (double)pts->imbalance;

1866:       for (p = 0; !vtxdist[p+1] && p < size; ++p);
1867:       if (vtxdist[p+1] == vtxdist[size]) {
1868:         if (rank == p) {
1869:           PTScotch_PartGraph_Seq(strat, imbal, nvtxs, xadj, adjncy, vwgt, adjwgt, nparts, assignment);
1870:         }
1871:       } else {
1872:         PTScotch_PartGraph_MPI(strat, imbal, vtxdist, xadj, adjncy, vwgt, adjwgt, nparts, assignment, comm);
1873:       }
1874:     }
1875:     PetscFree(vwgt);
1876:   }

1878:   /* Convert to PetscSection+IS */
1879:   PetscSectionSetChart(partSection, 0, nparts);
1880:   for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
1881:   PetscSectionSetUp(partSection);
1882:   PetscMalloc1(nvtxs, &points);
1883:   for (p = 0, i = 0; p < nparts; ++p) {
1884:     for (v = 0; v < nvtxs; ++v) {
1885:       if (assignment[v] == p) points[i++] = v;
1886:     }
1887:   }
1888:   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1889:   ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);

1891:   PetscFree2(vtxdist,assignment);
1892:   return(0);
1893: #else
1894:   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-ptscotch.");
1895: #endif
1896: }

1898: static PetscErrorCode PetscPartitionerInitialize_PTScotch(PetscPartitioner part)
1899: {
1901:   part->ops->view           = PetscPartitionerView_PTScotch;
1902:   part->ops->destroy        = PetscPartitionerDestroy_PTScotch;
1903:   part->ops->partition      = PetscPartitionerPartition_PTScotch;
1904:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_PTScotch;
1905:   return(0);
1906: }

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

1911:   Level: intermediate

1913: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1914: M*/

1916: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_PTScotch(PetscPartitioner part)
1917: {
1918:   PetscPartitioner_PTScotch *p;
1919:   PetscErrorCode          ierr;

1923:   PetscNewLog(part, &p);
1924:   part->data = p;

1926:   p->strategy  = 0;
1927:   p->imbalance = 0.01;

1929:   PetscPartitionerInitialize_PTScotch(part);
1930:   PetscCitationsRegister(PTScotchPartitionerCitation, &PTScotchPartitionercite);
1931:   return(0);
1932: }


1935: /*@
1936:   DMPlexGetPartitioner - Get the mesh partitioner

1938:   Not collective

1940:   Input Parameter:
1941: . dm - The DM

1943:   Output Parameter:
1944: . part - The PetscPartitioner

1946:   Level: developer

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

1950: .seealso DMPlexDistribute(), DMPlexSetPartitioner(), PetscPartitionerCreate()
1951: @*/
1952: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
1953: {
1954:   DM_Plex       *mesh = (DM_Plex *) dm->data;

1959:   *part = mesh->partitioner;
1960:   return(0);
1961: }

1963: /*@
1964:   DMPlexSetPartitioner - Set the mesh partitioner

1966:   logically collective on dm and part

1968:   Input Parameters:
1969: + dm - The DM
1970: - part - The partitioner

1972:   Level: developer

1974:   Note: Any existing PetscPartitioner will be destroyed.

1976: .seealso DMPlexDistribute(), DMPlexGetPartitioner(), PetscPartitionerCreate()
1977: @*/
1978: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
1979: {
1980:   DM_Plex       *mesh = (DM_Plex *) dm->data;

1986:   PetscObjectReference((PetscObject)part);
1987:   PetscPartitionerDestroy(&mesh->partitioner);
1988:   mesh->partitioner = part;
1989:   return(0);
1990: }

1992: static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
1993: {

1997:   if (up) {
1998:     PetscInt parent;

2000:     DMPlexGetTreeParent(dm,point,&parent,NULL);
2001:     if (parent != point) {
2002:       PetscInt closureSize, *closure = NULL, i;

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

2008:         PetscHSetIAdd(ht, cpoint);
2009:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
2010:       }
2011:       DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
2012:     }
2013:   }
2014:   if (down) {
2015:     PetscInt numChildren;
2016:     const PetscInt *children;

2018:     DMPlexGetTreeChildren(dm,point,&numChildren,&children);
2019:     if (numChildren) {
2020:       PetscInt i;

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

2025:         PetscHSetIAdd(ht, cpoint);
2026:         DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
2027:       }
2028:     }
2029:   }
2030:   return(0);
2031: }

2033: PETSC_INTERN PetscErrorCode DMPlexPartitionLabelClosure_Private(DM,DMLabel,PetscInt,PetscInt,const PetscInt[],IS*);

2035: PetscErrorCode DMPlexPartitionLabelClosure_Private(DM dm, DMLabel label, PetscInt rank, PetscInt numPoints, const PetscInt points[], IS *closureIS)
2036: {
2037:   DM_Plex        *mesh = (DM_Plex *)dm->data;
2038:   PetscBool      hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
2039:   PetscInt       *closure = NULL, closureSize, nelems, *elems, off = 0, p, c;
2040:   PetscHSetI     ht;

2044:   PetscHSetICreate(&ht);
2045:   PetscHSetIResize(ht, numPoints*16);
2046:   for (p = 0; p < numPoints; ++p) {
2047:     DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2048:     for (c = 0; c < closureSize*2; c += 2) {
2049:       PetscHSetIAdd(ht, closure[c]);
2050:       if (hasTree) {DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);}
2051:     }
2052:   }
2053:   if (closure) {DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);}
2054:   PetscHSetIGetSize(ht, &nelems);
2055:   PetscMalloc1(nelems, &elems);
2056:   PetscHSetIGetElems(ht, &off, elems);
2057:   PetscHSetIDestroy(&ht);
2058:   PetscSortInt(nelems, elems);
2059:   ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
2060:   return(0);
2061: }

2063: /*@
2064:   DMPlexPartitionLabelClosure - Add the closure of all points to the partition label

2066:   Input Parameters:
2067: + dm     - The DM
2068: - label  - DMLabel assinging ranks to remote roots

2070:   Level: developer

2072: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2073: @*/
2074: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
2075: {
2076:   IS              rankIS,   pointIS, closureIS;
2077:   const PetscInt *ranks,   *points;
2078:   PetscInt        numRanks, numPoints, r;
2079:   PetscErrorCode  ierr;

2082:   DMLabelGetValueIS(label, &rankIS);
2083:   ISGetLocalSize(rankIS, &numRanks);
2084:   ISGetIndices(rankIS, &ranks);
2085:   for (r = 0; r < numRanks; ++r) {
2086:     const PetscInt rank = ranks[r];
2087:     DMLabelGetStratumIS(label, rank, &pointIS);
2088:     ISGetLocalSize(pointIS, &numPoints);
2089:     ISGetIndices(pointIS, &points);
2090:     DMPlexPartitionLabelClosure_Private(dm, label, rank, numPoints, points, &closureIS);
2091:     ISRestoreIndices(pointIS, &points);
2092:     ISDestroy(&pointIS);
2093:     DMLabelSetStratumIS(label, rank, closureIS);
2094:     ISDestroy(&closureIS);
2095:   }
2096:   ISRestoreIndices(rankIS, &ranks);
2097:   ISDestroy(&rankIS);
2098:   return(0);
2099: }

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

2104:   Input Parameters:
2105: + dm     - The DM
2106: - label  - DMLabel assinging ranks to remote roots

2108:   Level: developer

2110: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2111: @*/
2112: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
2113: {
2114:   IS              rankIS,   pointIS;
2115:   const PetscInt *ranks,   *points;
2116:   PetscInt        numRanks, numPoints, r, p, a, adjSize;
2117:   PetscInt       *adj = NULL;
2118:   PetscErrorCode  ierr;

2121:   DMLabelGetValueIS(label, &rankIS);
2122:   ISGetLocalSize(rankIS, &numRanks);
2123:   ISGetIndices(rankIS, &ranks);
2124:   for (r = 0; r < numRanks; ++r) {
2125:     const PetscInt rank = ranks[r];

2127:     DMLabelGetStratumIS(label, rank, &pointIS);
2128:     ISGetLocalSize(pointIS, &numPoints);
2129:     ISGetIndices(pointIS, &points);
2130:     for (p = 0; p < numPoints; ++p) {
2131:       adjSize = PETSC_DETERMINE;
2132:       DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
2133:       for (a = 0; a < adjSize; ++a) {DMLabelSetValue(label, adj[a], rank);}
2134:     }
2135:     ISRestoreIndices(pointIS, &points);
2136:     ISDestroy(&pointIS);
2137:   }
2138:   ISRestoreIndices(rankIS, &ranks);
2139:   ISDestroy(&rankIS);
2140:   PetscFree(adj);
2141:   return(0);
2142: }

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

2147:   Input Parameters:
2148: + dm     - The DM
2149: - label  - DMLabel assinging ranks to remote roots

2151:   Level: developer

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

2156: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2157: @*/
2158: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
2159: {
2160:   MPI_Comm        comm;
2161:   PetscMPIInt     rank;
2162:   PetscSF         sfPoint;
2163:   DMLabel         lblRoots, lblLeaves;
2164:   IS              rankIS, pointIS;
2165:   const PetscInt *ranks;
2166:   PetscInt        numRanks, r;
2167:   PetscErrorCode  ierr;

2170:   PetscObjectGetComm((PetscObject) dm, &comm);
2171:   MPI_Comm_rank(comm, &rank);
2172:   DMGetPointSF(dm, &sfPoint);
2173:   /* Pull point contributions from remote leaves into local roots */
2174:   DMLabelGather(label, sfPoint, &lblLeaves);
2175:   DMLabelGetValueIS(lblLeaves, &rankIS);
2176:   ISGetLocalSize(rankIS, &numRanks);
2177:   ISGetIndices(rankIS, &ranks);
2178:   for (r = 0; r < numRanks; ++r) {
2179:     const PetscInt remoteRank = ranks[r];
2180:     if (remoteRank == rank) continue;
2181:     DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
2182:     DMLabelInsertIS(label, pointIS, remoteRank);
2183:     ISDestroy(&pointIS);
2184:   }
2185:   ISRestoreIndices(rankIS, &ranks);
2186:   ISDestroy(&rankIS);
2187:   DMLabelDestroy(&lblLeaves);
2188:   /* Push point contributions from roots into remote leaves */
2189:   DMLabelDistribute(label, sfPoint, &lblRoots);
2190:   DMLabelGetValueIS(lblRoots, &rankIS);
2191:   ISGetLocalSize(rankIS, &numRanks);
2192:   ISGetIndices(rankIS, &ranks);
2193:   for (r = 0; r < numRanks; ++r) {
2194:     const PetscInt remoteRank = ranks[r];
2195:     if (remoteRank == rank) continue;
2196:     DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
2197:     DMLabelInsertIS(label, pointIS, remoteRank);
2198:     ISDestroy(&pointIS);
2199:   }
2200:   ISRestoreIndices(rankIS, &ranks);
2201:   ISDestroy(&rankIS);
2202:   DMLabelDestroy(&lblRoots);
2203:   return(0);
2204: }

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

2209:   Input Parameters:
2210: + dm        - The DM
2211: . rootLabel - DMLabel assinging ranks to local roots
2212: . processSF - A star forest mapping into the local index on each remote rank

2214:   Output Parameter:
2215: - leafLabel - DMLabel assinging ranks to remote roots

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

2220:   Level: developer

2222: .seealso: DMPlexPartitionLabelCreateSF, DMPlexDistribute(), DMPlexCreateOverlap
2223: @*/
2224: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
2225: {
2226:   MPI_Comm           comm;
2227:   PetscMPIInt        rank, size, r;
2228:   PetscInt           p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
2229:   PetscSF            sfPoint;
2230:   PetscSection       rootSection;
2231:   PetscSFNode       *rootPoints, *leafPoints;
2232:   const PetscSFNode *remote;
2233:   const PetscInt    *local, *neighbors;
2234:   IS                 valueIS;
2235:   PetscBool          mpiOverflow = PETSC_FALSE;
2236:   PetscErrorCode     ierr;

2239:   PetscObjectGetComm((PetscObject) dm, &comm);
2240:   MPI_Comm_rank(comm, &rank);
2241:   MPI_Comm_size(comm, &size);
2242:   DMGetPointSF(dm, &sfPoint);

2244:   /* Convert to (point, rank) and use actual owners */
2245:   PetscSectionCreate(comm, &rootSection);
2246:   PetscSectionSetChart(rootSection, 0, size);
2247:   DMLabelGetValueIS(rootLabel, &valueIS);
2248:   ISGetLocalSize(valueIS, &numNeighbors);
2249:   ISGetIndices(valueIS, &neighbors);
2250:   for (n = 0; n < numNeighbors; ++n) {
2251:     DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
2252:     PetscSectionAddDof(rootSection, neighbors[n], numPoints);
2253:   }
2254:   PetscSectionSetUp(rootSection);
2255:   PetscSectionGetStorageSize(rootSection, &rootSize);
2256:   PetscMalloc1(rootSize, &rootPoints);
2257:   PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
2258:   for (n = 0; n < numNeighbors; ++n) {
2259:     IS              pointIS;
2260:     const PetscInt *points;

2262:     PetscSectionGetOffset(rootSection, neighbors[n], &off);
2263:     DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
2264:     ISGetLocalSize(pointIS, &numPoints);
2265:     ISGetIndices(pointIS, &points);
2266:     for (p = 0; p < numPoints; ++p) {
2267:       if (local) {PetscFindInt(points[p], nleaves, local, &l);}
2268:       else       {l = -1;}
2269:       if (l >= 0) {rootPoints[off+p] = remote[l];}
2270:       else        {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
2271:     }
2272:     ISRestoreIndices(pointIS, &points);
2273:     ISDestroy(&pointIS);
2274:   }

2276:   /* Try to communicate overlap using All-to-All */
2277:   if (!processSF) {
2278:     PetscInt64  counter = 0;
2279:     PetscBool   locOverflow = PETSC_FALSE;
2280:     PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;

2282:     PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
2283:     for (n = 0; n < numNeighbors; ++n) {
2284:       PetscSectionGetDof(rootSection, neighbors[n], &dof);
2285:       PetscSectionGetOffset(rootSection, neighbors[n], &off);
2286: #if defined(PETSC_USE_64BIT_INDICES)
2287:       if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2288:       if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2289: #endif
2290:       scounts[neighbors[n]] = (PetscMPIInt) dof;
2291:       sdispls[neighbors[n]] = (PetscMPIInt) off;
2292:     }
2293:     MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
2294:     for (r = 0; r < size; ++r) { rdispls[r] = (int)counter; counter += rcounts[r]; }
2295:     if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
2296:     MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
2297:     if (!mpiOverflow) {
2298:       leafSize = (PetscInt) counter;
2299:       PetscMalloc1(leafSize, &leafPoints);
2300:       MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
2301:     }
2302:     PetscFree4(scounts, sdispls, rcounts, rdispls);
2303:   }

2305:   /* Communicate overlap using process star forest */
2306:   if (processSF || mpiOverflow) {
2307:     PetscSF      procSF;
2308:     PetscSFNode  *remote;
2309:     PetscSection leafSection;

2311:     if (processSF) {
2312:       PetscObjectReference((PetscObject)processSF);
2313:       procSF = processSF;
2314:     } else {
2315:       PetscMalloc1(size, &remote);
2316:       for (r = 0; r < size; ++r) { remote[r].rank  = r; remote[r].index = rank; }
2317:       PetscSFCreate(comm, &procSF);
2318:       PetscSFSetGraph(procSF, size, size, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
2319:     }

2321:     PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
2322:     DMPlexDistributeData(dm, processSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
2323:     PetscSectionGetStorageSize(leafSection, &leafSize);
2324:     PetscSectionDestroy(&leafSection);
2325:     PetscSFDestroy(&procSF);
2326:   }

2328:   for (p = 0; p < leafSize; p++) {
2329:     DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
2330:   }

2332:   ISRestoreIndices(valueIS, &neighbors);
2333:   ISDestroy(&valueIS);
2334:   PetscSectionDestroy(&rootSection);
2335:   PetscFree(rootPoints);
2336:   PetscFree(leafPoints);
2337:   return(0);
2338: }

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

2343:   Input Parameters:
2344: + dm    - The DM
2345: . label - DMLabel assinging ranks to remote roots

2347:   Output Parameter:
2348: - sf    - The star forest communication context encapsulating the defined mapping

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

2352:   Level: developer

2354: .seealso: DMPlexDistribute(), DMPlexCreateOverlap
2355: @*/
2356: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
2357: {
2358:   PetscMPIInt     rank, size;
2359:   PetscInt        n, numRemote, p, numPoints, pStart, pEnd, idx = 0;
2360:   PetscSFNode    *remotePoints;
2361:   IS              remoteRootIS;
2362:   const PetscInt *remoteRoots;

2366:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2367:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);

2369:   for (numRemote = 0, n = 0; n < size; ++n) {
2370:     DMLabelGetStratumSize(label, n, &numPoints);
2371:     numRemote += numPoints;
2372:   }
2373:   PetscMalloc1(numRemote, &remotePoints);
2374:   /* Put owned points first */
2375:   DMLabelGetStratumSize(label, rank, &numPoints);
2376:   if (numPoints > 0) {
2377:     DMLabelGetStratumIS(label, rank, &remoteRootIS);
2378:     ISGetIndices(remoteRootIS, &remoteRoots);
2379:     for (p = 0; p < numPoints; p++) {
2380:       remotePoints[idx].index = remoteRoots[p];
2381:       remotePoints[idx].rank = rank;
2382:       idx++;
2383:     }
2384:     ISRestoreIndices(remoteRootIS, &remoteRoots);
2385:     ISDestroy(&remoteRootIS);
2386:   }
2387:   /* Now add remote points */
2388:   for (n = 0; n < size; ++n) {
2389:     DMLabelGetStratumSize(label, n, &numPoints);
2390:     if (numPoints <= 0 || n == rank) continue;
2391:     DMLabelGetStratumIS(label, n, &remoteRootIS);
2392:     ISGetIndices(remoteRootIS, &remoteRoots);
2393:     for (p = 0; p < numPoints; p++) {
2394:       remotePoints[idx].index = remoteRoots[p];
2395:       remotePoints[idx].rank = n;
2396:       idx++;
2397:     }
2398:     ISRestoreIndices(remoteRootIS, &remoteRoots);
2399:     ISDestroy(&remoteRootIS);
2400:   }
2401:   PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
2402:   DMPlexGetChart(dm, &pStart, &pEnd);
2403:   PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
2404:   return(0);
2405: }