Actual source code: plexpartition.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/partitionerimpl.h>
3: #include <petsc/private/hashseti.h>
5: const char *const DMPlexCSRAlgorithms[] = {"mat", "graph", "overlap", "DMPlexCSRAlgorithm", "DM_PLEX_CSR_", NULL};
7: static inline PetscInt DMPlex_GlobalID(PetscInt point)
8: {
9: return point >= 0 ? point : -(point + 1);
10: }
12: static PetscErrorCode DMPlexCreatePartitionerGraph_Overlap(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
13: {
14: DM ovdm;
15: PetscSF sfPoint;
16: IS cellNumbering;
17: const PetscInt *cellNum;
18: PetscInt *adj = NULL, *vOffsets = NULL, *vAdj = NULL;
19: PetscBool useCone, useClosure;
20: PetscInt dim, depth, overlap, cStart, cEnd, c, v;
21: PetscMPIInt rank, size;
23: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
24: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
25: DMGetDimension(dm, &dim);
26: DMPlexGetDepth(dm, &depth);
27: if (dim != depth) {
28: /* We do not handle the uninterpolated case here */
29: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
30: /* DMPlexCreateNeighborCSR does not make a numbering */
31: if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
32: /* Different behavior for empty graphs */
33: if (!*numVertices) {
34: PetscMalloc1(1, offsets);
35: (*offsets)[0] = 0;
36: }
37: /* Broken in parallel */
39: return 0;
40: }
41: /* Always use FVM adjacency to create partitioner graph */
42: DMGetBasicAdjacency(dm, &useCone, &useClosure);
43: DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
44: /* Need overlap >= 1 */
45: DMPlexGetOverlap(dm, &overlap);
46: if (size && overlap < 1) {
47: DMPlexDistributeOverlap(dm, 1, NULL, &ovdm);
48: } else {
49: PetscObjectReference((PetscObject)dm);
50: ovdm = dm;
51: }
52: DMGetPointSF(ovdm, &sfPoint);
53: DMPlexGetHeightStratum(ovdm, height, &cStart, &cEnd);
54: DMPlexCreateNumbering_Plex(ovdm, cStart, cEnd, 0, NULL, sfPoint, &cellNumbering);
55: if (globalNumbering) {
56: PetscObjectReference((PetscObject)cellNumbering);
57: *globalNumbering = cellNumbering;
58: }
59: ISGetIndices(cellNumbering, &cellNum);
60: /* Determine sizes */
61: for (*numVertices = 0, c = cStart; c < cEnd; ++c) {
62: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
63: if (cellNum[c - cStart] < 0) continue;
64: (*numVertices)++;
65: }
66: PetscCalloc1(*numVertices + 1, &vOffsets);
67: for (c = cStart, v = 0; c < cEnd; ++c) {
68: PetscInt adjSize = PETSC_DETERMINE, a, vsize = 0;
70: if (cellNum[c - cStart] < 0) continue;
71: DMPlexGetAdjacency(ovdm, c, &adjSize, &adj);
72: for (a = 0; a < adjSize; ++a) {
73: const PetscInt point = adj[a];
74: if (point != c && cStart <= point && point < cEnd) ++vsize;
75: }
76: vOffsets[v + 1] = vOffsets[v] + vsize;
77: ++v;
78: }
79: /* Determine adjacency */
80: PetscMalloc1(vOffsets[*numVertices], &vAdj);
81: for (c = cStart, v = 0; c < cEnd; ++c) {
82: PetscInt adjSize = PETSC_DETERMINE, a, off = vOffsets[v];
84: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
85: if (cellNum[c - cStart] < 0) continue;
86: DMPlexGetAdjacency(ovdm, c, &adjSize, &adj);
87: for (a = 0; a < adjSize; ++a) {
88: const PetscInt point = adj[a];
89: if (point != c && cStart <= point && point < cEnd) vAdj[off++] = DMPlex_GlobalID(cellNum[point - cStart]);
90: }
92: /* Sort adjacencies (not strictly necessary) */
93: PetscSortInt(off - vOffsets[v], &vAdj[vOffsets[v]]);
94: ++v;
95: }
96: PetscFree(adj);
97: ISRestoreIndices(cellNumbering, &cellNum);
98: ISDestroy(&cellNumbering);
99: DMSetBasicAdjacency(dm, useCone, useClosure);
100: DMDestroy(&ovdm);
101: if (offsets) {
102: *offsets = vOffsets;
103: } else PetscFree(vOffsets);
104: if (adjacency) {
105: *adjacency = vAdj;
106: } else PetscFree(vAdj);
107: return 0;
108: }
110: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
111: {
112: PetscInt dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
113: PetscInt *adj = NULL, *vOffsets = NULL, *graph = NULL;
114: IS cellNumbering;
115: const PetscInt *cellNum;
116: PetscBool useCone, useClosure;
117: PetscSection section;
118: PetscSegBuffer adjBuffer;
119: PetscSF sfPoint;
120: PetscInt *adjCells = NULL, *remoteCells = NULL;
121: const PetscInt *local;
122: PetscInt nroots, nleaves, l;
123: PetscMPIInt rank;
125: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
126: DMGetDimension(dm, &dim);
127: DMPlexGetDepth(dm, &depth);
128: if (dim != depth) {
129: /* We do not handle the uninterpolated case here */
130: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
131: /* DMPlexCreateNeighborCSR does not make a numbering */
132: if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
133: /* Different behavior for empty graphs */
134: if (!*numVertices) {
135: PetscMalloc1(1, offsets);
136: (*offsets)[0] = 0;
137: }
138: /* Broken in parallel */
140: return 0;
141: }
142: DMGetPointSF(dm, &sfPoint);
143: DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
144: /* Build adjacency graph via a section/segbuffer */
145: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
146: PetscSectionSetChart(section, pStart, pEnd);
147: PetscSegBufferCreate(sizeof(PetscInt), 1000, &adjBuffer);
148: /* Always use FVM adjacency to create partitioner graph */
149: DMGetBasicAdjacency(dm, &useCone, &useClosure);
150: DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
151: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
152: if (globalNumbering) {
153: PetscObjectReference((PetscObject)cellNumbering);
154: *globalNumbering = cellNumbering;
155: }
156: ISGetIndices(cellNumbering, &cellNum);
157: /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
158: Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
159: */
160: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
161: if (nroots >= 0) {
162: PetscInt fStart, fEnd, f;
164: PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
165: DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd);
166: for (l = 0; l < nroots; ++l) adjCells[l] = -3;
167: for (f = fStart; f < fEnd; ++f) {
168: const PetscInt *support;
169: PetscInt supportSize;
171: DMPlexGetSupport(dm, f, &support);
172: DMPlexGetSupportSize(dm, f, &supportSize);
173: if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
174: else if (supportSize == 2) {
175: PetscFindInt(support[0], nleaves, local, &p);
176: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
177: PetscFindInt(support[1], nleaves, local, &p);
178: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
179: }
180: /* Handle non-conforming meshes */
181: if (supportSize > 2) {
182: PetscInt numChildren, i;
183: const PetscInt *children;
185: DMPlexGetTreeChildren(dm, f, &numChildren, &children);
186: for (i = 0; i < numChildren; ++i) {
187: const PetscInt child = children[i];
188: if (fStart <= child && child < fEnd) {
189: DMPlexGetSupport(dm, child, &support);
190: DMPlexGetSupportSize(dm, child, &supportSize);
191: if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
192: else if (supportSize == 2) {
193: PetscFindInt(support[0], nleaves, local, &p);
194: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1] - pStart]);
195: PetscFindInt(support[1], nleaves, local, &p);
196: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0] - pStart]);
197: }
198: }
199: }
200: }
201: }
202: for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
203: PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE);
204: PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_REPLACE);
205: PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
206: PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
207: }
208: /* Combine local and global adjacencies */
209: for (*numVertices = 0, p = pStart; p < pEnd; p++) {
210: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
211: if (nroots > 0) {
212: if (cellNum[p - pStart] < 0) continue;
213: }
214: /* Add remote cells */
215: if (remoteCells) {
216: const PetscInt gp = DMPlex_GlobalID(cellNum[p - pStart]);
217: PetscInt coneSize, numChildren, c, i;
218: const PetscInt *cone, *children;
220: DMPlexGetCone(dm, p, &cone);
221: DMPlexGetConeSize(dm, p, &coneSize);
222: for (c = 0; c < coneSize; ++c) {
223: const PetscInt point = cone[c];
224: if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
225: PetscInt *PETSC_RESTRICT pBuf;
226: PetscSectionAddDof(section, p, 1);
227: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
228: *pBuf = remoteCells[point];
229: }
230: /* Handle non-conforming meshes */
231: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
232: for (i = 0; i < numChildren; ++i) {
233: const PetscInt child = children[i];
234: if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
235: PetscInt *PETSC_RESTRICT pBuf;
236: PetscSectionAddDof(section, p, 1);
237: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
238: *pBuf = remoteCells[child];
239: }
240: }
241: }
242: }
243: /* Add local cells */
244: adjSize = PETSC_DETERMINE;
245: DMPlexGetAdjacency(dm, p, &adjSize, &adj);
246: for (a = 0; a < adjSize; ++a) {
247: const PetscInt point = adj[a];
248: if (point != p && pStart <= point && point < pEnd) {
249: PetscInt *PETSC_RESTRICT pBuf;
250: PetscSectionAddDof(section, p, 1);
251: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
252: *pBuf = DMPlex_GlobalID(cellNum[point - pStart]);
253: }
254: }
255: (*numVertices)++;
256: }
257: PetscFree(adj);
258: PetscFree2(adjCells, remoteCells);
259: DMSetBasicAdjacency(dm, useCone, useClosure);
261: /* Derive CSR graph from section/segbuffer */
262: PetscSectionSetUp(section);
263: PetscSectionGetStorageSize(section, &size);
264: PetscMalloc1(*numVertices + 1, &vOffsets);
265: for (idx = 0, p = pStart; p < pEnd; p++) {
266: if (nroots > 0) {
267: if (cellNum[p - pStart] < 0) continue;
268: }
269: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
270: }
271: vOffsets[*numVertices] = size;
272: PetscSegBufferExtractAlloc(adjBuffer, &graph);
274: if (nroots >= 0) {
275: /* Filter out duplicate edges using section/segbuffer */
276: PetscSectionReset(section);
277: PetscSectionSetChart(section, 0, *numVertices);
278: for (p = 0; p < *numVertices; p++) {
279: PetscInt start = vOffsets[p], end = vOffsets[p + 1];
280: PetscInt numEdges = end - start, *PETSC_RESTRICT edges;
281: PetscSortRemoveDupsInt(&numEdges, &graph[start]);
282: PetscSectionSetDof(section, p, numEdges);
283: PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
284: PetscArraycpy(edges, &graph[start], numEdges);
285: }
286: PetscFree(vOffsets);
287: PetscFree(graph);
288: /* Derive CSR graph from section/segbuffer */
289: PetscSectionSetUp(section);
290: PetscSectionGetStorageSize(section, &size);
291: PetscMalloc1(*numVertices + 1, &vOffsets);
292: for (idx = 0, p = 0; p < *numVertices; p++) PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
293: vOffsets[*numVertices] = size;
294: PetscSegBufferExtractAlloc(adjBuffer, &graph);
295: } else {
296: /* Sort adjacencies (not strictly necessary) */
297: for (p = 0; p < *numVertices; p++) {
298: PetscInt start = vOffsets[p], end = vOffsets[p + 1];
299: PetscSortInt(end - start, &graph[start]);
300: }
301: }
303: if (offsets) *offsets = vOffsets;
304: if (adjacency) *adjacency = graph;
306: /* Cleanup */
307: PetscSegBufferDestroy(&adjBuffer);
308: PetscSectionDestroy(§ion);
309: ISRestoreIndices(cellNumbering, &cellNum);
310: ISDestroy(&cellNumbering);
311: return 0;
312: }
314: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
315: {
316: Mat conn, CSR;
317: IS fis, cis, cis_own;
318: PetscSF sfPoint;
319: const PetscInt *rows, *cols, *ii, *jj;
320: PetscInt *idxs, *idxs2;
321: PetscInt dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
322: PetscMPIInt rank;
323: PetscBool flg;
325: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
326: DMGetDimension(dm, &dim);
327: DMPlexGetDepth(dm, &depth);
328: if (dim != depth) {
329: /* We do not handle the uninterpolated case here */
330: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
331: /* DMPlexCreateNeighborCSR does not make a numbering */
332: if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
333: /* Different behavior for empty graphs */
334: if (!*numVertices) {
335: PetscMalloc1(1, offsets);
336: (*offsets)[0] = 0;
337: }
338: /* Broken in parallel */
340: return 0;
341: }
342: /* Interpolated and parallel case */
344: /* numbering */
345: DMGetPointSF(dm, &sfPoint);
346: DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
347: DMPlexGetHeightStratum(dm, height + 1, &fStart, &fEnd);
348: DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
349: DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
350: if (globalNumbering) ISDuplicate(cis, globalNumbering);
352: /* get positive global ids and local sizes for facets and cells */
353: ISGetLocalSize(fis, &m);
354: ISGetIndices(fis, &rows);
355: PetscMalloc1(m, &idxs);
356: for (i = 0, floc = 0; i < m; i++) {
357: const PetscInt p = rows[i];
359: if (p < 0) {
360: idxs[i] = -(p + 1);
361: } else {
362: idxs[i] = p;
363: floc += 1;
364: }
365: }
366: ISRestoreIndices(fis, &rows);
367: ISDestroy(&fis);
368: ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);
370: ISGetLocalSize(cis, &m);
371: ISGetIndices(cis, &cols);
372: PetscMalloc1(m, &idxs);
373: PetscMalloc1(m, &idxs2);
374: for (i = 0, cloc = 0; i < m; i++) {
375: const PetscInt p = cols[i];
377: if (p < 0) {
378: idxs[i] = -(p + 1);
379: } else {
380: idxs[i] = p;
381: idxs2[cloc++] = p;
382: }
383: }
384: ISRestoreIndices(cis, &cols);
385: ISDestroy(&cis);
386: ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
387: ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);
389: /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
390: MatCreate(PetscObjectComm((PetscObject)dm), &conn);
391: MatSetSizes(conn, floc, cloc, M, N);
392: MatSetType(conn, MATMPIAIJ);
393: DMPlexGetMaxSizes(dm, NULL, &lm);
394: MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
395: MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);
397: /* Assemble matrix */
398: ISGetIndices(fis, &rows);
399: ISGetIndices(cis, &cols);
400: for (c = cStart; c < cEnd; c++) {
401: const PetscInt *cone;
402: PetscInt coneSize, row, col, f;
404: col = cols[c - cStart];
405: DMPlexGetCone(dm, c, &cone);
406: DMPlexGetConeSize(dm, c, &coneSize);
407: for (f = 0; f < coneSize; f++) {
408: const PetscScalar v = 1.0;
409: const PetscInt *children;
410: PetscInt numChildren, ch;
412: row = rows[cone[f] - fStart];
413: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
415: /* non-conforming meshes */
416: DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
417: for (ch = 0; ch < numChildren; ch++) {
418: const PetscInt child = children[ch];
420: if (child < fStart || child >= fEnd) continue;
421: row = rows[child - fStart];
422: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
423: }
424: }
425: }
426: ISRestoreIndices(fis, &rows);
427: ISRestoreIndices(cis, &cols);
428: ISDestroy(&fis);
429: ISDestroy(&cis);
430: MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
431: MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);
433: /* Get parallel CSR by doing conn^T * conn */
434: MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
435: MatDestroy(&conn);
437: /* extract local part of the CSR */
438: MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
439: MatDestroy(&CSR);
440: MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
443: /* get back requested output */
444: if (numVertices) *numVertices = m;
445: if (offsets) {
446: PetscCalloc1(m + 1, &idxs);
447: for (i = 1; i < m + 1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
448: *offsets = idxs;
449: }
450: if (adjacency) {
451: PetscMalloc1(ii[m] - m, &idxs);
452: ISGetIndices(cis_own, &rows);
453: for (i = 0, c = 0; i < m; i++) {
454: PetscInt j, g = rows[i];
456: for (j = ii[i]; j < ii[i + 1]; j++) {
457: if (jj[j] == g) continue; /* again, self-connectivity */
458: idxs[c++] = jj[j];
459: }
460: }
462: ISRestoreIndices(cis_own, &rows);
463: *adjacency = idxs;
464: }
466: /* cleanup */
467: ISDestroy(&cis_own);
468: MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
470: MatDestroy(&conn);
471: return 0;
472: }
474: /*@C
475: DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner
477: Collective on dm
479: Input Parameters:
480: + dm - The mesh `DM`
481: - height - Height of the strata from which to construct the graph
483: Output Parameters:
484: + numVertices - Number of vertices in the graph
485: . offsets - Point offsets in the graph
486: . adjacency - Point connectivity in the graph
487: - 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.
489: Options Database Keys:
490: . -dm_plex_csr_alg <mat,graph,overlap> - Choose the algorithm for computing the CSR graph
492: Level: developer
494: Note:
495: The user can control the definition of adjacency for the mesh using `DMSetAdjacency()`. They should choose the combination appropriate for the function
496: representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().
498: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscPartitionerGetType()`, `PetscPartitionerCreate()`, `DMSetAdjacency()`
499: @*/
500: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
501: {
502: DMPlexCSRAlgorithm alg = DM_PLEX_CSR_GRAPH;
504: PetscOptionsGetEnum(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *)&alg, NULL);
505: switch (alg) {
506: case DM_PLEX_CSR_MAT:
507: DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
508: break;
509: case DM_PLEX_CSR_GRAPH:
510: DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
511: break;
512: case DM_PLEX_CSR_OVERLAP:
513: DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering);
514: break;
515: }
516: return 0;
517: }
519: /*@C
520: DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.
522: Collective on dm
524: Input Parameters:
525: + dm - The `DMPLEX`
526: - cellHeight - The height of mesh points to treat as cells (default should be 0)
528: Output Parameters:
529: + numVertices - The number of local vertices in the graph, or cells in the mesh.
530: . offsets - The offset to the adjacency list for each cell
531: - adjacency - The adjacency list for all cells
533: Level: advanced
535: Note:
536: This is suitable for input to a mesh partitioner like ParMetis.
538: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexCreate()`
539: @*/
540: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
541: {
542: const PetscInt maxFaceCases = 30;
543: PetscInt numFaceCases = 0;
544: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
545: PetscInt *off, *adj;
546: PetscInt *neighborCells = NULL;
547: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
549: /* For parallel partitioning, I think you have to communicate supports */
550: DMGetDimension(dm, &dim);
551: cellDim = dim - cellHeight;
552: DMPlexGetDepth(dm, &depth);
553: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
554: if (cEnd - cStart == 0) {
555: if (numVertices) *numVertices = 0;
556: if (offsets) *offsets = NULL;
557: if (adjacency) *adjacency = NULL;
558: return 0;
559: }
560: numCells = cEnd - cStart;
561: faceDepth = depth - cellHeight;
562: if (dim == depth) {
563: PetscInt f, fStart, fEnd;
565: PetscCalloc1(numCells + 1, &off);
566: /* Count neighboring cells */
567: DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd);
568: for (f = fStart; f < fEnd; ++f) {
569: const PetscInt *support;
570: PetscInt supportSize;
571: DMPlexGetSupportSize(dm, f, &supportSize);
572: DMPlexGetSupport(dm, f, &support);
573: if (supportSize == 2) {
574: PetscInt numChildren;
576: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
577: if (!numChildren) {
578: ++off[support[0] - cStart + 1];
579: ++off[support[1] - cStart + 1];
580: }
581: }
582: }
583: /* Prefix sum */
584: for (c = 1; c <= numCells; ++c) off[c] += off[c - 1];
585: if (adjacency) {
586: PetscInt *tmp;
588: PetscMalloc1(off[numCells], &adj);
589: PetscMalloc1(numCells + 1, &tmp);
590: PetscArraycpy(tmp, off, numCells + 1);
591: /* Get neighboring cells */
592: for (f = fStart; f < fEnd; ++f) {
593: const PetscInt *support;
594: PetscInt supportSize;
595: DMPlexGetSupportSize(dm, f, &supportSize);
596: DMPlexGetSupport(dm, f, &support);
597: if (supportSize == 2) {
598: PetscInt numChildren;
600: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
601: if (!numChildren) {
602: adj[tmp[support[0] - cStart]++] = support[1];
603: adj[tmp[support[1] - cStart]++] = support[0];
604: }
605: }
606: }
607: for (c = 0; c < cEnd - cStart; ++c) PetscAssert(tmp[c] == off[c + 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %" PetscInt_FMT " != %" PetscInt_FMT " for cell %" PetscInt_FMT, tmp[c], off[c], c + cStart);
608: PetscFree(tmp);
609: }
610: if (numVertices) *numVertices = numCells;
611: if (offsets) *offsets = off;
612: if (adjacency) *adjacency = adj;
613: return 0;
614: }
615: /* Setup face recognition */
616: if (faceDepth == 1) {
617: 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 */
619: for (c = cStart; c < cEnd; ++c) {
620: PetscInt corners;
622: DMPlexGetConeSize(dm, c, &corners);
623: if (!cornersSeen[corners]) {
624: PetscInt nFV;
627: cornersSeen[corners] = 1;
629: DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);
631: numFaceVertices[numFaceCases++] = nFV;
632: }
633: }
634: }
635: PetscCalloc1(numCells + 1, &off);
636: /* Count neighboring cells */
637: for (cell = cStart; cell < cEnd; ++cell) {
638: PetscInt numNeighbors = PETSC_DETERMINE, n;
640: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
641: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
642: for (n = 0; n < numNeighbors; ++n) {
643: PetscInt cellPair[2];
644: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
645: PetscInt meetSize = 0;
646: const PetscInt *meet = NULL;
648: cellPair[0] = cell;
649: cellPair[1] = neighborCells[n];
650: if (cellPair[0] == cellPair[1]) continue;
651: if (!found) {
652: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
653: if (meetSize) {
654: PetscInt f;
656: for (f = 0; f < numFaceCases; ++f) {
657: if (numFaceVertices[f] == meetSize) {
658: found = PETSC_TRUE;
659: break;
660: }
661: }
662: }
663: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
664: }
665: if (found) ++off[cell - cStart + 1];
666: }
667: }
668: /* Prefix sum */
669: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell - 1];
671: if (adjacency) {
672: PetscMalloc1(off[numCells], &adj);
673: /* Get neighboring cells */
674: for (cell = cStart; cell < cEnd; ++cell) {
675: PetscInt numNeighbors = PETSC_DETERMINE, n;
676: PetscInt cellOffset = 0;
678: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
679: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
680: for (n = 0; n < numNeighbors; ++n) {
681: PetscInt cellPair[2];
682: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
683: PetscInt meetSize = 0;
684: const PetscInt *meet = NULL;
686: cellPair[0] = cell;
687: cellPair[1] = neighborCells[n];
688: if (cellPair[0] == cellPair[1]) continue;
689: if (!found) {
690: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
691: if (meetSize) {
692: PetscInt f;
694: for (f = 0; f < numFaceCases; ++f) {
695: if (numFaceVertices[f] == meetSize) {
696: found = PETSC_TRUE;
697: break;
698: }
699: }
700: }
701: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
702: }
703: if (found) {
704: adj[off[cell - cStart] + cellOffset] = neighborCells[n];
705: ++cellOffset;
706: }
707: }
708: }
709: }
710: PetscFree(neighborCells);
711: if (numVertices) *numVertices = numCells;
712: if (offsets) *offsets = off;
713: if (adjacency) *adjacency = adj;
714: return 0;
715: }
717: /*@
718: PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh
720: Collective on part
722: Input Parameters:
723: + part - The `PetscPartitioner`
724: . targetSection - The `PetscSection` describing the absolute weight of each partition (can be NULL)
725: - dm - The mesh `DM`
727: Output Parameters:
728: + partSection - The `PetscSection` giving the division of points by partition
729: - partition - The list of points by partition
731: Level: developer
733: Note:
734: If the `DM` has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
735: by the section in the transitive closure of the point.
737: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `PetscPartitionerCreate()`, `PetscSectionCreate()`,
738: `PetscSectionSetChart()`, `PetscPartitionerPartition()`
739: @*/
740: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
741: {
742: PetscMPIInt size;
743: PetscBool isplex;
744: PetscSection vertSection = NULL;
751: PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex);
753: MPI_Comm_size(PetscObjectComm((PetscObject)part), &size);
754: if (size == 1) {
755: PetscInt *points;
756: PetscInt cStart, cEnd, c;
758: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
759: PetscSectionReset(partSection);
760: PetscSectionSetChart(partSection, 0, size);
761: PetscSectionSetDof(partSection, 0, cEnd - cStart);
762: PetscSectionSetUp(partSection);
763: PetscMalloc1(cEnd - cStart, &points);
764: for (c = cStart; c < cEnd; ++c) points[c] = c;
765: ISCreateGeneral(PetscObjectComm((PetscObject)part), cEnd - cStart, points, PETSC_OWN_POINTER, partition);
766: return 0;
767: }
768: if (part->height == 0) {
769: PetscInt numVertices = 0;
770: PetscInt *start = NULL;
771: PetscInt *adjacency = NULL;
772: IS globalNumbering;
774: if (!part->noGraph || part->viewGraph) {
775: DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
776: } else { /* only compute the number of owned local vertices */
777: const PetscInt *idxs;
778: PetscInt p, pStart, pEnd;
780: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
781: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
782: ISGetIndices(globalNumbering, &idxs);
783: for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
784: ISRestoreIndices(globalNumbering, &idxs);
785: }
786: if (part->usevwgt) {
787: PetscSection section = dm->localSection, clSection = NULL;
788: IS clPoints = NULL;
789: const PetscInt *gid, *clIdx;
790: PetscInt v, p, pStart, pEnd;
792: /* 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) */
793: /* We do this only if the local section has been set */
794: if (section) {
795: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
796: if (!clSection) DMPlexCreateClosureIndex(dm, NULL);
797: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
798: ISGetIndices(clPoints, &clIdx);
799: }
800: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
801: PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
802: PetscSectionSetChart(vertSection, 0, numVertices);
803: if (globalNumbering) {
804: ISGetIndices(globalNumbering, &gid);
805: } else gid = NULL;
806: for (p = pStart, v = 0; p < pEnd; ++p) {
807: PetscInt dof = 1;
809: /* skip cells in the overlap */
810: if (gid && gid[p - pStart] < 0) continue;
812: if (section) {
813: PetscInt cl, clSize, clOff;
815: dof = 0;
816: PetscSectionGetDof(clSection, p, &clSize);
817: PetscSectionGetOffset(clSection, p, &clOff);
818: for (cl = 0; cl < clSize; cl += 2) {
819: PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */
821: PetscSectionGetDof(section, clPoint, &clDof);
822: dof += clDof;
823: }
824: }
826: PetscSectionSetDof(vertSection, v, dof);
827: v++;
828: }
829: if (globalNumbering) ISRestoreIndices(globalNumbering, &gid);
830: if (clPoints) ISRestoreIndices(clPoints, &clIdx);
831: PetscSectionSetUp(vertSection);
832: }
833: PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
834: PetscFree(start);
835: PetscFree(adjacency);
836: if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
837: const PetscInt *globalNum;
838: const PetscInt *partIdx;
839: PetscInt *map, cStart, cEnd;
840: PetscInt *adjusted, i, localSize, offset;
841: IS newPartition;
843: ISGetLocalSize(*partition, &localSize);
844: PetscMalloc1(localSize, &adjusted);
845: ISGetIndices(globalNumbering, &globalNum);
846: ISGetIndices(*partition, &partIdx);
847: PetscMalloc1(localSize, &map);
848: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
849: for (i = cStart, offset = 0; i < cEnd; i++) {
850: if (globalNum[i - cStart] >= 0) map[offset++] = i;
851: }
852: for (i = 0; i < localSize; i++) adjusted[i] = map[partIdx[i]];
853: PetscFree(map);
854: ISRestoreIndices(*partition, &partIdx);
855: ISRestoreIndices(globalNumbering, &globalNum);
856: ISCreateGeneral(PETSC_COMM_SELF, localSize, adjusted, PETSC_OWN_POINTER, &newPartition);
857: ISDestroy(&globalNumbering);
858: ISDestroy(partition);
859: *partition = newPartition;
860: }
861: } else SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %" PetscInt_FMT " for points to partition", part->height);
862: PetscSectionDestroy(&vertSection);
863: return 0;
864: }
866: /*@
867: DMPlexGetPartitioner - Get the mesh partitioner
869: Not collective
871: Input Parameter:
872: . dm - The `DM`
874: Output Parameter:
875: . part - The `PetscPartitioner`
877: Level: developer
879: Note:
880: This gets a borrowed reference, so the user should not destroy this `PetscPartitioner`.
882: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`, `PetscSection`, `DMPlexDistribute()`, `DMPlexSetPartitioner()`, `PetscPartitionerDMPlexPartition()`, `PetscPartitionerCreate()`
883: @*/
884: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
885: {
886: DM_Plex *mesh = (DM_Plex *)dm->data;
890: *part = mesh->partitioner;
891: return 0;
892: }
894: /*@
895: DMPlexSetPartitioner - Set the mesh partitioner
897: logically collective on dm
899: Input Parameters:
900: + dm - The `DM`
901: - part - The partitioner
903: Level: developer
905: Note:
906: Any existing `PetscPartitioner` will be destroyed.
908: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `PetscPartitioner`,`DMPlexDistribute()`, `DMPlexGetPartitioner()`, `PetscPartitionerCreate()`
909: @*/
910: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
911: {
912: DM_Plex *mesh = (DM_Plex *)dm->data;
916: PetscObjectReference((PetscObject)part);
917: PetscPartitionerDestroy(&mesh->partitioner);
918: mesh->partitioner = part;
919: return 0;
920: }
922: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
923: {
924: const PetscInt *cone;
925: PetscInt coneSize, c;
926: PetscBool missing;
929: PetscHSetIQueryAdd(ht, point, &missing);
930: if (missing) {
931: DMPlexGetCone(dm, point, &cone);
932: DMPlexGetConeSize(dm, point, &coneSize);
933: for (c = 0; c < coneSize; c++) DMPlexAddClosure_Private(dm, ht, cone[c]);
934: }
935: return 0;
936: }
938: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
939: {
940: if (up) {
941: PetscInt parent;
943: DMPlexGetTreeParent(dm, point, &parent, NULL);
944: if (parent != point) {
945: PetscInt closureSize, *closure = NULL, i;
947: DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
948: for (i = 0; i < closureSize; i++) {
949: PetscInt cpoint = closure[2 * i];
951: PetscHSetIAdd(ht, cpoint);
952: DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_TRUE, PETSC_FALSE);
953: }
954: DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
955: }
956: }
957: if (down) {
958: PetscInt numChildren;
959: const PetscInt *children;
961: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
962: if (numChildren) {
963: PetscInt i;
965: for (i = 0; i < numChildren; i++) {
966: PetscInt cpoint = children[i];
968: PetscHSetIAdd(ht, cpoint);
969: DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_FALSE, PETSC_TRUE);
970: }
971: }
972: }
973: return 0;
974: }
976: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
977: {
978: PetscInt parent;
981: DMPlexGetTreeParent(dm, point, &parent, NULL);
982: if (point != parent) {
983: const PetscInt *cone;
984: PetscInt coneSize, c;
986: DMPlexAddClosureTree_Up_Private(dm, ht, parent);
987: DMPlexAddClosure_Private(dm, ht, parent);
988: DMPlexGetCone(dm, parent, &cone);
989: DMPlexGetConeSize(dm, parent, &coneSize);
990: for (c = 0; c < coneSize; c++) {
991: const PetscInt cp = cone[c];
993: DMPlexAddClosureTree_Up_Private(dm, ht, cp);
994: }
995: }
996: return 0;
997: }
999: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
1000: {
1001: PetscInt i, numChildren;
1002: const PetscInt *children;
1005: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
1006: for (i = 0; i < numChildren; i++) PetscHSetIAdd(ht, children[i]);
1007: return 0;
1008: }
1010: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
1011: {
1012: const PetscInt *cone;
1013: PetscInt coneSize, c;
1016: PetscHSetIAdd(ht, point);
1017: DMPlexAddClosureTree_Up_Private(dm, ht, point);
1018: DMPlexAddClosureTree_Down_Private(dm, ht, point);
1019: DMPlexGetCone(dm, point, &cone);
1020: DMPlexGetConeSize(dm, point, &coneSize);
1021: for (c = 0; c < coneSize; c++) DMPlexAddClosureTree_Private(dm, ht, cone[c]);
1022: return 0;
1023: }
1025: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
1026: {
1027: DM_Plex *mesh = (DM_Plex *)dm->data;
1028: const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
1029: PetscInt nelems, *elems, off = 0, p;
1030: PetscHSetI ht = NULL;
1032: PetscHSetICreate(&ht);
1033: PetscHSetIResize(ht, numPoints * 16);
1034: if (!hasTree) {
1035: for (p = 0; p < numPoints; ++p) DMPlexAddClosure_Private(dm, ht, points[p]);
1036: } else {
1037: #if 1
1038: for (p = 0; p < numPoints; ++p) DMPlexAddClosureTree_Private(dm, ht, points[p]);
1039: #else
1040: PetscInt *closure = NULL, closureSize, c;
1041: for (p = 0; p < numPoints; ++p) {
1042: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
1043: for (c = 0; c < closureSize * 2; c += 2) {
1044: PetscHSetIAdd(ht, closure[c]);
1045: if (hasTree) DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);
1046: }
1047: }
1048: if (closure) DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);
1049: #endif
1050: }
1051: PetscHSetIGetSize(ht, &nelems);
1052: PetscMalloc1(nelems, &elems);
1053: PetscHSetIGetElems(ht, &off, elems);
1054: PetscHSetIDestroy(&ht);
1055: PetscSortInt(nelems, elems);
1056: ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
1057: return 0;
1058: }
1060: /*@
1061: DMPlexPartitionLabelClosure - Add the closure of all points to the partition label
1063: Input Parameters:
1064: + dm - The `DM`
1065: - label - `DMLabel` assigning ranks to remote roots
1067: Level: developer
1069: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1070: @*/
1071: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
1072: {
1073: IS rankIS, pointIS, closureIS;
1074: const PetscInt *ranks, *points;
1075: PetscInt numRanks, numPoints, r;
1077: DMLabelGetValueIS(label, &rankIS);
1078: ISGetLocalSize(rankIS, &numRanks);
1079: ISGetIndices(rankIS, &ranks);
1080: for (r = 0; r < numRanks; ++r) {
1081: const PetscInt rank = ranks[r];
1082: DMLabelGetStratumIS(label, rank, &pointIS);
1083: ISGetLocalSize(pointIS, &numPoints);
1084: ISGetIndices(pointIS, &points);
1085: DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
1086: ISRestoreIndices(pointIS, &points);
1087: ISDestroy(&pointIS);
1088: DMLabelSetStratumIS(label, rank, closureIS);
1089: ISDestroy(&closureIS);
1090: }
1091: ISRestoreIndices(rankIS, &ranks);
1092: ISDestroy(&rankIS);
1093: return 0;
1094: }
1096: /*@
1097: DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label
1099: Input Parameters:
1100: + dm - The `DM`
1101: - label - `DMLabel` assigning ranks to remote roots
1103: Level: developer
1105: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1106: @*/
1107: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1108: {
1109: IS rankIS, pointIS;
1110: const PetscInt *ranks, *points;
1111: PetscInt numRanks, numPoints, r, p, a, adjSize;
1112: PetscInt *adj = NULL;
1114: DMLabelGetValueIS(label, &rankIS);
1115: ISGetLocalSize(rankIS, &numRanks);
1116: ISGetIndices(rankIS, &ranks);
1117: for (r = 0; r < numRanks; ++r) {
1118: const PetscInt rank = ranks[r];
1120: DMLabelGetStratumIS(label, rank, &pointIS);
1121: ISGetLocalSize(pointIS, &numPoints);
1122: ISGetIndices(pointIS, &points);
1123: for (p = 0; p < numPoints; ++p) {
1124: adjSize = PETSC_DETERMINE;
1125: DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
1126: for (a = 0; a < adjSize; ++a) DMLabelSetValue(label, adj[a], rank);
1127: }
1128: ISRestoreIndices(pointIS, &points);
1129: ISDestroy(&pointIS);
1130: }
1131: ISRestoreIndices(rankIS, &ranks);
1132: ISDestroy(&rankIS);
1133: PetscFree(adj);
1134: return 0;
1135: }
1137: /*@
1138: DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point `PetscSF`
1140: Input Parameters:
1141: + dm - The `DM`
1142: - label - `DMLabel` assigning ranks to remote roots
1144: Level: developer
1146: Note:
1147: This is required when generating multi-level overlaps to capture
1148: overlap points from non-neighbouring partitions.
1150: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1151: @*/
1152: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1153: {
1154: MPI_Comm comm;
1155: PetscMPIInt rank;
1156: PetscSF sfPoint;
1157: DMLabel lblRoots, lblLeaves;
1158: IS rankIS, pointIS;
1159: const PetscInt *ranks;
1160: PetscInt numRanks, r;
1162: PetscObjectGetComm((PetscObject)dm, &comm);
1163: MPI_Comm_rank(comm, &rank);
1164: DMGetPointSF(dm, &sfPoint);
1165: /* Pull point contributions from remote leaves into local roots */
1166: DMLabelGather(label, sfPoint, &lblLeaves);
1167: DMLabelGetValueIS(lblLeaves, &rankIS);
1168: ISGetLocalSize(rankIS, &numRanks);
1169: ISGetIndices(rankIS, &ranks);
1170: for (r = 0; r < numRanks; ++r) {
1171: const PetscInt remoteRank = ranks[r];
1172: if (remoteRank == rank) continue;
1173: DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
1174: DMLabelInsertIS(label, pointIS, remoteRank);
1175: ISDestroy(&pointIS);
1176: }
1177: ISRestoreIndices(rankIS, &ranks);
1178: ISDestroy(&rankIS);
1179: DMLabelDestroy(&lblLeaves);
1180: /* Push point contributions from roots into remote leaves */
1181: DMLabelDistribute(label, sfPoint, &lblRoots);
1182: DMLabelGetValueIS(lblRoots, &rankIS);
1183: ISGetLocalSize(rankIS, &numRanks);
1184: ISGetIndices(rankIS, &ranks);
1185: for (r = 0; r < numRanks; ++r) {
1186: const PetscInt remoteRank = ranks[r];
1187: if (remoteRank == rank) continue;
1188: DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
1189: DMLabelInsertIS(label, pointIS, remoteRank);
1190: ISDestroy(&pointIS);
1191: }
1192: ISRestoreIndices(rankIS, &ranks);
1193: ISDestroy(&rankIS);
1194: DMLabelDestroy(&lblRoots);
1195: return 0;
1196: }
1198: /*@
1199: DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label
1201: Input Parameters:
1202: + dm - The `DM`
1203: . rootLabel - `DMLabel` assigning ranks to local roots
1204: - processSF - A star forest mapping into the local index on each remote rank
1206: Output Parameter:
1207: . leafLabel - `DMLabel `assigning ranks to remote roots
1209: Level: developer
1211: Note:
1212: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
1213: resulting leafLabel is a receiver mapping of remote roots to their parent rank.
1215: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1216: @*/
1217: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1218: {
1219: MPI_Comm comm;
1220: PetscMPIInt rank, size, r;
1221: PetscInt p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1222: PetscSF sfPoint;
1223: PetscSection rootSection;
1224: PetscSFNode *rootPoints, *leafPoints;
1225: const PetscSFNode *remote;
1226: const PetscInt *local, *neighbors;
1227: IS valueIS;
1228: PetscBool mpiOverflow = PETSC_FALSE;
1230: PetscLogEventBegin(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1231: PetscObjectGetComm((PetscObject)dm, &comm);
1232: MPI_Comm_rank(comm, &rank);
1233: MPI_Comm_size(comm, &size);
1234: DMGetPointSF(dm, &sfPoint);
1236: /* Convert to (point, rank) and use actual owners */
1237: PetscSectionCreate(comm, &rootSection);
1238: PetscSectionSetChart(rootSection, 0, size);
1239: DMLabelGetValueIS(rootLabel, &valueIS);
1240: ISGetLocalSize(valueIS, &numNeighbors);
1241: ISGetIndices(valueIS, &neighbors);
1242: for (n = 0; n < numNeighbors; ++n) {
1243: DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
1244: PetscSectionAddDof(rootSection, neighbors[n], numPoints);
1245: }
1246: PetscSectionSetUp(rootSection);
1247: PetscSectionGetStorageSize(rootSection, &rootSize);
1248: PetscMalloc1(rootSize, &rootPoints);
1249: PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
1250: for (n = 0; n < numNeighbors; ++n) {
1251: IS pointIS;
1252: const PetscInt *points;
1254: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1255: DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
1256: ISGetLocalSize(pointIS, &numPoints);
1257: ISGetIndices(pointIS, &points);
1258: for (p = 0; p < numPoints; ++p) {
1259: if (local) PetscFindInt(points[p], nleaves, local, &l);
1260: else l = -1;
1261: if (l >= 0) {
1262: rootPoints[off + p] = remote[l];
1263: } else {
1264: rootPoints[off + p].index = points[p];
1265: rootPoints[off + p].rank = rank;
1266: }
1267: }
1268: ISRestoreIndices(pointIS, &points);
1269: ISDestroy(&pointIS);
1270: }
1272: /* Try to communicate overlap using All-to-All */
1273: if (!processSF) {
1274: PetscInt64 counter = 0;
1275: PetscBool locOverflow = PETSC_FALSE;
1276: PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;
1278: PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
1279: for (n = 0; n < numNeighbors; ++n) {
1280: PetscSectionGetDof(rootSection, neighbors[n], &dof);
1281: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1282: #if defined(PETSC_USE_64BIT_INDICES)
1283: if (dof > PETSC_MPI_INT_MAX) {
1284: locOverflow = PETSC_TRUE;
1285: break;
1286: }
1287: if (off > PETSC_MPI_INT_MAX) {
1288: locOverflow = PETSC_TRUE;
1289: break;
1290: }
1291: #endif
1292: scounts[neighbors[n]] = (PetscMPIInt)dof;
1293: sdispls[neighbors[n]] = (PetscMPIInt)off;
1294: }
1295: MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
1296: for (r = 0; r < size; ++r) {
1297: rdispls[r] = (int)counter;
1298: counter += rcounts[r];
1299: }
1300: if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
1301: MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
1302: if (!mpiOverflow) {
1303: PetscInfo(dm, "Using Alltoallv for mesh distribution\n");
1304: leafSize = (PetscInt)counter;
1305: PetscMalloc1(leafSize, &leafPoints);
1306: MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
1307: }
1308: PetscFree4(scounts, sdispls, rcounts, rdispls);
1309: }
1311: /* Communicate overlap using process star forest */
1312: if (processSF || mpiOverflow) {
1313: PetscSF procSF;
1314: PetscSection leafSection;
1316: if (processSF) {
1317: PetscInfo(dm, "Using processSF for mesh distribution\n");
1318: PetscObjectReference((PetscObject)processSF);
1319: procSF = processSF;
1320: } else {
1321: PetscInfo(dm, "Using processSF for mesh distribution (MPI overflow)\n");
1322: PetscSFCreate(comm, &procSF);
1323: PetscSFSetGraphWithPattern(procSF, NULL, PETSCSF_PATTERN_ALLTOALL);
1324: }
1326: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
1327: DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void **)&leafPoints);
1328: PetscSectionGetStorageSize(leafSection, &leafSize);
1329: PetscSectionDestroy(&leafSection);
1330: PetscSFDestroy(&procSF);
1331: }
1333: for (p = 0; p < leafSize; p++) DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
1335: ISRestoreIndices(valueIS, &neighbors);
1336: ISDestroy(&valueIS);
1337: PetscSectionDestroy(&rootSection);
1338: PetscFree(rootPoints);
1339: PetscFree(leafPoints);
1340: PetscLogEventEnd(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1341: return 0;
1342: }
1344: /*@
1345: DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points
1347: Input Parameters:
1348: + dm - The `DM`
1349: - label - `DMLabel` assigning ranks to remote roots
1351: Output Parameter:
1352: . sf - The star forest communication context encapsulating the defined mapping
1354: Level: developer
1356: Note:
1357: The incoming label is a receiver mapping of remote points to their parent rank.
1359: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1360: @*/
1361: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
1362: {
1363: PetscMPIInt rank;
1364: PetscInt n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1365: PetscSFNode *remotePoints;
1366: IS remoteRootIS, neighborsIS;
1367: const PetscInt *remoteRoots, *neighbors;
1369: PetscLogEventBegin(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1370: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1372: DMLabelGetValueIS(label, &neighborsIS);
1373: #if 0
1374: {
1375: IS is;
1376: ISDuplicate(neighborsIS, &is);
1377: ISSort(is);
1378: ISDestroy(&neighborsIS);
1379: neighborsIS = is;
1380: }
1381: #endif
1382: ISGetLocalSize(neighborsIS, &nNeighbors);
1383: ISGetIndices(neighborsIS, &neighbors);
1384: for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1385: DMLabelGetStratumSize(label, neighbors[n], &numPoints);
1386: numRemote += numPoints;
1387: }
1388: PetscMalloc1(numRemote, &remotePoints);
1389: /* Put owned points first */
1390: DMLabelGetStratumSize(label, rank, &numPoints);
1391: if (numPoints > 0) {
1392: DMLabelGetStratumIS(label, rank, &remoteRootIS);
1393: ISGetIndices(remoteRootIS, &remoteRoots);
1394: for (p = 0; p < numPoints; p++) {
1395: remotePoints[idx].index = remoteRoots[p];
1396: remotePoints[idx].rank = rank;
1397: idx++;
1398: }
1399: ISRestoreIndices(remoteRootIS, &remoteRoots);
1400: ISDestroy(&remoteRootIS);
1401: }
1402: /* Now add remote points */
1403: for (n = 0; n < nNeighbors; ++n) {
1404: const PetscInt nn = neighbors[n];
1406: DMLabelGetStratumSize(label, nn, &numPoints);
1407: if (nn == rank || numPoints <= 0) continue;
1408: DMLabelGetStratumIS(label, nn, &remoteRootIS);
1409: ISGetIndices(remoteRootIS, &remoteRoots);
1410: for (p = 0; p < numPoints; p++) {
1411: remotePoints[idx].index = remoteRoots[p];
1412: remotePoints[idx].rank = nn;
1413: idx++;
1414: }
1415: ISRestoreIndices(remoteRootIS, &remoteRoots);
1416: ISDestroy(&remoteRootIS);
1417: }
1418: PetscSFCreate(PetscObjectComm((PetscObject)dm), sf);
1419: DMPlexGetChart(dm, &pStart, &pEnd);
1420: PetscSFSetGraph(*sf, pEnd - pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1421: PetscSFSetUp(*sf);
1422: ISDestroy(&neighborsIS);
1423: PetscLogEventEnd(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1424: return 0;
1425: }
1427: #if defined(PETSC_HAVE_PARMETIS)
1428: #include <parmetis.h>
1429: #endif
1431: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
1432: * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
1433: * them out in that case. */
1434: #if defined(PETSC_HAVE_PARMETIS)
1435: /*@C
1437: DMPlexRewriteSF - Rewrites the ownership of the `PetsSF` of a `DM` (in place).
1439: Input parameters:
1440: + dm - The `DMPLEX` object.
1441: . n - The number of points.
1442: . pointsToRewrite - The points in the `PetscSF` whose ownership will change.
1443: . targetOwners - New owner for each element in pointsToRewrite.
1444: - degrees - Degrees of the points in the `PetscSF` as obtained by `PetscSFComputeDegreeBegin()`/`PetscSFComputeDegreeEnd()`.
1446: Level: developer
1448: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMLabel`, `PetscSF`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1449: @*/
1450: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1451: {
1452: PetscInt pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1453: PetscInt *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1454: PetscSFNode *leafLocationsNew;
1455: const PetscSFNode *iremote;
1456: const PetscInt *ilocal;
1457: PetscBool *isLeaf;
1458: PetscSF sf;
1459: MPI_Comm comm;
1460: PetscMPIInt rank, size;
1462: PetscObjectGetComm((PetscObject)dm, &comm);
1463: MPI_Comm_rank(comm, &rank);
1464: MPI_Comm_size(comm, &size);
1465: DMPlexGetChart(dm, &pStart, &pEnd);
1467: DMGetPointSF(dm, &sf);
1468: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1469: PetscMalloc1(pEnd - pStart, &isLeaf);
1470: for (i = 0; i < pEnd - pStart; i++) isLeaf[i] = PETSC_FALSE;
1471: for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;
1473: PetscMalloc1(pEnd - pStart + 1, &cumSumDegrees);
1474: cumSumDegrees[0] = 0;
1475: for (i = 1; i <= pEnd - pStart; i++) cumSumDegrees[i] = cumSumDegrees[i - 1] + degrees[i - 1];
1476: sumDegrees = cumSumDegrees[pEnd - pStart];
1477: /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */
1479: PetscMalloc1(sumDegrees, &locationsOfLeafs);
1480: PetscMalloc1(pEnd - pStart, &rankOnLeafs);
1481: for (i = 0; i < pEnd - pStart; i++) rankOnLeafs[i] = rank;
1482: PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1483: PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1484: PetscFree(rankOnLeafs);
1486: /* get the remote local points of my leaves */
1487: PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
1488: PetscMalloc1(pEnd - pStart, &points);
1489: for (i = 0; i < pEnd - pStart; i++) points[i] = pStart + i;
1490: PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1491: PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1492: PetscFree(points);
1493: /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1494: PetscMalloc1(pEnd - pStart, &newOwners);
1495: PetscMalloc1(pEnd - pStart, &newNumbers);
1496: for (i = 0; i < pEnd - pStart; i++) {
1497: newOwners[i] = -1;
1498: newNumbers[i] = -1;
1499: }
1500: {
1501: PetscInt oldNumber, newNumber, oldOwner, newOwner;
1502: for (i = 0; i < n; i++) {
1503: oldNumber = pointsToRewrite[i];
1504: newNumber = -1;
1505: oldOwner = rank;
1506: newOwner = targetOwners[i];
1507: if (oldOwner == newOwner) {
1508: newNumber = oldNumber;
1509: } else {
1510: for (j = 0; j < degrees[oldNumber]; j++) {
1511: if (locationsOfLeafs[cumSumDegrees[oldNumber] + j] == newOwner) {
1512: newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber] + j];
1513: break;
1514: }
1515: }
1516: }
1519: newOwners[oldNumber] = newOwner;
1520: newNumbers[oldNumber] = newNumber;
1521: }
1522: }
1523: PetscFree(cumSumDegrees);
1524: PetscFree(locationsOfLeafs);
1525: PetscFree(remoteLocalPointOfLeafs);
1527: PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1528: PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1529: PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);
1530: PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);
1532: /* Now count how many leafs we have on each processor. */
1533: leafCounter = 0;
1534: for (i = 0; i < pEnd - pStart; i++) {
1535: if (newOwners[i] >= 0) {
1536: if (newOwners[i] != rank) leafCounter++;
1537: } else {
1538: if (isLeaf[i]) leafCounter++;
1539: }
1540: }
1542: /* Now set up the new sf by creating the leaf arrays */
1543: PetscMalloc1(leafCounter, &leafsNew);
1544: PetscMalloc1(leafCounter, &leafLocationsNew);
1546: leafCounter = 0;
1547: counter = 0;
1548: for (i = 0; i < pEnd - pStart; i++) {
1549: if (newOwners[i] >= 0) {
1550: if (newOwners[i] != rank) {
1551: leafsNew[leafCounter] = i;
1552: leafLocationsNew[leafCounter].rank = newOwners[i];
1553: leafLocationsNew[leafCounter].index = newNumbers[i];
1554: leafCounter++;
1555: }
1556: } else {
1557: if (isLeaf[i]) {
1558: leafsNew[leafCounter] = i;
1559: leafLocationsNew[leafCounter].rank = iremote[counter].rank;
1560: leafLocationsNew[leafCounter].index = iremote[counter].index;
1561: leafCounter++;
1562: }
1563: }
1564: if (isLeaf[i]) counter++;
1565: }
1567: PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
1568: PetscFree(newOwners);
1569: PetscFree(newNumbers);
1570: PetscFree(isLeaf);
1571: return 0;
1572: }
1574: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1575: {
1576: PetscInt *distribution, min, max, sum;
1577: PetscMPIInt rank, size;
1579: MPI_Comm_size(comm, &size);
1580: MPI_Comm_rank(comm, &rank);
1581: PetscCalloc1(size, &distribution);
1582: for (PetscInt i = 0; i < n; i++) {
1583: if (part) distribution[part[i]] += vtxwgt[skip * i];
1584: else distribution[rank] += vtxwgt[skip * i];
1585: }
1586: MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
1587: min = distribution[0];
1588: max = distribution[0];
1589: sum = distribution[0];
1590: for (PetscInt i = 1; i < size; i++) {
1591: if (distribution[i] < min) min = distribution[i];
1592: if (distribution[i] > max) max = distribution[i];
1593: sum += distribution[i];
1594: }
1595: PetscViewerASCIIPrintf(viewer, "Min: %" PetscInt_FMT ", Avg: %" PetscInt_FMT ", Max: %" PetscInt_FMT ", Balance: %f\n", min, sum / size, max, (max * 1. * size) / sum);
1596: PetscFree(distribution);
1597: return 0;
1598: }
1600: #endif
1602: /*@
1603: DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the `PointSF` of the `DM` inplace.
1605: Input parameters:
1606: + dm - The `DMPLEX` object.
1607: . entityDepth - depth of the entity to balance (0 -> balance vertices).
1608: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
1609: - parallel - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.
1611: Output parameters:
1612: . success - whether the graph partitioning was successful or not, optional. Unsuccessful simply means no change to the partitioning
1614: Options Database Keys:
1615: + -dm_plex_rebalance_shared_points_parmetis - Use ParMetis instead of Metis for the partitioner
1616: . -dm_plex_rebalance_shared_points_use_initial_guess - Use current partition to bootstrap ParMetis partition
1617: . -dm_plex_rebalance_shared_points_use_mat_partitioning - Use the MatPartitioning object to perform the partition, the prefix for those operations is -dm_plex_rebalance_shared_points_
1618: - -dm_plex_rebalance_shared_points_monitor - Monitor the shared points rebalance process
1620: Level: intermediate
1622: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1623: @*/
1624: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1625: {
1626: #if defined(PETSC_HAVE_PARMETIS)
1627: PetscSF sf;
1628: PetscInt ierr, i, j, idx, jdx;
1629: PetscInt eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1630: const PetscInt *degrees, *ilocal;
1631: const PetscSFNode *iremote;
1632: PetscBool *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1633: PetscInt numExclusivelyOwned, numNonExclusivelyOwned;
1634: PetscMPIInt rank, size;
1635: PetscInt *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1636: const PetscInt *cumSumVertices;
1637: PetscInt offset, counter;
1638: PetscInt *vtxwgt;
1639: const PetscInt *xadj, *adjncy;
1640: PetscInt *part, *options;
1641: PetscInt nparts, wgtflag, numflag, ncon, edgecut;
1642: real_t *ubvec;
1643: PetscInt *firstVertices, *renumbering;
1644: PetscInt failed, failedGlobal;
1645: MPI_Comm comm;
1646: Mat A;
1647: PetscViewer viewer;
1648: PetscViewerFormat format;
1649: PetscLayout layout;
1650: real_t *tpwgts;
1651: PetscMPIInt *counts, *mpiCumSumVertices;
1652: PetscInt *pointsToRewrite;
1653: PetscInt numRows;
1654: PetscBool done, usematpartitioning = PETSC_FALSE;
1655: IS ispart = NULL;
1656: MatPartitioning mp;
1657: const char *prefix;
1659: PetscObjectGetComm((PetscObject)dm, &comm);
1660: MPI_Comm_size(comm, &size);
1661: if (size == 1) {
1662: if (success) *success = PETSC_TRUE;
1663: return 0;
1664: }
1665: if (success) *success = PETSC_FALSE;
1666: MPI_Comm_rank(comm, &rank);
1668: parallel = PETSC_FALSE;
1669: useInitialGuess = PETSC_FALSE;
1670: PetscObjectOptionsBegin((PetscObject)dm);
1671: PetscOptionsName("-dm_plex_rebalance_shared_points_parmetis", "Use ParMetis instead of Metis for the partitioner", "DMPlexRebalanceSharedPoints", ¶llel);
1672: PetscOptionsBool("-dm_plex_rebalance_shared_points_use_initial_guess", "Use current partition to bootstrap ParMetis partition", "DMPlexRebalanceSharedPoints", useInitialGuess, &useInitialGuess, NULL);
1673: PetscOptionsBool("-dm_plex_rebalance_shared_points_use_mat_partitioning", "Use the MatPartitioning object to partition", "DMPlexRebalanceSharedPoints", usematpartitioning, &usematpartitioning, NULL);
1674: PetscOptionsViewer("-dm_plex_rebalance_shared_points_monitor", "Monitor the shared points rebalance process", "DMPlexRebalanceSharedPoints", &viewer, &format, NULL);
1675: PetscOptionsEnd();
1676: if (viewer) PetscViewerPushFormat(viewer, format);
1678: PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1680: DMGetOptionsPrefix(dm, &prefix);
1681: PetscOptionsGetViewer(comm, ((PetscObject)dm)->options, prefix, "-dm_rebalance_partition_view", &viewer, &format, NULL);
1682: if (viewer) PetscViewerPushFormat(viewer, format);
1684: /* Figure out all points in the plex that we are interested in balancing. */
1685: DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
1686: DMPlexGetChart(dm, &pStart, &pEnd);
1687: PetscMalloc1(pEnd - pStart, &toBalance);
1688: for (i = 0; i < pEnd - pStart; i++) toBalance[i] = (PetscBool)(i >= eBegin && i < eEnd);
1690: /* There are three types of points:
1691: * exclusivelyOwned: points that are owned by this process and only seen by this process
1692: * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1693: * leaf: a point that is seen by this process but owned by a different process
1694: */
1695: DMGetPointSF(dm, &sf);
1696: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1697: PetscMalloc1(pEnd - pStart, &isLeaf);
1698: PetscMalloc1(pEnd - pStart, &isNonExclusivelyOwned);
1699: PetscMalloc1(pEnd - pStart, &isExclusivelyOwned);
1700: for (i = 0; i < pEnd - pStart; i++) {
1701: isNonExclusivelyOwned[i] = PETSC_FALSE;
1702: isExclusivelyOwned[i] = PETSC_FALSE;
1703: isLeaf[i] = PETSC_FALSE;
1704: }
1706: /* mark all the leafs */
1707: for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;
1709: /* for an owned point, we can figure out whether another processor sees it or
1710: * not by calculating its degree */
1711: PetscSFComputeDegreeBegin(sf, °rees);
1712: PetscSFComputeDegreeEnd(sf, °rees);
1713: numExclusivelyOwned = 0;
1714: numNonExclusivelyOwned = 0;
1715: for (i = 0; i < pEnd - pStart; i++) {
1716: if (toBalance[i]) {
1717: if (degrees[i] > 0) {
1718: isNonExclusivelyOwned[i] = PETSC_TRUE;
1719: numNonExclusivelyOwned += 1;
1720: } else {
1721: if (!isLeaf[i]) {
1722: isExclusivelyOwned[i] = PETSC_TRUE;
1723: numExclusivelyOwned += 1;
1724: }
1725: }
1726: }
1727: }
1729: /* Build a graph with one vertex per core representing the
1730: * exclusively owned points and then one vertex per nonExclusively owned
1731: * point. */
1732: PetscLayoutCreate(comm, &layout);
1733: PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
1734: PetscLayoutSetUp(layout);
1735: PetscLayoutGetRanges(layout, &cumSumVertices);
1736: PetscMalloc1(pEnd - pStart, &globalNumbersOfLocalOwnedVertices);
1737: for (i = 0; i < pEnd - pStart; i++) globalNumbersOfLocalOwnedVertices[i] = pStart - 1;
1738: offset = cumSumVertices[rank];
1739: counter = 0;
1740: for (i = 0; i < pEnd - pStart; i++) {
1741: if (toBalance[i]) {
1742: if (degrees[i] > 0) {
1743: globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1744: counter++;
1745: }
1746: }
1747: }
1749: /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
1750: PetscMalloc1(pEnd - pStart, &leafGlobalNumbers);
1751: PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE);
1752: PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE);
1754: /* Build the graph for partitioning */
1755: numRows = 1 + numNonExclusivelyOwned;
1756: PetscLogEventBegin(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);
1757: MatCreate(comm, &A);
1758: MatSetType(A, MATMPIADJ);
1759: MatSetSizes(A, numRows, numRows, cumSumVertices[size], cumSumVertices[size]);
1760: idx = cumSumVertices[rank];
1761: for (i = 0; i < pEnd - pStart; i++) {
1762: if (toBalance[i]) {
1763: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1764: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1765: else continue;
1766: MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
1767: MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
1768: }
1769: }
1770: PetscFree(globalNumbersOfLocalOwnedVertices);
1771: PetscFree(leafGlobalNumbers);
1772: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1773: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1774: PetscLogEventEnd(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);
1776: nparts = size;
1777: ncon = 1;
1778: PetscMalloc1(ncon * nparts, &tpwgts);
1779: for (i = 0; i < ncon * nparts; i++) tpwgts[i] = 1. / (nparts);
1780: PetscMalloc1(ncon, &ubvec);
1781: ubvec[0] = 1.05;
1782: ubvec[1] = 1.05;
1784: PetscMalloc1(ncon * (1 + numNonExclusivelyOwned), &vtxwgt);
1785: if (ncon == 2) {
1786: vtxwgt[0] = numExclusivelyOwned;
1787: vtxwgt[1] = 1;
1788: for (i = 0; i < numNonExclusivelyOwned; i++) {
1789: vtxwgt[ncon * (i + 1)] = 1;
1790: vtxwgt[ncon * (i + 1) + 1] = 0;
1791: }
1792: } else {
1793: PetscInt base, ms;
1794: MPI_Allreduce(&numExclusivelyOwned, &base, 1, MPIU_INT, MPIU_MAX, PetscObjectComm((PetscObject)dm));
1795: MatGetSize(A, &ms, NULL);
1796: ms -= size;
1797: base = PetscMax(base, ms);
1798: vtxwgt[0] = base + numExclusivelyOwned;
1799: for (i = 0; i < numNonExclusivelyOwned; i++) vtxwgt[i + 1] = 1;
1800: }
1802: if (viewer) {
1803: PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %" PetscInt_FMT " on interface of mesh distribution.\n", entityDepth);
1804: PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %" PetscInt_FMT "\n", cumSumVertices[size]);
1805: }
1806: /* TODO: Drop the parallel/sequential choice here and just use MatPartioner for much more flexibility */
1807: if (usematpartitioning) {
1808: const char *prefix;
1810: MatPartitioningCreate(PetscObjectComm((PetscObject)dm), &mp);
1811: PetscObjectSetOptionsPrefix((PetscObject)mp, "dm_plex_rebalance_shared_points_");
1812: PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
1813: PetscObjectPrependOptionsPrefix((PetscObject)mp, prefix);
1814: MatPartitioningSetAdjacency(mp, A);
1815: MatPartitioningSetNumberVertexWeights(mp, ncon);
1816: MatPartitioningSetVertexWeights(mp, vtxwgt);
1817: MatPartitioningSetFromOptions(mp);
1818: MatPartitioningApply(mp, &ispart);
1819: ISGetIndices(ispart, (const PetscInt **)&part);
1820: } else if (parallel) {
1821: if (viewer) PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n");
1822: PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1823: MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1825: PetscMalloc1(4, &options);
1826: options[0] = 1;
1827: options[1] = 0; /* Verbosity */
1828: if (viewer) options[1] = 1;
1829: options[2] = 0; /* Seed */
1830: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1831: wgtflag = 2;
1832: numflag = 0;
1833: if (useInitialGuess) {
1834: PetscViewerASCIIPrintf(viewer, "THIS DOES NOT WORK! I don't know why. Using current distribution of points as initial guess.\n");
1835: for (i = 0; i < numRows; i++) part[i] = rank;
1836: if (viewer) PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n");
1837: PetscStackPushExternal("ParMETIS_V3_RefineKway");
1838: PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1839: ParMETIS_V3_RefineKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1840: PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1841: PetscStackPop;
1843: } else {
1844: PetscStackPushExternal("ParMETIS_V3_PartKway");
1845: PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1846: ParMETIS_V3_PartKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1847: PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1848: PetscStackPop;
1850: }
1851: MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1852: PetscFree(options);
1853: } else {
1854: if (viewer) PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n");
1855: Mat As;
1856: PetscInt *partGlobal;
1857: PetscInt *numExclusivelyOwnedAll;
1859: PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1860: MatGetSize(A, &numRows, NULL);
1861: PetscLogEventBegin(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);
1862: MatMPIAdjToSeqRankZero(A, &As);
1863: PetscLogEventEnd(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);
1865: PetscMalloc1(size, &numExclusivelyOwnedAll);
1866: numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1867: MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, numExclusivelyOwnedAll, 1, MPIU_INT, comm);
1869: PetscMalloc1(numRows, &partGlobal);
1870: PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1871: if (rank == 0) {
1872: PetscInt *vtxwgt_g, numRows_g;
1874: MatGetRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1875: PetscMalloc1(2 * numRows_g, &vtxwgt_g);
1876: for (i = 0; i < size; i++) {
1877: vtxwgt_g[ncon * cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1878: if (ncon > 1) vtxwgt_g[ncon * cumSumVertices[i] + 1] = 1;
1879: for (j = cumSumVertices[i] + 1; j < cumSumVertices[i + 1]; j++) {
1880: vtxwgt_g[ncon * j] = 1;
1881: if (ncon > 1) vtxwgt_g[2 * j + 1] = 0;
1882: }
1883: }
1885: PetscMalloc1(64, &options);
1886: METIS_SetDefaultOptions(options); /* initialize all defaults */
1888: options[METIS_OPTION_CONTIG] = 1;
1889: PetscStackPushExternal("METIS_PartGraphKway");
1890: METIS_PartGraphKway(&numRows_g, &ncon, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1891: PetscStackPop;
1893: PetscFree(options);
1894: PetscFree(vtxwgt_g);
1895: MatRestoreRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1896: MatDestroy(&As);
1897: }
1898: PetscBarrier((PetscObject)dm);
1899: PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1900: PetscFree(numExclusivelyOwnedAll);
1902: /* scatter the partitioning information to ranks */
1903: PetscLogEventBegin(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1904: PetscMalloc1(size, &counts);
1905: PetscMalloc1(size + 1, &mpiCumSumVertices);
1906: for (i = 0; i < size; i++) PetscMPIIntCast(cumSumVertices[i + 1] - cumSumVertices[i], &(counts[i]));
1907: for (i = 0; i <= size; i++) PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
1908: MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
1909: PetscFree(counts);
1910: PetscFree(mpiCumSumVertices);
1911: PetscFree(partGlobal);
1912: PetscLogEventEnd(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1913: }
1914: PetscFree(ubvec);
1915: PetscFree(tpwgts);
1917: /* Rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1918: PetscMalloc2(size, &firstVertices, size, &renumbering);
1919: firstVertices[rank] = part[0];
1920: MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, firstVertices, 1, MPIU_INT, comm);
1921: for (i = 0; i < size; i++) renumbering[firstVertices[i]] = i;
1922: for (i = 0; i < cumSumVertices[rank + 1] - cumSumVertices[rank]; i++) part[i] = renumbering[part[i]];
1923: PetscFree2(firstVertices, renumbering);
1925: /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1926: failed = (PetscInt)(part[0] != rank);
1927: MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
1928: if (failedGlobal > 0) {
1930: PetscFree(vtxwgt);
1931: PetscFree(toBalance);
1932: PetscFree(isLeaf);
1933: PetscFree(isNonExclusivelyOwned);
1934: PetscFree(isExclusivelyOwned);
1935: if (usematpartitioning) {
1936: ISRestoreIndices(ispart, (const PetscInt **)&part);
1937: ISDestroy(&ispart);
1938: } else PetscFree(part);
1939: if (viewer) {
1940: PetscViewerPopFormat(viewer);
1941: PetscViewerDestroy(&viewer);
1942: }
1943: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1944: return 0;
1945: }
1947: /* Check how well we did distributing points*/
1948: if (viewer) {
1949: PetscViewerASCIIPrintf(viewer, "Number of owned entities of depth %" PetscInt_FMT ".\n", entityDepth);
1950: PetscViewerASCIIPrintf(viewer, "Initial ");
1951: DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
1952: PetscViewerASCIIPrintf(viewer, "Rebalanced ");
1953: DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1954: }
1956: /* Check that every vertex is owned by a process that it is actually connected to. */
1957: MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1958: for (i = 1; i <= numNonExclusivelyOwned; i++) {
1959: PetscInt loc = 0;
1960: PetscFindInt(cumSumVertices[part[i]], xadj[i + 1] - xadj[i], &adjncy[xadj[i]], &loc);
1961: /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1962: if (loc < 0) part[i] = rank;
1963: }
1964: MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1965: MatDestroy(&A);
1967: /* See how significant the influences of the previous fixing up step was.*/
1968: if (viewer) {
1969: PetscViewerASCIIPrintf(viewer, "After fix ");
1970: DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1971: }
1972: if (!usematpartitioning) PetscFree(vtxwgt);
1973: else MatPartitioningDestroy(&mp);
1975: PetscLayoutDestroy(&layout);
1977: PetscLogEventBegin(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);
1978: /* Rewrite the SF to reflect the new ownership. */
1979: PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
1980: counter = 0;
1981: for (i = 0; i < pEnd - pStart; i++) {
1982: if (toBalance[i]) {
1983: if (isNonExclusivelyOwned[i]) {
1984: pointsToRewrite[counter] = i + pStart;
1985: counter++;
1986: }
1987: }
1988: }
1989: DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part + 1, degrees);
1990: PetscFree(pointsToRewrite);
1991: PetscLogEventEnd(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);
1993: PetscFree(toBalance);
1994: PetscFree(isLeaf);
1995: PetscFree(isNonExclusivelyOwned);
1996: PetscFree(isExclusivelyOwned);
1997: if (usematpartitioning) {
1998: ISRestoreIndices(ispart, (const PetscInt **)&part);
1999: ISDestroy(&ispart);
2000: } else PetscFree(part);
2001: if (viewer) {
2002: PetscViewerPopFormat(viewer);
2003: PetscViewerDestroy(&viewer);
2004: }
2005: if (success) *success = PETSC_TRUE;
2006: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
2007: return 0;
2008: #else
2009: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2010: #endif
2011: }