Actual source code: plexpartition.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/hashseti.h>
4: PetscClassId PETSCPARTITIONER_CLASSID = 0;
6: PetscFunctionList PetscPartitionerList = NULL;
7: PetscBool PetscPartitionerRegisterAllCalled = PETSC_FALSE;
9: PetscBool ChacoPartitionercite = PETSC_FALSE;
10: const char ChacoPartitionerCitation[] = "@inproceedings{Chaco95,\n"
11: " author = {Bruce Hendrickson and Robert Leland},\n"
12: " title = {A multilevel algorithm for partitioning graphs},\n"
13: " booktitle = {Supercomputing '95: Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM)},"
14: " isbn = {0-89791-816-9},\n"
15: " pages = {28},\n"
16: " doi = {https://doi.acm.org/10.1145/224170.224228},\n"
17: " publisher = {ACM Press},\n"
18: " address = {New York},\n"
19: " year = {1995}\n}\n";
21: PetscBool ParMetisPartitionercite = PETSC_FALSE;
22: const char ParMetisPartitionerCitation[] = "@article{KarypisKumar98,\n"
23: " author = {George Karypis and Vipin Kumar},\n"
24: " title = {A Parallel Algorithm for Multilevel Graph Partitioning and Sparse Matrix Ordering},\n"
25: " journal = {Journal of Parallel and Distributed Computing},\n"
26: " volume = {48},\n"
27: " pages = {71--85},\n"
28: " year = {1998}\n}\n";
30: PetscBool PTScotchPartitionercite = PETSC_FALSE;
31: const char PTScotchPartitionerCitation[] =
32: "@article{PTSCOTCH,\n"
33: " author = {C. Chevalier and F. Pellegrini},\n"
34: " title = {{PT-SCOTCH}: a tool for efficient parallel graph ordering},\n"
35: " journal = {Parallel Computing},\n"
36: " volume = {34},\n"
37: " number = {6},\n"
38: " pages = {318--331},\n"
39: " year = {2008},\n"
40: " doi = {https://doi.org/10.1016/j.parco.2007.12.001}\n"
41: "}\n";
44: PETSC_STATIC_INLINE PetscInt DMPlex_GlobalID(PetscInt point) { return point >= 0 ? point : -(point+1); }
46: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
47: {
48: PetscInt dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
49: PetscInt *adj = NULL, *vOffsets = NULL, *graph = NULL;
50: IS cellNumbering;
51: const PetscInt *cellNum;
52: PetscBool useCone, useClosure;
53: PetscSection section;
54: PetscSegBuffer adjBuffer;
55: PetscSF sfPoint;
56: PetscInt *adjCells = NULL, *remoteCells = NULL;
57: const PetscInt *local;
58: PetscInt nroots, nleaves, l;
59: PetscMPIInt rank;
63: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
64: DMGetDimension(dm, &dim);
65: DMPlexGetDepth(dm, &depth);
66: if (dim != depth) {
67: /* We do not handle the uninterpolated case here */
68: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
69: /* DMPlexCreateNeighborCSR does not make a numbering */
70: if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
71: /* Different behavior for empty graphs */
72: if (!*numVertices) {
73: PetscMalloc1(1, offsets);
74: (*offsets)[0] = 0;
75: }
76: /* Broken in parallel */
77: if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
78: return(0);
79: }
80: DMGetPointSF(dm, &sfPoint);
81: DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
82: /* Build adjacency graph via a section/segbuffer */
83: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
84: PetscSectionSetChart(section, pStart, pEnd);
85: PetscSegBufferCreate(sizeof(PetscInt),1000,&adjBuffer);
86: /* Always use FVM adjacency to create partitioner graph */
87: DMGetBasicAdjacency(dm, &useCone, &useClosure);
88: DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
89: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
90: if (globalNumbering) {
91: PetscObjectReference((PetscObject)cellNumbering);
92: *globalNumbering = cellNumbering;
93: }
94: ISGetIndices(cellNumbering, &cellNum);
95: /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
96: Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
97: */
98: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
99: if (nroots >= 0) {
100: PetscInt fStart, fEnd, f;
102: PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
103: DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
104: for (l = 0; l < nroots; ++l) adjCells[l] = -3;
105: for (f = fStart; f < fEnd; ++f) {
106: const PetscInt *support;
107: PetscInt supportSize;
109: DMPlexGetSupport(dm, f, &support);
110: DMPlexGetSupportSize(dm, f, &supportSize);
111: if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
112: else if (supportSize == 2) {
113: PetscFindInt(support[0], nleaves, local, &p);
114: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1]]);
115: PetscFindInt(support[1], nleaves, local, &p);
116: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
117: }
118: /* Handle non-conforming meshes */
119: if (supportSize > 2) {
120: PetscInt numChildren, i;
121: const PetscInt *children;
123: DMPlexGetTreeChildren(dm, f, &numChildren, &children);
124: for (i = 0; i < numChildren; ++i) {
125: const PetscInt child = children[i];
126: if (fStart <= child && child < fEnd) {
127: DMPlexGetSupport(dm, child, &support);
128: DMPlexGetSupportSize(dm, child, &supportSize);
129: if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
130: else if (supportSize == 2) {
131: PetscFindInt(support[0], nleaves, local, &p);
132: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1]]);
133: PetscFindInt(support[1], nleaves, local, &p);
134: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
135: }
136: }
137: }
138: }
139: }
140: for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
141: PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells);
142: PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells);
143: PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
144: PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
145: }
146: /* Combine local and global adjacencies */
147: for (*numVertices = 0, p = pStart; p < pEnd; p++) {
148: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
149: if (nroots > 0) {if (cellNum[p] < 0) continue;}
150: /* Add remote cells */
151: if (remoteCells) {
152: const PetscInt gp = DMPlex_GlobalID(cellNum[p]);
153: PetscInt coneSize, numChildren, c, i;
154: const PetscInt *cone, *children;
156: DMPlexGetCone(dm, p, &cone);
157: DMPlexGetConeSize(dm, p, &coneSize);
158: for (c = 0; c < coneSize; ++c) {
159: const PetscInt point = cone[c];
160: if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
161: PetscInt *PETSC_RESTRICT pBuf;
162: PetscSectionAddDof(section, p, 1);
163: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
164: *pBuf = remoteCells[point];
165: }
166: /* Handle non-conforming meshes */
167: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
168: for (i = 0; i < numChildren; ++i) {
169: const PetscInt child = children[i];
170: if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
171: PetscInt *PETSC_RESTRICT pBuf;
172: PetscSectionAddDof(section, p, 1);
173: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
174: *pBuf = remoteCells[child];
175: }
176: }
177: }
178: }
179: /* Add local cells */
180: adjSize = PETSC_DETERMINE;
181: DMPlexGetAdjacency(dm, p, &adjSize, &adj);
182: for (a = 0; a < adjSize; ++a) {
183: const PetscInt point = adj[a];
184: if (point != p && pStart <= point && point < pEnd) {
185: PetscInt *PETSC_RESTRICT pBuf;
186: PetscSectionAddDof(section, p, 1);
187: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
188: *pBuf = DMPlex_GlobalID(cellNum[point]);
189: }
190: }
191: (*numVertices)++;
192: }
193: PetscFree(adj);
194: PetscFree2(adjCells, remoteCells);
195: DMSetBasicAdjacency(dm, useCone, useClosure);
197: /* Derive CSR graph from section/segbuffer */
198: PetscSectionSetUp(section);
199: PetscSectionGetStorageSize(section, &size);
200: PetscMalloc1(*numVertices+1, &vOffsets);
201: for (idx = 0, p = pStart; p < pEnd; p++) {
202: if (nroots > 0) {if (cellNum[p] < 0) continue;}
203: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
204: }
205: vOffsets[*numVertices] = size;
206: PetscSegBufferExtractAlloc(adjBuffer, &graph);
208: if (nroots >= 0) {
209: /* Filter out duplicate edges using section/segbuffer */
210: PetscSectionReset(section);
211: PetscSectionSetChart(section, 0, *numVertices);
212: for (p = 0; p < *numVertices; p++) {
213: PetscInt start = vOffsets[p], end = vOffsets[p+1];
214: PetscInt numEdges = end-start, *PETSC_RESTRICT edges;
215: PetscSortRemoveDupsInt(&numEdges, &graph[start]);
216: PetscSectionSetDof(section, p, numEdges);
217: PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
218: PetscArraycpy(edges, &graph[start], numEdges);
219: }
220: PetscFree(vOffsets);
221: PetscFree(graph);
222: /* Derive CSR graph from section/segbuffer */
223: PetscSectionSetUp(section);
224: PetscSectionGetStorageSize(section, &size);
225: PetscMalloc1(*numVertices+1, &vOffsets);
226: for (idx = 0, p = 0; p < *numVertices; p++) {
227: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
228: }
229: vOffsets[*numVertices] = size;
230: PetscSegBufferExtractAlloc(adjBuffer, &graph);
231: } else {
232: /* Sort adjacencies (not strictly necessary) */
233: for (p = 0; p < *numVertices; p++) {
234: PetscInt start = vOffsets[p], end = vOffsets[p+1];
235: PetscSortInt(end-start, &graph[start]);
236: }
237: }
239: if (offsets) *offsets = vOffsets;
240: if (adjacency) *adjacency = graph;
242: /* Cleanup */
243: PetscSegBufferDestroy(&adjBuffer);
244: PetscSectionDestroy(§ion);
245: ISRestoreIndices(cellNumbering, &cellNum);
246: ISDestroy(&cellNumbering);
247: return(0);
248: }
250: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
251: {
252: Mat conn, CSR;
253: IS fis, cis, cis_own;
254: PetscSF sfPoint;
255: const PetscInt *rows, *cols, *ii, *jj;
256: PetscInt *idxs,*idxs2;
257: PetscInt dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
258: PetscMPIInt rank;
259: PetscBool flg;
263: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
264: DMGetDimension(dm, &dim);
265: DMPlexGetDepth(dm, &depth);
266: if (dim != depth) {
267: /* We do not handle the uninterpolated case here */
268: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
269: /* DMPlexCreateNeighborCSR does not make a numbering */
270: if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
271: /* Different behavior for empty graphs */
272: if (!*numVertices) {
273: PetscMalloc1(1, offsets);
274: (*offsets)[0] = 0;
275: }
276: /* Broken in parallel */
277: if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
278: return(0);
279: }
280: /* Interpolated and parallel case */
282: /* numbering */
283: DMGetPointSF(dm, &sfPoint);
284: DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
285: DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
286: DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
287: DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
288: if (globalNumbering) {
289: ISDuplicate(cis, globalNumbering);
290: }
292: /* get positive global ids and local sizes for facets and cells */
293: ISGetLocalSize(fis, &m);
294: ISGetIndices(fis, &rows);
295: PetscMalloc1(m, &idxs);
296: for (i = 0, floc = 0; i < m; i++) {
297: const PetscInt p = rows[i];
299: if (p < 0) {
300: idxs[i] = -(p+1);
301: } else {
302: idxs[i] = p;
303: floc += 1;
304: }
305: }
306: ISRestoreIndices(fis, &rows);
307: ISDestroy(&fis);
308: ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);
310: ISGetLocalSize(cis, &m);
311: ISGetIndices(cis, &cols);
312: PetscMalloc1(m, &idxs);
313: PetscMalloc1(m, &idxs2);
314: for (i = 0, cloc = 0; i < m; i++) {
315: const PetscInt p = cols[i];
317: if (p < 0) {
318: idxs[i] = -(p+1);
319: } else {
320: idxs[i] = p;
321: idxs2[cloc++] = p;
322: }
323: }
324: ISRestoreIndices(cis, &cols);
325: ISDestroy(&cis);
326: ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
327: ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);
329: /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
330: MatCreate(PetscObjectComm((PetscObject)dm), &conn);
331: MatSetSizes(conn, floc, cloc, M, N);
332: MatSetType(conn, MATMPIAIJ);
333: DMPlexGetMaxSizes(dm, NULL, &lm);
334: MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dm));
335: MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);
337: /* Assemble matrix */
338: ISGetIndices(fis, &rows);
339: ISGetIndices(cis, &cols);
340: for (c = cStart; c < cEnd; c++) {
341: const PetscInt *cone;
342: PetscInt coneSize, row, col, f;
344: col = cols[c-cStart];
345: DMPlexGetCone(dm, c, &cone);
346: DMPlexGetConeSize(dm, c, &coneSize);
347: for (f = 0; f < coneSize; f++) {
348: const PetscScalar v = 1.0;
349: const PetscInt *children;
350: PetscInt numChildren, ch;
352: row = rows[cone[f]-fStart];
353: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
355: /* non-conforming meshes */
356: DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
357: for (ch = 0; ch < numChildren; ch++) {
358: const PetscInt child = children[ch];
360: if (child < fStart || child >= fEnd) continue;
361: row = rows[child-fStart];
362: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
363: }
364: }
365: }
366: ISRestoreIndices(fis, &rows);
367: ISRestoreIndices(cis, &cols);
368: ISDestroy(&fis);
369: ISDestroy(&cis);
370: MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
371: MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);
373: /* Get parallel CSR by doing conn^T * conn */
374: MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
375: MatDestroy(&conn);
377: /* extract local part of the CSR */
378: MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
379: MatDestroy(&CSR);
380: MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
381: if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
383: /* get back requested output */
384: if (numVertices) *numVertices = m;
385: if (offsets) {
386: PetscCalloc1(m+1, &idxs);
387: for (i = 1; i < m+1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
388: *offsets = idxs;
389: }
390: if (adjacency) {
391: PetscMalloc1(ii[m] - m, &idxs);
392: ISGetIndices(cis_own, &rows);
393: for (i = 0, c = 0; i < m; i++) {
394: PetscInt j, g = rows[i];
396: for (j = ii[i]; j < ii[i+1]; j++) {
397: if (jj[j] == g) continue; /* again, self-connectivity */
398: idxs[c++] = jj[j];
399: }
400: }
401: if (c != ii[m] - m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %D != %D",c,ii[m]-m);
402: ISRestoreIndices(cis_own, &rows);
403: *adjacency = idxs;
404: }
406: /* cleanup */
407: ISDestroy(&cis_own);
408: MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
409: if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
410: MatDestroy(&conn);
411: return(0);
412: }
414: /*@C
415: DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner
417: Input Parameters:
418: + dm - The mesh DM dm
419: - height - Height of the strata from which to construct the graph
421: Output Parameter:
422: + numVertices - Number of vertices in the graph
423: . offsets - Point offsets in the graph
424: . adjacency - Point connectivity in the graph
425: - globalNumbering - A map from the local cell numbering to the global numbering used in "adjacency". Negative indicates that the cell is a duplicate from another process.
427: The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function
428: representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().
430: Level: developer
432: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
433: @*/
434: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
435: {
437: PetscBool usemat = PETSC_FALSE;
440: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_csr_via_mat", &usemat, NULL);
441: if (usemat) {
442: DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
443: } else {
444: DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
445: }
446: return(0);
447: }
449: /*@C
450: DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.
452: Collective on DM
454: Input Arguments:
455: + dm - The DMPlex
456: - cellHeight - The height of mesh points to treat as cells (default should be 0)
458: Output Arguments:
459: + numVertices - The number of local vertices in the graph, or cells in the mesh.
460: . offsets - The offset to the adjacency list for each cell
461: - adjacency - The adjacency list for all cells
463: Note: This is suitable for input to a mesh partitioner like ParMetis.
465: Level: advanced
467: .seealso: DMPlexCreate()
468: @*/
469: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
470: {
471: const PetscInt maxFaceCases = 30;
472: PetscInt numFaceCases = 0;
473: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
474: PetscInt *off, *adj;
475: PetscInt *neighborCells = NULL;
476: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
480: /* For parallel partitioning, I think you have to communicate supports */
481: DMGetDimension(dm, &dim);
482: cellDim = dim - cellHeight;
483: DMPlexGetDepth(dm, &depth);
484: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
485: if (cEnd - cStart == 0) {
486: if (numVertices) *numVertices = 0;
487: if (offsets) *offsets = NULL;
488: if (adjacency) *adjacency = NULL;
489: return(0);
490: }
491: numCells = cEnd - cStart;
492: faceDepth = depth - cellHeight;
493: if (dim == depth) {
494: PetscInt f, fStart, fEnd;
496: PetscCalloc1(numCells+1, &off);
497: /* Count neighboring cells */
498: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
499: for (f = fStart; f < fEnd; ++f) {
500: const PetscInt *support;
501: PetscInt supportSize;
502: DMPlexGetSupportSize(dm, f, &supportSize);
503: DMPlexGetSupport(dm, f, &support);
504: if (supportSize == 2) {
505: PetscInt numChildren;
507: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
508: if (!numChildren) {
509: ++off[support[0]-cStart+1];
510: ++off[support[1]-cStart+1];
511: }
512: }
513: }
514: /* Prefix sum */
515: for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
516: if (adjacency) {
517: PetscInt *tmp;
519: PetscMalloc1(off[numCells], &adj);
520: PetscMalloc1(numCells+1, &tmp);
521: PetscArraycpy(tmp, off, numCells+1);
522: /* Get neighboring cells */
523: for (f = fStart; f < fEnd; ++f) {
524: const PetscInt *support;
525: PetscInt supportSize;
526: DMPlexGetSupportSize(dm, f, &supportSize);
527: DMPlexGetSupport(dm, f, &support);
528: if (supportSize == 2) {
529: PetscInt numChildren;
531: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
532: if (!numChildren) {
533: adj[tmp[support[0]-cStart]++] = support[1];
534: adj[tmp[support[1]-cStart]++] = support[0];
535: }
536: }
537: }
538: #if defined(PETSC_USE_DEBUG)
539: for (c = 0; c < cEnd-cStart; ++c) if (tmp[c] != off[c+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %d != %d for cell %d", tmp[c], off[c], c+cStart);
540: #endif
541: PetscFree(tmp);
542: }
543: if (numVertices) *numVertices = numCells;
544: if (offsets) *offsets = off;
545: if (adjacency) *adjacency = adj;
546: return(0);
547: }
548: /* Setup face recognition */
549: if (faceDepth == 1) {
550: PetscInt cornersSeen[30] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /* Could use PetscBT */
552: for (c = cStart; c < cEnd; ++c) {
553: PetscInt corners;
555: DMPlexGetConeSize(dm, c, &corners);
556: if (!cornersSeen[corners]) {
557: PetscInt nFV;
559: if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
560: cornersSeen[corners] = 1;
562: DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);
564: numFaceVertices[numFaceCases++] = nFV;
565: }
566: }
567: }
568: PetscCalloc1(numCells+1, &off);
569: /* Count neighboring cells */
570: for (cell = cStart; cell < cEnd; ++cell) {
571: PetscInt numNeighbors = PETSC_DETERMINE, n;
573: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
574: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
575: for (n = 0; n < numNeighbors; ++n) {
576: PetscInt cellPair[2];
577: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
578: PetscInt meetSize = 0;
579: const PetscInt *meet = NULL;
581: cellPair[0] = cell; cellPair[1] = neighborCells[n];
582: if (cellPair[0] == cellPair[1]) continue;
583: if (!found) {
584: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
585: if (meetSize) {
586: PetscInt f;
588: for (f = 0; f < numFaceCases; ++f) {
589: if (numFaceVertices[f] == meetSize) {
590: found = PETSC_TRUE;
591: break;
592: }
593: }
594: }
595: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
596: }
597: if (found) ++off[cell-cStart+1];
598: }
599: }
600: /* Prefix sum */
601: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];
603: if (adjacency) {
604: PetscMalloc1(off[numCells], &adj);
605: /* Get neighboring cells */
606: for (cell = cStart; cell < cEnd; ++cell) {
607: PetscInt numNeighbors = PETSC_DETERMINE, n;
608: PetscInt cellOffset = 0;
610: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
611: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
612: for (n = 0; n < numNeighbors; ++n) {
613: PetscInt cellPair[2];
614: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
615: PetscInt meetSize = 0;
616: const PetscInt *meet = NULL;
618: cellPair[0] = cell; cellPair[1] = neighborCells[n];
619: if (cellPair[0] == cellPair[1]) continue;
620: if (!found) {
621: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
622: if (meetSize) {
623: PetscInt f;
625: for (f = 0; f < numFaceCases; ++f) {
626: if (numFaceVertices[f] == meetSize) {
627: found = PETSC_TRUE;
628: break;
629: }
630: }
631: }
632: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
633: }
634: if (found) {
635: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
636: ++cellOffset;
637: }
638: }
639: }
640: }
641: PetscFree(neighborCells);
642: if (numVertices) *numVertices = numCells;
643: if (offsets) *offsets = off;
644: if (adjacency) *adjacency = adj;
645: return(0);
646: }
648: /*@C
649: PetscPartitionerRegister - Adds a new PetscPartitioner implementation
651: Not Collective
653: Input Parameters:
654: + name - The name of a new user-defined creation routine
655: - create_func - The creation routine itself
657: Notes:
658: PetscPartitionerRegister() may be called multiple times to add several user-defined PetscPartitioners
660: Sample usage:
661: .vb
662: PetscPartitionerRegister("my_part", MyPetscPartitionerCreate);
663: .ve
665: Then, your PetscPartitioner type can be chosen with the procedural interface via
666: .vb
667: PetscPartitionerCreate(MPI_Comm, PetscPartitioner *);
668: PetscPartitionerSetType(PetscPartitioner, "my_part");
669: .ve
670: or at runtime via the option
671: .vb
672: -petscpartitioner_type my_part
673: .ve
675: Level: advanced
677: .seealso: PetscPartitionerRegisterAll(), PetscPartitionerRegisterDestroy()
679: @*/
680: PetscErrorCode PetscPartitionerRegister(const char sname[], PetscErrorCode (*function)(PetscPartitioner))
681: {
685: PetscFunctionListAdd(&PetscPartitionerList, sname, function);
686: return(0);
687: }
689: /*@C
690: PetscPartitionerSetType - Builds a particular PetscPartitioner
692: Collective on PetscPartitioner
694: Input Parameters:
695: + part - The PetscPartitioner object
696: - name - The kind of partitioner
698: Options Database Key:
699: . -petscpartitioner_type <type> - Sets the PetscPartitioner type; use -help for a list of available types
701: Level: intermediate
703: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate()
704: @*/
705: PetscErrorCode PetscPartitionerSetType(PetscPartitioner part, PetscPartitionerType name)
706: {
707: PetscErrorCode (*r)(PetscPartitioner);
708: PetscBool match;
713: PetscObjectTypeCompare((PetscObject) part, name, &match);
714: if (match) return(0);
716: PetscFunctionListFind(PetscPartitionerList, name, &r);
717: if (!r) SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscPartitioner type: %s", name);
719: if (part->ops->destroy) {
720: (*part->ops->destroy)(part);
721: }
722: part->noGraph = PETSC_FALSE;
723: PetscMemzero(part->ops, sizeof(*part->ops));
724: PetscObjectChangeTypeName((PetscObject) part, name);
725: (*r)(part);
726: return(0);
727: }
729: /*@C
730: PetscPartitionerGetType - Gets the PetscPartitioner type name (as a string) from the object.
732: Not Collective
734: Input Parameter:
735: . part - The PetscPartitioner
737: Output Parameter:
738: . name - The PetscPartitioner type name
740: Level: intermediate
742: .seealso: PetscPartitionerSetType(), PetscPartitionerCreate()
743: @*/
744: PetscErrorCode PetscPartitionerGetType(PetscPartitioner part, PetscPartitionerType *name)
745: {
749: *name = ((PetscObject) part)->type_name;
750: return(0);
751: }
753: /*@C
754: PetscPartitionerViewFromOptions - View from Options
756: Collective on PetscPartitioner
758: Input Parameters:
759: + A - the PetscPartitioner object
760: . obj - Optional object
761: - name - command line option
763: Level: intermediate
764: .seealso: PetscPartitionerView(), PetscObjectViewFromOptions()
765: @*/
766: PetscErrorCode PetscPartitionerViewFromOptions(PetscPartitioner A,PetscObject obj,const char name[])
767: {
772: PetscObjectViewFromOptions((PetscObject)A,obj,name);
773: return(0);
774: }
776: /*@C
777: PetscPartitionerView - Views a PetscPartitioner
779: Collective on PetscPartitioner
781: Input Parameter:
782: + part - the PetscPartitioner object to view
783: - v - the viewer
785: Level: developer
787: .seealso: PetscPartitionerDestroy()
788: @*/
789: PetscErrorCode PetscPartitionerView(PetscPartitioner part, PetscViewer v)
790: {
791: PetscMPIInt size;
792: PetscBool isascii;
797: if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) part), &v);}
798: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &isascii);
799: if (isascii) {
800: MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
801: PetscViewerASCIIPrintf(v, "Graph Partitioner: %d MPI Process%s\n", size, size > 1 ? "es" : "");
802: PetscViewerASCIIPrintf(v, " type: %s\n", part->hdr.type_name);
803: PetscViewerASCIIPrintf(v, " edge cut: %D\n", part->edgeCut);
804: PetscViewerASCIIPrintf(v, " balance: %.2g\n", part->balance);
805: PetscViewerASCIIPrintf(v, " use vertex weights: %d\n", part->usevwgt);
806: }
807: if (part->ops->view) {(*part->ops->view)(part, v);}
808: return(0);
809: }
811: static PetscErrorCode PetscPartitionerGetDefaultType(const char *currentType, const char **defaultType)
812: {
814: if (!currentType) {
815: #if defined(PETSC_HAVE_PARMETIS)
816: *defaultType = PETSCPARTITIONERPARMETIS;
817: #elif defined(PETSC_HAVE_PTSCOTCH)
818: *defaultType = PETSCPARTITIONERPTSCOTCH;
819: #elif defined(PETSC_HAVE_CHACO)
820: *defaultType = PETSCPARTITIONERCHACO;
821: #else
822: *defaultType = PETSCPARTITIONERSIMPLE;
823: #endif
824: } else {
825: *defaultType = currentType;
826: }
827: return(0);
828: }
830: /*@
831: PetscPartitionerSetFromOptions - sets parameters in a PetscPartitioner from the options database
833: Collective on PetscPartitioner
835: Input Parameter:
836: . part - the PetscPartitioner object to set options for
838: Options Database Keys:
839: + -petscpartitioner_type <type> - Sets the PetscPartitioner type; use -help for a list of available types
840: . -petscpartitioner_use_vertex_weights - Uses weights associated with the graph vertices
841: - -petscpartitioner_view_graph - View the graph each time PetscPartitionerPartition is called. Viewer can be customized, see PetscOptionsGetViewer()
843: Level: developer
845: .seealso: PetscPartitionerView(), PetscPartitionerSetType(), PetscPartitionerPartition()
846: @*/
847: PetscErrorCode PetscPartitionerSetFromOptions(PetscPartitioner part)
848: {
849: const char *defaultType;
850: char name[256];
851: PetscBool flg;
856: PetscPartitionerGetDefaultType(((PetscObject) part)->type_name,&defaultType);
857: PetscObjectOptionsBegin((PetscObject) part);
858: PetscOptionsFList("-petscpartitioner_type", "Graph partitioner", "PetscPartitionerSetType", PetscPartitionerList, defaultType, name, sizeof(name), &flg);
859: if (flg) {
860: PetscPartitionerSetType(part, name);
861: } else if (!((PetscObject) part)->type_name) {
862: PetscPartitionerSetType(part, defaultType);
863: }
864: PetscOptionsBool("-petscpartitioner_use_vertex_weights","Use vertex weights","",part->usevwgt,&part->usevwgt,NULL);
865: if (part->ops->setfromoptions) {
866: (*part->ops->setfromoptions)(PetscOptionsObject,part);
867: }
868: PetscViewerDestroy(&part->viewer);
869: PetscViewerDestroy(&part->viewerGraph);
870: PetscOptionsGetViewer(((PetscObject) part)->comm, ((PetscObject) part)->options, ((PetscObject) part)->prefix, "-petscpartitioner_view", &part->viewer, NULL, NULL);
871: PetscOptionsGetViewer(((PetscObject) part)->comm, ((PetscObject) part)->options, ((PetscObject) part)->prefix, "-petscpartitioner_view_graph", &part->viewerGraph, NULL, &part->viewGraph);
872: /* process any options handlers added with PetscObjectAddOptionsHandler() */
873: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) part);
874: PetscOptionsEnd();
875: return(0);
876: }
878: /*@C
879: PetscPartitionerSetUp - Construct data structures for the PetscPartitioner
881: Collective on PetscPartitioner
883: Input Parameter:
884: . part - the PetscPartitioner object to setup
886: Level: developer
888: .seealso: PetscPartitionerView(), PetscPartitionerDestroy()
889: @*/
890: PetscErrorCode PetscPartitionerSetUp(PetscPartitioner part)
891: {
896: if (part->ops->setup) {(*part->ops->setup)(part);}
897: return(0);
898: }
900: /*@
901: PetscPartitionerDestroy - Destroys a PetscPartitioner object
903: Collective on PetscPartitioner
905: Input Parameter:
906: . part - the PetscPartitioner object to destroy
908: Level: developer
910: .seealso: PetscPartitionerView()
911: @*/
912: PetscErrorCode PetscPartitionerDestroy(PetscPartitioner *part)
913: {
917: if (!*part) return(0);
920: if (--((PetscObject)(*part))->refct > 0) {*part = 0; return(0);}
921: ((PetscObject) (*part))->refct = 0;
923: PetscViewerDestroy(&(*part)->viewer);
924: PetscViewerDestroy(&(*part)->viewerGraph);
925: if ((*part)->ops->destroy) {(*(*part)->ops->destroy)(*part);}
926: PetscHeaderDestroy(part);
927: return(0);
928: }
930: /*@
931: PetscPartitionerPartition - Partition a graph
933: Collective on PetscPartitioner
935: Input Parameters:
936: + part - The PetscPartitioner
937: . nparts - Number of partitions
938: . numVertices - Number of vertices in the local part of the graph
939: . start - row pointers for the local part of the graph (CSR style)
940: . adjacency - adjacency list (CSR style)
941: . vertexSection - PetscSection describing the absolute weight of each local vertex (can be NULL)
942: - targetSection - PetscSection describing the absolute weight of each partition (can be NULL)
944: Output Parameters:
945: + partSection - The PetscSection giving the division of points by partition
946: - partition - The list of points by partition
948: Options Database:
949: + -petscpartitioner_view - View the partitioner information
950: - -petscpartitioner_view_graph - View the graph we are partitioning
952: Notes:
953: The chart of the vertexSection (if present) must contain [0,numVertices), with the number of dofs in the section specifying the absolute weight for each vertex.
954: The chart of the targetSection (if present) must contain [0,nparts), with the number of dofs in the section specifying the absolute weight for each partition. This information must be the same across processes, PETSc does not check it.
956: Level: developer
958: .seealso PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscSectionSetDof()
959: @*/
960: PetscErrorCode PetscPartitionerPartition(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertexSection, PetscSection targetSection, PetscSection partSection, IS *partition)
961: {
967: if (nparts <= 0) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Number of parts must be positive");
968: if (numVertices < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices must be non-negative");
969: if (numVertices && !part->noGraph) {
973: }
974: if (vertexSection) {
975: PetscInt s,e;
978: PetscSectionGetChart(vertexSection, &s, &e);
979: if (s > 0 || e < numVertices) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid vertexSection chart [%D,%D)",s,e);
980: }
981: if (targetSection) {
982: PetscInt s,e;
985: PetscSectionGetChart(targetSection, &s, &e);
986: if (s > 0 || e < nparts) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid targetSection chart [%D,%D)",s,e);
987: }
991: PetscSectionReset(partSection);
992: PetscSectionSetChart(partSection, 0, nparts);
993: if (nparts == 1) { /* quick */
994: PetscSectionSetDof(partSection, 0, numVertices);
995: ISCreateStride(PetscObjectComm((PetscObject)part),numVertices,0,1,partition);
996: } else {
997: if (!part->ops->partition) SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "PetscPartitioner %s has no partitioning method", ((PetscObject)part)->type_name);
998: (*part->ops->partition)(part, nparts, numVertices, start, adjacency, vertexSection, targetSection, partSection, partition);
999: }
1000: PetscSectionSetUp(partSection);
1001: if (part->viewerGraph) {
1002: PetscViewer viewer = part->viewerGraph;
1003: PetscBool isascii;
1004: PetscInt v, i;
1005: PetscMPIInt rank;
1007: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
1008: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
1009: if (isascii) {
1010: PetscViewerASCIIPushSynchronized(viewer);
1011: PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Nv: %D\n", rank, numVertices);
1012: for (v = 0; v < numVertices; ++v) {
1013: const PetscInt s = start[v];
1014: const PetscInt e = start[v+1];
1016: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] ", rank);
1017: for (i = s; i < e; ++i) {PetscViewerASCIISynchronizedPrintf(viewer, "%D ", adjacency[i]);}
1018: PetscViewerASCIISynchronizedPrintf(viewer, "[%D-%D)\n", s, e);
1019: }
1020: PetscViewerFlush(viewer);
1021: PetscViewerASCIIPopSynchronized(viewer);
1022: }
1023: }
1024: if (part->viewer) {
1025: PetscPartitionerView(part,part->viewer);
1026: }
1027: return(0);
1028: }
1030: /*@
1031: PetscPartitionerCreate - Creates an empty PetscPartitioner object. The type can then be set with PetscPartitionerSetType().
1033: Collective
1035: Input Parameter:
1036: . comm - The communicator for the PetscPartitioner object
1038: Output Parameter:
1039: . part - The PetscPartitioner object
1041: Level: beginner
1043: .seealso: PetscPartitionerSetType(), PETSCPARTITIONERCHACO, PETSCPARTITIONERPARMETIS, PETSCPARTITIONERSHELL, PETSCPARTITIONERSIMPLE, PETSCPARTITIONERGATHER
1044: @*/
1045: PetscErrorCode PetscPartitionerCreate(MPI_Comm comm, PetscPartitioner *part)
1046: {
1047: PetscPartitioner p;
1048: const char *partitionerType = NULL;
1049: PetscErrorCode ierr;
1053: *part = NULL;
1054: DMInitializePackage();
1056: PetscHeaderCreate(p, PETSCPARTITIONER_CLASSID, "PetscPartitioner", "Graph Partitioner", "PetscPartitioner", comm, PetscPartitionerDestroy, PetscPartitionerView);
1057: PetscPartitionerGetDefaultType(NULL,&partitionerType);
1058: PetscPartitionerSetType(p,partitionerType);
1060: p->edgeCut = 0;
1061: p->balance = 0.0;
1062: p->usevwgt = PETSC_TRUE;
1064: *part = p;
1065: return(0);
1066: }
1068: /*@
1069: PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh
1071: Collective on PetscPartitioner
1073: Input Parameters:
1074: + part - The PetscPartitioner
1075: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
1076: - dm - The mesh DM
1078: Output Parameters:
1079: + partSection - The PetscSection giving the division of points by partition
1080: - partition - The list of points by partition
1082: Notes:
1083: If the DM has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
1084: by the section in the transitive closure of the point.
1086: Level: developer
1088: .seealso DMPlexDistribute(), PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscPartitionerPartition()
1089: @*/
1090: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
1091: {
1092: PetscMPIInt size;
1093: PetscBool isplex;
1095: PetscSection vertSection = NULL;
1103: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isplex);
1104: if (!isplex) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not for type %s",((PetscObject)dm)->type_name);
1105: MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
1106: if (size == 1) {
1107: PetscInt *points;
1108: PetscInt cStart, cEnd, c;
1110: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
1111: PetscSectionReset(partSection);
1112: PetscSectionSetChart(partSection, 0, size);
1113: PetscSectionSetDof(partSection, 0, cEnd-cStart);
1114: PetscSectionSetUp(partSection);
1115: PetscMalloc1(cEnd-cStart, &points);
1116: for (c = cStart; c < cEnd; ++c) points[c] = c;
1117: ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
1118: return(0);
1119: }
1120: if (part->height == 0) {
1121: PetscInt numVertices = 0;
1122: PetscInt *start = NULL;
1123: PetscInt *adjacency = NULL;
1124: IS globalNumbering;
1126: if (!part->noGraph || part->viewGraph) {
1127: DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
1128: } else { /* only compute the number of owned local vertices */
1129: const PetscInt *idxs;
1130: PetscInt p, pStart, pEnd;
1132: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
1133: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
1134: ISGetIndices(globalNumbering, &idxs);
1135: for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
1136: ISRestoreIndices(globalNumbering, &idxs);
1137: }
1138: if (part->usevwgt) {
1139: PetscSection section = dm->localSection, clSection = NULL;
1140: IS clPoints = NULL;
1141: const PetscInt *gid,*clIdx;
1142: PetscInt v, p, pStart, pEnd;
1144: /* dm->localSection encodes degrees of freedom per point, not per cell. We need to get the closure index to properly specify cell weights (aka dofs) */
1145: /* We do this only if the local section has been set */
1146: if (section) {
1147: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
1148: if (!clSection) {
1149: DMPlexCreateClosureIndex(dm,NULL);
1150: }
1151: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
1152: ISGetIndices(clPoints,&clIdx);
1153: }
1154: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
1155: PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
1156: PetscSectionSetChart(vertSection, 0, numVertices);
1157: if (globalNumbering) {
1158: ISGetIndices(globalNumbering,&gid);
1159: } else gid = NULL;
1160: for (p = pStart, v = 0; p < pEnd; ++p) {
1161: PetscInt dof = 1;
1163: /* skip cells in the overlap */
1164: if (gid && gid[p-pStart] < 0) continue;
1166: if (section) {
1167: PetscInt cl, clSize, clOff;
1169: dof = 0;
1170: PetscSectionGetDof(clSection, p, &clSize);
1171: PetscSectionGetOffset(clSection, p, &clOff);
1172: for (cl = 0; cl < clSize; cl+=2) {
1173: PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */
1175: PetscSectionGetDof(section, clPoint, &clDof);
1176: dof += clDof;
1177: }
1178: }
1179: if (!dof) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of dofs for point %D in the local section should be positive",p);
1180: PetscSectionSetDof(vertSection, v, dof);
1181: v++;
1182: }
1183: if (globalNumbering) {
1184: ISRestoreIndices(globalNumbering,&gid);
1185: }
1186: if (clPoints) {
1187: ISRestoreIndices(clPoints,&clIdx);
1188: }
1189: PetscSectionSetUp(vertSection);
1190: }
1191: PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
1192: PetscFree(start);
1193: PetscFree(adjacency);
1194: if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
1195: const PetscInt *globalNum;
1196: const PetscInt *partIdx;
1197: PetscInt *map, cStart, cEnd;
1198: PetscInt *adjusted, i, localSize, offset;
1199: IS newPartition;
1201: ISGetLocalSize(*partition,&localSize);
1202: PetscMalloc1(localSize,&adjusted);
1203: ISGetIndices(globalNumbering,&globalNum);
1204: ISGetIndices(*partition,&partIdx);
1205: PetscMalloc1(localSize,&map);
1206: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
1207: for (i = cStart, offset = 0; i < cEnd; i++) {
1208: if (globalNum[i - cStart] >= 0) map[offset++] = i;
1209: }
1210: for (i = 0; i < localSize; i++) {
1211: adjusted[i] = map[partIdx[i]];
1212: }
1213: PetscFree(map);
1214: ISRestoreIndices(*partition,&partIdx);
1215: ISRestoreIndices(globalNumbering,&globalNum);
1216: ISCreateGeneral(PETSC_COMM_SELF,localSize,adjusted,PETSC_OWN_POINTER,&newPartition);
1217: ISDestroy(&globalNumbering);
1218: ISDestroy(partition);
1219: *partition = newPartition;
1220: }
1221: } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
1222: PetscSectionDestroy(&vertSection);
1223: return(0);
1224: }
1226: static PetscErrorCode PetscPartitionerDestroy_Shell(PetscPartitioner part)
1227: {
1228: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1229: PetscErrorCode ierr;
1232: PetscSectionDestroy(&p->section);
1233: ISDestroy(&p->partition);
1234: PetscFree(p);
1235: return(0);
1236: }
1238: static PetscErrorCode PetscPartitionerView_Shell_Ascii(PetscPartitioner part, PetscViewer viewer)
1239: {
1240: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1241: PetscErrorCode ierr;
1244: if (p->random) {
1245: PetscViewerASCIIPushTab(viewer);
1246: PetscViewerASCIIPrintf(viewer, "using random partition\n");
1247: PetscViewerASCIIPopTab(viewer);
1248: }
1249: return(0);
1250: }
1252: static PetscErrorCode PetscPartitionerView_Shell(PetscPartitioner part, PetscViewer viewer)
1253: {
1254: PetscBool iascii;
1260: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1261: if (iascii) {PetscPartitionerView_Shell_Ascii(part, viewer);}
1262: return(0);
1263: }
1265: static PetscErrorCode PetscPartitionerSetFromOptions_Shell(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1266: {
1267: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1268: PetscErrorCode ierr;
1271: PetscOptionsHead(PetscOptionsObject, "PetscPartitioner Shell Options");
1272: PetscOptionsBool("-petscpartitioner_shell_random", "Use a random partition", "PetscPartitionerView", PETSC_FALSE, &p->random, NULL);
1273: PetscOptionsTail();
1274: return(0);
1275: }
1277: static PetscErrorCode PetscPartitionerPartition_Shell(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1278: {
1279: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1280: PetscInt np;
1281: PetscErrorCode ierr;
1284: if (p->random) {
1285: PetscRandom r;
1286: PetscInt *sizes, *points, v, p;
1287: PetscMPIInt rank;
1289: MPI_Comm_rank(PetscObjectComm((PetscObject) part), &rank);
1290: PetscRandomCreate(PETSC_COMM_SELF, &r);
1291: PetscRandomSetInterval(r, 0.0, (PetscScalar) nparts);
1292: PetscRandomSetFromOptions(r);
1293: PetscCalloc2(nparts, &sizes, numVertices, &points);
1294: for (v = 0; v < numVertices; ++v) {points[v] = v;}
1295: for (p = 0; p < nparts; ++p) {sizes[p] = numVertices/nparts + (PetscInt) (p < numVertices % nparts);}
1296: for (v = numVertices-1; v > 0; --v) {
1297: PetscReal val;
1298: PetscInt w, tmp;
1300: PetscRandomSetInterval(r, 0.0, (PetscScalar) (v+1));
1301: PetscRandomGetValueReal(r, &val);
1302: w = PetscFloorReal(val);
1303: tmp = points[v];
1304: points[v] = points[w];
1305: points[w] = tmp;
1306: }
1307: PetscRandomDestroy(&r);
1308: PetscPartitionerShellSetPartition(part, nparts, sizes, points);
1309: PetscFree2(sizes, points);
1310: }
1311: if (!p->section) SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_WRONG, "Shell partitioner information not provided. Please call PetscPartitionerShellSetPartition()");
1312: PetscSectionGetChart(p->section, NULL, &np);
1313: if (nparts != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of requested partitions %d != configured partitions %d", nparts, np);
1314: ISGetLocalSize(p->partition, &np);
1315: if (numVertices != np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of input vertices %d != configured vertices %d", numVertices, np);
1316: PetscSectionCopy(p->section, partSection);
1317: *partition = p->partition;
1318: PetscObjectReference((PetscObject) p->partition);
1319: return(0);
1320: }
1322: static PetscErrorCode PetscPartitionerInitialize_Shell(PetscPartitioner part)
1323: {
1325: part->noGraph = PETSC_TRUE; /* PetscPartitionerShell cannot overload the partition call, so it is safe for now */
1326: part->ops->view = PetscPartitionerView_Shell;
1327: part->ops->setfromoptions = PetscPartitionerSetFromOptions_Shell;
1328: part->ops->destroy = PetscPartitionerDestroy_Shell;
1329: part->ops->partition = PetscPartitionerPartition_Shell;
1330: return(0);
1331: }
1333: /*MC
1334: PETSCPARTITIONERSHELL = "shell" - A PetscPartitioner object
1336: Level: intermediate
1338: Options Database Keys:
1339: . -petscpartitioner_shell_random - Use a random partition
1341: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1342: M*/
1344: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Shell(PetscPartitioner part)
1345: {
1346: PetscPartitioner_Shell *p;
1347: PetscErrorCode ierr;
1351: PetscNewLog(part, &p);
1352: part->data = p;
1354: PetscPartitionerInitialize_Shell(part);
1355: p->random = PETSC_FALSE;
1356: return(0);
1357: }
1359: /*@C
1360: PetscPartitionerShellSetPartition - Set an artifical partition for a mesh
1362: Collective on PetscPartitioner
1364: Input Parameters:
1365: + part - The PetscPartitioner
1366: . size - The number of partitions
1367: . sizes - array of length size (or NULL) providing the number of points in each partition
1368: - points - array of length sum(sizes) (may be NULL iff sizes is NULL), a permutation of the points that groups those assigned to each partition in order (i.e., partition 0 first, partition 1 next, etc.)
1370: Level: developer
1372: Notes:
1373: It is safe to free the sizes and points arrays after use in this routine.
1375: .seealso DMPlexDistribute(), PetscPartitionerCreate()
1376: @*/
1377: PetscErrorCode PetscPartitionerShellSetPartition(PetscPartitioner part, PetscInt size, const PetscInt sizes[], const PetscInt points[])
1378: {
1379: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1380: PetscInt proc, numPoints;
1381: PetscErrorCode ierr;
1387: PetscSectionDestroy(&p->section);
1388: ISDestroy(&p->partition);
1389: PetscSectionCreate(PetscObjectComm((PetscObject) part), &p->section);
1390: PetscSectionSetChart(p->section, 0, size);
1391: if (sizes) {
1392: for (proc = 0; proc < size; ++proc) {
1393: PetscSectionSetDof(p->section, proc, sizes[proc]);
1394: }
1395: }
1396: PetscSectionSetUp(p->section);
1397: PetscSectionGetStorageSize(p->section, &numPoints);
1398: ISCreateGeneral(PetscObjectComm((PetscObject) part), numPoints, points, PETSC_COPY_VALUES, &p->partition);
1399: return(0);
1400: }
1402: /*@
1403: PetscPartitionerShellSetRandom - Set the flag to use a random partition
1405: Collective on PetscPartitioner
1407: Input Parameters:
1408: + part - The PetscPartitioner
1409: - random - The flag to use a random partition
1411: Level: intermediate
1413: .seealso PetscPartitionerShellGetRandom(), PetscPartitionerCreate()
1414: @*/
1415: PetscErrorCode PetscPartitionerShellSetRandom(PetscPartitioner part, PetscBool random)
1416: {
1417: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1421: p->random = random;
1422: return(0);
1423: }
1425: /*@
1426: PetscPartitionerShellGetRandom - get the flag to use a random partition
1428: Collective on PetscPartitioner
1430: Input Parameter:
1431: . part - The PetscPartitioner
1433: Output Parameter:
1434: . random - The flag to use a random partition
1436: Level: intermediate
1438: .seealso PetscPartitionerShellSetRandom(), PetscPartitionerCreate()
1439: @*/
1440: PetscErrorCode PetscPartitionerShellGetRandom(PetscPartitioner part, PetscBool *random)
1441: {
1442: PetscPartitioner_Shell *p = (PetscPartitioner_Shell *) part->data;
1447: *random = p->random;
1448: return(0);
1449: }
1451: static PetscErrorCode PetscPartitionerDestroy_Simple(PetscPartitioner part)
1452: {
1453: PetscPartitioner_Simple *p = (PetscPartitioner_Simple *) part->data;
1454: PetscErrorCode ierr;
1457: PetscFree(p);
1458: return(0);
1459: }
1461: static PetscErrorCode PetscPartitionerView_Simple_Ascii(PetscPartitioner part, PetscViewer viewer)
1462: {
1464: return(0);
1465: }
1467: static PetscErrorCode PetscPartitionerView_Simple(PetscPartitioner part, PetscViewer viewer)
1468: {
1469: PetscBool iascii;
1475: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1476: if (iascii) {PetscPartitionerView_Simple_Ascii(part, viewer);}
1477: return(0);
1478: }
1480: static PetscErrorCode PetscPartitionerPartition_Simple(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1481: {
1482: MPI_Comm comm;
1483: PetscInt np, *tpwgts = NULL, sumw = 0, numVerticesGlobal = 0;
1484: PetscMPIInt size;
1488: if (vertSection) { PetscInfo(part,"PETSCPARTITIONERSIMPLE ignores vertex weights\n"); }
1489: comm = PetscObjectComm((PetscObject)part);
1490: MPI_Comm_size(comm,&size);
1491: if (targetSection) {
1492: MPIU_Allreduce(&numVertices, &numVerticesGlobal, 1, MPIU_INT, MPI_SUM, comm);
1493: PetscCalloc1(nparts,&tpwgts);
1494: for (np = 0; np < nparts; ++np) {
1495: PetscSectionGetDof(targetSection,np,&tpwgts[np]);
1496: sumw += tpwgts[np];
1497: }
1498: if (!sumw) {
1499: PetscFree(tpwgts);
1500: } else {
1501: PetscInt m,mp;
1502: for (np = 0; np < nparts; ++np) tpwgts[np] = (tpwgts[np]*numVerticesGlobal)/sumw;
1503: for (np = 0, m = -1, mp = 0, sumw = 0; np < nparts; ++np) {
1504: if (m < tpwgts[np]) { m = tpwgts[np]; mp = np; }
1505: sumw += tpwgts[np];
1506: }
1507: if (sumw != numVerticesGlobal) tpwgts[mp] += numVerticesGlobal - sumw;
1508: }
1509: }
1511: ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1512: if (size == 1) {
1513: if (tpwgts) {
1514: for (np = 0; np < nparts; ++np) {
1515: PetscSectionSetDof(partSection, np, tpwgts[np]);
1516: }
1517: } else {
1518: for (np = 0; np < nparts; ++np) {
1519: PetscSectionSetDof(partSection, np, numVertices/nparts + ((numVertices % nparts) > np));
1520: }
1521: }
1522: } else {
1523: if (tpwgts) {
1524: Vec v;
1525: PetscScalar *array;
1526: PetscInt st,j;
1527: PetscMPIInt rank;
1529: VecCreate(comm,&v);
1530: VecSetSizes(v,numVertices,numVerticesGlobal);
1531: VecSetType(v,VECSTANDARD);
1532: MPI_Comm_rank(comm,&rank);
1533: for (np = 0,st = 0; np < nparts; ++np) {
1534: if (rank == np || (rank == size-1 && size < nparts && np >= size)) {
1535: for (j = 0; j < tpwgts[np]; j++) {
1536: VecSetValue(v,st+j,np,INSERT_VALUES);
1537: }
1538: }
1539: st += tpwgts[np];
1540: }
1541: VecAssemblyBegin(v);
1542: VecAssemblyEnd(v);
1543: VecGetArray(v,&array);
1544: for (j = 0; j < numVertices; ++j) {
1545: PetscSectionAddDof(partSection,PetscRealPart(array[j]),1);
1546: }
1547: VecRestoreArray(v,&array);
1548: VecDestroy(&v);
1549: } else {
1550: PetscMPIInt rank;
1551: PetscInt nvGlobal, *offsets, myFirst, myLast;
1553: PetscMalloc1(size+1,&offsets);
1554: offsets[0] = 0;
1555: MPI_Allgather(&numVertices,1,MPIU_INT,&offsets[1],1,MPIU_INT,comm);
1556: for (np = 2; np <= size; np++) {
1557: offsets[np] += offsets[np-1];
1558: }
1559: nvGlobal = offsets[size];
1560: MPI_Comm_rank(comm,&rank);
1561: myFirst = offsets[rank];
1562: myLast = offsets[rank + 1] - 1;
1563: PetscFree(offsets);
1564: if (numVertices) {
1565: PetscInt firstPart = 0, firstLargePart = 0;
1566: PetscInt lastPart = 0, lastLargePart = 0;
1567: PetscInt rem = nvGlobal % nparts;
1568: PetscInt pSmall = nvGlobal/nparts;
1569: PetscInt pBig = nvGlobal/nparts + 1;
1571: if (rem) {
1572: firstLargePart = myFirst / pBig;
1573: lastLargePart = myLast / pBig;
1575: if (firstLargePart < rem) {
1576: firstPart = firstLargePart;
1577: } else {
1578: firstPart = rem + (myFirst - (rem * pBig)) / pSmall;
1579: }
1580: if (lastLargePart < rem) {
1581: lastPart = lastLargePart;
1582: } else {
1583: lastPart = rem + (myLast - (rem * pBig)) / pSmall;
1584: }
1585: } else {
1586: firstPart = myFirst / (nvGlobal/nparts);
1587: lastPart = myLast / (nvGlobal/nparts);
1588: }
1590: for (np = firstPart; np <= lastPart; np++) {
1591: PetscInt PartStart = np * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np);
1592: PetscInt PartEnd = (np+1) * (nvGlobal/nparts) + PetscMin(nvGlobal % nparts,np+1);
1594: PartStart = PetscMax(PartStart,myFirst);
1595: PartEnd = PetscMin(PartEnd,myLast+1);
1596: PetscSectionSetDof(partSection,np,PartEnd-PartStart);
1597: }
1598: }
1599: }
1600: }
1601: PetscFree(tpwgts);
1602: return(0);
1603: }
1605: static PetscErrorCode PetscPartitionerInitialize_Simple(PetscPartitioner part)
1606: {
1608: part->noGraph = PETSC_TRUE;
1609: part->ops->view = PetscPartitionerView_Simple;
1610: part->ops->destroy = PetscPartitionerDestroy_Simple;
1611: part->ops->partition = PetscPartitionerPartition_Simple;
1612: return(0);
1613: }
1615: /*MC
1616: PETSCPARTITIONERSIMPLE = "simple" - A PetscPartitioner object
1618: Level: intermediate
1620: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1621: M*/
1623: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Simple(PetscPartitioner part)
1624: {
1625: PetscPartitioner_Simple *p;
1626: PetscErrorCode ierr;
1630: PetscNewLog(part, &p);
1631: part->data = p;
1633: PetscPartitionerInitialize_Simple(part);
1634: return(0);
1635: }
1637: static PetscErrorCode PetscPartitionerDestroy_Gather(PetscPartitioner part)
1638: {
1639: PetscPartitioner_Gather *p = (PetscPartitioner_Gather *) part->data;
1640: PetscErrorCode ierr;
1643: PetscFree(p);
1644: return(0);
1645: }
1647: static PetscErrorCode PetscPartitionerView_Gather_Ascii(PetscPartitioner part, PetscViewer viewer)
1648: {
1650: return(0);
1651: }
1653: static PetscErrorCode PetscPartitionerView_Gather(PetscPartitioner part, PetscViewer viewer)
1654: {
1655: PetscBool iascii;
1661: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1662: if (iascii) {PetscPartitionerView_Gather_Ascii(part, viewer);}
1663: return(0);
1664: }
1666: static PetscErrorCode PetscPartitionerPartition_Gather(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1667: {
1668: PetscInt np;
1672: ISCreateStride(PETSC_COMM_SELF, numVertices, 0, 1, partition);
1673: PetscSectionSetDof(partSection,0,numVertices);
1674: for (np = 1; np < nparts; ++np) {PetscSectionSetDof(partSection, np, 0);}
1675: return(0);
1676: }
1678: static PetscErrorCode PetscPartitionerInitialize_Gather(PetscPartitioner part)
1679: {
1681: part->noGraph = PETSC_TRUE;
1682: part->ops->view = PetscPartitionerView_Gather;
1683: part->ops->destroy = PetscPartitionerDestroy_Gather;
1684: part->ops->partition = PetscPartitionerPartition_Gather;
1685: return(0);
1686: }
1688: /*MC
1689: PETSCPARTITIONERGATHER = "gather" - A PetscPartitioner object
1691: Level: intermediate
1693: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1694: M*/
1696: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Gather(PetscPartitioner part)
1697: {
1698: PetscPartitioner_Gather *p;
1699: PetscErrorCode ierr;
1703: PetscNewLog(part, &p);
1704: part->data = p;
1706: PetscPartitionerInitialize_Gather(part);
1707: return(0);
1708: }
1711: static PetscErrorCode PetscPartitionerDestroy_Chaco(PetscPartitioner part)
1712: {
1713: PetscPartitioner_Chaco *p = (PetscPartitioner_Chaco *) part->data;
1714: PetscErrorCode ierr;
1717: PetscFree(p);
1718: return(0);
1719: }
1721: static PetscErrorCode PetscPartitionerView_Chaco_Ascii(PetscPartitioner part, PetscViewer viewer)
1722: {
1724: return(0);
1725: }
1727: static PetscErrorCode PetscPartitionerView_Chaco(PetscPartitioner part, PetscViewer viewer)
1728: {
1729: PetscBool iascii;
1735: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1736: if (iascii) {PetscPartitionerView_Chaco_Ascii(part, viewer);}
1737: return(0);
1738: }
1740: #if defined(PETSC_HAVE_CHACO)
1741: #if defined(PETSC_HAVE_UNISTD_H)
1742: #include <unistd.h>
1743: #endif
1744: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1745: #include <chaco.h>
1746: #else
1747: /* Older versions of Chaco do not have an include file */
1748: PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts,
1749: float *ewgts, float *x, float *y, float *z, char *outassignname,
1750: char *outfilename, short *assignment, int architecture, int ndims_tot,
1751: int mesh_dims[3], double *goal, int global_method, int local_method,
1752: int rqi_flag, int vmax, int ndims, double eigtol, long seed);
1753: #endif
1754: extern int FREE_GRAPH;
1755: #endif
1757: static PetscErrorCode PetscPartitionerPartition_Chaco(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1758: {
1759: #if defined(PETSC_HAVE_CHACO)
1760: enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
1761: MPI_Comm comm;
1762: int nvtxs = numVertices; /* number of vertices in full graph */
1763: int *vwgts = NULL; /* weights for all vertices */
1764: float *ewgts = NULL; /* weights for all edges */
1765: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
1766: char *outassignname = NULL; /* name of assignment output file */
1767: char *outfilename = NULL; /* output file name */
1768: int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */
1769: int ndims_tot = 0; /* total number of cube dimensions to divide */
1770: int mesh_dims[3]; /* dimensions of mesh of processors */
1771: double *goal = NULL; /* desired set sizes for each set */
1772: int global_method = 1; /* global partitioning algorithm */
1773: int local_method = 1; /* local partitioning algorithm */
1774: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
1775: int vmax = 200; /* how many vertices to coarsen down to? */
1776: int ndims = 1; /* number of eigenvectors (2^d sets) */
1777: double eigtol = 0.001; /* tolerance on eigenvectors */
1778: long seed = 123636512; /* for random graph mutations */
1779: #if defined(PETSC_HAVE_CHACO_INT_ASSIGNMENT)
1780: int *assignment; /* Output partition */
1781: #else
1782: short int *assignment; /* Output partition */
1783: #endif
1784: int fd_stdout, fd_pipe[2];
1785: PetscInt *points;
1786: int i, v, p;
1790: PetscObjectGetComm((PetscObject)part,&comm);
1791: #if defined (PETSC_USE_DEBUG)
1792: {
1793: int ival,isum;
1794: PetscBool distributed;
1796: ival = (numVertices > 0);
1797: MPI_Allreduce(&ival, &isum, 1, MPI_INT, MPI_SUM, comm);
1798: distributed = (isum > 1) ? PETSC_TRUE : PETSC_FALSE;
1799: if (distributed) SETERRQ(comm, PETSC_ERR_SUP, "Chaco cannot partition a distributed graph");
1800: }
1801: #endif
1802: if (!numVertices) { /* distributed case, return if not holding the graph */
1803: ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
1804: return(0);
1805: }
1806: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
1807: for (i = 0; i < start[numVertices]; ++i) ++adjacency[i];
1809: if (global_method == INERTIAL_METHOD) {
1810: /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
1811: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
1812: }
1813: mesh_dims[0] = nparts;
1814: mesh_dims[1] = 1;
1815: mesh_dims[2] = 1;
1816: PetscMalloc1(nvtxs, &assignment);
1817: /* Chaco outputs to stdout. We redirect this to a buffer. */
1818: /* TODO: check error codes for UNIX calls */
1819: #if defined(PETSC_HAVE_UNISTD_H)
1820: {
1821: int piperet;
1822: piperet = pipe(fd_pipe);
1823: if (piperet) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"Could not create pipe");
1824: fd_stdout = dup(1);
1825: close(1);
1826: dup2(fd_pipe[1], 1);
1827: }
1828: #endif
1829: if (part->usevwgt) { PetscInfo(part,"PETSCPARTITIONERCHACO ignores vertex weights\n"); }
1830: interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
1831: assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
1832: vmax, ndims, eigtol, seed);
1833: #if defined(PETSC_HAVE_UNISTD_H)
1834: {
1835: char msgLog[10000];
1836: int count;
1838: fflush(stdout);
1839: count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
1840: if (count < 0) count = 0;
1841: msgLog[count] = 0;
1842: close(1);
1843: dup2(fd_stdout, 1);
1844: close(fd_stdout);
1845: close(fd_pipe[0]);
1846: close(fd_pipe[1]);
1847: if (ierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
1848: }
1849: #else
1850: if (ierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in Chaco library: %s", "error in stdout");
1851: #endif
1852: /* Convert to PetscSection+IS */
1853: for (v = 0; v < nvtxs; ++v) {
1854: PetscSectionAddDof(partSection, assignment[v], 1);
1855: }
1856: PetscMalloc1(nvtxs, &points);
1857: for (p = 0, i = 0; p < nparts; ++p) {
1858: for (v = 0; v < nvtxs; ++v) {
1859: if (assignment[v] == p) points[i++] = v;
1860: }
1861: }
1862: if (i != nvtxs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
1863: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
1864: if (global_method == INERTIAL_METHOD) {
1865: /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
1866: }
1867: PetscFree(assignment);
1868: for (i = 0; i < start[numVertices]; ++i) --adjacency[i];
1869: return(0);
1870: #else
1871: SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco.");
1872: #endif
1873: }
1875: static PetscErrorCode PetscPartitionerInitialize_Chaco(PetscPartitioner part)
1876: {
1878: part->noGraph = PETSC_FALSE;
1879: part->ops->view = PetscPartitionerView_Chaco;
1880: part->ops->destroy = PetscPartitionerDestroy_Chaco;
1881: part->ops->partition = PetscPartitionerPartition_Chaco;
1882: return(0);
1883: }
1885: /*MC
1886: PETSCPARTITIONERCHACO = "chaco" - A PetscPartitioner object using the Chaco library
1888: Level: intermediate
1890: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
1891: M*/
1893: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_Chaco(PetscPartitioner part)
1894: {
1895: PetscPartitioner_Chaco *p;
1896: PetscErrorCode ierr;
1900: PetscNewLog(part, &p);
1901: part->data = p;
1903: PetscPartitionerInitialize_Chaco(part);
1904: PetscCitationsRegister(ChacoPartitionerCitation, &ChacoPartitionercite);
1905: return(0);
1906: }
1908: static const char *ptypes[] = {"kway", "rb"};
1910: static PetscErrorCode PetscPartitionerDestroy_ParMetis(PetscPartitioner part)
1911: {
1912: PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1913: PetscErrorCode ierr;
1916: MPI_Comm_free(&p->pcomm);
1917: PetscFree(p);
1918: return(0);
1919: }
1921: static PetscErrorCode PetscPartitionerView_ParMetis_Ascii(PetscPartitioner part, PetscViewer viewer)
1922: {
1923: PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1924: PetscErrorCode ierr;
1927: PetscViewerASCIIPushTab(viewer);
1928: PetscViewerASCIIPrintf(viewer, "ParMetis type: %s\n", ptypes[p->ptype]);
1929: PetscViewerASCIIPrintf(viewer, "load imbalance ratio %g\n", (double) p->imbalanceRatio);
1930: PetscViewerASCIIPrintf(viewer, "debug flag %D\n", p->debugFlag);
1931: PetscViewerASCIIPrintf(viewer, "random seed %D\n", p->randomSeed);
1932: PetscViewerASCIIPopTab(viewer);
1933: return(0);
1934: }
1936: static PetscErrorCode PetscPartitionerView_ParMetis(PetscPartitioner part, PetscViewer viewer)
1937: {
1938: PetscBool iascii;
1944: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1945: if (iascii) {PetscPartitionerView_ParMetis_Ascii(part, viewer);}
1946: return(0);
1947: }
1949: static PetscErrorCode PetscPartitionerSetFromOptions_ParMetis(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
1950: {
1951: PetscPartitioner_ParMetis *p = (PetscPartitioner_ParMetis *) part->data;
1952: PetscErrorCode ierr;
1955: PetscOptionsHead(PetscOptionsObject, "PetscPartitioner ParMetis Options");
1956: PetscOptionsEList("-petscpartitioner_parmetis_type", "Partitioning method", "", ptypes, 2, ptypes[p->ptype], &p->ptype, NULL);
1957: PetscOptionsReal("-petscpartitioner_parmetis_imbalance_ratio", "Load imbalance ratio limit", "", p->imbalanceRatio, &p->imbalanceRatio, NULL);
1958: PetscOptionsInt("-petscpartitioner_parmetis_debug", "Debugging flag", "", p->debugFlag, &p->debugFlag, NULL);
1959: PetscOptionsInt("-petscpartitioner_parmetis_seed", "Random seed", "", p->randomSeed, &p->randomSeed, NULL);
1960: PetscOptionsTail();
1961: return(0);
1962: }
1964: #if defined(PETSC_HAVE_PARMETIS)
1965: #include <parmetis.h>
1966: #endif
1968: static PetscErrorCode PetscPartitionerPartition_ParMetis(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
1969: {
1970: #if defined(PETSC_HAVE_PARMETIS)
1971: PetscPartitioner_ParMetis *pm = (PetscPartitioner_ParMetis *) part->data;
1972: MPI_Comm comm;
1973: PetscInt nvtxs = numVertices; /* The number of vertices in full graph */
1974: PetscInt *vtxdist; /* Distribution of vertices across processes */
1975: PetscInt *xadj = start; /* Start of edge list for each vertex */
1976: PetscInt *adjncy = adjacency; /* Edge lists for all vertices */
1977: PetscInt *vwgt = NULL; /* Vertex weights */
1978: PetscInt *adjwgt = NULL; /* Edge weights */
1979: PetscInt wgtflag = 0; /* Indicates which weights are present */
1980: PetscInt numflag = 0; /* Indicates initial offset (0 or 1) */
1981: PetscInt ncon = 1; /* The number of weights per vertex */
1982: PetscInt metis_ptype = pm->ptype; /* kway or recursive bisection */
1983: real_t *tpwgts; /* The fraction of vertex weights assigned to each partition */
1984: real_t *ubvec; /* The balance intolerance for vertex weights */
1985: PetscInt options[64]; /* Options */
1986: PetscInt v, i, *assignment, *points;
1987: PetscMPIInt p, size, rank;
1988: PetscBool hasempty = PETSC_FALSE;
1992: PetscObjectGetComm((PetscObject) part, &comm);
1993: MPI_Comm_size(comm, &size);
1994: MPI_Comm_rank(comm, &rank);
1995: /* Calculate vertex distribution */
1996: PetscMalloc4(size+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment);
1997: vtxdist[0] = 0;
1998: MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
1999: for (p = 2; p <= size; ++p) {
2000: hasempty = (PetscBool)(hasempty || !vtxdist[p-1] || !vtxdist[p]);
2001: vtxdist[p] += vtxdist[p-1];
2002: }
2003: /* null graph */
2004: if (vtxdist[size] == 0) {
2005: PetscFree4(vtxdist,tpwgts,ubvec,assignment);
2006: ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
2007: return(0);
2008: }
2009: /* Calculate partition weights */
2010: if (targetSection) {
2011: PetscInt p;
2012: real_t sumt = 0.0;
2014: for (p = 0; p < nparts; ++p) {
2015: PetscInt tpd;
2017: PetscSectionGetDof(targetSection,p,&tpd);
2018: sumt += tpd;
2019: tpwgts[p] = tpd;
2020: }
2021: if (sumt) { /* METIS/ParMETIS do not like exactly zero weight */
2022: for (p = 0, sumt = 0.0; p < nparts; ++p) {
2023: tpwgts[p] = PetscMax(tpwgts[p],PETSC_SMALL);
2024: sumt += tpwgts[p];
2025: }
2026: for (p = 0; p < nparts; ++p) tpwgts[p] /= sumt;
2027: for (p = 0, sumt = 0.0; p < nparts-1; ++p) sumt += tpwgts[p];
2028: tpwgts[nparts - 1] = 1. - sumt;
2029: }
2030: } else {
2031: for (p = 0; p < nparts; ++p) tpwgts[p] = 1.0/nparts;
2032: }
2033: ubvec[0] = pm->imbalanceRatio;
2035: /* Weight cells */
2036: if (vertSection) {
2037: PetscMalloc1(nvtxs,&vwgt);
2038: for (v = 0; v < nvtxs; ++v) {
2039: PetscSectionGetDof(vertSection, v, &vwgt[v]);
2040: }
2041: wgtflag |= 2; /* have weights on graph vertices */
2042: }
2044: for (p = 0; !vtxdist[p+1] && p < size; ++p);
2045: if (vtxdist[p+1] == vtxdist[size]) {
2046: if (rank == p) {
2047: METIS_SetDefaultOptions(options); /* initialize all defaults */
2048: options[METIS_OPTION_DBGLVL] = pm->debugFlag;
2049: options[METIS_OPTION_SEED] = pm->randomSeed;
2050: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
2051: if (metis_ptype == 1) {
2052: PetscStackPush("METIS_PartGraphRecursive");
2053: METIS_PartGraphRecursive(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
2054: PetscStackPop;
2055: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphRecursive()");
2056: } else {
2057: /*
2058: It would be nice to activate the two options below, but they would need some actual testing.
2059: - Turning on these options may exercise path of the METIS code that have bugs and may break production runs.
2060: - 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.
2061: */
2062: /* options[METIS_OPTION_CONTIG] = 1; */ /* try to produce partitions that are contiguous */
2063: /* options[METIS_OPTION_MINCONN] = 1; */ /* minimize the maximum degree of the subdomain graph */
2064: PetscStackPush("METIS_PartGraphKway");
2065: METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment);
2066: PetscStackPop;
2067: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
2068: }
2069: }
2070: } else {
2071: MPI_Comm pcomm;
2073: options[0] = 1; /*use options */
2074: options[1] = pm->debugFlag;
2075: options[2] = (pm->randomSeed == -1) ? 15 : pm->randomSeed; /* default is GLOBAL_SEED=15 from `libparmetis/defs.h` */
2077: if (hasempty) { /* parmetis does not support empty graphs on some of the processes */
2078: PetscInt cnt;
2080: MPI_Comm_split(pm->pcomm,!!nvtxs,rank,&pcomm);
2081: for (p=0,cnt=0;p<size;p++) {
2082: if (vtxdist[p+1] != vtxdist[p]) {
2083: vtxdist[cnt+1] = vtxdist[p+1];
2084: cnt++;
2085: }
2086: }
2087: } else pcomm = pm->pcomm;
2088: if (nvtxs) {
2089: PetscStackPush("ParMETIS_V3_PartKway");
2090: ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &part->edgeCut, assignment, &pcomm);
2091: PetscStackPop;
2092: if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error %d in ParMETIS_V3_PartKway()", ierr);
2093: }
2094: if (hasempty) {
2095: MPI_Comm_free(&pcomm);
2096: }
2097: }
2099: /* Convert to PetscSection+IS */
2100: for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
2101: PetscMalloc1(nvtxs, &points);
2102: for (p = 0, i = 0; p < nparts; ++p) {
2103: for (v = 0; v < nvtxs; ++v) {
2104: if (assignment[v] == p) points[i++] = v;
2105: }
2106: }
2107: if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
2108: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
2109: PetscFree4(vtxdist,tpwgts,ubvec,assignment);
2110: PetscFree(vwgt);
2111: return(0);
2112: #else
2113: SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2114: #endif
2115: }
2117: static PetscErrorCode PetscPartitionerInitialize_ParMetis(PetscPartitioner part)
2118: {
2120: part->noGraph = PETSC_FALSE;
2121: part->ops->view = PetscPartitionerView_ParMetis;
2122: part->ops->setfromoptions = PetscPartitionerSetFromOptions_ParMetis;
2123: part->ops->destroy = PetscPartitionerDestroy_ParMetis;
2124: part->ops->partition = PetscPartitionerPartition_ParMetis;
2125: return(0);
2126: }
2128: /*MC
2129: PETSCPARTITIONERPARMETIS = "parmetis" - A PetscPartitioner object using the ParMETIS library
2131: Level: intermediate
2133: Options Database Keys:
2134: + -petscpartitioner_parmetis_type <string> - ParMETIS partitioning type. Either "kway" or "rb" (recursive bisection)
2135: . -petscpartitioner_parmetis_imbalance_ratio <value> - Load imbalance ratio limit
2136: . -petscpartitioner_parmetis_debug <int> - Debugging flag passed to ParMETIS/METIS routines
2137: - -petscpartitioner_parmetis_seed <int> - Random seed
2139: Notes: when the graph is on a single process, this partitioner actually calls METIS and not ParMETIS
2141: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
2142: M*/
2144: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_ParMetis(PetscPartitioner part)
2145: {
2146: PetscPartitioner_ParMetis *p;
2147: PetscErrorCode ierr;
2151: PetscNewLog(part, &p);
2152: part->data = p;
2154: MPI_Comm_dup(PetscObjectComm((PetscObject)part),&p->pcomm);
2155: p->ptype = 0;
2156: p->imbalanceRatio = 1.05;
2157: p->debugFlag = 0;
2158: p->randomSeed = -1; /* defaults to GLOBAL_SEED=15 from `libparmetis/defs.h` */
2160: PetscPartitionerInitialize_ParMetis(part);
2161: PetscCitationsRegister(ParMetisPartitionerCitation, &ParMetisPartitionercite);
2162: return(0);
2163: }
2165: #if defined(PETSC_HAVE_PTSCOTCH)
2167: EXTERN_C_BEGIN
2168: #include <ptscotch.h>
2169: EXTERN_C_END
2171: #define CHKERRPTSCOTCH(ierr) do { if (ierr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error calling PT-Scotch library"); } while(0)
2173: static int PTScotch_Strategy(PetscInt strategy)
2174: {
2175: switch (strategy) {
2176: case 0: return SCOTCH_STRATDEFAULT;
2177: case 1: return SCOTCH_STRATQUALITY;
2178: case 2: return SCOTCH_STRATSPEED;
2179: case 3: return SCOTCH_STRATBALANCE;
2180: case 4: return SCOTCH_STRATSAFETY;
2181: case 5: return SCOTCH_STRATSCALABILITY;
2182: case 6: return SCOTCH_STRATRECURSIVE;
2183: case 7: return SCOTCH_STRATREMAP;
2184: default: return SCOTCH_STRATDEFAULT;
2185: }
2186: }
2188: static PetscErrorCode PTScotch_PartGraph_Seq(SCOTCH_Num strategy, double imbalance, SCOTCH_Num n, SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
2189: SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num tpart[], SCOTCH_Num part[])
2190: {
2191: SCOTCH_Graph grafdat;
2192: SCOTCH_Strat stradat;
2193: SCOTCH_Num vertnbr = n;
2194: SCOTCH_Num edgenbr = xadj[n];
2195: SCOTCH_Num* velotab = vtxwgt;
2196: SCOTCH_Num* edlotab = adjwgt;
2197: SCOTCH_Num flagval = strategy;
2198: double kbalval = imbalance;
2202: {
2203: PetscBool flg = PETSC_TRUE;
2204: PetscOptionsDeprecatedNoObject("-petscpartititoner_ptscotch_vertex_weight",NULL,"3.13","Use -petscpartitioner_use_vertex_weights");
2205: PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
2206: if (!flg) velotab = NULL;
2207: }
2208: SCOTCH_graphInit(&grafdat);CHKERRPTSCOTCH(ierr);
2209: SCOTCH_graphBuild(&grafdat, 0, vertnbr, xadj, xadj + 1, velotab, NULL, edgenbr, adjncy, edlotab);CHKERRPTSCOTCH(ierr);
2210: SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
2211: SCOTCH_stratGraphMapBuild(&stradat, flagval, nparts, kbalval);CHKERRPTSCOTCH(ierr);
2212: if (tpart) {
2213: SCOTCH_Arch archdat;
2214: SCOTCH_archInit(&archdat);CHKERRPTSCOTCH(ierr);
2215: SCOTCH_archCmpltw(&archdat, nparts, tpart);CHKERRPTSCOTCH(ierr);
2216: SCOTCH_graphMap(&grafdat, &archdat, &stradat, part);CHKERRPTSCOTCH(ierr);
2217: SCOTCH_archExit(&archdat);
2218: } else {
2219: SCOTCH_graphPart(&grafdat, nparts, &stradat, part);CHKERRPTSCOTCH(ierr);
2220: }
2221: SCOTCH_stratExit(&stradat);
2222: SCOTCH_graphExit(&grafdat);
2223: return(0);
2224: }
2226: static PetscErrorCode PTScotch_PartGraph_MPI(SCOTCH_Num strategy, double imbalance, SCOTCH_Num vtxdist[], SCOTCH_Num xadj[], SCOTCH_Num adjncy[],
2227: SCOTCH_Num vtxwgt[], SCOTCH_Num adjwgt[], SCOTCH_Num nparts, SCOTCH_Num tpart[], SCOTCH_Num part[], MPI_Comm comm)
2228: {
2229: PetscMPIInt procglbnbr;
2230: PetscMPIInt proclocnum;
2231: SCOTCH_Arch archdat;
2232: SCOTCH_Dgraph grafdat;
2233: SCOTCH_Dmapping mappdat;
2234: SCOTCH_Strat stradat;
2235: SCOTCH_Num vertlocnbr;
2236: SCOTCH_Num edgelocnbr;
2237: SCOTCH_Num* veloloctab = vtxwgt;
2238: SCOTCH_Num* edloloctab = adjwgt;
2239: SCOTCH_Num flagval = strategy;
2240: double kbalval = imbalance;
2241: PetscErrorCode ierr;
2244: {
2245: PetscBool flg = PETSC_TRUE;
2246: PetscOptionsDeprecatedNoObject("-petscpartititoner_ptscotch_vertex_weight",NULL,"3.13","Use -petscpartitioner_use_vertex_weights");
2247: PetscOptionsGetBool(NULL, NULL, "-petscpartititoner_ptscotch_vertex_weight", &flg, NULL);
2248: if (!flg) veloloctab = NULL;
2249: }
2250: MPI_Comm_size(comm, &procglbnbr);
2251: MPI_Comm_rank(comm, &proclocnum);
2252: vertlocnbr = vtxdist[proclocnum + 1] - vtxdist[proclocnum];
2253: edgelocnbr = xadj[vertlocnbr];
2255: SCOTCH_dgraphInit(&grafdat, comm);CHKERRPTSCOTCH(ierr);
2256: SCOTCH_dgraphBuild(&grafdat, 0, vertlocnbr, vertlocnbr, xadj, xadj + 1, veloloctab, NULL, edgelocnbr, edgelocnbr, adjncy, NULL, edloloctab);CHKERRPTSCOTCH(ierr);
2257: SCOTCH_stratInit(&stradat);CHKERRPTSCOTCH(ierr);
2258: SCOTCH_stratDgraphMapBuild(&stradat, flagval, procglbnbr, nparts, kbalval);
2259: SCOTCH_archInit(&archdat);CHKERRPTSCOTCH(ierr);
2260: if (tpart) { /* target partition weights */
2261: SCOTCH_archCmpltw(&archdat, nparts, tpart);CHKERRPTSCOTCH(ierr);
2262: } else {
2263: SCOTCH_archCmplt(&archdat, nparts);CHKERRPTSCOTCH(ierr);
2264: }
2265: SCOTCH_dgraphMapInit(&grafdat, &mappdat, &archdat, part);CHKERRPTSCOTCH(ierr);
2267: SCOTCH_dgraphMapCompute(&grafdat, &mappdat, &stradat);CHKERRPTSCOTCH(ierr);
2268: SCOTCH_dgraphMapExit(&grafdat, &mappdat);
2269: SCOTCH_archExit(&archdat);
2270: SCOTCH_stratExit(&stradat);
2271: SCOTCH_dgraphExit(&grafdat);
2272: return(0);
2273: }
2275: #endif /* PETSC_HAVE_PTSCOTCH */
2277: static PetscErrorCode PetscPartitionerDestroy_PTScotch(PetscPartitioner part)
2278: {
2279: PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2280: PetscErrorCode ierr;
2283: MPI_Comm_free(&p->pcomm);
2284: PetscFree(p);
2285: return(0);
2286: }
2288: static PetscErrorCode PetscPartitionerView_PTScotch_Ascii(PetscPartitioner part, PetscViewer viewer)
2289: {
2290: PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2291: PetscErrorCode ierr;
2294: PetscViewerASCIIPushTab(viewer);
2295: PetscViewerASCIIPrintf(viewer,"using partitioning strategy %s\n",PTScotchStrategyList[p->strategy]);
2296: PetscViewerASCIIPrintf(viewer,"using load imbalance ratio %g\n",(double)p->imbalance);
2297: PetscViewerASCIIPopTab(viewer);
2298: return(0);
2299: }
2301: static PetscErrorCode PetscPartitionerView_PTScotch(PetscPartitioner part, PetscViewer viewer)
2302: {
2303: PetscBool iascii;
2309: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2310: if (iascii) {PetscPartitionerView_PTScotch_Ascii(part, viewer);}
2311: return(0);
2312: }
2314: static PetscErrorCode PetscPartitionerSetFromOptions_PTScotch(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
2315: {
2316: PetscPartitioner_PTScotch *p = (PetscPartitioner_PTScotch *) part->data;
2317: const char *const *slist = PTScotchStrategyList;
2318: PetscInt nlist = (PetscInt)(sizeof(PTScotchStrategyList)/sizeof(PTScotchStrategyList[0]));
2319: PetscBool flag;
2320: PetscErrorCode ierr;
2323: PetscOptionsHead(PetscOptionsObject, "PetscPartitioner PTScotch Options");
2324: PetscOptionsEList("-petscpartitioner_ptscotch_strategy","Partitioning strategy","",slist,nlist,slist[p->strategy],&p->strategy,&flag);
2325: PetscOptionsReal("-petscpartitioner_ptscotch_imbalance","Load imbalance ratio","",p->imbalance,&p->imbalance,&flag);
2326: PetscOptionsTail();
2327: return(0);
2328: }
2330: static PetscErrorCode PetscPartitionerPartition_PTScotch(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *partition)
2331: {
2332: #if defined(PETSC_HAVE_PTSCOTCH)
2333: MPI_Comm comm;
2334: PetscInt nvtxs = numVertices; /* The number of vertices in full graph */
2335: PetscInt *vtxdist; /* Distribution of vertices across processes */
2336: PetscInt *xadj = start; /* Start of edge list for each vertex */
2337: PetscInt *adjncy = adjacency; /* Edge lists for all vertices */
2338: PetscInt *vwgt = NULL; /* Vertex weights */
2339: PetscInt *adjwgt = NULL; /* Edge weights */
2340: PetscInt v, i, *assignment, *points;
2341: PetscMPIInt size, rank, p;
2342: PetscBool hasempty = PETSC_FALSE;
2343: PetscInt *tpwgts = NULL;
2347: PetscObjectGetComm((PetscObject)part,&comm);
2348: MPI_Comm_size(comm, &size);
2349: MPI_Comm_rank(comm, &rank);
2350: PetscMalloc2(size+1,&vtxdist,PetscMax(nvtxs,1),&assignment);
2351: /* Calculate vertex distribution */
2352: vtxdist[0] = 0;
2353: MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);
2354: for (p = 2; p <= size; ++p) {
2355: hasempty = (PetscBool)(hasempty || !vtxdist[p-1] || !vtxdist[p]);
2356: vtxdist[p] += vtxdist[p-1];
2357: }
2358: /* null graph */
2359: if (vtxdist[size] == 0) {
2360: PetscFree2(vtxdist, assignment);
2361: ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);
2362: return(0);
2363: }
2365: /* Calculate vertex weights */
2366: if (vertSection) {
2367: PetscMalloc1(nvtxs,&vwgt);
2368: for (v = 0; v < nvtxs; ++v) {
2369: PetscSectionGetDof(vertSection, v, &vwgt[v]);
2370: }
2371: }
2373: /* Calculate partition weights */
2374: if (targetSection) {
2375: PetscInt sumw;
2377: PetscCalloc1(nparts,&tpwgts);
2378: for (p = 0, sumw = 0; p < nparts; ++p) {
2379: PetscSectionGetDof(targetSection,p,&tpwgts[p]);
2380: sumw += tpwgts[p];
2381: }
2382: if (!sumw) {
2383: PetscFree(tpwgts);
2384: }
2385: }
2387: {
2388: PetscPartitioner_PTScotch *pts = (PetscPartitioner_PTScotch *) part->data;
2389: int strat = PTScotch_Strategy(pts->strategy);
2390: double imbal = (double)pts->imbalance;
2392: for (p = 0; !vtxdist[p+1] && p < size; ++p);
2393: if (vtxdist[p+1] == vtxdist[size]) {
2394: if (rank == p) {
2395: PTScotch_PartGraph_Seq(strat, imbal, nvtxs, xadj, adjncy, vwgt, adjwgt, nparts, tpwgts, assignment);
2396: }
2397: } else {
2398: PetscInt cnt;
2399: MPI_Comm pcomm;
2401: if (hasempty) {
2402: MPI_Comm_split(pts->pcomm,!!nvtxs,rank,&pcomm);
2403: for (p=0,cnt=0;p<size;p++) {
2404: if (vtxdist[p+1] != vtxdist[p]) {
2405: vtxdist[cnt+1] = vtxdist[p+1];
2406: cnt++;
2407: }
2408: }
2409: } else pcomm = pts->pcomm;
2410: if (nvtxs) {
2411: PTScotch_PartGraph_MPI(strat, imbal, vtxdist, xadj, adjncy, vwgt, adjwgt, nparts, tpwgts, assignment, pcomm);
2412: }
2413: if (hasempty) {
2414: MPI_Comm_free(&pcomm);
2415: }
2416: }
2417: }
2418: PetscFree(vwgt);
2419: PetscFree(tpwgts);
2421: /* Convert to PetscSection+IS */
2422: for (v = 0; v < nvtxs; ++v) {PetscSectionAddDof(partSection, assignment[v], 1);}
2423: PetscMalloc1(nvtxs, &points);
2424: for (p = 0, i = 0; p < nparts; ++p) {
2425: for (v = 0; v < nvtxs; ++v) {
2426: if (assignment[v] == p) points[i++] = v;
2427: }
2428: }
2429: if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
2430: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
2432: PetscFree2(vtxdist,assignment);
2433: return(0);
2434: #else
2435: SETERRQ(PetscObjectComm((PetscObject) part), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-ptscotch.");
2436: #endif
2437: }
2439: static PetscErrorCode PetscPartitionerInitialize_PTScotch(PetscPartitioner part)
2440: {
2442: part->noGraph = PETSC_FALSE;
2443: part->ops->view = PetscPartitionerView_PTScotch;
2444: part->ops->destroy = PetscPartitionerDestroy_PTScotch;
2445: part->ops->partition = PetscPartitionerPartition_PTScotch;
2446: part->ops->setfromoptions = PetscPartitionerSetFromOptions_PTScotch;
2447: return(0);
2448: }
2450: /*MC
2451: PETSCPARTITIONERPTSCOTCH = "ptscotch" - A PetscPartitioner object using the PT-Scotch library
2453: Level: intermediate
2455: Options Database Keys:
2456: + -petscpartitioner_ptscotch_strategy <string> - PT-Scotch strategy. Choose one of default quality speed balance safety scalability recursive remap
2457: - -petscpartitioner_ptscotch_imbalance <val> - Load imbalance ratio
2459: Notes: when the graph is on a single process, this partitioner actually uses Scotch and not PT-Scotch
2461: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
2462: M*/
2464: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_PTScotch(PetscPartitioner part)
2465: {
2466: PetscPartitioner_PTScotch *p;
2467: PetscErrorCode ierr;
2471: PetscNewLog(part, &p);
2472: part->data = p;
2474: MPI_Comm_dup(PetscObjectComm((PetscObject)part),&p->pcomm);
2475: p->strategy = 0;
2476: p->imbalance = 0.01;
2478: PetscPartitionerInitialize_PTScotch(part);
2479: PetscCitationsRegister(PTScotchPartitionerCitation, &PTScotchPartitionercite);
2480: return(0);
2481: }
2483: /*@
2484: DMPlexGetPartitioner - Get the mesh partitioner
2486: Not collective
2488: Input Parameter:
2489: . dm - The DM
2491: Output Parameter:
2492: . part - The PetscPartitioner
2494: Level: developer
2496: Note: This gets a borrowed reference, so the user should not destroy this PetscPartitioner.
2498: .seealso DMPlexDistribute(), DMPlexSetPartitioner(), PetscPartitionerDMPlexPartition(), PetscPartitionerCreate()
2499: @*/
2500: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
2501: {
2502: DM_Plex *mesh = (DM_Plex *) dm->data;
2507: *part = mesh->partitioner;
2508: return(0);
2509: }
2511: /*@
2512: DMPlexSetPartitioner - Set the mesh partitioner
2514: logically collective on DM
2516: Input Parameters:
2517: + dm - The DM
2518: - part - The partitioner
2520: Level: developer
2522: Note: Any existing PetscPartitioner will be destroyed.
2524: .seealso DMPlexDistribute(), DMPlexGetPartitioner(), PetscPartitionerCreate()
2525: @*/
2526: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
2527: {
2528: DM_Plex *mesh = (DM_Plex *) dm->data;
2534: PetscObjectReference((PetscObject)part);
2535: PetscPartitionerDestroy(&mesh->partitioner);
2536: mesh->partitioner = part;
2537: return(0);
2538: }
2540: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
2541: {
2542: const PetscInt *cone;
2543: PetscInt coneSize, c;
2544: PetscBool missing;
2548: PetscHSetIQueryAdd(ht, point, &missing);
2549: if (missing) {
2550: DMPlexGetCone(dm, point, &cone);
2551: DMPlexGetConeSize(dm, point, &coneSize);
2552: for (c = 0; c < coneSize; c++) {
2553: DMPlexAddClosure_Private(dm, ht, cone[c]);
2554: }
2555: }
2556: return(0);
2557: }
2559: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
2560: {
2564: if (up) {
2565: PetscInt parent;
2567: DMPlexGetTreeParent(dm,point,&parent,NULL);
2568: if (parent != point) {
2569: PetscInt closureSize, *closure = NULL, i;
2571: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
2572: for (i = 0; i < closureSize; i++) {
2573: PetscInt cpoint = closure[2*i];
2575: PetscHSetIAdd(ht, cpoint);
2576: DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
2577: }
2578: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
2579: }
2580: }
2581: if (down) {
2582: PetscInt numChildren;
2583: const PetscInt *children;
2585: DMPlexGetTreeChildren(dm,point,&numChildren,&children);
2586: if (numChildren) {
2587: PetscInt i;
2589: for (i = 0; i < numChildren; i++) {
2590: PetscInt cpoint = children[i];
2592: PetscHSetIAdd(ht, cpoint);
2593: DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
2594: }
2595: }
2596: }
2597: return(0);
2598: }
2600: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
2601: {
2602: PetscInt parent;
2606: DMPlexGetTreeParent(dm, point, &parent,NULL);
2607: if (point != parent) {
2608: const PetscInt *cone;
2609: PetscInt coneSize, c;
2611: DMPlexAddClosureTree_Up_Private(dm, ht, parent);
2612: DMPlexAddClosure_Private(dm, ht, parent);
2613: DMPlexGetCone(dm, parent, &cone);
2614: DMPlexGetConeSize(dm, parent, &coneSize);
2615: for (c = 0; c < coneSize; c++) {
2616: const PetscInt cp = cone[c];
2618: DMPlexAddClosureTree_Up_Private(dm, ht, cp);
2619: }
2620: }
2621: return(0);
2622: }
2624: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
2625: {
2626: PetscInt i, numChildren;
2627: const PetscInt *children;
2631: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
2632: for (i = 0; i < numChildren; i++) {
2633: PetscHSetIAdd(ht, children[i]);
2634: }
2635: return(0);
2636: }
2638: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
2639: {
2640: const PetscInt *cone;
2641: PetscInt coneSize, c;
2645: PetscHSetIAdd(ht, point);
2646: DMPlexAddClosureTree_Up_Private(dm, ht, point);
2647: DMPlexAddClosureTree_Down_Private(dm, ht, point);
2648: DMPlexGetCone(dm, point, &cone);
2649: DMPlexGetConeSize(dm, point, &coneSize);
2650: for (c = 0; c < coneSize; c++) {
2651: DMPlexAddClosureTree_Private(dm, ht, cone[c]);
2652: }
2653: return(0);
2654: }
2656: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
2657: {
2658: DM_Plex *mesh = (DM_Plex *)dm->data;
2659: const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
2660: PetscInt nelems, *elems, off = 0, p;
2661: PetscHSetI ht;
2662: PetscErrorCode ierr;
2665: PetscHSetICreate(&ht);
2666: PetscHSetIResize(ht, numPoints*16);
2667: if (!hasTree) {
2668: for (p = 0; p < numPoints; ++p) {
2669: DMPlexAddClosure_Private(dm, ht, points[p]);
2670: }
2671: } else {
2672: #if 1
2673: for (p = 0; p < numPoints; ++p) {
2674: DMPlexAddClosureTree_Private(dm, ht, points[p]);
2675: }
2676: #else
2677: PetscInt *closure = NULL, closureSize, c;
2678: for (p = 0; p < numPoints; ++p) {
2679: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
2680: for (c = 0; c < closureSize*2; c += 2) {
2681: PetscHSetIAdd(ht, closure[c]);
2682: if (hasTree) {DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);}
2683: }
2684: }
2685: if (closure) {DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);}
2686: #endif
2687: }
2688: PetscHSetIGetSize(ht, &nelems);
2689: PetscMalloc1(nelems, &elems);
2690: PetscHSetIGetElems(ht, &off, elems);
2691: PetscHSetIDestroy(&ht);
2692: PetscSortInt(nelems, elems);
2693: ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
2694: return(0);
2695: }
2697: /*@
2698: DMPlexPartitionLabelClosure - Add the closure of all points to the partition label
2700: Input Parameters:
2701: + dm - The DM
2702: - label - DMLabel assinging ranks to remote roots
2704: Level: developer
2706: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
2707: @*/
2708: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
2709: {
2710: IS rankIS, pointIS, closureIS;
2711: const PetscInt *ranks, *points;
2712: PetscInt numRanks, numPoints, r;
2713: PetscErrorCode ierr;
2716: DMLabelGetValueIS(label, &rankIS);
2717: ISGetLocalSize(rankIS, &numRanks);
2718: ISGetIndices(rankIS, &ranks);
2719: for (r = 0; r < numRanks; ++r) {
2720: const PetscInt rank = ranks[r];
2721: DMLabelGetStratumIS(label, rank, &pointIS);
2722: ISGetLocalSize(pointIS, &numPoints);
2723: ISGetIndices(pointIS, &points);
2724: DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
2725: ISRestoreIndices(pointIS, &points);
2726: ISDestroy(&pointIS);
2727: DMLabelSetStratumIS(label, rank, closureIS);
2728: ISDestroy(&closureIS);
2729: }
2730: ISRestoreIndices(rankIS, &ranks);
2731: ISDestroy(&rankIS);
2732: return(0);
2733: }
2735: /*@
2736: DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label
2738: Input Parameters:
2739: + dm - The DM
2740: - label - DMLabel assinging ranks to remote roots
2742: Level: developer
2744: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
2745: @*/
2746: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
2747: {
2748: IS rankIS, pointIS;
2749: const PetscInt *ranks, *points;
2750: PetscInt numRanks, numPoints, r, p, a, adjSize;
2751: PetscInt *adj = NULL;
2752: PetscErrorCode ierr;
2755: DMLabelGetValueIS(label, &rankIS);
2756: ISGetLocalSize(rankIS, &numRanks);
2757: ISGetIndices(rankIS, &ranks);
2758: for (r = 0; r < numRanks; ++r) {
2759: const PetscInt rank = ranks[r];
2761: DMLabelGetStratumIS(label, rank, &pointIS);
2762: ISGetLocalSize(pointIS, &numPoints);
2763: ISGetIndices(pointIS, &points);
2764: for (p = 0; p < numPoints; ++p) {
2765: adjSize = PETSC_DETERMINE;
2766: DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
2767: for (a = 0; a < adjSize; ++a) {DMLabelSetValue(label, adj[a], rank);}
2768: }
2769: ISRestoreIndices(pointIS, &points);
2770: ISDestroy(&pointIS);
2771: }
2772: ISRestoreIndices(rankIS, &ranks);
2773: ISDestroy(&rankIS);
2774: PetscFree(adj);
2775: return(0);
2776: }
2778: /*@
2779: DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point SF
2781: Input Parameters:
2782: + dm - The DM
2783: - label - DMLabel assinging ranks to remote roots
2785: Level: developer
2787: Note: This is required when generating multi-level overlaps to capture
2788: overlap points from non-neighbouring partitions.
2790: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
2791: @*/
2792: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
2793: {
2794: MPI_Comm comm;
2795: PetscMPIInt rank;
2796: PetscSF sfPoint;
2797: DMLabel lblRoots, lblLeaves;
2798: IS rankIS, pointIS;
2799: const PetscInt *ranks;
2800: PetscInt numRanks, r;
2801: PetscErrorCode ierr;
2804: PetscObjectGetComm((PetscObject) dm, &comm);
2805: MPI_Comm_rank(comm, &rank);
2806: DMGetPointSF(dm, &sfPoint);
2807: /* Pull point contributions from remote leaves into local roots */
2808: DMLabelGather(label, sfPoint, &lblLeaves);
2809: DMLabelGetValueIS(lblLeaves, &rankIS);
2810: ISGetLocalSize(rankIS, &numRanks);
2811: ISGetIndices(rankIS, &ranks);
2812: for (r = 0; r < numRanks; ++r) {
2813: const PetscInt remoteRank = ranks[r];
2814: if (remoteRank == rank) continue;
2815: DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
2816: DMLabelInsertIS(label, pointIS, remoteRank);
2817: ISDestroy(&pointIS);
2818: }
2819: ISRestoreIndices(rankIS, &ranks);
2820: ISDestroy(&rankIS);
2821: DMLabelDestroy(&lblLeaves);
2822: /* Push point contributions from roots into remote leaves */
2823: DMLabelDistribute(label, sfPoint, &lblRoots);
2824: DMLabelGetValueIS(lblRoots, &rankIS);
2825: ISGetLocalSize(rankIS, &numRanks);
2826: ISGetIndices(rankIS, &ranks);
2827: for (r = 0; r < numRanks; ++r) {
2828: const PetscInt remoteRank = ranks[r];
2829: if (remoteRank == rank) continue;
2830: DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
2831: DMLabelInsertIS(label, pointIS, remoteRank);
2832: ISDestroy(&pointIS);
2833: }
2834: ISRestoreIndices(rankIS, &ranks);
2835: ISDestroy(&rankIS);
2836: DMLabelDestroy(&lblRoots);
2837: return(0);
2838: }
2840: /*@
2841: DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label
2843: Input Parameters:
2844: + dm - The DM
2845: . rootLabel - DMLabel assinging ranks to local roots
2846: - processSF - A star forest mapping into the local index on each remote rank
2848: Output Parameter:
2849: . leafLabel - DMLabel assinging ranks to remote roots
2851: Note: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
2852: resulting leafLabel is a receiver mapping of remote roots to their parent rank.
2854: Level: developer
2856: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
2857: @*/
2858: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
2859: {
2860: MPI_Comm comm;
2861: PetscMPIInt rank, size, r;
2862: PetscInt p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
2863: PetscSF sfPoint;
2864: PetscSection rootSection;
2865: PetscSFNode *rootPoints, *leafPoints;
2866: const PetscSFNode *remote;
2867: const PetscInt *local, *neighbors;
2868: IS valueIS;
2869: PetscBool mpiOverflow = PETSC_FALSE;
2870: PetscErrorCode ierr;
2873: PetscLogEventBegin(DMPLEX_PartLabelInvert,dm,0,0,0);
2874: PetscObjectGetComm((PetscObject) dm, &comm);
2875: MPI_Comm_rank(comm, &rank);
2876: MPI_Comm_size(comm, &size);
2877: DMGetPointSF(dm, &sfPoint);
2879: /* Convert to (point, rank) and use actual owners */
2880: PetscSectionCreate(comm, &rootSection);
2881: PetscSectionSetChart(rootSection, 0, size);
2882: DMLabelGetValueIS(rootLabel, &valueIS);
2883: ISGetLocalSize(valueIS, &numNeighbors);
2884: ISGetIndices(valueIS, &neighbors);
2885: for (n = 0; n < numNeighbors; ++n) {
2886: DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
2887: PetscSectionAddDof(rootSection, neighbors[n], numPoints);
2888: }
2889: PetscSectionSetUp(rootSection);
2890: PetscSectionGetStorageSize(rootSection, &rootSize);
2891: PetscMalloc1(rootSize, &rootPoints);
2892: PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
2893: for (n = 0; n < numNeighbors; ++n) {
2894: IS pointIS;
2895: const PetscInt *points;
2897: PetscSectionGetOffset(rootSection, neighbors[n], &off);
2898: DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
2899: ISGetLocalSize(pointIS, &numPoints);
2900: ISGetIndices(pointIS, &points);
2901: for (p = 0; p < numPoints; ++p) {
2902: if (local) {PetscFindInt(points[p], nleaves, local, &l);}
2903: else {l = -1;}
2904: if (l >= 0) {rootPoints[off+p] = remote[l];}
2905: else {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
2906: }
2907: ISRestoreIndices(pointIS, &points);
2908: ISDestroy(&pointIS);
2909: }
2911: /* Try to communicate overlap using All-to-All */
2912: if (!processSF) {
2913: PetscInt64 counter = 0;
2914: PetscBool locOverflow = PETSC_FALSE;
2915: PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;
2917: PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
2918: for (n = 0; n < numNeighbors; ++n) {
2919: PetscSectionGetDof(rootSection, neighbors[n], &dof);
2920: PetscSectionGetOffset(rootSection, neighbors[n], &off);
2921: #if defined(PETSC_USE_64BIT_INDICES)
2922: if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2923: if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
2924: #endif
2925: scounts[neighbors[n]] = (PetscMPIInt) dof;
2926: sdispls[neighbors[n]] = (PetscMPIInt) off;
2927: }
2928: MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
2929: for (r = 0; r < size; ++r) { rdispls[r] = (int)counter; counter += rcounts[r]; }
2930: if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
2931: MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
2932: if (!mpiOverflow) {
2933: PetscInfo(dm,"Using Alltoallv for mesh distribution\n");
2934: leafSize = (PetscInt) counter;
2935: PetscMalloc1(leafSize, &leafPoints);
2936: MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
2937: }
2938: PetscFree4(scounts, sdispls, rcounts, rdispls);
2939: }
2941: /* Communicate overlap using process star forest */
2942: if (processSF || mpiOverflow) {
2943: PetscSF procSF;
2944: PetscSection leafSection;
2946: if (processSF) {
2947: PetscInfo(dm,"Using processSF for mesh distribution\n");
2948: PetscObjectReference((PetscObject)processSF);
2949: procSF = processSF;
2950: } else {
2951: PetscInfo(dm,"Using processSF for mesh distribution (MPI overflow)\n");
2952: PetscSFCreate(comm,&procSF);
2953: PetscSFSetGraphWithPattern(procSF,NULL,PETSCSF_PATTERN_ALLTOALL);
2954: }
2956: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
2957: DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
2958: PetscSectionGetStorageSize(leafSection, &leafSize);
2959: PetscSectionDestroy(&leafSection);
2960: PetscSFDestroy(&procSF);
2961: }
2963: for (p = 0; p < leafSize; p++) {
2964: DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
2965: }
2967: ISRestoreIndices(valueIS, &neighbors);
2968: ISDestroy(&valueIS);
2969: PetscSectionDestroy(&rootSection);
2970: PetscFree(rootPoints);
2971: PetscFree(leafPoints);
2972: PetscLogEventEnd(DMPLEX_PartLabelInvert,dm,0,0,0);
2973: return(0);
2974: }
2976: /*@
2977: DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points
2979: Input Parameters:
2980: + dm - The DM
2981: - label - DMLabel assinging ranks to remote roots
2983: Output Parameter:
2984: . sf - The star forest communication context encapsulating the defined mapping
2986: Note: The incoming label is a receiver mapping of remote points to their parent rank.
2988: Level: developer
2990: .seealso: DMPlexDistribute(), DMPlexCreateOverlap()
2991: @*/
2992: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
2993: {
2994: PetscMPIInt rank;
2995: PetscInt n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
2996: PetscSFNode *remotePoints;
2997: IS remoteRootIS, neighborsIS;
2998: const PetscInt *remoteRoots, *neighbors;
3002: PetscLogEventBegin(DMPLEX_PartLabelCreateSF,dm,0,0,0);
3003: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3005: DMLabelGetValueIS(label, &neighborsIS);
3006: #if 0
3007: {
3008: IS is;
3009: ISDuplicate(neighborsIS, &is);
3010: ISSort(is);
3011: ISDestroy(&neighborsIS);
3012: neighborsIS = is;
3013: }
3014: #endif
3015: ISGetLocalSize(neighborsIS, &nNeighbors);
3016: ISGetIndices(neighborsIS, &neighbors);
3017: for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
3018: DMLabelGetStratumSize(label, neighbors[n], &numPoints);
3019: numRemote += numPoints;
3020: }
3021: PetscMalloc1(numRemote, &remotePoints);
3022: /* Put owned points first */
3023: DMLabelGetStratumSize(label, rank, &numPoints);
3024: if (numPoints > 0) {
3025: DMLabelGetStratumIS(label, rank, &remoteRootIS);
3026: ISGetIndices(remoteRootIS, &remoteRoots);
3027: for (p = 0; p < numPoints; p++) {
3028: remotePoints[idx].index = remoteRoots[p];
3029: remotePoints[idx].rank = rank;
3030: idx++;
3031: }
3032: ISRestoreIndices(remoteRootIS, &remoteRoots);
3033: ISDestroy(&remoteRootIS);
3034: }
3035: /* Now add remote points */
3036: for (n = 0; n < nNeighbors; ++n) {
3037: const PetscInt nn = neighbors[n];
3039: DMLabelGetStratumSize(label, nn, &numPoints);
3040: if (nn == rank || numPoints <= 0) continue;
3041: DMLabelGetStratumIS(label, nn, &remoteRootIS);
3042: ISGetIndices(remoteRootIS, &remoteRoots);
3043: for (p = 0; p < numPoints; p++) {
3044: remotePoints[idx].index = remoteRoots[p];
3045: remotePoints[idx].rank = nn;
3046: idx++;
3047: }
3048: ISRestoreIndices(remoteRootIS, &remoteRoots);
3049: ISDestroy(&remoteRootIS);
3050: }
3051: PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
3052: DMPlexGetChart(dm, &pStart, &pEnd);
3053: PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
3054: PetscSFSetUp(*sf);
3055: ISDestroy(&neighborsIS);
3056: PetscLogEventEnd(DMPLEX_PartLabelCreateSF,dm,0,0,0);
3057: return(0);
3058: }
3060: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
3061: * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
3062: * them out in that case. */
3063: #if defined(PETSC_HAVE_PARMETIS)
3064: /*@C
3066: DMPlexRewriteSF - Rewrites the ownership of the SF of a DM (in place).
3068: Input parameters:
3069: + dm - The DMPlex object.
3070: . n - The number of points.
3071: . pointsToRewrite - The points in the SF whose ownership will change.
3072: . targetOwners - New owner for each element in pointsToRewrite.
3073: - degrees - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.
3075: Level: developer
3077: @*/
3078: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
3079: {
3080: PetscInt ierr, pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
3081: PetscInt *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
3082: PetscSFNode *leafLocationsNew;
3083: const PetscSFNode *iremote;
3084: const PetscInt *ilocal;
3085: PetscBool *isLeaf;
3086: PetscSF sf;
3087: MPI_Comm comm;
3088: PetscMPIInt rank, size;
3091: PetscObjectGetComm((PetscObject) dm, &comm);
3092: MPI_Comm_rank(comm, &rank);
3093: MPI_Comm_size(comm, &size);
3094: DMPlexGetChart(dm, &pStart, &pEnd);
3096: DMGetPointSF(dm, &sf);
3097: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
3098: PetscMalloc1(pEnd-pStart, &isLeaf);
3099: for (i=0; i<pEnd-pStart; i++) {
3100: isLeaf[i] = PETSC_FALSE;
3101: }
3102: for (i=0; i<nleafs; i++) {
3103: isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
3104: }
3106: PetscMalloc1(pEnd-pStart+1, &cumSumDegrees);
3107: cumSumDegrees[0] = 0;
3108: for (i=1; i<=pEnd-pStart; i++) {
3109: cumSumDegrees[i] = cumSumDegrees[i-1] + degrees[i-1];
3110: }
3111: sumDegrees = cumSumDegrees[pEnd-pStart];
3112: /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */
3114: PetscMalloc1(sumDegrees, &locationsOfLeafs);
3115: PetscMalloc1(pEnd-pStart, &rankOnLeafs);
3116: for (i=0; i<pEnd-pStart; i++) {
3117: rankOnLeafs[i] = rank;
3118: }
3119: PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
3120: PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
3121: PetscFree(rankOnLeafs);
3123: /* get the remote local points of my leaves */
3124: PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
3125: PetscMalloc1(pEnd-pStart, &points);
3126: for (i=0; i<pEnd-pStart; i++) {
3127: points[i] = pStart+i;
3128: }
3129: PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
3130: PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
3131: PetscFree(points);
3132: /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
3133: PetscMalloc1(pEnd-pStart, &newOwners);
3134: PetscMalloc1(pEnd-pStart, &newNumbers);
3135: for (i=0; i<pEnd-pStart; i++) {
3136: newOwners[i] = -1;
3137: newNumbers[i] = -1;
3138: }
3139: {
3140: PetscInt oldNumber, newNumber, oldOwner, newOwner;
3141: for (i=0; i<n; i++) {
3142: oldNumber = pointsToRewrite[i];
3143: newNumber = -1;
3144: oldOwner = rank;
3145: newOwner = targetOwners[i];
3146: if (oldOwner == newOwner) {
3147: newNumber = oldNumber;
3148: } else {
3149: for (j=0; j<degrees[oldNumber]; j++) {
3150: if (locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) {
3151: newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j];
3152: break;
3153: }
3154: }
3155: }
3156: if (newNumber == -1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");
3158: newOwners[oldNumber] = newOwner;
3159: newNumbers[oldNumber] = newNumber;
3160: }
3161: }
3162: PetscFree(cumSumDegrees);
3163: PetscFree(locationsOfLeafs);
3164: PetscFree(remoteLocalPointOfLeafs);
3166: PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners);
3167: PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners);
3168: PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers);
3169: PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers);
3171: /* Now count how many leafs we have on each processor. */
3172: leafCounter=0;
3173: for (i=0; i<pEnd-pStart; i++) {
3174: if (newOwners[i] >= 0) {
3175: if (newOwners[i] != rank) {
3176: leafCounter++;
3177: }
3178: } else {
3179: if (isLeaf[i]) {
3180: leafCounter++;
3181: }
3182: }
3183: }
3185: /* Now set up the new sf by creating the leaf arrays */
3186: PetscMalloc1(leafCounter, &leafsNew);
3187: PetscMalloc1(leafCounter, &leafLocationsNew);
3189: leafCounter = 0;
3190: counter = 0;
3191: for (i=0; i<pEnd-pStart; i++) {
3192: if (newOwners[i] >= 0) {
3193: if (newOwners[i] != rank) {
3194: leafsNew[leafCounter] = i;
3195: leafLocationsNew[leafCounter].rank = newOwners[i];
3196: leafLocationsNew[leafCounter].index = newNumbers[i];
3197: leafCounter++;
3198: }
3199: } else {
3200: if (isLeaf[i]) {
3201: leafsNew[leafCounter] = i;
3202: leafLocationsNew[leafCounter].rank = iremote[counter].rank;
3203: leafLocationsNew[leafCounter].index = iremote[counter].index;
3204: leafCounter++;
3205: }
3206: }
3207: if (isLeaf[i]) {
3208: counter++;
3209: }
3210: }
3212: PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
3213: PetscFree(newOwners);
3214: PetscFree(newNumbers);
3215: PetscFree(isLeaf);
3216: return(0);
3217: }
3219: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
3220: {
3221: PetscInt *distribution, min, max, sum, i, ierr;
3222: PetscMPIInt rank, size;
3224: MPI_Comm_size(comm, &size);
3225: MPI_Comm_rank(comm, &rank);
3226: PetscCalloc1(size, &distribution);
3227: for (i=0; i<n; i++) {
3228: if (part) distribution[part[i]] += vtxwgt[skip*i];
3229: else distribution[rank] += vtxwgt[skip*i];
3230: }
3231: MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
3232: min = distribution[0];
3233: max = distribution[0];
3234: sum = distribution[0];
3235: for (i=1; i<size; i++) {
3236: if (distribution[i]<min) min=distribution[i];
3237: if (distribution[i]>max) max=distribution[i];
3238: sum += distribution[i];
3239: }
3240: PetscViewerASCIIPrintf(viewer, "Min: %D, Avg: %D, Max: %D, Balance: %f\n", min, sum/size, max, (max*1.*size)/sum);
3241: PetscFree(distribution);
3242: return(0);
3243: }
3245: #endif
3247: /*@
3248: DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the PointSF of the DM inplace.
3250: Input parameters:
3251: + dm - The DMPlex object.
3252: . entityDepth - depth of the entity to balance (0 -> balance vertices).
3253: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
3254: - parallel - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.
3256: Output parameters:
3257: . success - whether the graph partitioning was successful or not. If not, try useInitialGuess=True and parallel=True.
3259: Level: intermediate
3261: @*/
3263: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
3264: {
3265: #if defined(PETSC_HAVE_PARMETIS)
3266: PetscSF sf;
3267: PetscInt ierr, i, j, idx, jdx;
3268: PetscInt eBegin, eEnd, nroots, nleafs, pStart, pEnd;
3269: const PetscInt *degrees, *ilocal;
3270: const PetscSFNode *iremote;
3271: PetscBool *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
3272: PetscInt numExclusivelyOwned, numNonExclusivelyOwned;
3273: PetscMPIInt rank, size;
3274: PetscInt *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
3275: const PetscInt *cumSumVertices;
3276: PetscInt offset, counter;
3277: PetscInt lenadjncy;
3278: PetscInt *xadj, *adjncy, *vtxwgt;
3279: PetscInt lenxadj;
3280: PetscInt *adjwgt = NULL;
3281: PetscInt *part, *options;
3282: PetscInt nparts, wgtflag, numflag, ncon, edgecut;
3283: real_t *ubvec;
3284: PetscInt *firstVertices, *renumbering;
3285: PetscInt failed, failedGlobal;
3286: MPI_Comm comm;
3287: Mat A, Apre;
3288: const char *prefix = NULL;
3289: PetscViewer viewer;
3290: PetscViewerFormat format;
3291: PetscLayout layout;
3294: if (success) *success = PETSC_FALSE;
3295: PetscObjectGetComm((PetscObject) dm, &comm);
3296: MPI_Comm_rank(comm, &rank);
3297: MPI_Comm_size(comm, &size);
3298: if (size==1) return(0);
3300: PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3302: PetscOptionsGetViewer(comm,((PetscObject)dm)->options, prefix,"-dm_rebalance_partition_view",&viewer,&format,NULL);
3303: if (viewer) {
3304: PetscViewerPushFormat(viewer,format);
3305: }
3307: /* Figure out all points in the plex that we are interested in balancing. */
3308: DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
3309: DMPlexGetChart(dm, &pStart, &pEnd);
3310: PetscMalloc1(pEnd-pStart, &toBalance);
3312: for (i=0; i<pEnd-pStart; i++) {
3313: toBalance[i] = (PetscBool)(i-pStart>=eBegin && i-pStart<eEnd);
3314: }
3316: /* There are three types of points:
3317: * exclusivelyOwned: points that are owned by this process and only seen by this process
3318: * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
3319: * leaf: a point that is seen by this process but owned by a different process
3320: */
3321: DMGetPointSF(dm, &sf);
3322: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
3323: PetscMalloc1(pEnd-pStart, &isLeaf);
3324: PetscMalloc1(pEnd-pStart, &isNonExclusivelyOwned);
3325: PetscMalloc1(pEnd-pStart, &isExclusivelyOwned);
3326: for (i=0; i<pEnd-pStart; i++) {
3327: isNonExclusivelyOwned[i] = PETSC_FALSE;
3328: isExclusivelyOwned[i] = PETSC_FALSE;
3329: isLeaf[i] = PETSC_FALSE;
3330: }
3332: /* start by marking all the leafs */
3333: for (i=0; i<nleafs; i++) {
3334: isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
3335: }
3337: /* for an owned point, we can figure out whether another processor sees it or
3338: * not by calculating its degree */
3339: PetscSFComputeDegreeBegin(sf, °rees);
3340: PetscSFComputeDegreeEnd(sf, °rees);
3342: numExclusivelyOwned = 0;
3343: numNonExclusivelyOwned = 0;
3344: for (i=0; i<pEnd-pStart; i++) {
3345: if (toBalance[i]) {
3346: if (degrees[i] > 0) {
3347: isNonExclusivelyOwned[i] = PETSC_TRUE;
3348: numNonExclusivelyOwned += 1;
3349: } else {
3350: if (!isLeaf[i]) {
3351: isExclusivelyOwned[i] = PETSC_TRUE;
3352: numExclusivelyOwned += 1;
3353: }
3354: }
3355: }
3356: }
3358: /* We are going to build a graph with one vertex per core representing the
3359: * exclusively owned points and then one vertex per nonExclusively owned
3360: * point. */
3362: PetscLayoutCreate(comm, &layout);
3363: PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
3364: PetscLayoutSetUp(layout);
3365: PetscLayoutGetRanges(layout, &cumSumVertices);
3367: PetscMalloc1(pEnd-pStart, &globalNumbersOfLocalOwnedVertices);
3368: for (i=0; i<pEnd-pStart; i++) {globalNumbersOfLocalOwnedVertices[i] = pStart - 1;}
3369: offset = cumSumVertices[rank];
3370: counter = 0;
3371: for (i=0; i<pEnd-pStart; i++) {
3372: if (toBalance[i]) {
3373: if (degrees[i] > 0) {
3374: globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
3375: counter++;
3376: }
3377: }
3378: }
3380: /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
3381: PetscMalloc1(pEnd-pStart, &leafGlobalNumbers);
3382: PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);
3383: PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);
3385: /* Now start building the data structure for ParMETIS */
3387: MatCreate(comm, &Apre);
3388: MatSetType(Apre, MATPREALLOCATOR);
3389: MatSetSizes(Apre, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
3390: MatSetUp(Apre);
3392: for (i=0; i<pEnd-pStart; i++) {
3393: if (toBalance[i]) {
3394: idx = cumSumVertices[rank];
3395: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
3396: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
3397: else continue;
3398: MatSetValue(Apre, idx, jdx, 1, INSERT_VALUES);
3399: MatSetValue(Apre, jdx, idx, 1, INSERT_VALUES);
3400: }
3401: }
3403: MatAssemblyBegin(Apre, MAT_FINAL_ASSEMBLY);
3404: MatAssemblyEnd(Apre, MAT_FINAL_ASSEMBLY);
3406: MatCreate(comm, &A);
3407: MatSetType(A, MATMPIAIJ);
3408: MatSetSizes(A, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
3409: MatPreallocatorPreallocate(Apre, PETSC_FALSE, A);
3410: MatDestroy(&Apre);
3412: for (i=0; i<pEnd-pStart; i++) {
3413: if (toBalance[i]) {
3414: idx = cumSumVertices[rank];
3415: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
3416: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
3417: else continue;
3418: MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
3419: MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
3420: }
3421: }
3423: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
3424: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
3425: PetscFree(leafGlobalNumbers);
3426: PetscFree(globalNumbersOfLocalOwnedVertices);
3428: nparts = size;
3429: wgtflag = 2;
3430: numflag = 0;
3431: ncon = 2;
3432: real_t *tpwgts;
3433: PetscMalloc1(ncon * nparts, &tpwgts);
3434: for (i=0; i<ncon*nparts; i++) {
3435: tpwgts[i] = 1./(nparts);
3436: }
3438: PetscMalloc1(ncon, &ubvec);
3439: ubvec[0] = 1.01;
3440: ubvec[1] = 1.01;
3441: lenadjncy = 0;
3442: for (i=0; i<1+numNonExclusivelyOwned; i++) {
3443: PetscInt temp=0;
3444: MatGetRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
3445: lenadjncy += temp;
3446: MatRestoreRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
3447: }
3448: PetscMalloc1(lenadjncy, &adjncy);
3449: lenxadj = 2 + numNonExclusivelyOwned;
3450: PetscMalloc1(lenxadj, &xadj);
3451: xadj[0] = 0;
3452: counter = 0;
3453: for (i=0; i<1+numNonExclusivelyOwned; i++) {
3454: PetscInt temp=0;
3455: const PetscInt *cols;
3456: MatGetRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
3457: PetscArraycpy(&adjncy[counter], cols, temp);
3458: counter += temp;
3459: xadj[i+1] = counter;
3460: MatRestoreRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
3461: }
3463: PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);
3464: PetscMalloc1(ncon*(1 + numNonExclusivelyOwned), &vtxwgt);
3465: vtxwgt[0] = numExclusivelyOwned;
3466: if (ncon>1) vtxwgt[1] = 1;
3467: for (i=0; i<numNonExclusivelyOwned; i++) {
3468: vtxwgt[ncon*(i+1)] = 1;
3469: if (ncon>1) vtxwgt[ncon*(i+1)+1] = 0;
3470: }
3472: if (viewer) {
3473: PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %D on interface of mesh distribution.\n", entityDepth);
3474: PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %D\n", cumSumVertices[size]);
3475: }
3476: if (parallel) {
3477: PetscMalloc1(4, &options);
3478: options[0] = 1;
3479: options[1] = 0; /* Verbosity */
3480: options[2] = 0; /* Seed */
3481: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
3482: if (viewer) { PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"); }
3483: if (useInitialGuess) {
3484: if (viewer) { PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"); }
3485: PetscStackPush("ParMETIS_V3_RefineKway");
3486: ParMETIS_V3_RefineKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
3487: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_RefineKway()");
3488: PetscStackPop;
3489: } else {
3490: PetscStackPush("ParMETIS_V3_PartKway");
3491: ParMETIS_V3_PartKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
3492: PetscStackPop;
3493: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
3494: }
3495: PetscFree(options);
3496: } else {
3497: if (viewer) { PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"); }
3498: Mat As;
3499: PetscInt numRows;
3500: PetscInt *partGlobal;
3501: MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &As);
3503: PetscInt *numExclusivelyOwnedAll;
3504: PetscMalloc1(size, &numExclusivelyOwnedAll);
3505: numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
3506: MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,numExclusivelyOwnedAll,1,MPIU_INT,comm);
3508: MatGetSize(As, &numRows, NULL);
3509: PetscMalloc1(numRows, &partGlobal);
3510: if (!rank) {
3511: PetscInt *adjncy_g, *xadj_g, *vtxwgt_g;
3512: lenadjncy = 0;
3514: for (i=0; i<numRows; i++) {
3515: PetscInt temp=0;
3516: MatGetRow(As, i, &temp, NULL, NULL);
3517: lenadjncy += temp;
3518: MatRestoreRow(As, i, &temp, NULL, NULL);
3519: }
3520: PetscMalloc1(lenadjncy, &adjncy_g);
3521: lenxadj = 1 + numRows;
3522: PetscMalloc1(lenxadj, &xadj_g);
3523: xadj_g[0] = 0;
3524: counter = 0;
3525: for (i=0; i<numRows; i++) {
3526: PetscInt temp=0;
3527: const PetscInt *cols;
3528: MatGetRow(As, i, &temp, &cols, NULL);
3529: PetscArraycpy(&adjncy_g[counter], cols, temp);
3530: counter += temp;
3531: xadj_g[i+1] = counter;
3532: MatRestoreRow(As, i, &temp, &cols, NULL);
3533: }
3534: PetscMalloc1(2*numRows, &vtxwgt_g);
3535: for (i=0; i<size; i++){
3536: vtxwgt_g[ncon*cumSumVertices[i]] = numExclusivelyOwnedAll[i];
3537: if (ncon>1) vtxwgt_g[ncon*cumSumVertices[i]+1] = 1;
3538: for (j=cumSumVertices[i]+1; j<cumSumVertices[i+1]; j++) {
3539: vtxwgt_g[ncon*j] = 1;
3540: if (ncon>1) vtxwgt_g[2*j+1] = 0;
3541: }
3542: }
3543: PetscMalloc1(64, &options);
3544: METIS_SetDefaultOptions(options); /* initialize all defaults */
3545: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
3546: options[METIS_OPTION_CONTIG] = 1;
3547: PetscStackPush("METIS_PartGraphKway");
3548: METIS_PartGraphKway(&numRows, &ncon, xadj_g, adjncy_g, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
3549: PetscStackPop;
3550: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
3551: PetscFree(options);
3552: PetscFree(xadj_g);
3553: PetscFree(adjncy_g);
3554: PetscFree(vtxwgt_g);
3555: }
3556: PetscFree(numExclusivelyOwnedAll);
3558: /* Now scatter the parts array. */
3559: {
3560: PetscMPIInt *counts, *mpiCumSumVertices;
3561: PetscMalloc1(size, &counts);
3562: PetscMalloc1(size+1, &mpiCumSumVertices);
3563: for(i=0; i<size; i++) {
3564: PetscMPIIntCast(cumSumVertices[i+1] - cumSumVertices[i], &(counts[i]));
3565: }
3566: for(i=0; i<=size; i++) {
3567: PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
3568: }
3569: MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
3570: PetscFree(counts);
3571: PetscFree(mpiCumSumVertices);
3572: }
3574: PetscFree(partGlobal);
3575: MatDestroy(&As);
3576: }
3578: MatDestroy(&A);
3579: PetscFree(ubvec);
3580: PetscFree(tpwgts);
3582: /* Now rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
3584: PetscMalloc1(size, &firstVertices);
3585: PetscMalloc1(size, &renumbering);
3586: firstVertices[rank] = part[0];
3587: MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);
3588: for (i=0; i<size; i++) {
3589: renumbering[firstVertices[i]] = i;
3590: }
3591: for (i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) {
3592: part[i] = renumbering[part[i]];
3593: }
3594: /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
3595: failed = (PetscInt)(part[0] != rank);
3596: MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
3598: PetscFree(firstVertices);
3599: PetscFree(renumbering);
3601: if (failedGlobal > 0) {
3602: PetscLayoutDestroy(&layout);
3603: PetscFree(xadj);
3604: PetscFree(adjncy);
3605: PetscFree(vtxwgt);
3606: PetscFree(toBalance);
3607: PetscFree(isLeaf);
3608: PetscFree(isNonExclusivelyOwned);
3609: PetscFree(isExclusivelyOwned);
3610: PetscFree(part);
3611: if (viewer) {
3612: PetscViewerPopFormat(viewer);
3613: PetscViewerDestroy(&viewer);
3614: }
3615: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3616: return(0);
3617: }
3619: /*Let's check how well we did distributing points*/
3620: if (viewer) {
3621: PetscViewerASCIIPrintf(viewer, "Comparing number of owned entities of depth %D on each process before rebalancing, after rebalancing, and after consistency checks.\n", entityDepth);
3622: PetscViewerASCIIPrintf(viewer, "Initial. ");
3623: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
3624: PetscViewerASCIIPrintf(viewer, "Rebalanced. ");
3625: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
3626: }
3628: /* Now check that every vertex is owned by a process that it is actually connected to. */
3629: for (i=1; i<=numNonExclusivelyOwned; i++) {
3630: PetscInt loc = 0;
3631: PetscFindInt(cumSumVertices[part[i]], xadj[i+1]-xadj[i], &adjncy[xadj[i]], &loc);
3632: /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
3633: if (loc<0) {
3634: part[i] = rank;
3635: }
3636: }
3638: /* Let's see how significant the influences of the previous fixing up step was.*/
3639: if (viewer) {
3640: PetscViewerASCIIPrintf(viewer, "After. ");
3641: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
3642: }
3644: PetscLayoutDestroy(&layout);
3645: PetscFree(xadj);
3646: PetscFree(adjncy);
3647: PetscFree(vtxwgt);
3649: /* Almost done, now rewrite the SF to reflect the new ownership. */
3650: {
3651: PetscInt *pointsToRewrite;
3652: PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
3653: counter = 0;
3654: for(i=0; i<pEnd-pStart; i++) {
3655: if (toBalance[i]) {
3656: if (isNonExclusivelyOwned[i]) {
3657: pointsToRewrite[counter] = i + pStart;
3658: counter++;
3659: }
3660: }
3661: }
3662: DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part+1, degrees);
3663: PetscFree(pointsToRewrite);
3664: }
3666: PetscFree(toBalance);
3667: PetscFree(isLeaf);
3668: PetscFree(isNonExclusivelyOwned);
3669: PetscFree(isExclusivelyOwned);
3670: PetscFree(part);
3671: if (viewer) {
3672: PetscViewerPopFormat(viewer);
3673: PetscViewerDestroy(&viewer);
3674: }
3675: if (success) *success = PETSC_TRUE;
3676: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
3677: #else
3678: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
3679: #endif
3680: return(0);
3681: }