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: Input Parameters:
478: + dm - The mesh DM dm
479: - height - Height of the strata from which to construct the graph
481: Output Parameters:
482: + numVertices - Number of vertices in the graph
483: . offsets - Point offsets in the graph
484: . adjacency - Point connectivity in the graph
485: - 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.
487: The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function
488: representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().
490: Options Database Keys:
491: . -dm_plex_csr_alg <mat,graph,overlap> - Choose the algorithm for computing the CSR graph
493: Level: developer
495: .seealso: `PetscPartitionerGetType()`, `PetscPartitionerCreate()`, `DMSetAdjacency()`
496: @*/
497: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
498: {
499: DMPlexCSRAlgorithm alg = DM_PLEX_CSR_GRAPH;
501: PetscOptionsGetEnum(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *)&alg, NULL);
502: switch (alg) {
503: case DM_PLEX_CSR_MAT:
504: DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
505: break;
506: case DM_PLEX_CSR_GRAPH:
507: DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
508: break;
509: case DM_PLEX_CSR_OVERLAP:
510: DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering);
511: break;
512: }
513: return 0;
514: }
516: /*@C
517: DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.
519: Collective on DM
521: Input Parameters:
522: + dm - The DMPlex
523: - cellHeight - The height of mesh points to treat as cells (default should be 0)
525: Output Parameters:
526: + numVertices - The number of local vertices in the graph, or cells in the mesh.
527: . offsets - The offset to the adjacency list for each cell
528: - adjacency - The adjacency list for all cells
530: Note: This is suitable for input to a mesh partitioner like ParMetis.
532: Level: advanced
534: .seealso: `DMPlexCreate()`
535: @*/
536: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
537: {
538: const PetscInt maxFaceCases = 30;
539: PetscInt numFaceCases = 0;
540: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
541: PetscInt *off, *adj;
542: PetscInt *neighborCells = NULL;
543: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
545: /* For parallel partitioning, I think you have to communicate supports */
546: DMGetDimension(dm, &dim);
547: cellDim = dim - cellHeight;
548: DMPlexGetDepth(dm, &depth);
549: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
550: if (cEnd - cStart == 0) {
551: if (numVertices) *numVertices = 0;
552: if (offsets) *offsets = NULL;
553: if (adjacency) *adjacency = NULL;
554: return 0;
555: }
556: numCells = cEnd - cStart;
557: faceDepth = depth - cellHeight;
558: if (dim == depth) {
559: PetscInt f, fStart, fEnd;
561: PetscCalloc1(numCells + 1, &off);
562: /* Count neighboring cells */
563: DMPlexGetHeightStratum(dm, cellHeight + 1, &fStart, &fEnd);
564: for (f = fStart; f < fEnd; ++f) {
565: const PetscInt *support;
566: PetscInt supportSize;
567: DMPlexGetSupportSize(dm, f, &supportSize);
568: DMPlexGetSupport(dm, f, &support);
569: if (supportSize == 2) {
570: PetscInt numChildren;
572: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
573: if (!numChildren) {
574: ++off[support[0] - cStart + 1];
575: ++off[support[1] - cStart + 1];
576: }
577: }
578: }
579: /* Prefix sum */
580: for (c = 1; c <= numCells; ++c) off[c] += off[c - 1];
581: if (adjacency) {
582: PetscInt *tmp;
584: PetscMalloc1(off[numCells], &adj);
585: PetscMalloc1(numCells + 1, &tmp);
586: PetscArraycpy(tmp, off, numCells + 1);
587: /* Get neighboring cells */
588: for (f = fStart; f < fEnd; ++f) {
589: const PetscInt *support;
590: PetscInt supportSize;
591: DMPlexGetSupportSize(dm, f, &supportSize);
592: DMPlexGetSupport(dm, f, &support);
593: if (supportSize == 2) {
594: PetscInt numChildren;
596: DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
597: if (!numChildren) {
598: adj[tmp[support[0] - cStart]++] = support[1];
599: adj[tmp[support[1] - cStart]++] = support[0];
600: }
601: }
602: }
603: 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);
604: PetscFree(tmp);
605: }
606: if (numVertices) *numVertices = numCells;
607: if (offsets) *offsets = off;
608: if (adjacency) *adjacency = adj;
609: return 0;
610: }
611: /* Setup face recognition */
612: if (faceDepth == 1) {
613: 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 */
615: for (c = cStart; c < cEnd; ++c) {
616: PetscInt corners;
618: DMPlexGetConeSize(dm, c, &corners);
619: if (!cornersSeen[corners]) {
620: PetscInt nFV;
623: cornersSeen[corners] = 1;
625: DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);
627: numFaceVertices[numFaceCases++] = nFV;
628: }
629: }
630: }
631: PetscCalloc1(numCells + 1, &off);
632: /* Count neighboring cells */
633: for (cell = cStart; cell < cEnd; ++cell) {
634: PetscInt numNeighbors = PETSC_DETERMINE, n;
636: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
637: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
638: for (n = 0; n < numNeighbors; ++n) {
639: PetscInt cellPair[2];
640: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
641: PetscInt meetSize = 0;
642: const PetscInt *meet = NULL;
644: cellPair[0] = cell;
645: cellPair[1] = neighborCells[n];
646: if (cellPair[0] == cellPair[1]) continue;
647: if (!found) {
648: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
649: if (meetSize) {
650: PetscInt f;
652: for (f = 0; f < numFaceCases; ++f) {
653: if (numFaceVertices[f] == meetSize) {
654: found = PETSC_TRUE;
655: break;
656: }
657: }
658: }
659: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
660: }
661: if (found) ++off[cell - cStart + 1];
662: }
663: }
664: /* Prefix sum */
665: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell - 1];
667: if (adjacency) {
668: PetscMalloc1(off[numCells], &adj);
669: /* Get neighboring cells */
670: for (cell = cStart; cell < cEnd; ++cell) {
671: PetscInt numNeighbors = PETSC_DETERMINE, n;
672: PetscInt cellOffset = 0;
674: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
675: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
676: for (n = 0; n < numNeighbors; ++n) {
677: PetscInt cellPair[2];
678: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
679: PetscInt meetSize = 0;
680: const PetscInt *meet = NULL;
682: cellPair[0] = cell;
683: cellPair[1] = neighborCells[n];
684: if (cellPair[0] == cellPair[1]) continue;
685: if (!found) {
686: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
687: if (meetSize) {
688: PetscInt f;
690: for (f = 0; f < numFaceCases; ++f) {
691: if (numFaceVertices[f] == meetSize) {
692: found = PETSC_TRUE;
693: break;
694: }
695: }
696: }
697: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
698: }
699: if (found) {
700: adj[off[cell - cStart] + cellOffset] = neighborCells[n];
701: ++cellOffset;
702: }
703: }
704: }
705: }
706: PetscFree(neighborCells);
707: if (numVertices) *numVertices = numCells;
708: if (offsets) *offsets = off;
709: if (adjacency) *adjacency = adj;
710: return 0;
711: }
713: /*@
714: PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh
716: Collective on PetscPartitioner
718: Input Parameters:
719: + part - The PetscPartitioner
720: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
721: - dm - The mesh DM
723: Output Parameters:
724: + partSection - The PetscSection giving the division of points by partition
725: - partition - The list of points by partition
727: Notes:
728: If the DM has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
729: by the section in the transitive closure of the point.
731: Level: developer
733: .seealso `DMPlexDistribute()`, `PetscPartitionerCreate()`, `PetscSectionCreate()`, `PetscSectionSetChart()`, `PetscPartitionerPartition()`
734: @*/
735: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
736: {
737: PetscMPIInt size;
738: PetscBool isplex;
739: PetscSection vertSection = NULL;
746: PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex);
748: MPI_Comm_size(PetscObjectComm((PetscObject)part), &size);
749: if (size == 1) {
750: PetscInt *points;
751: PetscInt cStart, cEnd, c;
753: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
754: PetscSectionReset(partSection);
755: PetscSectionSetChart(partSection, 0, size);
756: PetscSectionSetDof(partSection, 0, cEnd - cStart);
757: PetscSectionSetUp(partSection);
758: PetscMalloc1(cEnd - cStart, &points);
759: for (c = cStart; c < cEnd; ++c) points[c] = c;
760: ISCreateGeneral(PetscObjectComm((PetscObject)part), cEnd - cStart, points, PETSC_OWN_POINTER, partition);
761: return 0;
762: }
763: if (part->height == 0) {
764: PetscInt numVertices = 0;
765: PetscInt *start = NULL;
766: PetscInt *adjacency = NULL;
767: IS globalNumbering;
769: if (!part->noGraph || part->viewGraph) {
770: DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
771: } else { /* only compute the number of owned local vertices */
772: const PetscInt *idxs;
773: PetscInt p, pStart, pEnd;
775: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
776: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
777: ISGetIndices(globalNumbering, &idxs);
778: for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
779: ISRestoreIndices(globalNumbering, &idxs);
780: }
781: if (part->usevwgt) {
782: PetscSection section = dm->localSection, clSection = NULL;
783: IS clPoints = NULL;
784: const PetscInt *gid, *clIdx;
785: PetscInt v, p, pStart, pEnd;
787: /* 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) */
788: /* We do this only if the local section has been set */
789: if (section) {
790: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
791: if (!clSection) DMPlexCreateClosureIndex(dm, NULL);
792: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
793: ISGetIndices(clPoints, &clIdx);
794: }
795: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
796: PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
797: PetscSectionSetChart(vertSection, 0, numVertices);
798: if (globalNumbering) {
799: ISGetIndices(globalNumbering, &gid);
800: } else gid = NULL;
801: for (p = pStart, v = 0; p < pEnd; ++p) {
802: PetscInt dof = 1;
804: /* skip cells in the overlap */
805: if (gid && gid[p - pStart] < 0) continue;
807: if (section) {
808: PetscInt cl, clSize, clOff;
810: dof = 0;
811: PetscSectionGetDof(clSection, p, &clSize);
812: PetscSectionGetOffset(clSection, p, &clOff);
813: for (cl = 0; cl < clSize; cl += 2) {
814: PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */
816: PetscSectionGetDof(section, clPoint, &clDof);
817: dof += clDof;
818: }
819: }
821: PetscSectionSetDof(vertSection, v, dof);
822: v++;
823: }
824: if (globalNumbering) ISRestoreIndices(globalNumbering, &gid);
825: if (clPoints) ISRestoreIndices(clPoints, &clIdx);
826: PetscSectionSetUp(vertSection);
827: }
828: PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
829: PetscFree(start);
830: PetscFree(adjacency);
831: if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
832: const PetscInt *globalNum;
833: const PetscInt *partIdx;
834: PetscInt *map, cStart, cEnd;
835: PetscInt *adjusted, i, localSize, offset;
836: IS newPartition;
838: ISGetLocalSize(*partition, &localSize);
839: PetscMalloc1(localSize, &adjusted);
840: ISGetIndices(globalNumbering, &globalNum);
841: ISGetIndices(*partition, &partIdx);
842: PetscMalloc1(localSize, &map);
843: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
844: for (i = cStart, offset = 0; i < cEnd; i++) {
845: if (globalNum[i - cStart] >= 0) map[offset++] = i;
846: }
847: for (i = 0; i < localSize; i++) adjusted[i] = map[partIdx[i]];
848: PetscFree(map);
849: ISRestoreIndices(*partition, &partIdx);
850: ISRestoreIndices(globalNumbering, &globalNum);
851: ISCreateGeneral(PETSC_COMM_SELF, localSize, adjusted, PETSC_OWN_POINTER, &newPartition);
852: ISDestroy(&globalNumbering);
853: ISDestroy(partition);
854: *partition = newPartition;
855: }
856: } else SETERRQ(PetscObjectComm((PetscObject)part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %" PetscInt_FMT " for points to partition", part->height);
857: PetscSectionDestroy(&vertSection);
858: return 0;
859: }
861: /*@
862: DMPlexGetPartitioner - Get the mesh partitioner
864: Not collective
866: Input Parameter:
867: . dm - The DM
869: Output Parameter:
870: . part - The PetscPartitioner
872: Level: developer
874: Note: This gets a borrowed reference, so the user should not destroy this PetscPartitioner.
876: .seealso `DMPlexDistribute()`, `DMPlexSetPartitioner()`, `PetscPartitionerDMPlexPartition()`, `PetscPartitionerCreate()`
877: @*/
878: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
879: {
880: DM_Plex *mesh = (DM_Plex *)dm->data;
884: *part = mesh->partitioner;
885: return 0;
886: }
888: /*@
889: DMPlexSetPartitioner - Set the mesh partitioner
891: logically collective on DM
893: Input Parameters:
894: + dm - The DM
895: - part - The partitioner
897: Level: developer
899: Note: Any existing PetscPartitioner will be destroyed.
901: .seealso `DMPlexDistribute()`, `DMPlexGetPartitioner()`, `PetscPartitionerCreate()`
902: @*/
903: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
904: {
905: DM_Plex *mesh = (DM_Plex *)dm->data;
909: PetscObjectReference((PetscObject)part);
910: PetscPartitionerDestroy(&mesh->partitioner);
911: mesh->partitioner = part;
912: return 0;
913: }
915: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
916: {
917: const PetscInt *cone;
918: PetscInt coneSize, c;
919: PetscBool missing;
922: PetscHSetIQueryAdd(ht, point, &missing);
923: if (missing) {
924: DMPlexGetCone(dm, point, &cone);
925: DMPlexGetConeSize(dm, point, &coneSize);
926: for (c = 0; c < coneSize; c++) DMPlexAddClosure_Private(dm, ht, cone[c]);
927: }
928: return 0;
929: }
931: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
932: {
933: if (up) {
934: PetscInt parent;
936: DMPlexGetTreeParent(dm, point, &parent, NULL);
937: if (parent != point) {
938: PetscInt closureSize, *closure = NULL, i;
940: DMPlexGetTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
941: for (i = 0; i < closureSize; i++) {
942: PetscInt cpoint = closure[2 * i];
944: PetscHSetIAdd(ht, cpoint);
945: DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_TRUE, PETSC_FALSE);
946: }
947: DMPlexRestoreTransitiveClosure(dm, parent, PETSC_TRUE, &closureSize, &closure);
948: }
949: }
950: if (down) {
951: PetscInt numChildren;
952: const PetscInt *children;
954: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
955: if (numChildren) {
956: PetscInt i;
958: for (i = 0; i < numChildren; i++) {
959: PetscInt cpoint = children[i];
961: PetscHSetIAdd(ht, cpoint);
962: DMPlexAddClosure_Tree(dm, ht, cpoint, PETSC_FALSE, PETSC_TRUE);
963: }
964: }
965: }
966: return 0;
967: }
969: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
970: {
971: PetscInt parent;
974: DMPlexGetTreeParent(dm, point, &parent, NULL);
975: if (point != parent) {
976: const PetscInt *cone;
977: PetscInt coneSize, c;
979: DMPlexAddClosureTree_Up_Private(dm, ht, parent);
980: DMPlexAddClosure_Private(dm, ht, parent);
981: DMPlexGetCone(dm, parent, &cone);
982: DMPlexGetConeSize(dm, parent, &coneSize);
983: for (c = 0; c < coneSize; c++) {
984: const PetscInt cp = cone[c];
986: DMPlexAddClosureTree_Up_Private(dm, ht, cp);
987: }
988: }
989: return 0;
990: }
992: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
993: {
994: PetscInt i, numChildren;
995: const PetscInt *children;
998: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
999: for (i = 0; i < numChildren; i++) PetscHSetIAdd(ht, children[i]);
1000: return 0;
1001: }
1003: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
1004: {
1005: const PetscInt *cone;
1006: PetscInt coneSize, c;
1009: PetscHSetIAdd(ht, point);
1010: DMPlexAddClosureTree_Up_Private(dm, ht, point);
1011: DMPlexAddClosureTree_Down_Private(dm, ht, point);
1012: DMPlexGetCone(dm, point, &cone);
1013: DMPlexGetConeSize(dm, point, &coneSize);
1014: for (c = 0; c < coneSize; c++) DMPlexAddClosureTree_Private(dm, ht, cone[c]);
1015: return 0;
1016: }
1018: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
1019: {
1020: DM_Plex *mesh = (DM_Plex *)dm->data;
1021: const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
1022: PetscInt nelems, *elems, off = 0, p;
1023: PetscHSetI ht = NULL;
1025: PetscHSetICreate(&ht);
1026: PetscHSetIResize(ht, numPoints * 16);
1027: if (!hasTree) {
1028: for (p = 0; p < numPoints; ++p) DMPlexAddClosure_Private(dm, ht, points[p]);
1029: } else {
1030: #if 1
1031: for (p = 0; p < numPoints; ++p) DMPlexAddClosureTree_Private(dm, ht, points[p]);
1032: #else
1033: PetscInt *closure = NULL, closureSize, c;
1034: for (p = 0; p < numPoints; ++p) {
1035: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
1036: for (c = 0; c < closureSize * 2; c += 2) {
1037: PetscHSetIAdd(ht, closure[c]);
1038: if (hasTree) DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);
1039: }
1040: }
1041: if (closure) DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);
1042: #endif
1043: }
1044: PetscHSetIGetSize(ht, &nelems);
1045: PetscMalloc1(nelems, &elems);
1046: PetscHSetIGetElems(ht, &off, elems);
1047: PetscHSetIDestroy(&ht);
1048: PetscSortInt(nelems, elems);
1049: ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
1050: return 0;
1051: }
1053: /*@
1054: DMPlexPartitionLabelClosure - Add the closure of all points to the partition label
1056: Input Parameters:
1057: + dm - The DM
1058: - label - DMLabel assigning ranks to remote roots
1060: Level: developer
1062: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1063: @*/
1064: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
1065: {
1066: IS rankIS, pointIS, closureIS;
1067: const PetscInt *ranks, *points;
1068: PetscInt numRanks, numPoints, r;
1070: DMLabelGetValueIS(label, &rankIS);
1071: ISGetLocalSize(rankIS, &numRanks);
1072: ISGetIndices(rankIS, &ranks);
1073: for (r = 0; r < numRanks; ++r) {
1074: const PetscInt rank = ranks[r];
1075: DMLabelGetStratumIS(label, rank, &pointIS);
1076: ISGetLocalSize(pointIS, &numPoints);
1077: ISGetIndices(pointIS, &points);
1078: DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
1079: ISRestoreIndices(pointIS, &points);
1080: ISDestroy(&pointIS);
1081: DMLabelSetStratumIS(label, rank, closureIS);
1082: ISDestroy(&closureIS);
1083: }
1084: ISRestoreIndices(rankIS, &ranks);
1085: ISDestroy(&rankIS);
1086: return 0;
1087: }
1089: /*@
1090: DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label
1092: Input Parameters:
1093: + dm - The DM
1094: - label - DMLabel assigning ranks to remote roots
1096: Level: developer
1098: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1099: @*/
1100: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1101: {
1102: IS rankIS, pointIS;
1103: const PetscInt *ranks, *points;
1104: PetscInt numRanks, numPoints, r, p, a, adjSize;
1105: PetscInt *adj = NULL;
1107: DMLabelGetValueIS(label, &rankIS);
1108: ISGetLocalSize(rankIS, &numRanks);
1109: ISGetIndices(rankIS, &ranks);
1110: for (r = 0; r < numRanks; ++r) {
1111: const PetscInt rank = ranks[r];
1113: DMLabelGetStratumIS(label, rank, &pointIS);
1114: ISGetLocalSize(pointIS, &numPoints);
1115: ISGetIndices(pointIS, &points);
1116: for (p = 0; p < numPoints; ++p) {
1117: adjSize = PETSC_DETERMINE;
1118: DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
1119: for (a = 0; a < adjSize; ++a) DMLabelSetValue(label, adj[a], rank);
1120: }
1121: ISRestoreIndices(pointIS, &points);
1122: ISDestroy(&pointIS);
1123: }
1124: ISRestoreIndices(rankIS, &ranks);
1125: ISDestroy(&rankIS);
1126: PetscFree(adj);
1127: return 0;
1128: }
1130: /*@
1131: DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point SF
1133: Input Parameters:
1134: + dm - The DM
1135: - label - DMLabel assigning ranks to remote roots
1137: Level: developer
1139: Note: This is required when generating multi-level overlaps to capture
1140: overlap points from non-neighbouring partitions.
1142: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1143: @*/
1144: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1145: {
1146: MPI_Comm comm;
1147: PetscMPIInt rank;
1148: PetscSF sfPoint;
1149: DMLabel lblRoots, lblLeaves;
1150: IS rankIS, pointIS;
1151: const PetscInt *ranks;
1152: PetscInt numRanks, r;
1154: PetscObjectGetComm((PetscObject)dm, &comm);
1155: MPI_Comm_rank(comm, &rank);
1156: DMGetPointSF(dm, &sfPoint);
1157: /* Pull point contributions from remote leaves into local roots */
1158: DMLabelGather(label, sfPoint, &lblLeaves);
1159: DMLabelGetValueIS(lblLeaves, &rankIS);
1160: ISGetLocalSize(rankIS, &numRanks);
1161: ISGetIndices(rankIS, &ranks);
1162: for (r = 0; r < numRanks; ++r) {
1163: const PetscInt remoteRank = ranks[r];
1164: if (remoteRank == rank) continue;
1165: DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
1166: DMLabelInsertIS(label, pointIS, remoteRank);
1167: ISDestroy(&pointIS);
1168: }
1169: ISRestoreIndices(rankIS, &ranks);
1170: ISDestroy(&rankIS);
1171: DMLabelDestroy(&lblLeaves);
1172: /* Push point contributions from roots into remote leaves */
1173: DMLabelDistribute(label, sfPoint, &lblRoots);
1174: DMLabelGetValueIS(lblRoots, &rankIS);
1175: ISGetLocalSize(rankIS, &numRanks);
1176: ISGetIndices(rankIS, &ranks);
1177: for (r = 0; r < numRanks; ++r) {
1178: const PetscInt remoteRank = ranks[r];
1179: if (remoteRank == rank) continue;
1180: DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
1181: DMLabelInsertIS(label, pointIS, remoteRank);
1182: ISDestroy(&pointIS);
1183: }
1184: ISRestoreIndices(rankIS, &ranks);
1185: ISDestroy(&rankIS);
1186: DMLabelDestroy(&lblRoots);
1187: return 0;
1188: }
1190: /*@
1191: DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label
1193: Input Parameters:
1194: + dm - The DM
1195: . rootLabel - DMLabel assigning ranks to local roots
1196: - processSF - A star forest mapping into the local index on each remote rank
1198: Output Parameter:
1199: . leafLabel - DMLabel assigning ranks to remote roots
1201: Note: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
1202: resulting leafLabel is a receiver mapping of remote roots to their parent rank.
1204: Level: developer
1206: .seealso: `DMPlexPartitionLabelCreateSF()`, `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1207: @*/
1208: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1209: {
1210: MPI_Comm comm;
1211: PetscMPIInt rank, size, r;
1212: PetscInt p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1213: PetscSF sfPoint;
1214: PetscSection rootSection;
1215: PetscSFNode *rootPoints, *leafPoints;
1216: const PetscSFNode *remote;
1217: const PetscInt *local, *neighbors;
1218: IS valueIS;
1219: PetscBool mpiOverflow = PETSC_FALSE;
1221: PetscLogEventBegin(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1222: PetscObjectGetComm((PetscObject)dm, &comm);
1223: MPI_Comm_rank(comm, &rank);
1224: MPI_Comm_size(comm, &size);
1225: DMGetPointSF(dm, &sfPoint);
1227: /* Convert to (point, rank) and use actual owners */
1228: PetscSectionCreate(comm, &rootSection);
1229: PetscSectionSetChart(rootSection, 0, size);
1230: DMLabelGetValueIS(rootLabel, &valueIS);
1231: ISGetLocalSize(valueIS, &numNeighbors);
1232: ISGetIndices(valueIS, &neighbors);
1233: for (n = 0; n < numNeighbors; ++n) {
1234: DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
1235: PetscSectionAddDof(rootSection, neighbors[n], numPoints);
1236: }
1237: PetscSectionSetUp(rootSection);
1238: PetscSectionGetStorageSize(rootSection, &rootSize);
1239: PetscMalloc1(rootSize, &rootPoints);
1240: PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
1241: for (n = 0; n < numNeighbors; ++n) {
1242: IS pointIS;
1243: const PetscInt *points;
1245: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1246: DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
1247: ISGetLocalSize(pointIS, &numPoints);
1248: ISGetIndices(pointIS, &points);
1249: for (p = 0; p < numPoints; ++p) {
1250: if (local) PetscFindInt(points[p], nleaves, local, &l);
1251: else l = -1;
1252: if (l >= 0) {
1253: rootPoints[off + p] = remote[l];
1254: } else {
1255: rootPoints[off + p].index = points[p];
1256: rootPoints[off + p].rank = rank;
1257: }
1258: }
1259: ISRestoreIndices(pointIS, &points);
1260: ISDestroy(&pointIS);
1261: }
1263: /* Try to communicate overlap using All-to-All */
1264: if (!processSF) {
1265: PetscInt64 counter = 0;
1266: PetscBool locOverflow = PETSC_FALSE;
1267: PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;
1269: PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
1270: for (n = 0; n < numNeighbors; ++n) {
1271: PetscSectionGetDof(rootSection, neighbors[n], &dof);
1272: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1273: #if defined(PETSC_USE_64BIT_INDICES)
1274: if (dof > PETSC_MPI_INT_MAX) {
1275: locOverflow = PETSC_TRUE;
1276: break;
1277: }
1278: if (off > PETSC_MPI_INT_MAX) {
1279: locOverflow = PETSC_TRUE;
1280: break;
1281: }
1282: #endif
1283: scounts[neighbors[n]] = (PetscMPIInt)dof;
1284: sdispls[neighbors[n]] = (PetscMPIInt)off;
1285: }
1286: MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
1287: for (r = 0; r < size; ++r) {
1288: rdispls[r] = (int)counter;
1289: counter += rcounts[r];
1290: }
1291: if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
1292: MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
1293: if (!mpiOverflow) {
1294: PetscInfo(dm, "Using Alltoallv for mesh distribution\n");
1295: leafSize = (PetscInt)counter;
1296: PetscMalloc1(leafSize, &leafPoints);
1297: MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
1298: }
1299: PetscFree4(scounts, sdispls, rcounts, rdispls);
1300: }
1302: /* Communicate overlap using process star forest */
1303: if (processSF || mpiOverflow) {
1304: PetscSF procSF;
1305: PetscSection leafSection;
1307: if (processSF) {
1308: PetscInfo(dm, "Using processSF for mesh distribution\n");
1309: PetscObjectReference((PetscObject)processSF);
1310: procSF = processSF;
1311: } else {
1312: PetscInfo(dm, "Using processSF for mesh distribution (MPI overflow)\n");
1313: PetscSFCreate(comm, &procSF);
1314: PetscSFSetGraphWithPattern(procSF, NULL, PETSCSF_PATTERN_ALLTOALL);
1315: }
1317: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
1318: DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void **)&leafPoints);
1319: PetscSectionGetStorageSize(leafSection, &leafSize);
1320: PetscSectionDestroy(&leafSection);
1321: PetscSFDestroy(&procSF);
1322: }
1324: for (p = 0; p < leafSize; p++) DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
1326: ISRestoreIndices(valueIS, &neighbors);
1327: ISDestroy(&valueIS);
1328: PetscSectionDestroy(&rootSection);
1329: PetscFree(rootPoints);
1330: PetscFree(leafPoints);
1331: PetscLogEventEnd(DMPLEX_PartLabelInvert, dm, 0, 0, 0);
1332: return 0;
1333: }
1335: /*@
1336: DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points
1338: Input Parameters:
1339: + dm - The DM
1340: - label - DMLabel assigning ranks to remote roots
1342: Output Parameter:
1343: . sf - The star forest communication context encapsulating the defined mapping
1345: Note: The incoming label is a receiver mapping of remote points to their parent rank.
1347: Level: developer
1349: .seealso: `DMPlexDistribute()`, `DMPlexCreateOverlap()`
1350: @*/
1351: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
1352: {
1353: PetscMPIInt rank;
1354: PetscInt n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1355: PetscSFNode *remotePoints;
1356: IS remoteRootIS, neighborsIS;
1357: const PetscInt *remoteRoots, *neighbors;
1359: PetscLogEventBegin(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1360: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1362: DMLabelGetValueIS(label, &neighborsIS);
1363: #if 0
1364: {
1365: IS is;
1366: ISDuplicate(neighborsIS, &is);
1367: ISSort(is);
1368: ISDestroy(&neighborsIS);
1369: neighborsIS = is;
1370: }
1371: #endif
1372: ISGetLocalSize(neighborsIS, &nNeighbors);
1373: ISGetIndices(neighborsIS, &neighbors);
1374: for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1375: DMLabelGetStratumSize(label, neighbors[n], &numPoints);
1376: numRemote += numPoints;
1377: }
1378: PetscMalloc1(numRemote, &remotePoints);
1379: /* Put owned points first */
1380: DMLabelGetStratumSize(label, rank, &numPoints);
1381: if (numPoints > 0) {
1382: DMLabelGetStratumIS(label, rank, &remoteRootIS);
1383: ISGetIndices(remoteRootIS, &remoteRoots);
1384: for (p = 0; p < numPoints; p++) {
1385: remotePoints[idx].index = remoteRoots[p];
1386: remotePoints[idx].rank = rank;
1387: idx++;
1388: }
1389: ISRestoreIndices(remoteRootIS, &remoteRoots);
1390: ISDestroy(&remoteRootIS);
1391: }
1392: /* Now add remote points */
1393: for (n = 0; n < nNeighbors; ++n) {
1394: const PetscInt nn = neighbors[n];
1396: DMLabelGetStratumSize(label, nn, &numPoints);
1397: if (nn == rank || numPoints <= 0) continue;
1398: DMLabelGetStratumIS(label, nn, &remoteRootIS);
1399: ISGetIndices(remoteRootIS, &remoteRoots);
1400: for (p = 0; p < numPoints; p++) {
1401: remotePoints[idx].index = remoteRoots[p];
1402: remotePoints[idx].rank = nn;
1403: idx++;
1404: }
1405: ISRestoreIndices(remoteRootIS, &remoteRoots);
1406: ISDestroy(&remoteRootIS);
1407: }
1408: PetscSFCreate(PetscObjectComm((PetscObject)dm), sf);
1409: DMPlexGetChart(dm, &pStart, &pEnd);
1410: PetscSFSetGraph(*sf, pEnd - pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1411: PetscSFSetUp(*sf);
1412: ISDestroy(&neighborsIS);
1413: PetscLogEventEnd(DMPLEX_PartLabelCreateSF, dm, 0, 0, 0);
1414: return 0;
1415: }
1417: #if defined(PETSC_HAVE_PARMETIS)
1418: #include <parmetis.h>
1419: #endif
1421: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
1422: * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
1423: * them out in that case. */
1424: #if defined(PETSC_HAVE_PARMETIS)
1425: /*@C
1427: DMPlexRewriteSF - Rewrites the ownership of the SF of a DM (in place).
1429: Input parameters:
1430: + dm - The DMPlex object.
1431: . n - The number of points.
1432: . pointsToRewrite - The points in the SF whose ownership will change.
1433: . targetOwners - New owner for each element in pointsToRewrite.
1434: - degrees - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.
1436: Level: developer
1438: @*/
1439: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1440: {
1441: PetscInt pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1442: PetscInt *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1443: PetscSFNode *leafLocationsNew;
1444: const PetscSFNode *iremote;
1445: const PetscInt *ilocal;
1446: PetscBool *isLeaf;
1447: PetscSF sf;
1448: MPI_Comm comm;
1449: PetscMPIInt rank, size;
1451: PetscObjectGetComm((PetscObject)dm, &comm);
1452: MPI_Comm_rank(comm, &rank);
1453: MPI_Comm_size(comm, &size);
1454: DMPlexGetChart(dm, &pStart, &pEnd);
1456: DMGetPointSF(dm, &sf);
1457: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1458: PetscMalloc1(pEnd - pStart, &isLeaf);
1459: for (i = 0; i < pEnd - pStart; i++) isLeaf[i] = PETSC_FALSE;
1460: for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;
1462: PetscMalloc1(pEnd - pStart + 1, &cumSumDegrees);
1463: cumSumDegrees[0] = 0;
1464: for (i = 1; i <= pEnd - pStart; i++) cumSumDegrees[i] = cumSumDegrees[i - 1] + degrees[i - 1];
1465: sumDegrees = cumSumDegrees[pEnd - pStart];
1466: /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */
1468: PetscMalloc1(sumDegrees, &locationsOfLeafs);
1469: PetscMalloc1(pEnd - pStart, &rankOnLeafs);
1470: for (i = 0; i < pEnd - pStart; i++) rankOnLeafs[i] = rank;
1471: PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1472: PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1473: PetscFree(rankOnLeafs);
1475: /* get the remote local points of my leaves */
1476: PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
1477: PetscMalloc1(pEnd - pStart, &points);
1478: for (i = 0; i < pEnd - pStart; i++) points[i] = pStart + i;
1479: PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1480: PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1481: PetscFree(points);
1482: /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1483: PetscMalloc1(pEnd - pStart, &newOwners);
1484: PetscMalloc1(pEnd - pStart, &newNumbers);
1485: for (i = 0; i < pEnd - pStart; i++) {
1486: newOwners[i] = -1;
1487: newNumbers[i] = -1;
1488: }
1489: {
1490: PetscInt oldNumber, newNumber, oldOwner, newOwner;
1491: for (i = 0; i < n; i++) {
1492: oldNumber = pointsToRewrite[i];
1493: newNumber = -1;
1494: oldOwner = rank;
1495: newOwner = targetOwners[i];
1496: if (oldOwner == newOwner) {
1497: newNumber = oldNumber;
1498: } else {
1499: for (j = 0; j < degrees[oldNumber]; j++) {
1500: if (locationsOfLeafs[cumSumDegrees[oldNumber] + j] == newOwner) {
1501: newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber] + j];
1502: break;
1503: }
1504: }
1505: }
1508: newOwners[oldNumber] = newOwner;
1509: newNumbers[oldNumber] = newNumber;
1510: }
1511: }
1512: PetscFree(cumSumDegrees);
1513: PetscFree(locationsOfLeafs);
1514: PetscFree(remoteLocalPointOfLeafs);
1516: PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1517: PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners, MPI_REPLACE);
1518: PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);
1519: PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers, MPI_REPLACE);
1521: /* Now count how many leafs we have on each processor. */
1522: leafCounter = 0;
1523: for (i = 0; i < pEnd - pStart; i++) {
1524: if (newOwners[i] >= 0) {
1525: if (newOwners[i] != rank) leafCounter++;
1526: } else {
1527: if (isLeaf[i]) leafCounter++;
1528: }
1529: }
1531: /* Now set up the new sf by creating the leaf arrays */
1532: PetscMalloc1(leafCounter, &leafsNew);
1533: PetscMalloc1(leafCounter, &leafLocationsNew);
1535: leafCounter = 0;
1536: counter = 0;
1537: for (i = 0; i < pEnd - pStart; i++) {
1538: if (newOwners[i] >= 0) {
1539: if (newOwners[i] != rank) {
1540: leafsNew[leafCounter] = i;
1541: leafLocationsNew[leafCounter].rank = newOwners[i];
1542: leafLocationsNew[leafCounter].index = newNumbers[i];
1543: leafCounter++;
1544: }
1545: } else {
1546: if (isLeaf[i]) {
1547: leafsNew[leafCounter] = i;
1548: leafLocationsNew[leafCounter].rank = iremote[counter].rank;
1549: leafLocationsNew[leafCounter].index = iremote[counter].index;
1550: leafCounter++;
1551: }
1552: }
1553: if (isLeaf[i]) counter++;
1554: }
1556: PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
1557: PetscFree(newOwners);
1558: PetscFree(newNumbers);
1559: PetscFree(isLeaf);
1560: return 0;
1561: }
1563: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1564: {
1565: PetscInt *distribution, min, max, sum;
1566: PetscMPIInt rank, size;
1568: MPI_Comm_size(comm, &size);
1569: MPI_Comm_rank(comm, &rank);
1570: PetscCalloc1(size, &distribution);
1571: for (PetscInt i = 0; i < n; i++) {
1572: if (part) distribution[part[i]] += vtxwgt[skip * i];
1573: else distribution[rank] += vtxwgt[skip * i];
1574: }
1575: MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
1576: min = distribution[0];
1577: max = distribution[0];
1578: sum = distribution[0];
1579: for (PetscInt i = 1; i < size; i++) {
1580: if (distribution[i] < min) min = distribution[i];
1581: if (distribution[i] > max) max = distribution[i];
1582: sum += distribution[i];
1583: }
1584: PetscViewerASCIIPrintf(viewer, "Min: %" PetscInt_FMT ", Avg: %" PetscInt_FMT ", Max: %" PetscInt_FMT ", Balance: %f\n", min, sum / size, max, (max * 1. * size) / sum);
1585: PetscFree(distribution);
1586: return 0;
1587: }
1589: #endif
1591: /*@
1592: DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the PointSF of the DM inplace.
1594: Input parameters:
1595: + dm - The DMPlex object.
1596: . entityDepth - depth of the entity to balance (0 -> balance vertices).
1597: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
1598: - parallel - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.
1600: Output parameters:
1601: . success - whether the graph partitioning was successful or not, optional. Unsuccessful simply means no change to the partitioning
1603: Options Database:
1604: + -dm_plex_rebalance_shared_points_parmetis - Use ParMetis instead of Metis for the partitioner
1605: . -dm_plex_rebalance_shared_points_use_initial_guess - Use current partition to bootstrap ParMetis partition
1606: . -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_
1607: - -dm_plex_rebalance_shared_points_monitor - Monitor the shared points rebalance process
1609: Developer Notes:
1610: This should use MatPartitioning to allow the use of any partitioner and not be hardwired to use ParMetis
1612: Level: intermediate
1613: @*/
1614: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1615: {
1616: #if defined(PETSC_HAVE_PARMETIS)
1617: PetscSF sf;
1618: PetscInt ierr, i, j, idx, jdx;
1619: PetscInt eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1620: const PetscInt *degrees, *ilocal;
1621: const PetscSFNode *iremote;
1622: PetscBool *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1623: PetscInt numExclusivelyOwned, numNonExclusivelyOwned;
1624: PetscMPIInt rank, size;
1625: PetscInt *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1626: const PetscInt *cumSumVertices;
1627: PetscInt offset, counter;
1628: PetscInt *vtxwgt;
1629: const PetscInt *xadj, *adjncy;
1630: PetscInt *part, *options;
1631: PetscInt nparts, wgtflag, numflag, ncon, edgecut;
1632: real_t *ubvec;
1633: PetscInt *firstVertices, *renumbering;
1634: PetscInt failed, failedGlobal;
1635: MPI_Comm comm;
1636: Mat A;
1637: PetscViewer viewer;
1638: PetscViewerFormat format;
1639: PetscLayout layout;
1640: real_t *tpwgts;
1641: PetscMPIInt *counts, *mpiCumSumVertices;
1642: PetscInt *pointsToRewrite;
1643: PetscInt numRows;
1644: PetscBool done, usematpartitioning = PETSC_FALSE;
1645: IS ispart = NULL;
1646: MatPartitioning mp;
1647: const char *prefix;
1649: PetscObjectGetComm((PetscObject)dm, &comm);
1650: MPI_Comm_size(comm, &size);
1651: if (size == 1) {
1652: if (success) *success = PETSC_TRUE;
1653: return 0;
1654: }
1655: if (success) *success = PETSC_FALSE;
1656: MPI_Comm_rank(comm, &rank);
1658: parallel = PETSC_FALSE;
1659: useInitialGuess = PETSC_FALSE;
1660: PetscObjectOptionsBegin((PetscObject)dm);
1661: PetscOptionsName("-dm_plex_rebalance_shared_points_parmetis", "Use ParMetis instead of Metis for the partitioner", "DMPlexRebalanceSharedPoints", ¶llel);
1662: PetscOptionsBool("-dm_plex_rebalance_shared_points_use_initial_guess", "Use current partition to bootstrap ParMetis partition", "DMPlexRebalanceSharedPoints", useInitialGuess, &useInitialGuess, NULL);
1663: PetscOptionsBool("-dm_plex_rebalance_shared_points_use_mat_partitioning", "Use the MatPartitioning object to partition", "DMPlexRebalanceSharedPoints", usematpartitioning, &usematpartitioning, NULL);
1664: PetscOptionsViewer("-dm_plex_rebalance_shared_points_monitor", "Monitor the shared points rebalance process", "DMPlexRebalanceSharedPoints", &viewer, &format, NULL);
1665: PetscOptionsEnd();
1666: if (viewer) PetscViewerPushFormat(viewer, format);
1668: PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1670: DMGetOptionsPrefix(dm, &prefix);
1671: PetscOptionsGetViewer(comm, ((PetscObject)dm)->options, prefix, "-dm_rebalance_partition_view", &viewer, &format, NULL);
1672: if (viewer) PetscViewerPushFormat(viewer, format);
1674: /* Figure out all points in the plex that we are interested in balancing. */
1675: DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
1676: DMPlexGetChart(dm, &pStart, &pEnd);
1677: PetscMalloc1(pEnd - pStart, &toBalance);
1678: for (i = 0; i < pEnd - pStart; i++) toBalance[i] = (PetscBool)(i >= eBegin && i < eEnd);
1680: /* There are three types of points:
1681: * exclusivelyOwned: points that are owned by this process and only seen by this process
1682: * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1683: * leaf: a point that is seen by this process but owned by a different process
1684: */
1685: DMGetPointSF(dm, &sf);
1686: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1687: PetscMalloc1(pEnd - pStart, &isLeaf);
1688: PetscMalloc1(pEnd - pStart, &isNonExclusivelyOwned);
1689: PetscMalloc1(pEnd - pStart, &isExclusivelyOwned);
1690: for (i = 0; i < pEnd - pStart; i++) {
1691: isNonExclusivelyOwned[i] = PETSC_FALSE;
1692: isExclusivelyOwned[i] = PETSC_FALSE;
1693: isLeaf[i] = PETSC_FALSE;
1694: }
1696: /* mark all the leafs */
1697: for (i = 0; i < nleafs; i++) isLeaf[ilocal[i] - pStart] = PETSC_TRUE;
1699: /* for an owned point, we can figure out whether another processor sees it or
1700: * not by calculating its degree */
1701: PetscSFComputeDegreeBegin(sf, °rees);
1702: PetscSFComputeDegreeEnd(sf, °rees);
1703: numExclusivelyOwned = 0;
1704: numNonExclusivelyOwned = 0;
1705: for (i = 0; i < pEnd - pStart; i++) {
1706: if (toBalance[i]) {
1707: if (degrees[i] > 0) {
1708: isNonExclusivelyOwned[i] = PETSC_TRUE;
1709: numNonExclusivelyOwned += 1;
1710: } else {
1711: if (!isLeaf[i]) {
1712: isExclusivelyOwned[i] = PETSC_TRUE;
1713: numExclusivelyOwned += 1;
1714: }
1715: }
1716: }
1717: }
1719: /* Build a graph with one vertex per core representing the
1720: * exclusively owned points and then one vertex per nonExclusively owned
1721: * point. */
1722: PetscLayoutCreate(comm, &layout);
1723: PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
1724: PetscLayoutSetUp(layout);
1725: PetscLayoutGetRanges(layout, &cumSumVertices);
1726: PetscMalloc1(pEnd - pStart, &globalNumbersOfLocalOwnedVertices);
1727: for (i = 0; i < pEnd - pStart; i++) globalNumbersOfLocalOwnedVertices[i] = pStart - 1;
1728: offset = cumSumVertices[rank];
1729: counter = 0;
1730: for (i = 0; i < pEnd - pStart; i++) {
1731: if (toBalance[i]) {
1732: if (degrees[i] > 0) {
1733: globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1734: counter++;
1735: }
1736: }
1737: }
1739: /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
1740: PetscMalloc1(pEnd - pStart, &leafGlobalNumbers);
1741: PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE);
1742: PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers, MPI_REPLACE);
1744: /* Build the graph for partitioning */
1745: numRows = 1 + numNonExclusivelyOwned;
1746: PetscLogEventBegin(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);
1747: MatCreate(comm, &A);
1748: MatSetType(A, MATMPIADJ);
1749: MatSetSizes(A, numRows, numRows, cumSumVertices[size], cumSumVertices[size]);
1750: idx = cumSumVertices[rank];
1751: for (i = 0; i < pEnd - pStart; i++) {
1752: if (toBalance[i]) {
1753: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1754: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1755: else continue;
1756: MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
1757: MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
1758: }
1759: }
1760: PetscFree(globalNumbersOfLocalOwnedVertices);
1761: PetscFree(leafGlobalNumbers);
1762: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1763: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1764: PetscLogEventEnd(DMPLEX_RebalBuildGraph, dm, 0, 0, 0);
1766: nparts = size;
1767: ncon = 1;
1768: PetscMalloc1(ncon * nparts, &tpwgts);
1769: for (i = 0; i < ncon * nparts; i++) tpwgts[i] = 1. / (nparts);
1770: PetscMalloc1(ncon, &ubvec);
1771: ubvec[0] = 1.05;
1772: ubvec[1] = 1.05;
1774: PetscMalloc1(ncon * (1 + numNonExclusivelyOwned), &vtxwgt);
1775: if (ncon == 2) {
1776: vtxwgt[0] = numExclusivelyOwned;
1777: vtxwgt[1] = 1;
1778: for (i = 0; i < numNonExclusivelyOwned; i++) {
1779: vtxwgt[ncon * (i + 1)] = 1;
1780: vtxwgt[ncon * (i + 1) + 1] = 0;
1781: }
1782: } else {
1783: PetscInt base, ms;
1784: MPI_Allreduce(&numExclusivelyOwned, &base, 1, MPIU_INT, MPIU_MAX, PetscObjectComm((PetscObject)dm));
1785: MatGetSize(A, &ms, NULL);
1786: ms -= size;
1787: base = PetscMax(base, ms);
1788: vtxwgt[0] = base + numExclusivelyOwned;
1789: for (i = 0; i < numNonExclusivelyOwned; i++) vtxwgt[i + 1] = 1;
1790: }
1792: if (viewer) {
1793: PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %" PetscInt_FMT " on interface of mesh distribution.\n", entityDepth);
1794: PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %" PetscInt_FMT "\n", cumSumVertices[size]);
1795: }
1796: /* TODO: Drop the parallel/sequential choice here and just use MatPartioner for much more flexibility */
1797: if (usematpartitioning) {
1798: const char *prefix;
1800: MatPartitioningCreate(PetscObjectComm((PetscObject)dm), &mp);
1801: PetscObjectSetOptionsPrefix((PetscObject)mp, "dm_plex_rebalance_shared_points_");
1802: PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
1803: PetscObjectPrependOptionsPrefix((PetscObject)mp, prefix);
1804: MatPartitioningSetAdjacency(mp, A);
1805: MatPartitioningSetNumberVertexWeights(mp, ncon);
1806: MatPartitioningSetVertexWeights(mp, vtxwgt);
1807: MatPartitioningSetFromOptions(mp);
1808: MatPartitioningApply(mp, &ispart);
1809: ISGetIndices(ispart, (const PetscInt **)&part);
1810: } else if (parallel) {
1811: if (viewer) PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n");
1812: PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1813: MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1815: PetscMalloc1(4, &options);
1816: options[0] = 1;
1817: options[1] = 0; /* Verbosity */
1818: if (viewer) options[1] = 1;
1819: options[2] = 0; /* Seed */
1820: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1821: wgtflag = 2;
1822: numflag = 0;
1823: if (useInitialGuess) {
1824: PetscViewerASCIIPrintf(viewer, "THIS DOES NOT WORK! I don't know why. Using current distribution of points as initial guess.\n");
1825: for (i = 0; i < numRows; i++) part[i] = rank;
1826: if (viewer) PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n");
1827: PetscStackPushExternal("ParMETIS_V3_RefineKway");
1828: PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1829: ParMETIS_V3_RefineKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1830: PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1831: PetscStackPop;
1833: } else {
1834: PetscStackPushExternal("ParMETIS_V3_PartKway");
1835: PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1836: ParMETIS_V3_PartKway((PetscInt *)cumSumVertices, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt, NULL, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1837: PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1838: PetscStackPop;
1840: }
1841: MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, &xadj, &adjncy, &done);
1842: PetscFree(options);
1843: } else {
1844: if (viewer) PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n");
1845: Mat As;
1846: PetscInt *partGlobal;
1847: PetscInt *numExclusivelyOwnedAll;
1849: PetscMalloc1(cumSumVertices[rank + 1] - cumSumVertices[rank], &part);
1850: MatGetSize(A, &numRows, NULL);
1851: PetscLogEventBegin(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);
1852: MatMPIAdjToSeqRankZero(A, &As);
1853: PetscLogEventEnd(DMPLEX_RebalGatherGraph, dm, 0, 0, 0);
1855: PetscMalloc1(size, &numExclusivelyOwnedAll);
1856: numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1857: MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, numExclusivelyOwnedAll, 1, MPIU_INT, comm);
1859: PetscMalloc1(numRows, &partGlobal);
1860: PetscLogEventBegin(DMPLEX_RebalPartition, 0, 0, 0, 0);
1861: if (rank == 0) {
1862: PetscInt *vtxwgt_g, numRows_g;
1864: MatGetRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1865: PetscMalloc1(2 * numRows_g, &vtxwgt_g);
1866: for (i = 0; i < size; i++) {
1867: vtxwgt_g[ncon * cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1868: if (ncon > 1) vtxwgt_g[ncon * cumSumVertices[i] + 1] = 1;
1869: for (j = cumSumVertices[i] + 1; j < cumSumVertices[i + 1]; j++) {
1870: vtxwgt_g[ncon * j] = 1;
1871: if (ncon > 1) vtxwgt_g[2 * j + 1] = 0;
1872: }
1873: }
1875: PetscMalloc1(64, &options);
1876: METIS_SetDefaultOptions(options); /* initialize all defaults */
1878: options[METIS_OPTION_CONTIG] = 1;
1879: PetscStackPushExternal("METIS_PartGraphKway");
1880: METIS_PartGraphKway(&numRows_g, &ncon, (idx_t *)xadj, (idx_t *)adjncy, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1881: PetscStackPop;
1883: PetscFree(options);
1884: PetscFree(vtxwgt_g);
1885: MatRestoreRowIJ(As, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows_g, &xadj, &adjncy, &done);
1886: MatDestroy(&As);
1887: }
1888: PetscBarrier((PetscObject)dm);
1889: PetscLogEventEnd(DMPLEX_RebalPartition, 0, 0, 0, 0);
1890: PetscFree(numExclusivelyOwnedAll);
1892: /* scatter the partitioning information to ranks */
1893: PetscLogEventBegin(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1894: PetscMalloc1(size, &counts);
1895: PetscMalloc1(size + 1, &mpiCumSumVertices);
1896: for (i = 0; i < size; i++) PetscMPIIntCast(cumSumVertices[i + 1] - cumSumVertices[i], &(counts[i]));
1897: for (i = 0; i <= size; i++) PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
1898: MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
1899: PetscFree(counts);
1900: PetscFree(mpiCumSumVertices);
1901: PetscFree(partGlobal);
1902: PetscLogEventEnd(DMPLEX_RebalScatterPart, 0, 0, 0, 0);
1903: }
1904: PetscFree(ubvec);
1905: PetscFree(tpwgts);
1907: /* Rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1908: PetscMalloc2(size, &firstVertices, size, &renumbering);
1909: firstVertices[rank] = part[0];
1910: MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, firstVertices, 1, MPIU_INT, comm);
1911: for (i = 0; i < size; i++) renumbering[firstVertices[i]] = i;
1912: for (i = 0; i < cumSumVertices[rank + 1] - cumSumVertices[rank]; i++) part[i] = renumbering[part[i]];
1913: PetscFree2(firstVertices, renumbering);
1915: /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1916: failed = (PetscInt)(part[0] != rank);
1917: MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
1918: if (failedGlobal > 0) {
1920: PetscFree(vtxwgt);
1921: PetscFree(toBalance);
1922: PetscFree(isLeaf);
1923: PetscFree(isNonExclusivelyOwned);
1924: PetscFree(isExclusivelyOwned);
1925: if (usematpartitioning) {
1926: ISRestoreIndices(ispart, (const PetscInt **)&part);
1927: ISDestroy(&ispart);
1928: } else PetscFree(part);
1929: if (viewer) {
1930: PetscViewerPopFormat(viewer);
1931: PetscViewerDestroy(&viewer);
1932: }
1933: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1934: return 0;
1935: }
1937: /* Check how well we did distributing points*/
1938: if (viewer) {
1939: PetscViewerASCIIPrintf(viewer, "Number of owned entities of depth %" PetscInt_FMT ".\n", entityDepth);
1940: PetscViewerASCIIPrintf(viewer, "Initial ");
1941: DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
1942: PetscViewerASCIIPrintf(viewer, "Rebalanced ");
1943: DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1944: }
1946: /* Check that every vertex is owned by a process that it is actually connected to. */
1947: MatGetRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1948: for (i = 1; i <= numNonExclusivelyOwned; i++) {
1949: PetscInt loc = 0;
1950: PetscFindInt(cumSumVertices[part[i]], xadj[i + 1] - xadj[i], &adjncy[xadj[i]], &loc);
1951: /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1952: if (loc < 0) part[i] = rank;
1953: }
1954: MatRestoreRowIJ(A, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, &numRows, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done);
1955: MatDestroy(&A);
1957: /* See how significant the influences of the previous fixing up step was.*/
1958: if (viewer) {
1959: PetscViewerASCIIPrintf(viewer, "After fix ");
1960: DMPlexViewDistribution(comm, cumSumVertices[rank + 1] - cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1961: }
1962: if (!usematpartitioning) PetscFree(vtxwgt);
1963: else MatPartitioningDestroy(&mp);
1965: PetscLayoutDestroy(&layout);
1967: PetscLogEventBegin(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);
1968: /* Rewrite the SF to reflect the new ownership. */
1969: PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
1970: counter = 0;
1971: for (i = 0; i < pEnd - pStart; i++) {
1972: if (toBalance[i]) {
1973: if (isNonExclusivelyOwned[i]) {
1974: pointsToRewrite[counter] = i + pStart;
1975: counter++;
1976: }
1977: }
1978: }
1979: DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part + 1, degrees);
1980: PetscFree(pointsToRewrite);
1981: PetscLogEventEnd(DMPLEX_RebalRewriteSF, dm, 0, 0, 0);
1983: PetscFree(toBalance);
1984: PetscFree(isLeaf);
1985: PetscFree(isNonExclusivelyOwned);
1986: PetscFree(isExclusivelyOwned);
1987: if (usematpartitioning) {
1988: ISRestoreIndices(ispart, (const PetscInt **)&part);
1989: ISDestroy(&ispart);
1990: } else PetscFree(part);
1991: if (viewer) {
1992: PetscViewerPopFormat(viewer);
1993: PetscViewerDestroy(&viewer);
1994: }
1995: if (success) *success = PETSC_TRUE;
1996: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1997: return 0;
1998: #else
1999: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2000: #endif
2001: }