xref: /petsc/src/dm/impls/plex/plexpartition.c (revision 4b15ede2af8a3e0bc78f86cb79f2c0c764c0a2d7) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 PetscClassId PETSCPARTITIONER_CLASSID = 0;
4 
5 PetscFunctionList PetscPartitionerList              = NULL;
6 PetscBool         PetscPartitionerRegisterAllCalled = PETSC_FALSE;
7 
8 PetscBool ChacoPartitionercite = PETSC_FALSE;
9 const char ChacoPartitionerCitation[] = "@inproceedings{Chaco95,\n"
10                                "  author    = {Bruce Hendrickson and Robert Leland},\n"
11                                "  title     = {A multilevel algorithm for partitioning graphs},\n"
12                                "  booktitle = {Supercomputing '95: Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM)},"
13                                "  isbn      = {0-89791-816-9},\n"
14                                "  pages     = {28},\n"
15                                "  doi       = {http://doi.acm.org/10.1145/224170.224228},\n"
16                                "  publisher = {ACM Press},\n"
17                                "  address   = {New York},\n"
18                                "  year      = {1995}\n}\n";
19 
20 PetscBool ParMetisPartitionercite = PETSC_FALSE;
21 const char ParMetisPartitionerCitation[] = "@article{KarypisKumar98,\n"
22                                "  author  = {George Karypis and Vipin Kumar},\n"
23                                "  title   = {A Parallel Algorithm for Multilevel Graph Partitioning and Sparse Matrix Ordering},\n"
24                                "  journal = {Journal of Parallel and Distributed Computing},\n"
25                                "  volume  = {48},\n"
26                                "  pages   = {71--85},\n"
27                                "  year    = {1998}\n}\n";
28 
29 #undef __FUNCT__
30 #define __FUNCT__ "DMPlexCreateNeighborCSR"
31 PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
32 {
33   const PetscInt maxFaceCases = 30;
34   PetscInt       numFaceCases = 0;
35   PetscInt       numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
36   PetscInt      *off, *adj;
37   PetscInt      *neighborCells = NULL;
38   PetscInt       dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
39   PetscErrorCode ierr;
40 
41   PetscFunctionBegin;
42   /* For parallel partitioning, I think you have to communicate supports */
43   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
44   cellDim = dim - cellHeight;
45   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
46   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
47   if (cEnd - cStart == 0) {
48     if (numVertices) *numVertices = 0;
49     if (offsets)   *offsets   = NULL;
50     if (adjacency) *adjacency = NULL;
51     PetscFunctionReturn(0);
52   }
53   numCells  = cEnd - cStart;
54   faceDepth = depth - cellHeight;
55   if (dim == depth) {
56     PetscInt f, fStart, fEnd;
57 
58     ierr = PetscCalloc1(numCells+1, &off);CHKERRQ(ierr);
59     /* Count neighboring cells */
60     ierr = DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);CHKERRQ(ierr);
61     for (f = fStart; f < fEnd; ++f) {
62       const PetscInt *support;
63       PetscInt        supportSize;
64       ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr);
65       ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr);
66       if (supportSize == 2) {
67         PetscInt numChildren;
68 
69         ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr);
70         if (!numChildren) {
71           ++off[support[0]-cStart+1];
72           ++off[support[1]-cStart+1];
73         }
74       }
75     }
76     /* Prefix sum */
77     for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
78     if (adjacency) {
79       PetscInt *tmp;
80 
81       ierr = PetscMalloc1(off[numCells], &adj);CHKERRQ(ierr);
82       ierr = PetscMalloc1((numCells+1), &tmp);CHKERRQ(ierr);
83       ierr = PetscMemcpy(tmp, off, (numCells+1) * sizeof(PetscInt));CHKERRQ(ierr);
84       /* Get neighboring cells */
85       for (f = fStart; f < fEnd; ++f) {
86         const PetscInt *support;
87         PetscInt        supportSize;
88         ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr);
89         ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr);
90         if (supportSize == 2) {
91           PetscInt numChildren;
92 
93           ierr = DMPlexGetTreeChildren(dm,f,&numChildren,NULL);CHKERRQ(ierr);
94           if (!numChildren) {
95             adj[tmp[support[0]-cStart]++] = support[1];
96             adj[tmp[support[1]-cStart]++] = support[0];
97           }
98         }
99       }
100 #if defined(PETSC_USE_DEBUG)
101       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);
102 #endif
103       ierr = PetscFree(tmp);CHKERRQ(ierr);
104     }
105     if (numVertices) *numVertices = numCells;
106     if (offsets)   *offsets   = off;
107     if (adjacency) *adjacency = adj;
108     PetscFunctionReturn(0);
109   }
110   /* Setup face recognition */
111   if (faceDepth == 1) {
112     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 */
113 
114     for (c = cStart; c < cEnd; ++c) {
115       PetscInt corners;
116 
117       ierr = DMPlexGetConeSize(dm, c, &corners);CHKERRQ(ierr);
118       if (!cornersSeen[corners]) {
119         PetscInt nFV;
120 
121         if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
122         cornersSeen[corners] = 1;
123 
124         ierr = DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);CHKERRQ(ierr);
125 
126         numFaceVertices[numFaceCases++] = nFV;
127       }
128     }
129   }
130   ierr = PetscCalloc1(numCells+1, &off);CHKERRQ(ierr);
131   /* Count neighboring cells */
132   for (cell = cStart; cell < cEnd; ++cell) {
133     PetscInt numNeighbors = PETSC_DETERMINE, n;
134 
135     ierr = DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);CHKERRQ(ierr);
136     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
137     for (n = 0; n < numNeighbors; ++n) {
138       PetscInt        cellPair[2];
139       PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
140       PetscInt        meetSize = 0;
141       const PetscInt *meet    = NULL;
142 
143       cellPair[0] = cell; cellPair[1] = neighborCells[n];
144       if (cellPair[0] == cellPair[1]) continue;
145       if (!found) {
146         ierr = DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr);
147         if (meetSize) {
148           PetscInt f;
149 
150           for (f = 0; f < numFaceCases; ++f) {
151             if (numFaceVertices[f] == meetSize) {
152               found = PETSC_TRUE;
153               break;
154             }
155           }
156         }
157         ierr = DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr);
158       }
159       if (found) ++off[cell-cStart+1];
160     }
161   }
162   /* Prefix sum */
163   for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];
164 
165   if (adjacency) {
166     ierr = PetscMalloc1(off[numCells], &adj);CHKERRQ(ierr);
167     /* Get neighboring cells */
168     for (cell = cStart; cell < cEnd; ++cell) {
169       PetscInt numNeighbors = PETSC_DETERMINE, n;
170       PetscInt cellOffset   = 0;
171 
172       ierr = DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);CHKERRQ(ierr);
173       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
174       for (n = 0; n < numNeighbors; ++n) {
175         PetscInt        cellPair[2];
176         PetscBool       found    = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
177         PetscInt        meetSize = 0;
178         const PetscInt *meet    = NULL;
179 
180         cellPair[0] = cell; cellPair[1] = neighborCells[n];
181         if (cellPair[0] == cellPair[1]) continue;
182         if (!found) {
183           ierr = DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr);
184           if (meetSize) {
185             PetscInt f;
186 
187             for (f = 0; f < numFaceCases; ++f) {
188               if (numFaceVertices[f] == meetSize) {
189                 found = PETSC_TRUE;
190                 break;
191               }
192             }
193           }
194           ierr = DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr);
195         }
196         if (found) {
197           adj[off[cell-cStart]+cellOffset] = neighborCells[n];
198           ++cellOffset;
199         }
200       }
201     }
202   }
203   ierr = PetscFree(neighborCells);CHKERRQ(ierr);
204   if (numVertices) *numVertices = numCells;
205   if (offsets)   *offsets   = off;
206   if (adjacency) *adjacency = adj;
207   PetscFunctionReturn(0);
208 }
209 
210 #undef __FUNCT__
211 #define __FUNCT__ "DMPlexEnlargePartition"
212 /* Expand the partition by BFS on the adjacency graph */
213 PetscErrorCode DMPlexEnlargePartition(DM dm, const PetscInt start[], const PetscInt adjacency[], PetscSection origPartSection, IS origPartition, PetscSection partSection, IS *partition)
214 {
215   PetscHashI      h;
216   const PetscInt *points;
217   PetscInt      **tmpPoints, *newPoints, totPoints = 0;
218   PetscInt        pStart, pEnd, part, q;
219   PetscBool       useCone;
220   PetscErrorCode  ierr;
221 
222   PetscFunctionBegin;
223   PetscHashICreate(h);
224   ierr = PetscSectionGetChart(origPartSection, &pStart, &pEnd);CHKERRQ(ierr);
225   ierr = PetscSectionSetChart(partSection, pStart, pEnd);CHKERRQ(ierr);
226   ierr = ISGetIndices(origPartition, &points);CHKERRQ(ierr);
227   ierr = PetscMalloc1((pEnd - pStart), &tmpPoints);CHKERRQ(ierr);
228   ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr);
229   ierr = DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);CHKERRQ(ierr);
230   for (part = pStart; part < pEnd; ++part) {
231     PetscInt *adj = NULL;
232     PetscInt  numPoints, nP, numNewPoints, off, p, n = 0;
233 
234     PetscHashIClear(h);
235     ierr = PetscSectionGetDof(origPartSection, part, &numPoints);CHKERRQ(ierr);
236     ierr = PetscSectionGetOffset(origPartSection, part, &off);CHKERRQ(ierr);
237     /* Add all existing points to h */
238     for (p = 0; p < numPoints; ++p) {
239       const PetscInt point = points[off+p];
240       PetscHashIAdd(h, point, 1);
241     }
242     PetscHashISize(h, nP);
243     if (nP != numPoints) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid partition has %d points, but only %d were unique", numPoints, nP);
244     /* Add all points in next BFS level */
245     for (p = 0; p < numPoints; ++p) {
246       const PetscInt point   = points[off+p];
247       PetscInt       adjSize = PETSC_DETERMINE, a;
248 
249       ierr = DMPlexGetAdjacency(dm, point, &adjSize, &adj);CHKERRQ(ierr);
250       for (a = 0; a < adjSize; ++a) PetscHashIAdd(h, adj[a], 1);
251     }
252     PetscHashISize(h, numNewPoints);
253     ierr = PetscSectionSetDof(partSection, part, numNewPoints);CHKERRQ(ierr);
254     ierr = PetscMalloc1(numNewPoints, &tmpPoints[part]);CHKERRQ(ierr);
255     ierr = PetscHashIGetKeys(h, &n, tmpPoints[part]);CHKERRQ(ierr);
256     ierr = PetscFree(adj);CHKERRQ(ierr);
257     totPoints += numNewPoints;
258   }
259   ierr = DMPlexSetAdjacencyUseCone(dm, useCone);CHKERRQ(ierr);
260   ierr = ISRestoreIndices(origPartition, &points);CHKERRQ(ierr);
261   PetscHashIDestroy(h);
262   ierr = PetscSectionSetUp(partSection);CHKERRQ(ierr);
263   ierr = PetscMalloc1(totPoints, &newPoints);CHKERRQ(ierr);
264   for (part = pStart, q = 0; part < pEnd; ++part) {
265     PetscInt numPoints, p;
266 
267     ierr = PetscSectionGetDof(partSection, part, &numPoints);CHKERRQ(ierr);
268     for (p = 0; p < numPoints; ++p, ++q) newPoints[q] = tmpPoints[part][p];
269     ierr = PetscFree(tmpPoints[part]);CHKERRQ(ierr);
270   }
271   ierr = PetscFree(tmpPoints);CHKERRQ(ierr);
272   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), totPoints, newPoints, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
273   PetscFunctionReturn(0);
274 }
275 
276 #undef __FUNCT__
277 #define __FUNCT__ "PetscPartitionerRegister"
278 /*@C
279   PetscPartitionerRegister - Adds a new PetscPartitioner implementation
280 
281   Not Collective
282 
283   Input Parameters:
284 + name        - The name of a new user-defined creation routine
285 - create_func - The creation routine itself
286 
287   Notes:
288   PetscPartitionerRegister() may be called multiple times to add several user-defined PetscPartitioners
289 
290   Sample usage:
291 .vb
292     PetscPartitionerRegister("my_part", MyPetscPartitionerCreate);
293 .ve
294 
295   Then, your PetscPartitioner type can be chosen with the procedural interface via
296 .vb
297     PetscPartitionerCreate(MPI_Comm, PetscPartitioner *);
298     PetscPartitionerSetType(PetscPartitioner, "my_part");
299 .ve
300    or at runtime via the option
301 .vb
302     -petscpartitioner_type my_part
303 .ve
304 
305   Level: advanced
306 
307 .keywords: PetscPartitioner, register
308 .seealso: PetscPartitionerRegisterAll(), PetscPartitionerRegisterDestroy()
309 
310 @*/
311 PetscErrorCode PetscPartitionerRegister(const char sname[], PetscErrorCode (*function)(PetscPartitioner))
312 {
313   PetscErrorCode ierr;
314 
315   PetscFunctionBegin;
316   ierr = PetscFunctionListAdd(&PetscPartitionerList, sname, function);CHKERRQ(ierr);
317   PetscFunctionReturn(0);
318 }
319 
320 #undef __FUNCT__
321 #define __FUNCT__ "PetscPartitionerSetType"
322 /*@C
323   PetscPartitionerSetType - Builds a particular PetscPartitioner
324 
325   Collective on PetscPartitioner
326 
327   Input Parameters:
328 + part - The PetscPartitioner object
329 - name - The kind of partitioner
330 
331   Options Database Key:
332 . -petscpartitioner_type <type> - Sets the PetscPartitioner type; use -help for a list of available types
333 
334   Level: intermediate
335 
336 .keywords: PetscPartitioner, set, type
337 .seealso: PetscPartitionerGetType(), PetscPartitionerCreate()
338 @*/
339 PetscErrorCode PetscPartitionerSetType(PetscPartitioner part, PetscPartitionerType name)
340 {
341   PetscErrorCode (*r)(PetscPartitioner);
342   PetscBool      match;
343   PetscErrorCode ierr;
344 
345   PetscFunctionBegin;
346   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
347   ierr = PetscObjectTypeCompare((PetscObject) part, name, &match);CHKERRQ(ierr);
348   if (match) PetscFunctionReturn(0);
349 
350   if (!PetscPartitionerRegisterAllCalled) {ierr = PetscPartitionerRegisterAll();CHKERRQ(ierr);}
351   ierr = PetscFunctionListFind(PetscPartitionerList, name, &r);CHKERRQ(ierr);
352   if (!r) SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscPartitioner type: %s", name);
353 
354   if (part->ops->destroy) {
355     ierr              = (*part->ops->destroy)(part);CHKERRQ(ierr);
356     part->ops->destroy = NULL;
357   }
358   ierr = (*r)(part);CHKERRQ(ierr);
359   ierr = PetscObjectChangeTypeName((PetscObject) part, name);CHKERRQ(ierr);
360   PetscFunctionReturn(0);
361 }
362 
363 #undef __FUNCT__
364 #define __FUNCT__ "PetscPartitionerGetType"
365 /*@C
366   PetscPartitionerGetType - Gets the PetscPartitioner type name (as a string) from the object.
367 
368   Not Collective
369 
370   Input Parameter:
371 . part - The PetscPartitioner
372 
373   Output Parameter:
374 . name - The PetscPartitioner type name
375 
376   Level: intermediate
377 
378 .keywords: PetscPartitioner, get, type, name
379 .seealso: PetscPartitionerSetType(), PetscPartitionerCreate()
380 @*/
381 PetscErrorCode PetscPartitionerGetType(PetscPartitioner part, PetscPartitionerType *name)
382 {
383   PetscErrorCode ierr;
384 
385   PetscFunctionBegin;
386   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
387   PetscValidPointer(name, 2);
388   if (!PetscPartitionerRegisterAllCalled) {ierr = PetscPartitionerRegisterAll();CHKERRQ(ierr);}
389   *name = ((PetscObject) part)->type_name;
390   PetscFunctionReturn(0);
391 }
392 
393 #undef __FUNCT__
394 #define __FUNCT__ "PetscPartitionerView"
395 /*@C
396   PetscPartitionerView - Views a PetscPartitioner
397 
398   Collective on PetscPartitioner
399 
400   Input Parameter:
401 + part - the PetscPartitioner object to view
402 - v    - the viewer
403 
404   Level: developer
405 
406 .seealso: PetscPartitionerDestroy()
407 @*/
408 PetscErrorCode PetscPartitionerView(PetscPartitioner part, PetscViewer v)
409 {
410   PetscErrorCode ierr;
411 
412   PetscFunctionBegin;
413   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
414   if (!v) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) part), &v);CHKERRQ(ierr);}
415   if (part->ops->view) {ierr = (*part->ops->view)(part, v);CHKERRQ(ierr);}
416   PetscFunctionReturn(0);
417 }
418 
419 #undef __FUNCT__
420 #define __FUNCT__ "PetscPartitionerViewFromOptions"
421 /*
422   PetscPartitionerViewFromOptions - Processes command line options to determine if/how a PetscPartitioner is to be viewed.
423 
424   Collective on PetscPartitioner
425 
426   Input Parameters:
427 + part   - the PetscPartitioner
428 . prefix - prefix to use for viewing, or NULL
429 - optionname - option to activate viewing
430 
431   Level: intermediate
432 
433 .keywords: PetscPartitioner, view, options, database
434 .seealso: VecViewFromOptions(), MatViewFromOptions()
435 */
436 PetscErrorCode PetscPartitionerViewFromOptions(PetscPartitioner part, const char prefix[], const char optionname[])
437 {
438   PetscViewer       viewer;
439   PetscViewerFormat format;
440   PetscBool         flg;
441   PetscErrorCode    ierr;
442 
443   PetscFunctionBegin;
444   if (prefix) {ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) part), prefix,                       optionname, &viewer, &format, &flg);CHKERRQ(ierr);}
445   else        {ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject) part), ((PetscObject) part)->prefix, optionname, &viewer, &format, &flg);CHKERRQ(ierr);}
446   if (flg) {
447     ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr);
448     ierr = PetscPartitionerView(part, viewer);CHKERRQ(ierr);
449     ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
450     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
451   }
452   PetscFunctionReturn(0);
453 }
454 #undef __FUNCT__
455 #define __FUNCT__ "PetscPartitionerSetTypeFromOptions_Internal"
456 PetscErrorCode PetscPartitionerSetTypeFromOptions_Internal(PetscPartitioner part)
457 {
458   const char    *defaultType;
459   char           name[256];
460   PetscBool      flg;
461   PetscErrorCode ierr;
462 
463   PetscFunctionBegin;
464   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
465   if (!((PetscObject) part)->type_name) defaultType = PETSCPARTITIONERCHACO;
466   else                                  defaultType = ((PetscObject) part)->type_name;
467   if (!PetscPartitionerRegisterAllCalled) {ierr = PetscPartitionerRegisterAll();CHKERRQ(ierr);}
468 
469   ierr = PetscObjectOptionsBegin((PetscObject) part);CHKERRQ(ierr);
470   ierr = PetscOptionsFList("-petscpartitioner_type", "Graph partitioner", "PetscPartitionerSetType", PetscPartitionerList, defaultType, name, 256, &flg);CHKERRQ(ierr);
471   if (flg) {
472     ierr = PetscPartitionerSetType(part, name);CHKERRQ(ierr);
473   } else if (!((PetscObject) part)->type_name) {
474     ierr = PetscPartitionerSetType(part, defaultType);CHKERRQ(ierr);
475   }
476   ierr = PetscOptionsEnd();CHKERRQ(ierr);
477   PetscFunctionReturn(0);
478 }
479 
480 #undef __FUNCT__
481 #define __FUNCT__ "PetscPartitionerSetFromOptions"
482 /*@
483   PetscPartitionerSetFromOptions - sets parameters in a PetscPartitioner from the options database
484 
485   Collective on PetscPartitioner
486 
487   Input Parameter:
488 . part - the PetscPartitioner object to set options for
489 
490   Level: developer
491 
492 .seealso: PetscPartitionerView()
493 @*/
494 PetscErrorCode PetscPartitionerSetFromOptions(PetscPartitioner part)
495 {
496   PetscErrorCode ierr;
497 
498   PetscFunctionBegin;
499   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
500   ierr = PetscPartitionerSetTypeFromOptions_Internal(part);CHKERRQ(ierr);
501 
502   ierr = PetscObjectOptionsBegin((PetscObject) part);CHKERRQ(ierr);
503   if (part->ops->setfromoptions) {ierr = (*part->ops->setfromoptions)(part);CHKERRQ(ierr);}
504   /* process any options handlers added with PetscObjectAddOptionsHandler() */
505   ierr = PetscObjectProcessOptionsHandlers((PetscObject) part);CHKERRQ(ierr);
506   ierr = PetscOptionsEnd();CHKERRQ(ierr);
507   ierr = PetscPartitionerViewFromOptions(part, NULL, "-petscpartitioner_view");CHKERRQ(ierr);
508   PetscFunctionReturn(0);
509 }
510 
511 #undef __FUNCT__
512 #define __FUNCT__ "PetscPartitionerSetUp"
513 /*@C
514   PetscPartitionerSetUp - Construct data structures for the PetscPartitioner
515 
516   Collective on PetscPartitioner
517 
518   Input Parameter:
519 . part - the PetscPartitioner object to setup
520 
521   Level: developer
522 
523 .seealso: PetscPartitionerView(), PetscPartitionerDestroy()
524 @*/
525 PetscErrorCode PetscPartitionerSetUp(PetscPartitioner part)
526 {
527   PetscErrorCode ierr;
528 
529   PetscFunctionBegin;
530   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
531   if (part->ops->setup) {ierr = (*part->ops->setup)(part);CHKERRQ(ierr);}
532   PetscFunctionReturn(0);
533 }
534 
535 #undef __FUNCT__
536 #define __FUNCT__ "PetscPartitionerDestroy"
537 /*@
538   PetscPartitionerDestroy - Destroys a PetscPartitioner object
539 
540   Collective on PetscPartitioner
541 
542   Input Parameter:
543 . part - the PetscPartitioner object to destroy
544 
545   Level: developer
546 
547 .seealso: PetscPartitionerView()
548 @*/
549 PetscErrorCode PetscPartitionerDestroy(PetscPartitioner *part)
550 {
551   PetscErrorCode ierr;
552 
553   PetscFunctionBegin;
554   if (!*part) PetscFunctionReturn(0);
555   PetscValidHeaderSpecific((*part), PETSCPARTITIONER_CLASSID, 1);
556 
557   if (--((PetscObject)(*part))->refct > 0) {*part = 0; PetscFunctionReturn(0);}
558   ((PetscObject) (*part))->refct = 0;
559 
560   if ((*part)->ops->destroy) {ierr = (*(*part)->ops->destroy)(*part);CHKERRQ(ierr);}
561   ierr = PetscHeaderDestroy(part);CHKERRQ(ierr);
562   PetscFunctionReturn(0);
563 }
564 
565 #undef __FUNCT__
566 #define __FUNCT__ "PetscPartitionerCreate"
567 /*@
568   PetscPartitionerCreate - Creates an empty PetscPartitioner object. The type can then be set with PetscPartitionerSetType().
569 
570   Collective on MPI_Comm
571 
572   Input Parameter:
573 . comm - The communicator for the PetscPartitioner object
574 
575   Output Parameter:
576 . part - The PetscPartitioner object
577 
578   Level: beginner
579 
580 .seealso: PetscPartitionerSetType(), PETSCPARTITIONERCHACO, PETSCPARTITIONERPARMETIS, PETSCPARTITIONERSHELL
581 @*/
582 PetscErrorCode PetscPartitionerCreate(MPI_Comm comm, PetscPartitioner *part)
583 {
584   PetscPartitioner p;
585   PetscErrorCode   ierr;
586 
587   PetscFunctionBegin;
588   PetscValidPointer(part, 2);
589   *part = NULL;
590   ierr = PetscFVInitializePackage();CHKERRQ(ierr);
591 
592   ierr = PetscHeaderCreate(p, _p_PetscPartitioner, struct _PetscPartitionerOps, PETSCPARTITIONER_CLASSID, "PetscPartitioner", "Graph Partitioner", "PetscPartitioner", comm, PetscPartitionerDestroy, PetscPartitionerView);CHKERRQ(ierr);
593   ierr = PetscMemzero(p->ops, sizeof(struct _PetscPartitionerOps));CHKERRQ(ierr);
594 
595   *part = p;
596   PetscFunctionReturn(0);
597 }
598 
599 #undef __FUNCT__
600 #define __FUNCT__ "PetscPartitionerPartition"
601 /*@
602   PetscPartitionerPartition - Create a non-overlapping partition of the cells in the mesh
603 
604   Collective on DM
605 
606   Input Parameters:
607 + part    - The PetscPartitioner
608 . dm      - The mesh DM
609 - enlarge - Expand each partition with neighbors
610 
611   Output Parameters:
612 + partSection     - The PetscSection giving the division of points by partition
613 . partition       - The list of points by partition
614 . origPartSection - If enlarge is true, the PetscSection giving the division of points before enlarging by partition, otherwise NULL
615 - origPartition   - If enlarge is true, the list of points before enlarging by partition, otherwise NULL
616 
617   Note: Instead of cells, points at a given height can be partitioned by calling PetscPartitionerSetPointHeight()
618 
619   Level: developer
620 
621 .seealso DMPlexDistribute(), PetscPartitionerSetPointHeight(), PetscPartitionerCreate()
622 @*/
623 PetscErrorCode PetscPartitionerPartition(PetscPartitioner part, DM dm, PetscBool enlarge, PetscSection partSection, IS *partition, PetscSection origPartSection, IS *origPartition)
624 {
625   PetscMPIInt    size;
626   PetscErrorCode ierr;
627 
628   PetscFunctionBegin;
629   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
630   PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
631   PetscValidHeaderSpecific(partSection, PETSC_SECTION_CLASSID, 4);
632   PetscValidPointer(partition, 5);
633   if (enlarge) {PetscValidHeaderSpecific(origPartSection, PETSC_SECTION_CLASSID, 6);}
634   if (enlarge) {PetscValidPointer(origPartition, 7);}
635   ierr = MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);CHKERRQ(ierr);
636   *origPartition = NULL;
637   if (size == 1) {
638     PetscInt *points;
639     PetscInt  cStart, cEnd, c;
640 
641     ierr = DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);CHKERRQ(ierr);
642     ierr = PetscSectionSetChart(partSection, 0, size);CHKERRQ(ierr);
643     ierr = PetscSectionSetDof(partSection, 0, cEnd-cStart);CHKERRQ(ierr);
644     ierr = PetscSectionSetUp(partSection);CHKERRQ(ierr);
645     ierr = PetscMalloc1(cEnd-cStart, &points);CHKERRQ(ierr);
646     for (c = cStart; c < cEnd; ++c) points[c] = c;
647     ierr = ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
648     PetscFunctionReturn(0);
649   }
650   if (part->height == 0) {
651     PetscInt  numVertices;
652     PetscInt *start     = NULL;
653     PetscInt *adjacency = NULL;
654 
655     ierr = DMPlexCreateNeighborCSR(dm, 0, &numVertices, &start, &adjacency);CHKERRQ(ierr);
656     if (!part->ops->partition) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_WRONGSTATE, "PetscPartitioner has no type");
657     ierr = (*part->ops->partition)(part, dm, size, numVertices, start, adjacency, partSection, partition);CHKERRQ(ierr);
658     if (enlarge) {
659       origPartSection = partSection;
660       *origPartition  = *partition;
661 
662       ierr = DMPlexEnlargePartition(dm, start, adjacency, origPartSection, *origPartition, partSection, partition);CHKERRQ(ierr);
663     }
664     ierr = PetscFree(start);CHKERRQ(ierr);
665     ierr = PetscFree(adjacency);CHKERRQ(ierr);
666   } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
667   PetscFunctionReturn(0);
668 }
669 
670 #undef __FUNCT__
671 #define __FUNCT__ "PetscPartitionerDestroy_Shell"
672 PetscErrorCode PetscPartitionerDestroy_Shell(PetscPartitioner part)
673 {
674   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
675   PetscErrorCode          ierr;
676 
677   PetscFunctionBegin;
678   ierr = PetscSectionDestroy(&p->section);CHKERRQ(ierr);
679   ierr = ISDestroy(&p->partition);CHKERRQ(ierr);
680   ierr = PetscFree(p);CHKERRQ(ierr);
681   PetscFunctionReturn(0);
682 }
683 
684 #undef __FUNCT__
685 #define __FUNCT__ "PetscPartitionerView_Shell_Ascii"
686 PetscErrorCode PetscPartitionerView_Shell_Ascii(PetscPartitioner part, PetscViewer viewer)
687 {
688   PetscViewerFormat format;
689   PetscErrorCode    ierr;
690 
691   PetscFunctionBegin;
692   ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
693   ierr = PetscViewerASCIIPrintf(viewer, "Shell Graph Partitioner:\n");CHKERRQ(ierr);
694   PetscFunctionReturn(0);
695 }
696 
697 #undef __FUNCT__
698 #define __FUNCT__ "PetscPartitionerView_Shell"
699 PetscErrorCode PetscPartitionerView_Shell(PetscPartitioner part, PetscViewer viewer)
700 {
701   PetscBool      iascii;
702   PetscErrorCode ierr;
703 
704   PetscFunctionBegin;
705   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
706   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
707   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
708   if (iascii) {ierr = PetscPartitionerView_Shell_Ascii(part, viewer);CHKERRQ(ierr);}
709   PetscFunctionReturn(0);
710 }
711 
712 #undef __FUNCT__
713 #define __FUNCT__ "PetscPartitionerPartition_Shell"
714 PetscErrorCode PetscPartitionerPartition_Shell(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
715 {
716   PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
717   PetscInt                np;
718   PetscErrorCode          ierr;
719 
720   PetscFunctionBegin;
721   if (!numVertices) {
722     ierr = PetscSectionSetChart(partSection, 0, nparts);CHKERRQ(ierr);
723     ierr = PetscSectionSetUp(partSection);CHKERRQ(ierr);
724     ierr = ISCreateGeneral(PetscObjectComm((PetscObject) part), 0, NULL, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
725     PetscFunctionReturn(0);
726   }
727   ierr = PetscSectionGetChart(p->section, NULL, &np);CHKERRQ(ierr);
728   if (nparts != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of requested partitions %d != configured partitions %d", nparts, np);
729   ierr = ISGetLocalSize(p->partition, &np);CHKERRQ(ierr);
730   if (numVertices != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input vertices %d != configured vertices %d", numVertices, np);
731   partSection = p->section;
732   ierr = PetscObjectReference((PetscObject) p->section);CHKERRQ(ierr);
733   *partition  = p->partition;
734   ierr = PetscObjectReference((PetscObject) p->partition);CHKERRQ(ierr);
735   PetscFunctionReturn(0);
736 }
737 
738 #undef __FUNCT__
739 #define __FUNCT__ "PetscPartitionerInitialize_Shell"
740 PetscErrorCode PetscPartitionerInitialize_Shell(PetscPartitioner part)
741 {
742   PetscFunctionBegin;
743   part->ops->view      = PetscPartitionerView_Shell;
744   part->ops->destroy   = PetscPartitionerDestroy_Shell;
745   part->ops->partition = PetscPartitionerPartition_Shell;
746   PetscFunctionReturn(0);
747 }
748 
749 /*MC
750   PETSCPARTITIONERSHELL = "shell" - A PetscPartitioner object
751 
752   Level: intermediate
753 
754 .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
755 M*/
756 
757 #undef __FUNCT__
758 #define __FUNCT__ "PetscPartitionerCreate_Shell"
759 PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Shell(PetscPartitioner part)
760 {
761   PetscPartitioner_Shell *p;
762   PetscErrorCode          ierr;
763 
764   PetscFunctionBegin;
765   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
766   ierr       = PetscNewLog(part, &p);CHKERRQ(ierr);
767   part->data = p;
768 
769   ierr = PetscPartitionerInitialize_Shell(part);CHKERRQ(ierr);
770   PetscFunctionReturn(0);
771 }
772 
773 #undef __FUNCT__
774 #define __FUNCT__ "PetscPartitionerDestroy_Chaco"
775 PetscErrorCode PetscPartitionerDestroy_Chaco(PetscPartitioner part)
776 {
777   PetscPartitioner_Chaco *p = (PetscPartitioner_Chaco *) part->data;
778   PetscErrorCode          ierr;
779 
780   PetscFunctionBegin;
781   ierr = PetscFree(p);CHKERRQ(ierr);
782   PetscFunctionReturn(0);
783 }
784 
785 #undef __FUNCT__
786 #define __FUNCT__ "PetscPartitionerView_Chaco_Ascii"
787 PetscErrorCode PetscPartitionerView_Chaco_Ascii(PetscPartitioner part, PetscViewer viewer)
788 {
789   PetscViewerFormat format;
790   PetscErrorCode    ierr;
791 
792   PetscFunctionBegin;
793   ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
794   ierr = PetscViewerASCIIPrintf(viewer, "Chaco Graph Partitioner:\n");CHKERRQ(ierr);
795   PetscFunctionReturn(0);
796 }
797 
798 #undef __FUNCT__
799 #define __FUNCT__ "PetscPartitionerView_Chaco"
800 PetscErrorCode PetscPartitionerView_Chaco(PetscPartitioner part, PetscViewer viewer)
801 {
802   PetscBool      iascii;
803   PetscErrorCode ierr;
804 
805   PetscFunctionBegin;
806   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
807   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
808   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
809   if (iascii) {ierr = PetscPartitionerView_Chaco_Ascii(part, viewer);CHKERRQ(ierr);}
810   PetscFunctionReturn(0);
811 }
812 
813 #if defined(PETSC_HAVE_CHACO)
814 #if defined(PETSC_HAVE_UNISTD_H)
815 #include <unistd.h>
816 #endif
817 /* Chaco does not have an include file */
818 PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
819                        float *ewgts, float *x, float *y, float *z, char *outassignname,
820                        char *outfilename, short *assignment, int architecture, int ndims_tot,
821                        int mesh_dims[3], double *goal, int global_method, int local_method,
822                        int rqi_flag, int vmax, int ndims, double eigtol, long seed);
823 
824 extern int FREE_GRAPH;
825 #endif
826 
827 #undef __FUNCT__
828 #define __FUNCT__ "PetscPartitionerPartition_Chaco"
829 PetscErrorCode PetscPartitionerPartition_Chaco(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
830 {
831 #if defined(PETSC_HAVE_CHACO)
832   enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
833   MPI_Comm       comm;
834   int            nvtxs          = numVertices; /* number of vertices in full graph */
835   int           *vwgts          = NULL;   /* weights for all vertices */
836   float         *ewgts          = NULL;   /* weights for all edges */
837   float         *x              = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
838   char          *outassignname  = NULL;   /*  name of assignment output file */
839   char          *outfilename    = NULL;   /* output file name */
840   int            architecture   = 1;      /* 0 => hypercube, d => d-dimensional mesh */
841   int            ndims_tot      = 0;      /* total number of cube dimensions to divide */
842   int            mesh_dims[3];            /* dimensions of mesh of processors */
843   double        *goal          = NULL;    /* desired set sizes for each set */
844   int            global_method = 1;       /* global partitioning algorithm */
845   int            local_method  = 1;       /* local partitioning algorithm */
846   int            rqi_flag      = 0;       /* should I use RQI/Symmlq eigensolver? */
847   int            vmax          = 200;     /* how many vertices to coarsen down to? */
848   int            ndims         = 1;       /* number of eigenvectors (2^d sets) */
849   double         eigtol        = 0.001;   /* tolerance on eigenvectors */
850   long           seed          = 123636512; /* for random graph mutations */
851   short int     *assignment;              /* Output partition */
852   int            fd_stdout, fd_pipe[2];
853   PetscInt      *points;
854   int            i, v, p;
855   PetscErrorCode ierr;
856 
857   PetscFunctionBegin;
858   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
859   if (!numVertices) {
860     ierr = PetscSectionSetChart(partSection, 0, nparts);CHKERRQ(ierr);
861     ierr = PetscSectionSetUp(partSection);CHKERRQ(ierr);
862     ierr = ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
863     PetscFunctionReturn(0);
864   }
865   FREE_GRAPH = 0;                         /* Do not let Chaco free my memory */
866   for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];
867 
868   if (global_method == INERTIAL_METHOD) {
869     /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
870     SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
871   }
872   mesh_dims[0] = nparts;
873   mesh_dims[1] = 1;
874   mesh_dims[2] = 1;
875   ierr = PetscMalloc1(nvtxs, &assignment);CHKERRQ(ierr);
876   /* Chaco outputs to stdout. We redirect this to a buffer. */
877   /* TODO: check error codes for UNIX calls */
878 #if defined(PETSC_HAVE_UNISTD_H)
879   {
880     int piperet;
881     piperet = pipe(fd_pipe);
882     if (piperet) SETERRQ(comm,PETSC_ERR_SYS,"Could not create pipe");
883     fd_stdout = dup(1);
884     close(1);
885     dup2(fd_pipe[1], 1);
886   }
887 #endif
888   ierr = interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
889                    assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
890                    vmax, ndims, eigtol, seed);
891 #if defined(PETSC_HAVE_UNISTD_H)
892   {
893     char msgLog[10000];
894     int  count;
895 
896     fflush(stdout);
897     count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
898     if (count < 0) count = 0;
899     msgLog[count] = 0;
900     close(1);
901     dup2(fd_stdout, 1);
902     close(fd_stdout);
903     close(fd_pipe[0]);
904     close(fd_pipe[1]);
905     if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
906   }
907 #endif
908   /* Convert to PetscSection+IS */
909   ierr = PetscSectionSetChart(partSection, 0, nparts);CHKERRQ(ierr);
910   for (v = 0; v < nvtxs; ++v) {
911     ierr = PetscSectionAddDof(partSection, assignment[v], 1);CHKERRQ(ierr);
912   }
913   ierr = PetscSectionSetUp(partSection);CHKERRQ(ierr);
914   ierr = PetscMalloc1(nvtxs, &points);CHKERRQ(ierr);
915   for (p = 0, i = 0; p < nparts; ++p) {
916     for (v = 0; v < nvtxs; ++v) {
917       if (assignment[v] == p) points[i++] = v;
918     }
919   }
920   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
921   ierr = ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
922   if (global_method == INERTIAL_METHOD) {
923     /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
924   }
925   ierr = PetscFree(assignment);CHKERRQ(ierr);
926   for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
927   PetscFunctionReturn(0);
928 #else
929   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco.");
930 #endif
931 }
932 
933 #undef __FUNCT__
934 #define __FUNCT__ "PetscPartitionerInitialize_Chaco"
935 PetscErrorCode PetscPartitionerInitialize_Chaco(PetscPartitioner part)
936 {
937   PetscFunctionBegin;
938   part->ops->view      = PetscPartitionerView_Chaco;
939   part->ops->destroy   = PetscPartitionerDestroy_Chaco;
940   part->ops->partition = PetscPartitionerPartition_Chaco;
941   PetscFunctionReturn(0);
942 }
943 
944 /*MC
945   PETSCPARTITIONERCHACO = "chaco" - A PetscPartitioner object using the Chaco library
946 
947   Level: intermediate
948 
949 .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
950 M*/
951 
952 #undef __FUNCT__
953 #define __FUNCT__ "PetscPartitionerCreate_Chaco"
954 PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Chaco(PetscPartitioner part)
955 {
956   PetscPartitioner_Chaco *p;
957   PetscErrorCode          ierr;
958 
959   PetscFunctionBegin;
960   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
961   ierr       = PetscNewLog(part, &p);CHKERRQ(ierr);
962   part->data = p;
963 
964   ierr = PetscPartitionerInitialize_Chaco(part);CHKERRQ(ierr);
965   ierr = PetscCitationsRegister(ChacoPartitionerCitation, &ChacoPartitionercite);CHKERRQ(ierr);
966   PetscFunctionReturn(0);
967 }
968 
969 #undef __FUNCT__
970 #define __FUNCT__ "PetscPartitionerDestroy_ParMetis"
971 PetscErrorCode PetscPartitionerDestroy_ParMetis(PetscPartitioner part)
972 {
973   PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
974   PetscErrorCode             ierr;
975 
976   PetscFunctionBegin;
977   ierr = PetscFree(p);CHKERRQ(ierr);
978   PetscFunctionReturn(0);
979 }
980 
981 #undef __FUNCT__
982 #define __FUNCT__ "PetscPartitionerView_ParMetis_Ascii"
983 PetscErrorCode PetscPartitionerView_ParMetis_Ascii(PetscPartitioner part, PetscViewer viewer)
984 {
985   PetscViewerFormat format;
986   PetscErrorCode    ierr;
987 
988   PetscFunctionBegin;
989   ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
990   ierr = PetscViewerASCIIPrintf(viewer, "ParMetis Graph Partitioner:\n");CHKERRQ(ierr);
991   PetscFunctionReturn(0);
992 }
993 
994 #undef __FUNCT__
995 #define __FUNCT__ "PetscPartitionerView_ParMetis"
996 PetscErrorCode PetscPartitionerView_ParMetis(PetscPartitioner part, PetscViewer viewer)
997 {
998   PetscBool      iascii;
999   PetscErrorCode ierr;
1000 
1001   PetscFunctionBegin;
1002   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
1003   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
1004   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
1005   if (iascii) {ierr = PetscPartitionerView_ParMetis_Ascii(part, viewer);CHKERRQ(ierr);}
1006   PetscFunctionReturn(0);
1007 }
1008 
1009 #if defined(PETSC_HAVE_PARMETIS)
1010 #include <parmetis.h>
1011 #endif
1012 
1013 #undef __FUNCT__
1014 #define __FUNCT__ "PetscPartitionerPartition_ParMetis"
1015 PetscErrorCode PetscPartitionerPartition_ParMetis(PetscPartitioner part, DM dm, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection partSection, IS *partition)
1016 {
1017 #if defined(PETSC_HAVE_PARMETIS)
1018   MPI_Comm       comm;
1019   PetscInt       nvtxs      = numVertices; /* The number of vertices in full graph */
1020   PetscInt      *vtxdist;                  /* Distribution of vertices across processes */
1021   PetscInt      *xadj       = start;       /* Start of edge list for each vertex */
1022   PetscInt      *adjncy     = adjacency;   /* Edge lists for all vertices */
1023   PetscInt      *vwgt       = NULL;        /* Vertex weights */
1024   PetscInt      *adjwgt     = NULL;        /* Edge weights */
1025   PetscInt       wgtflag    = 0;           /* Indicates which weights are present */
1026   PetscInt       numflag    = 0;           /* Indicates initial offset (0 or 1) */
1027   PetscInt       ncon       = 1;           /* The number of weights per vertex */
1028   PetscReal     *tpwgts;                   /* The fraction of vertex weights assigned to each partition */
1029   PetscReal     *ubvec;                    /* The balance intolerance for vertex weights */
1030   PetscInt       options[5];               /* Options */
1031   /* Outputs */
1032   PetscInt       edgeCut;                  /* The number of edges cut by the partition */
1033   PetscInt      *assignment, *points;
1034   PetscMPIInt    rank, p, v, i;
1035   PetscErrorCode ierr;
1036 
1037   PetscFunctionBegin;
1038   ierr = PetscObjectGetComm((PetscObject) part, &comm);CHKERRQ(ierr);
1039   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1040   options[0] = 0; /* Use all defaults */
1041   /* Calculate vertex distribution */
1042   ierr = PetscMalloc4(nparts+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment);CHKERRQ(ierr);
1043   vtxdist[0] = 0;
1044   ierr = MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);CHKERRQ(ierr);
1045   for (p = 2; p <= nparts; ++p) {
1046     vtxdist[p] += vtxdist[p-1];
1047   }
1048   /* Calculate weights */
1049   for (p = 0; p < nparts; ++p) {
1050     tpwgts[p] = 1.0/nparts;
1051   }
1052   ubvec[0] = 1.05;
1053 
1054   if (nparts == 1) {
1055     ierr = PetscMemzero(assignment, nvtxs * sizeof(PetscInt));
1056   } else {
1057     if (vtxdist[1] == vtxdist[nparts]) {
1058       if (!rank) {
1059         PetscStackPush("METIS_PartGraphKway");
1060         ierr = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, NULL, &edgeCut, assignment);
1061         PetscStackPop;
1062         if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
1063       }
1064     } else {
1065       PetscStackPush("ParMETIS_V3_PartKway");
1066       ierr = ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm);
1067       PetscStackPop;
1068       if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
1069     }
1070   }
1071   /* Convert to PetscSection+IS */
1072   ierr = PetscSectionSetChart(partSection, 0, nparts);CHKERRQ(ierr);
1073   for (v = 0; v < nvtxs; ++v) {ierr = PetscSectionAddDof(partSection, assignment[v], 1);CHKERRQ(ierr);}
1074   ierr = PetscSectionSetUp(partSection);CHKERRQ(ierr);
1075   ierr = PetscMalloc1(nvtxs, &points);CHKERRQ(ierr);
1076   for (p = 0, i = 0; p < nparts; ++p) {
1077     for (v = 0; v < nvtxs; ++v) {
1078       if (assignment[v] == p) points[i++] = v;
1079     }
1080   }
1081   if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1082   ierr = ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
1083   ierr = PetscFree4(vtxdist,tpwgts,ubvec,assignment);CHKERRQ(ierr);
1084 #else
1085   SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
1086 #endif
1087   PetscFunctionReturn(0);
1088 }
1089 
1090 #undef __FUNCT__
1091 #define __FUNCT__ "PetscPartitionerInitialize_ParMetis"
1092 PetscErrorCode PetscPartitionerInitialize_ParMetis(PetscPartitioner part)
1093 {
1094   PetscFunctionBegin;
1095   part->ops->view      = PetscPartitionerView_ParMetis;
1096   part->ops->destroy   = PetscPartitionerDestroy_ParMetis;
1097   part->ops->partition = PetscPartitionerPartition_ParMetis;
1098   PetscFunctionReturn(0);
1099 }
1100 
1101 /*MC
1102   PETSCPARTITIONERPARMETIS = "parmetis" - A PetscPartitioner object using the ParMetis library
1103 
1104   Level: intermediate
1105 
1106 .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1107 M*/
1108 
1109 #undef __FUNCT__
1110 #define __FUNCT__ "PetscPartitionerCreate_ParMetis"
1111 PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMetis(PetscPartitioner part)
1112 {
1113   PetscPartitioner_ParMetis *p;
1114   PetscErrorCode          ierr;
1115 
1116   PetscFunctionBegin;
1117   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
1118   ierr       = PetscNewLog(part, &p);CHKERRQ(ierr);
1119   part->data = p;
1120 
1121   ierr = PetscPartitionerInitialize_ParMetis(part);CHKERRQ(ierr);
1122   ierr = PetscCitationsRegister(ParMetisPartitionerCitation, &ParMetisPartitionercite);CHKERRQ(ierr);
1123   PetscFunctionReturn(0);
1124 }
1125 
1126 #undef __FUNCT__
1127 #define __FUNCT__ "DMPlexMarkTreeClosure"
1128 static PetscErrorCode DMPlexMarkTreeClosure(DM dm, PetscSegBuffer segpart, PetscBT bt, PetscInt point, PetscInt *partSize)
1129 {
1130   PetscInt       parent, closureSize, *closure = NULL, i, pStart, pEnd;
1131   PetscErrorCode ierr;
1132 
1133   PetscFunctionBegin;
1134   ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr);
1135   if (parent == point) PetscFunctionReturn(0);
1136   ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
1137   ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1138   for (i = 0; i < closureSize; i++) {
1139     PetscInt cpoint = closure[2*i];
1140 
1141     if (!PetscBTLookupSet(bt,cpoint-pStart)) {
1142       PetscInt *PETSC_RESTRICT pt;
1143       (*partSize)++;
1144       ierr = PetscSegBufferGetInts(segpart,1,&pt);CHKERRQ(ierr);
1145       *pt = cpoint;
1146     }
1147     ierr = DMPlexMarkTreeClosure(dm,segpart,bt,cpoint,partSize);CHKERRQ(ierr);
1148   }
1149   ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
1150   PetscFunctionReturn(0);
1151 }
1152 
1153 #undef __FUNCT__
1154 #define __FUNCT__ "DMPlexCreatePartitionClosure"
1155 PetscErrorCode DMPlexCreatePartitionClosure(DM dm, PetscSection pointSection, IS pointPartition, PetscSection *section, IS *partition)
1156 {
1157   /* const PetscInt  height = 0; */
1158   const PetscInt *partArray;
1159   PetscInt       *allPoints, *packPoints;
1160   PetscInt        rStart, rEnd, rank, pStart, pEnd, newSize;
1161   PetscErrorCode  ierr;
1162   PetscBT         bt;
1163   PetscSegBuffer  segpack,segpart;
1164 
1165   PetscFunctionBegin;
1166   ierr = PetscSectionGetChart(pointSection, &rStart, &rEnd);CHKERRQ(ierr);
1167   ierr = ISGetIndices(pointPartition, &partArray);CHKERRQ(ierr);
1168   ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);CHKERRQ(ierr);
1169   ierr = PetscSectionSetChart(*section, rStart, rEnd);CHKERRQ(ierr);
1170   ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
1171   ierr = PetscBTCreate(pEnd-pStart,&bt);CHKERRQ(ierr);
1172   ierr = PetscSegBufferCreate(sizeof(PetscInt),1000,&segpack);CHKERRQ(ierr);
1173   ierr = PetscSegBufferCreate(sizeof(PetscInt),1000,&segpart);CHKERRQ(ierr);
1174   for (rank = rStart; rank < rEnd; ++rank) {
1175     PetscInt partSize = 0, numPoints, offset, p, *PETSC_RESTRICT placePoints;
1176 
1177     ierr = PetscSectionGetDof(pointSection, rank, &numPoints);CHKERRQ(ierr);
1178     ierr = PetscSectionGetOffset(pointSection, rank, &offset);CHKERRQ(ierr);
1179     for (p = 0; p < numPoints; ++p) {
1180       PetscInt  point   = partArray[offset+p], closureSize, c;
1181       PetscInt *closure = NULL;
1182 
1183       /* TODO Include support for height > 0 case */
1184       ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1185       for (c=0; c<closureSize; c++) {
1186         PetscInt cpoint = closure[c*2];
1187         if (!PetscBTLookupSet(bt,cpoint-pStart)) {
1188           PetscInt *PETSC_RESTRICT pt;
1189           partSize++;
1190           ierr = PetscSegBufferGetInts(segpart,1,&pt);CHKERRQ(ierr);
1191           *pt = cpoint;
1192         }
1193         ierr = DMPlexMarkTreeClosure(dm,segpart,bt,cpoint,&partSize);CHKERRQ(ierr);
1194       }
1195       ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1196     }
1197     ierr = PetscSectionSetDof(*section, rank, partSize);CHKERRQ(ierr);
1198     ierr = PetscSegBufferGetInts(segpack,partSize,&placePoints);CHKERRQ(ierr);
1199     ierr = PetscSegBufferExtractTo(segpart,placePoints);CHKERRQ(ierr);
1200     ierr = PetscSortInt(partSize,placePoints);CHKERRQ(ierr);
1201     for (p=0; p<partSize; p++) {ierr = PetscBTClear(bt,placePoints[p]-pStart);CHKERRQ(ierr);}
1202   }
1203   ierr = PetscBTDestroy(&bt);CHKERRQ(ierr);
1204   ierr = PetscSegBufferDestroy(&segpart);CHKERRQ(ierr);
1205 
1206   ierr = PetscSectionSetUp(*section);CHKERRQ(ierr);
1207   ierr = PetscSectionGetStorageSize(*section, &newSize);CHKERRQ(ierr);
1208   ierr = PetscMalloc1(newSize, &allPoints);CHKERRQ(ierr);
1209 
1210   ierr = PetscSegBufferExtractInPlace(segpack,&packPoints);CHKERRQ(ierr);
1211   for (rank = rStart; rank < rEnd; ++rank) {
1212     PetscInt numPoints, offset;
1213 
1214     ierr = PetscSectionGetDof(*section, rank, &numPoints);CHKERRQ(ierr);
1215     ierr = PetscSectionGetOffset(*section, rank, &offset);CHKERRQ(ierr);
1216     ierr = PetscMemcpy(&allPoints[offset], packPoints, numPoints * sizeof(PetscInt));CHKERRQ(ierr);
1217     packPoints += numPoints;
1218   }
1219 
1220   ierr = PetscSegBufferDestroy(&segpack);CHKERRQ(ierr);
1221   ierr = ISRestoreIndices(pointPartition, &partArray);CHKERRQ(ierr);
1222   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), newSize, allPoints, PETSC_OWN_POINTER, partition);CHKERRQ(ierr);
1223   PetscFunctionReturn(0);
1224 }
1225