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) { return point >= 0 ? point : -(point+1); }
9: static PetscErrorCode DMPlexCreatePartitionerGraph_Overlap(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
10: {
11: DM ovdm;
12: PetscSF sfPoint;
13: IS cellNumbering;
14: const PetscInt *cellNum;
15: PetscInt *adj = NULL, *vOffsets = NULL, *vAdj = NULL;
16: PetscBool useCone, useClosure;
17: PetscInt dim, depth, overlap, cStart, cEnd, c, v;
18: PetscMPIInt rank, size;
20: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
21: MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);
22: DMGetDimension(dm, &dim);
23: DMPlexGetDepth(dm, &depth);
24: if (dim != depth) {
25: /* We do not handle the uninterpolated case here */
26: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
27: /* DMPlexCreateNeighborCSR does not make a numbering */
28: if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
29: /* Different behavior for empty graphs */
30: if (!*numVertices) {
31: PetscMalloc1(1, offsets);
32: (*offsets)[0] = 0;
33: }
34: /* Broken in parallel */
36: return 0;
37: }
38: /* Always use FVM adjacency to create partitioner graph */
39: DMGetBasicAdjacency(dm, &useCone, &useClosure);
40: DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
41: /* Need overlap >= 1 */
42: DMPlexGetOverlap(dm, &overlap);
43: if (size && overlap < 1) {
44: DMPlexDistributeOverlap(dm, 1, NULL, &ovdm);
45: } else {
46: PetscObjectReference((PetscObject) dm);
47: ovdm = dm;
48: }
49: DMGetPointSF(ovdm, &sfPoint);
50: DMPlexGetHeightStratum(ovdm, height, &cStart, &cEnd);
51: DMPlexCreateNumbering_Plex(ovdm, cStart, cEnd, 0, NULL, sfPoint, &cellNumbering);
52: if (globalNumbering) {
53: PetscObjectReference((PetscObject) cellNumbering);
54: *globalNumbering = cellNumbering;
55: }
56: ISGetIndices(cellNumbering, &cellNum);
57: /* Determine sizes */
58: for (*numVertices = 0, c = cStart; c < cEnd; ++c) {
59: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
60: if (cellNum[c-cStart] < 0) continue;
61: (*numVertices)++;
62: }
63: PetscCalloc1(*numVertices+1, &vOffsets);
64: for (c = cStart, v = 0; c < cEnd; ++c) {
65: PetscInt adjSize = PETSC_DETERMINE, a, vsize = 0;
67: if (cellNum[c-cStart] < 0) continue;
68: DMPlexGetAdjacency(ovdm, c, &adjSize, &adj);
69: for (a = 0; a < adjSize; ++a) {
70: const PetscInt point = adj[a];
71: if (point != c && cStart <= point && point < cEnd) ++vsize;
72: }
73: vOffsets[v+1] = vOffsets[v] + vsize;
74: ++v;
75: }
76: /* Determine adjacency */
77: PetscMalloc1(vOffsets[*numVertices], &vAdj);
78: for (c = cStart, v = 0; c < cEnd; ++c) {
79: PetscInt adjSize = PETSC_DETERMINE, a, off = vOffsets[v];
81: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
82: if (cellNum[c-cStart] < 0) continue;
83: DMPlexGetAdjacency(ovdm, c, &adjSize, &adj);
84: for (a = 0; a < adjSize; ++a) {
85: const PetscInt point = adj[a];
86: if (point != c && cStart <= point && point < cEnd) {
87: vAdj[off++] = DMPlex_GlobalID(cellNum[point-cStart]);
88: }
89: }
91: /* Sort adjacencies (not strictly necessary) */
92: PetscSortInt(off-vOffsets[v], &vAdj[vOffsets[v]]);
93: ++v;
94: }
95: PetscFree(adj);
96: ISRestoreIndices(cellNumbering, &cellNum);
97: ISDestroy(&cellNumbering);
98: DMSetBasicAdjacency(dm, useCone, useClosure);
99: DMDestroy(&ovdm);
100: if (offsets) {*offsets = vOffsets;}
101: else PetscFree(vOffsets);
102: if (adjacency) {*adjacency = vAdj;}
103: else PetscFree(vAdj);
104: return 0;
105: }
107: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
108: {
109: PetscInt dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
110: PetscInt *adj = NULL, *vOffsets = NULL, *graph = NULL;
111: IS cellNumbering;
112: const PetscInt *cellNum;
113: PetscBool useCone, useClosure;
114: PetscSection section;
115: PetscSegBuffer adjBuffer;
116: PetscSF sfPoint;
117: PetscInt *adjCells = NULL, *remoteCells = NULL;
118: const PetscInt *local;
119: PetscInt nroots, nleaves, l;
120: PetscMPIInt rank;
122: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
123: DMGetDimension(dm, &dim);
124: DMPlexGetDepth(dm, &depth);
125: if (dim != depth) {
126: /* We do not handle the uninterpolated case here */
127: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
128: /* DMPlexCreateNeighborCSR does not make a numbering */
129: if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
130: /* Different behavior for empty graphs */
131: if (!*numVertices) {
132: PetscMalloc1(1, offsets);
133: (*offsets)[0] = 0;
134: }
135: /* Broken in parallel */
137: return 0;
138: }
139: DMGetPointSF(dm, &sfPoint);
140: DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
141: /* Build adjacency graph via a section/segbuffer */
142: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
143: PetscSectionSetChart(section, pStart, pEnd);
144: PetscSegBufferCreate(sizeof(PetscInt),1000,&adjBuffer);
145: /* Always use FVM adjacency to create partitioner graph */
146: DMGetBasicAdjacency(dm, &useCone, &useClosure);
147: DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
148: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
149: if (globalNumbering) {
150: PetscObjectReference((PetscObject)cellNumbering);
151: *globalNumbering = cellNumbering;
152: }
153: ISGetIndices(cellNumbering, &cellNum);
154: /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
155: Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
156: */
157: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
158: if (nroots >= 0) {
159: PetscInt fStart, fEnd, f;
161: PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
162: DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
163: for (l = 0; l < nroots; ++l) adjCells[l] = -3;
164: for (f = fStart; f < fEnd; ++f) {
165: const PetscInt *support;
166: PetscInt supportSize;
168: DMPlexGetSupport(dm, f, &support);
169: DMPlexGetSupportSize(dm, f, &supportSize);
170: if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]-pStart]);
171: else if (supportSize == 2) {
172: PetscFindInt(support[0], nleaves, local, &p);
173: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1]-pStart]);
174: PetscFindInt(support[1], nleaves, local, &p);
175: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]-pStart]);
176: }
177: /* Handle non-conforming meshes */
178: if (supportSize > 2) {
179: PetscInt numChildren, i;
180: const PetscInt *children;
182: DMPlexGetTreeChildren(dm, f, &numChildren, &children);
183: for (i = 0; i < numChildren; ++i) {
184: const PetscInt child = children[i];
185: if (fStart <= child && child < fEnd) {
186: DMPlexGetSupport(dm, child, &support);
187: DMPlexGetSupportSize(dm, child, &supportSize);
188: if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]-pStart]);
189: else if (supportSize == 2) {
190: PetscFindInt(support[0], nleaves, local, &p);
191: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1]-pStart]);
192: PetscFindInt(support[1], nleaves, local, &p);
193: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]-pStart]);
194: }
195: }
196: }
197: }
198: }
199: for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
200: PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells,MPI_REPLACE);
201: PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells,MPI_REPLACE);
202: PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
203: PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
204: }
205: /* Combine local and global adjacencies */
206: for (*numVertices = 0, p = pStart; p < pEnd; p++) {
207: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
208: if (nroots > 0) {if (cellNum[p-pStart] < 0) continue;}
209: /* Add remote cells */
210: if (remoteCells) {
211: const PetscInt gp = DMPlex_GlobalID(cellNum[p-pStart]);
212: PetscInt coneSize, numChildren, c, i;
213: const PetscInt *cone, *children;
215: DMPlexGetCone(dm, p, &cone);
216: DMPlexGetConeSize(dm, p, &coneSize);
217: for (c = 0; c < coneSize; ++c) {
218: const PetscInt point = cone[c];
219: if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
220: PetscInt *PETSC_RESTRICT pBuf;
221: PetscSectionAddDof(section, p, 1);
222: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
223: *pBuf = remoteCells[point];
224: }
225: /* Handle non-conforming meshes */
226: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
227: for (i = 0; i < numChildren; ++i) {
228: const PetscInt child = children[i];
229: if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
230: PetscInt *PETSC_RESTRICT pBuf;
231: PetscSectionAddDof(section, p, 1);
232: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
233: *pBuf = remoteCells[child];
234: }
235: }
236: }
237: }
238: /* Add local cells */
239: adjSize = PETSC_DETERMINE;
240: DMPlexGetAdjacency(dm, p, &adjSize, &adj);
241: for (a = 0; a < adjSize; ++a) {
242: const PetscInt point = adj[a];
243: if (point != p && pStart <= point && point < pEnd) {
244: PetscInt *PETSC_RESTRICT pBuf;
245: PetscSectionAddDof(section, p, 1);
246: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
247: *pBuf = DMPlex_GlobalID(cellNum[point-pStart]);
248: }
249: }
250: (*numVertices)++;
251: }
252: PetscFree(adj);
253: PetscFree2(adjCells, remoteCells);
254: DMSetBasicAdjacency(dm, useCone, useClosure);
256: /* Derive CSR graph from section/segbuffer */
257: PetscSectionSetUp(section);
258: PetscSectionGetStorageSize(section, &size);
259: PetscMalloc1(*numVertices+1, &vOffsets);
260: for (idx = 0, p = pStart; p < pEnd; p++) {
261: if (nroots > 0) {if (cellNum[p-pStart] < 0) continue;}
262: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
263: }
264: vOffsets[*numVertices] = size;
265: PetscSegBufferExtractAlloc(adjBuffer, &graph);
267: if (nroots >= 0) {
268: /* Filter out duplicate edges using section/segbuffer */
269: PetscSectionReset(section);
270: PetscSectionSetChart(section, 0, *numVertices);
271: for (p = 0; p < *numVertices; p++) {
272: PetscInt start = vOffsets[p], end = vOffsets[p+1];
273: PetscInt numEdges = end-start, *PETSC_RESTRICT edges;
274: PetscSortRemoveDupsInt(&numEdges, &graph[start]);
275: PetscSectionSetDof(section, p, numEdges);
276: PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
277: PetscArraycpy(edges, &graph[start], numEdges);
278: }
279: PetscFree(vOffsets);
280: PetscFree(graph);
281: /* Derive CSR graph from section/segbuffer */
282: PetscSectionSetUp(section);
283: PetscSectionGetStorageSize(section, &size);
284: PetscMalloc1(*numVertices+1, &vOffsets);
285: for (idx = 0, p = 0; p < *numVertices; p++) {
286: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
287: }
288: vOffsets[*numVertices] = size;
289: PetscSegBufferExtractAlloc(adjBuffer, &graph);
290: } else {
291: /* Sort adjacencies (not strictly necessary) */
292: for (p = 0; p < *numVertices; p++) {
293: PetscInt start = vOffsets[p], end = vOffsets[p+1];
294: PetscSortInt(end-start, &graph[start]);
295: }
296: }
298: if (offsets) *offsets = vOffsets;
299: if (adjacency) *adjacency = graph;
301: /* Cleanup */
302: PetscSegBufferDestroy(&adjBuffer);
303: PetscSectionDestroy(§ion);
304: ISRestoreIndices(cellNumbering, &cellNum);
305: ISDestroy(&cellNumbering);
306: return 0;
307: }
309: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
310: {
311: Mat conn, CSR;
312: IS fis, cis, cis_own;
313: PetscSF sfPoint;
314: const PetscInt *rows, *cols, *ii, *jj;
315: PetscInt *idxs,*idxs2;
316: PetscInt dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
317: PetscMPIInt rank;
318: PetscBool flg;
320: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
321: DMGetDimension(dm, &dim);
322: DMPlexGetDepth(dm, &depth);
323: if (dim != depth) {
324: /* We do not handle the uninterpolated case here */
325: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
326: /* DMPlexCreateNeighborCSR does not make a numbering */
327: if (globalNumbering) DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);
328: /* Different behavior for empty graphs */
329: if (!*numVertices) {
330: PetscMalloc1(1, offsets);
331: (*offsets)[0] = 0;
332: }
333: /* Broken in parallel */
335: return 0;
336: }
337: /* Interpolated and parallel case */
339: /* numbering */
340: DMGetPointSF(dm, &sfPoint);
341: DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
342: DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
343: DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
344: DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
345: if (globalNumbering) {
346: ISDuplicate(cis, globalNumbering);
347: }
349: /* get positive global ids and local sizes for facets and cells */
350: ISGetLocalSize(fis, &m);
351: ISGetIndices(fis, &rows);
352: PetscMalloc1(m, &idxs);
353: for (i = 0, floc = 0; i < m; i++) {
354: const PetscInt p = rows[i];
356: if (p < 0) {
357: idxs[i] = -(p+1);
358: } else {
359: idxs[i] = p;
360: floc += 1;
361: }
362: }
363: ISRestoreIndices(fis, &rows);
364: ISDestroy(&fis);
365: ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);
367: ISGetLocalSize(cis, &m);
368: ISGetIndices(cis, &cols);
369: PetscMalloc1(m, &idxs);
370: PetscMalloc1(m, &idxs2);
371: for (i = 0, cloc = 0; i < m; i++) {
372: const PetscInt p = cols[i];
374: if (p < 0) {
375: idxs[i] = -(p+1);
376: } else {
377: idxs[i] = p;
378: idxs2[cloc++] = p;
379: }
380: }
381: ISRestoreIndices(cis, &cols);
382: ISDestroy(&cis);
383: ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
384: ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);
386: /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
387: MatCreate(PetscObjectComm((PetscObject)dm), &conn);
388: MatSetSizes(conn, floc, cloc, M, N);
389: MatSetType(conn, MATMPIAIJ);
390: DMPlexGetMaxSizes(dm, NULL, &lm);
391: MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dm));
392: MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);
394: /* Assemble matrix */
395: ISGetIndices(fis, &rows);
396: ISGetIndices(cis, &cols);
397: for (c = cStart; c < cEnd; c++) {
398: const PetscInt *cone;
399: PetscInt coneSize, row, col, f;
401: col = cols[c-cStart];
402: DMPlexGetCone(dm, c, &cone);
403: DMPlexGetConeSize(dm, c, &coneSize);
404: for (f = 0; f < coneSize; f++) {
405: const PetscScalar v = 1.0;
406: const PetscInt *children;
407: PetscInt numChildren, ch;
409: row = rows[cone[f]-fStart];
410: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
412: /* non-conforming meshes */
413: DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
414: for (ch = 0; ch < numChildren; ch++) {
415: const PetscInt child = children[ch];
417: if (child < fStart || child >= fEnd) continue;
418: row = rows[child-fStart];
419: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
420: }
421: }
422: }
423: ISRestoreIndices(fis, &rows);
424: ISRestoreIndices(cis, &cols);
425: ISDestroy(&fis);
426: ISDestroy(&cis);
427: MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
428: MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);
430: /* Get parallel CSR by doing conn^T * conn */
431: MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
432: MatDestroy(&conn);
434: /* extract local part of the CSR */
435: MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
436: MatDestroy(&CSR);
437: MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
440: /* get back requested output */
441: if (numVertices) *numVertices = m;
442: if (offsets) {
443: PetscCalloc1(m+1, &idxs);
444: for (i = 1; i < m+1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
445: *offsets = idxs;
446: }
447: if (adjacency) {
448: PetscMalloc1(ii[m] - m, &idxs);
449: ISGetIndices(cis_own, &rows);
450: for (i = 0, c = 0; i < m; i++) {
451: PetscInt j, g = rows[i];
453: for (j = ii[i]; j < ii[i+1]; j++) {
454: if (jj[j] == g) continue; /* again, self-connectivity */
455: idxs[c++] = jj[j];
456: }
457: }
459: ISRestoreIndices(cis_own, &rows);
460: *adjacency = idxs;
461: }
463: /* cleanup */
464: ISDestroy(&cis_own);
465: MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
467: MatDestroy(&conn);
468: return 0;
469: }
471: /*@C
472: DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner
474: Input Parameters:
475: + dm - The mesh DM dm
476: - height - Height of the strata from which to construct the graph
478: Output Parameters:
479: + numVertices - Number of vertices in the graph
480: . offsets - Point offsets in the graph
481: . adjacency - Point connectivity in the graph
482: - 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.
484: The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function
485: representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().
487: Options Database Keys:
488: . -dm_plex_csr_alg <mat,graph,overlap> - Choose the algorithm for computing the CSR graph
490: Level: developer
492: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
493: @*/
494: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
495: {
496: DMPlexCSRAlgorithm alg = DM_PLEX_CSR_GRAPH;
498: PetscOptionsGetEnum(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_csr_alg", DMPlexCSRAlgorithms, (PetscEnum *) &alg, NULL);
499: switch (alg) {
500: case DM_PLEX_CSR_MAT:
501: DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);break;
502: case DM_PLEX_CSR_GRAPH:
503: DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);break;
504: case DM_PLEX_CSR_OVERLAP:
505: DMPlexCreatePartitionerGraph_Overlap(dm, height, numVertices, offsets, adjacency, globalNumbering);break;
506: }
507: return 0;
508: }
510: /*@C
511: DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.
513: Collective on DM
515: Input Parameters:
516: + dm - The DMPlex
517: - cellHeight - The height of mesh points to treat as cells (default should be 0)
519: Output Parameters:
520: + numVertices - The number of local vertices in the graph, or cells in the mesh.
521: . offsets - The offset to the adjacency list for each cell
522: - adjacency - The adjacency list for all cells
524: Note: This is suitable for input to a mesh partitioner like ParMetis.
526: Level: advanced
528: .seealso: DMPlexCreate()
529: @*/
530: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
531: {
532: const PetscInt maxFaceCases = 30;
533: PetscInt numFaceCases = 0;
534: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
535: PetscInt *off, *adj;
536: PetscInt *neighborCells = NULL;
537: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
539: /* For parallel partitioning, I think you have to communicate supports */
540: DMGetDimension(dm, &dim);
541: cellDim = dim - cellHeight;
542: DMPlexGetDepth(dm, &depth);
543: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
544: if (cEnd - cStart == 0) {
545: if (numVertices) *numVertices = 0;
546: if (offsets) *offsets = NULL;
547: if (adjacency) *adjacency = NULL;
548: return 0;
549: }
550: numCells = cEnd - cStart;
551: faceDepth = depth - cellHeight;
552: if (dim == depth) {
553: PetscInt f, fStart, fEnd;
555: PetscCalloc1(numCells+1, &off);
556: /* Count neighboring cells */
557: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
558: for (f = fStart; f < fEnd; ++f) {
559: const PetscInt *support;
560: PetscInt supportSize;
561: DMPlexGetSupportSize(dm, f, &supportSize);
562: DMPlexGetSupport(dm, f, &support);
563: if (supportSize == 2) {
564: PetscInt numChildren;
566: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
567: if (!numChildren) {
568: ++off[support[0]-cStart+1];
569: ++off[support[1]-cStart+1];
570: }
571: }
572: }
573: /* Prefix sum */
574: for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
575: if (adjacency) {
576: PetscInt *tmp;
578: PetscMalloc1(off[numCells], &adj);
579: PetscMalloc1(numCells+1, &tmp);
580: PetscArraycpy(tmp, off, numCells+1);
581: /* Get neighboring cells */
582: for (f = fStart; f < fEnd; ++f) {
583: const PetscInt *support;
584: PetscInt supportSize;
585: DMPlexGetSupportSize(dm, f, &supportSize);
586: DMPlexGetSupport(dm, f, &support);
587: if (supportSize == 2) {
588: PetscInt numChildren;
590: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
591: if (!numChildren) {
592: adj[tmp[support[0]-cStart]++] = support[1];
593: adj[tmp[support[1]-cStart]++] = support[0];
594: }
595: }
596: }
597: for (c = 0; c < cEnd-cStart; ++c) PetscAssert(tmp[c] == off[c+1],PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %d != %d for cell %d", tmp[c], off[c], c+cStart);
598: PetscFree(tmp);
599: }
600: if (numVertices) *numVertices = numCells;
601: if (offsets) *offsets = off;
602: if (adjacency) *adjacency = adj;
603: return 0;
604: }
605: /* Setup face recognition */
606: if (faceDepth == 1) {
607: 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 */
609: for (c = cStart; c < cEnd; ++c) {
610: PetscInt corners;
612: DMPlexGetConeSize(dm, c, &corners);
613: if (!cornersSeen[corners]) {
614: PetscInt nFV;
617: cornersSeen[corners] = 1;
619: DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);
621: numFaceVertices[numFaceCases++] = nFV;
622: }
623: }
624: }
625: PetscCalloc1(numCells+1, &off);
626: /* Count neighboring cells */
627: for (cell = cStart; cell < cEnd; ++cell) {
628: PetscInt numNeighbors = PETSC_DETERMINE, n;
630: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
631: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
632: for (n = 0; n < numNeighbors; ++n) {
633: PetscInt cellPair[2];
634: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
635: PetscInt meetSize = 0;
636: const PetscInt *meet = NULL;
638: cellPair[0] = cell; cellPair[1] = neighborCells[n];
639: if (cellPair[0] == cellPair[1]) continue;
640: if (!found) {
641: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
642: if (meetSize) {
643: PetscInt f;
645: for (f = 0; f < numFaceCases; ++f) {
646: if (numFaceVertices[f] == meetSize) {
647: found = PETSC_TRUE;
648: break;
649: }
650: }
651: }
652: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
653: }
654: if (found) ++off[cell-cStart+1];
655: }
656: }
657: /* Prefix sum */
658: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];
660: if (adjacency) {
661: PetscMalloc1(off[numCells], &adj);
662: /* Get neighboring cells */
663: for (cell = cStart; cell < cEnd; ++cell) {
664: PetscInt numNeighbors = PETSC_DETERMINE, n;
665: PetscInt cellOffset = 0;
667: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
668: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
669: for (n = 0; n < numNeighbors; ++n) {
670: PetscInt cellPair[2];
671: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
672: PetscInt meetSize = 0;
673: const PetscInt *meet = NULL;
675: cellPair[0] = cell; cellPair[1] = neighborCells[n];
676: if (cellPair[0] == cellPair[1]) continue;
677: if (!found) {
678: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
679: if (meetSize) {
680: PetscInt f;
682: for (f = 0; f < numFaceCases; ++f) {
683: if (numFaceVertices[f] == meetSize) {
684: found = PETSC_TRUE;
685: break;
686: }
687: }
688: }
689: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
690: }
691: if (found) {
692: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
693: ++cellOffset;
694: }
695: }
696: }
697: }
698: PetscFree(neighborCells);
699: if (numVertices) *numVertices = numCells;
700: if (offsets) *offsets = off;
701: if (adjacency) *adjacency = adj;
702: return 0;
703: }
705: /*@
706: PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh
708: Collective on PetscPartitioner
710: Input Parameters:
711: + part - The PetscPartitioner
712: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
713: - dm - The mesh DM
715: Output Parameters:
716: + partSection - The PetscSection giving the division of points by partition
717: - partition - The list of points by partition
719: Notes:
720: If the DM has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
721: by the section in the transitive closure of the point.
723: Level: developer
725: .seealso DMPlexDistribute(), PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscPartitionerPartition()
726: @*/
727: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
728: {
729: PetscMPIInt size;
730: PetscBool isplex;
731: PetscSection vertSection = NULL;
738: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isplex);
740: MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
741: if (size == 1) {
742: PetscInt *points;
743: PetscInt cStart, cEnd, c;
745: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
746: PetscSectionReset(partSection);
747: PetscSectionSetChart(partSection, 0, size);
748: PetscSectionSetDof(partSection, 0, cEnd-cStart);
749: PetscSectionSetUp(partSection);
750: PetscMalloc1(cEnd-cStart, &points);
751: for (c = cStart; c < cEnd; ++c) points[c] = c;
752: ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
753: return 0;
754: }
755: if (part->height == 0) {
756: PetscInt numVertices = 0;
757: PetscInt *start = NULL;
758: PetscInt *adjacency = NULL;
759: IS globalNumbering;
761: if (!part->noGraph || part->viewGraph) {
762: DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
763: } else { /* only compute the number of owned local vertices */
764: const PetscInt *idxs;
765: PetscInt p, pStart, pEnd;
767: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
768: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
769: ISGetIndices(globalNumbering, &idxs);
770: for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
771: ISRestoreIndices(globalNumbering, &idxs);
772: }
773: if (part->usevwgt) {
774: PetscSection section = dm->localSection, clSection = NULL;
775: IS clPoints = NULL;
776: const PetscInt *gid,*clIdx;
777: PetscInt v, p, pStart, pEnd;
779: /* 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) */
780: /* We do this only if the local section has been set */
781: if (section) {
782: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
783: if (!clSection) {
784: DMPlexCreateClosureIndex(dm,NULL);
785: }
786: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
787: ISGetIndices(clPoints,&clIdx);
788: }
789: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
790: PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
791: PetscSectionSetChart(vertSection, 0, numVertices);
792: if (globalNumbering) {
793: ISGetIndices(globalNumbering,&gid);
794: } else gid = NULL;
795: for (p = pStart, v = 0; p < pEnd; ++p) {
796: PetscInt dof = 1;
798: /* skip cells in the overlap */
799: if (gid && gid[p-pStart] < 0) continue;
801: if (section) {
802: PetscInt cl, clSize, clOff;
804: dof = 0;
805: PetscSectionGetDof(clSection, p, &clSize);
806: PetscSectionGetOffset(clSection, p, &clOff);
807: for (cl = 0; cl < clSize; cl+=2) {
808: PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */
810: PetscSectionGetDof(section, clPoint, &clDof);
811: dof += clDof;
812: }
813: }
815: PetscSectionSetDof(vertSection, v, dof);
816: v++;
817: }
818: if (globalNumbering) {
819: ISRestoreIndices(globalNumbering,&gid);
820: }
821: if (clPoints) {
822: ISRestoreIndices(clPoints,&clIdx);
823: }
824: PetscSectionSetUp(vertSection);
825: }
826: PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
827: PetscFree(start);
828: PetscFree(adjacency);
829: if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
830: const PetscInt *globalNum;
831: const PetscInt *partIdx;
832: PetscInt *map, cStart, cEnd;
833: PetscInt *adjusted, i, localSize, offset;
834: IS newPartition;
836: ISGetLocalSize(*partition,&localSize);
837: PetscMalloc1(localSize,&adjusted);
838: ISGetIndices(globalNumbering,&globalNum);
839: ISGetIndices(*partition,&partIdx);
840: PetscMalloc1(localSize,&map);
841: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
842: for (i = cStart, offset = 0; i < cEnd; i++) {
843: if (globalNum[i - cStart] >= 0) map[offset++] = i;
844: }
845: for (i = 0; i < localSize; i++) {
846: adjusted[i] = map[partIdx[i]];
847: }
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 %D 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++) {
927: DMPlexAddClosure_Private(dm, ht, cone[c]);
928: }
929: }
930: return 0;
931: }
933: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
934: {
935: if (up) {
936: PetscInt parent;
938: DMPlexGetTreeParent(dm,point,&parent,NULL);
939: if (parent != point) {
940: PetscInt closureSize, *closure = NULL, i;
942: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
943: for (i = 0; i < closureSize; i++) {
944: PetscInt cpoint = closure[2*i];
946: PetscHSetIAdd(ht, cpoint);
947: DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
948: }
949: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
950: }
951: }
952: if (down) {
953: PetscInt numChildren;
954: const PetscInt *children;
956: DMPlexGetTreeChildren(dm,point,&numChildren,&children);
957: if (numChildren) {
958: PetscInt i;
960: for (i = 0; i < numChildren; i++) {
961: PetscInt cpoint = children[i];
963: PetscHSetIAdd(ht, cpoint);
964: DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
965: }
966: }
967: }
968: return 0;
969: }
971: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
972: {
973: PetscInt parent;
976: DMPlexGetTreeParent(dm, point, &parent,NULL);
977: if (point != parent) {
978: const PetscInt *cone;
979: PetscInt coneSize, c;
981: DMPlexAddClosureTree_Up_Private(dm, ht, parent);
982: DMPlexAddClosure_Private(dm, ht, parent);
983: DMPlexGetCone(dm, parent, &cone);
984: DMPlexGetConeSize(dm, parent, &coneSize);
985: for (c = 0; c < coneSize; c++) {
986: const PetscInt cp = cone[c];
988: DMPlexAddClosureTree_Up_Private(dm, ht, cp);
989: }
990: }
991: return 0;
992: }
994: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
995: {
996: PetscInt i, numChildren;
997: const PetscInt *children;
1000: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
1001: for (i = 0; i < numChildren; i++) {
1002: PetscHSetIAdd(ht, children[i]);
1003: }
1004: return 0;
1005: }
1007: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
1008: {
1009: const PetscInt *cone;
1010: PetscInt coneSize, c;
1013: PetscHSetIAdd(ht, point);
1014: DMPlexAddClosureTree_Up_Private(dm, ht, point);
1015: DMPlexAddClosureTree_Down_Private(dm, ht, point);
1016: DMPlexGetCone(dm, point, &cone);
1017: DMPlexGetConeSize(dm, point, &coneSize);
1018: for (c = 0; c < coneSize; c++) {
1019: DMPlexAddClosureTree_Private(dm, ht, cone[c]);
1020: }
1021: return 0;
1022: }
1024: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
1025: {
1026: DM_Plex *mesh = (DM_Plex *)dm->data;
1027: const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
1028: PetscInt nelems, *elems, off = 0, p;
1029: PetscHSetI ht = NULL;
1031: PetscHSetICreate(&ht);
1032: PetscHSetIResize(ht, numPoints*16);
1033: if (!hasTree) {
1034: for (p = 0; p < numPoints; ++p) {
1035: DMPlexAddClosure_Private(dm, ht, points[p]);
1036: }
1037: } else {
1038: #if 1
1039: for (p = 0; p < numPoints; ++p) {
1040: DMPlexAddClosureTree_Private(dm, ht, points[p]);
1041: }
1042: #else
1043: PetscInt *closure = NULL, closureSize, c;
1044: for (p = 0; p < numPoints; ++p) {
1045: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
1046: for (c = 0; c < closureSize*2; c += 2) {
1047: PetscHSetIAdd(ht, closure[c]);
1048: if (hasTree) DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);
1049: }
1050: }
1051: if (closure) DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);
1052: #endif
1053: }
1054: PetscHSetIGetSize(ht, &nelems);
1055: PetscMalloc1(nelems, &elems);
1056: PetscHSetIGetElems(ht, &off, elems);
1057: PetscHSetIDestroy(&ht);
1058: PetscSortInt(nelems, elems);
1059: ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
1060: return 0;
1061: }
1063: /*@
1064: DMPlexPartitionLabelClosure - Add the closure of all points to the partition label
1066: Input Parameters:
1067: + dm - The DM
1068: - label - DMLabel assigning ranks to remote roots
1070: Level: developer
1072: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1073: @*/
1074: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
1075: {
1076: IS rankIS, pointIS, closureIS;
1077: const PetscInt *ranks, *points;
1078: PetscInt numRanks, numPoints, r;
1080: DMLabelGetValueIS(label, &rankIS);
1081: ISGetLocalSize(rankIS, &numRanks);
1082: ISGetIndices(rankIS, &ranks);
1083: for (r = 0; r < numRanks; ++r) {
1084: const PetscInt rank = ranks[r];
1085: DMLabelGetStratumIS(label, rank, &pointIS);
1086: ISGetLocalSize(pointIS, &numPoints);
1087: ISGetIndices(pointIS, &points);
1088: DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
1089: ISRestoreIndices(pointIS, &points);
1090: ISDestroy(&pointIS);
1091: DMLabelSetStratumIS(label, rank, closureIS);
1092: ISDestroy(&closureIS);
1093: }
1094: ISRestoreIndices(rankIS, &ranks);
1095: ISDestroy(&rankIS);
1096: return 0;
1097: }
1099: /*@
1100: DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label
1102: Input Parameters:
1103: + dm - The DM
1104: - label - DMLabel assigning ranks to remote roots
1106: Level: developer
1108: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1109: @*/
1110: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1111: {
1112: IS rankIS, pointIS;
1113: const PetscInt *ranks, *points;
1114: PetscInt numRanks, numPoints, r, p, a, adjSize;
1115: PetscInt *adj = NULL;
1117: DMLabelGetValueIS(label, &rankIS);
1118: ISGetLocalSize(rankIS, &numRanks);
1119: ISGetIndices(rankIS, &ranks);
1120: for (r = 0; r < numRanks; ++r) {
1121: const PetscInt rank = ranks[r];
1123: DMLabelGetStratumIS(label, rank, &pointIS);
1124: ISGetLocalSize(pointIS, &numPoints);
1125: ISGetIndices(pointIS, &points);
1126: for (p = 0; p < numPoints; ++p) {
1127: adjSize = PETSC_DETERMINE;
1128: DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
1129: for (a = 0; a < adjSize; ++a) DMLabelSetValue(label, adj[a], rank);
1130: }
1131: ISRestoreIndices(pointIS, &points);
1132: ISDestroy(&pointIS);
1133: }
1134: ISRestoreIndices(rankIS, &ranks);
1135: ISDestroy(&rankIS);
1136: PetscFree(adj);
1137: return 0;
1138: }
1140: /*@
1141: DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point SF
1143: Input Parameters:
1144: + dm - The DM
1145: - label - DMLabel assigning ranks to remote roots
1147: Level: developer
1149: Note: This is required when generating multi-level overlaps to capture
1150: overlap points from non-neighbouring partitions.
1152: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1153: @*/
1154: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1155: {
1156: MPI_Comm comm;
1157: PetscMPIInt rank;
1158: PetscSF sfPoint;
1159: DMLabel lblRoots, lblLeaves;
1160: IS rankIS, pointIS;
1161: const PetscInt *ranks;
1162: PetscInt numRanks, r;
1164: PetscObjectGetComm((PetscObject) dm, &comm);
1165: MPI_Comm_rank(comm, &rank);
1166: DMGetPointSF(dm, &sfPoint);
1167: /* Pull point contributions from remote leaves into local roots */
1168: DMLabelGather(label, sfPoint, &lblLeaves);
1169: DMLabelGetValueIS(lblLeaves, &rankIS);
1170: ISGetLocalSize(rankIS, &numRanks);
1171: ISGetIndices(rankIS, &ranks);
1172: for (r = 0; r < numRanks; ++r) {
1173: const PetscInt remoteRank = ranks[r];
1174: if (remoteRank == rank) continue;
1175: DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
1176: DMLabelInsertIS(label, pointIS, remoteRank);
1177: ISDestroy(&pointIS);
1178: }
1179: ISRestoreIndices(rankIS, &ranks);
1180: ISDestroy(&rankIS);
1181: DMLabelDestroy(&lblLeaves);
1182: /* Push point contributions from roots into remote leaves */
1183: DMLabelDistribute(label, sfPoint, &lblRoots);
1184: DMLabelGetValueIS(lblRoots, &rankIS);
1185: ISGetLocalSize(rankIS, &numRanks);
1186: ISGetIndices(rankIS, &ranks);
1187: for (r = 0; r < numRanks; ++r) {
1188: const PetscInt remoteRank = ranks[r];
1189: if (remoteRank == rank) continue;
1190: DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
1191: DMLabelInsertIS(label, pointIS, remoteRank);
1192: ISDestroy(&pointIS);
1193: }
1194: ISRestoreIndices(rankIS, &ranks);
1195: ISDestroy(&rankIS);
1196: DMLabelDestroy(&lblRoots);
1197: return 0;
1198: }
1200: /*@
1201: DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label
1203: Input Parameters:
1204: + dm - The DM
1205: . rootLabel - DMLabel assigning ranks to local roots
1206: - processSF - A star forest mapping into the local index on each remote rank
1208: Output Parameter:
1209: . leafLabel - DMLabel assigning ranks to remote roots
1211: Note: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
1212: resulting leafLabel is a receiver mapping of remote roots to their parent rank.
1214: Level: developer
1216: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1217: @*/
1218: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1219: {
1220: MPI_Comm comm;
1221: PetscMPIInt rank, size, r;
1222: PetscInt p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1223: PetscSF sfPoint;
1224: PetscSection rootSection;
1225: PetscSFNode *rootPoints, *leafPoints;
1226: const PetscSFNode *remote;
1227: const PetscInt *local, *neighbors;
1228: IS valueIS;
1229: PetscBool mpiOverflow = PETSC_FALSE;
1231: PetscLogEventBegin(DMPLEX_PartLabelInvert,dm,0,0,0);
1232: PetscObjectGetComm((PetscObject) dm, &comm);
1233: MPI_Comm_rank(comm, &rank);
1234: MPI_Comm_size(comm, &size);
1235: DMGetPointSF(dm, &sfPoint);
1237: /* Convert to (point, rank) and use actual owners */
1238: PetscSectionCreate(comm, &rootSection);
1239: PetscSectionSetChart(rootSection, 0, size);
1240: DMLabelGetValueIS(rootLabel, &valueIS);
1241: ISGetLocalSize(valueIS, &numNeighbors);
1242: ISGetIndices(valueIS, &neighbors);
1243: for (n = 0; n < numNeighbors; ++n) {
1244: DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
1245: PetscSectionAddDof(rootSection, neighbors[n], numPoints);
1246: }
1247: PetscSectionSetUp(rootSection);
1248: PetscSectionGetStorageSize(rootSection, &rootSize);
1249: PetscMalloc1(rootSize, &rootPoints);
1250: PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
1251: for (n = 0; n < numNeighbors; ++n) {
1252: IS pointIS;
1253: const PetscInt *points;
1255: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1256: DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
1257: ISGetLocalSize(pointIS, &numPoints);
1258: ISGetIndices(pointIS, &points);
1259: for (p = 0; p < numPoints; ++p) {
1260: if (local) PetscFindInt(points[p], nleaves, local, &l);
1261: else {l = -1;}
1262: if (l >= 0) {rootPoints[off+p] = remote[l];}
1263: else {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
1264: }
1265: ISRestoreIndices(pointIS, &points);
1266: ISDestroy(&pointIS);
1267: }
1269: /* Try to communicate overlap using All-to-All */
1270: if (!processSF) {
1271: PetscInt64 counter = 0;
1272: PetscBool locOverflow = PETSC_FALSE;
1273: PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;
1275: PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
1276: for (n = 0; n < numNeighbors; ++n) {
1277: PetscSectionGetDof(rootSection, neighbors[n], &dof);
1278: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1279: #if defined(PETSC_USE_64BIT_INDICES)
1280: if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
1281: if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
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) { rdispls[r] = (int)counter; counter += rcounts[r]; }
1288: if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
1289: MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
1290: if (!mpiOverflow) {
1291: PetscInfo(dm,"Using Alltoallv for mesh distribution\n");
1292: leafSize = (PetscInt) counter;
1293: PetscMalloc1(leafSize, &leafPoints);
1294: MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
1295: }
1296: PetscFree4(scounts, sdispls, rcounts, rdispls);
1297: }
1299: /* Communicate overlap using process star forest */
1300: if (processSF || mpiOverflow) {
1301: PetscSF procSF;
1302: PetscSection leafSection;
1304: if (processSF) {
1305: PetscInfo(dm,"Using processSF for mesh distribution\n");
1306: PetscObjectReference((PetscObject)processSF);
1307: procSF = processSF;
1308: } else {
1309: PetscInfo(dm,"Using processSF for mesh distribution (MPI overflow)\n");
1310: PetscSFCreate(comm,&procSF);
1311: PetscSFSetGraphWithPattern(procSF,NULL,PETSCSF_PATTERN_ALLTOALL);
1312: }
1314: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
1315: DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
1316: PetscSectionGetStorageSize(leafSection, &leafSize);
1317: PetscSectionDestroy(&leafSection);
1318: PetscSFDestroy(&procSF);
1319: }
1321: for (p = 0; p < leafSize; p++) {
1322: DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
1323: }
1325: ISRestoreIndices(valueIS, &neighbors);
1326: ISDestroy(&valueIS);
1327: PetscSectionDestroy(&rootSection);
1328: PetscFree(rootPoints);
1329: PetscFree(leafPoints);
1330: PetscLogEventEnd(DMPLEX_PartLabelInvert,dm,0,0,0);
1331: return 0;
1332: }
1334: /*@
1335: DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points
1337: Input Parameters:
1338: + dm - The DM
1339: - label - DMLabel assigning ranks to remote roots
1341: Output Parameter:
1342: . sf - The star forest communication context encapsulating the defined mapping
1344: Note: The incoming label is a receiver mapping of remote points to their parent rank.
1346: Level: developer
1348: .seealso: DMPlexDistribute(), DMPlexCreateOverlap()
1349: @*/
1350: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
1351: {
1352: PetscMPIInt rank;
1353: PetscInt n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1354: PetscSFNode *remotePoints;
1355: IS remoteRootIS, neighborsIS;
1356: const PetscInt *remoteRoots, *neighbors;
1358: PetscLogEventBegin(DMPLEX_PartLabelCreateSF,dm,0,0,0);
1359: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1361: DMLabelGetValueIS(label, &neighborsIS);
1362: #if 0
1363: {
1364: IS is;
1365: ISDuplicate(neighborsIS, &is);
1366: ISSort(is);
1367: ISDestroy(&neighborsIS);
1368: neighborsIS = is;
1369: }
1370: #endif
1371: ISGetLocalSize(neighborsIS, &nNeighbors);
1372: ISGetIndices(neighborsIS, &neighbors);
1373: for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1374: DMLabelGetStratumSize(label, neighbors[n], &numPoints);
1375: numRemote += numPoints;
1376: }
1377: PetscMalloc1(numRemote, &remotePoints);
1378: /* Put owned points first */
1379: DMLabelGetStratumSize(label, rank, &numPoints);
1380: if (numPoints > 0) {
1381: DMLabelGetStratumIS(label, rank, &remoteRootIS);
1382: ISGetIndices(remoteRootIS, &remoteRoots);
1383: for (p = 0; p < numPoints; p++) {
1384: remotePoints[idx].index = remoteRoots[p];
1385: remotePoints[idx].rank = rank;
1386: idx++;
1387: }
1388: ISRestoreIndices(remoteRootIS, &remoteRoots);
1389: ISDestroy(&remoteRootIS);
1390: }
1391: /* Now add remote points */
1392: for (n = 0; n < nNeighbors; ++n) {
1393: const PetscInt nn = neighbors[n];
1395: DMLabelGetStratumSize(label, nn, &numPoints);
1396: if (nn == rank || numPoints <= 0) continue;
1397: DMLabelGetStratumIS(label, nn, &remoteRootIS);
1398: ISGetIndices(remoteRootIS, &remoteRoots);
1399: for (p = 0; p < numPoints; p++) {
1400: remotePoints[idx].index = remoteRoots[p];
1401: remotePoints[idx].rank = nn;
1402: idx++;
1403: }
1404: ISRestoreIndices(remoteRootIS, &remoteRoots);
1405: ISDestroy(&remoteRootIS);
1406: }
1407: PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
1408: DMPlexGetChart(dm, &pStart, &pEnd);
1409: PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1410: PetscSFSetUp(*sf);
1411: ISDestroy(&neighborsIS);
1412: PetscLogEventEnd(DMPLEX_PartLabelCreateSF,dm,0,0,0);
1413: return 0;
1414: }
1416: #if defined(PETSC_HAVE_PARMETIS)
1417: #include <parmetis.h>
1418: #endif
1420: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
1421: * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
1422: * them out in that case. */
1423: #if defined(PETSC_HAVE_PARMETIS)
1424: /*@C
1426: DMPlexRewriteSF - Rewrites the ownership of the SF of a DM (in place).
1428: Input parameters:
1429: + dm - The DMPlex object.
1430: . n - The number of points.
1431: . pointsToRewrite - The points in the SF whose ownership will change.
1432: . targetOwners - New owner for each element in pointsToRewrite.
1433: - degrees - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.
1435: Level: developer
1437: @*/
1438: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1439: {
1440: PetscInt pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1441: PetscInt *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1442: PetscSFNode *leafLocationsNew;
1443: const PetscSFNode *iremote;
1444: const PetscInt *ilocal;
1445: PetscBool *isLeaf;
1446: PetscSF sf;
1447: MPI_Comm comm;
1448: PetscMPIInt rank, size;
1450: PetscObjectGetComm((PetscObject) dm, &comm);
1451: MPI_Comm_rank(comm, &rank);
1452: MPI_Comm_size(comm, &size);
1453: DMPlexGetChart(dm, &pStart, &pEnd);
1455: DMGetPointSF(dm, &sf);
1456: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1457: PetscMalloc1(pEnd-pStart, &isLeaf);
1458: for (i=0; i<pEnd-pStart; i++) {
1459: isLeaf[i] = PETSC_FALSE;
1460: }
1461: for (i=0; i<nleafs; i++) {
1462: isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
1463: }
1465: PetscMalloc1(pEnd-pStart+1, &cumSumDegrees);
1466: cumSumDegrees[0] = 0;
1467: for (i=1; i<=pEnd-pStart; i++) {
1468: cumSumDegrees[i] = cumSumDegrees[i-1] + degrees[i-1];
1469: }
1470: sumDegrees = cumSumDegrees[pEnd-pStart];
1471: /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */
1473: PetscMalloc1(sumDegrees, &locationsOfLeafs);
1474: PetscMalloc1(pEnd-pStart, &rankOnLeafs);
1475: for (i=0; i<pEnd-pStart; i++) {
1476: rankOnLeafs[i] = rank;
1477: }
1478: PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1479: PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1480: PetscFree(rankOnLeafs);
1482: /* get the remote local points of my leaves */
1483: PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
1484: PetscMalloc1(pEnd-pStart, &points);
1485: for (i=0; i<pEnd-pStart; i++) {
1486: points[i] = pStart+i;
1487: }
1488: PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1489: PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1490: PetscFree(points);
1491: /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1492: PetscMalloc1(pEnd-pStart, &newOwners);
1493: PetscMalloc1(pEnd-pStart, &newNumbers);
1494: for (i=0; i<pEnd-pStart; i++) {
1495: newOwners[i] = -1;
1496: newNumbers[i] = -1;
1497: }
1498: {
1499: PetscInt oldNumber, newNumber, oldOwner, newOwner;
1500: for (i=0; i<n; i++) {
1501: oldNumber = pointsToRewrite[i];
1502: newNumber = -1;
1503: oldOwner = rank;
1504: newOwner = targetOwners[i];
1505: if (oldOwner == newOwner) {
1506: newNumber = oldNumber;
1507: } else {
1508: for (j=0; j<degrees[oldNumber]; j++) {
1509: if (locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) {
1510: newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j];
1511: break;
1512: }
1513: }
1514: }
1517: newOwners[oldNumber] = newOwner;
1518: newNumbers[oldNumber] = newNumber;
1519: }
1520: }
1521: PetscFree(cumSumDegrees);
1522: PetscFree(locationsOfLeafs);
1523: PetscFree(remoteLocalPointOfLeafs);
1525: PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners,MPI_REPLACE);
1526: PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners,MPI_REPLACE);
1527: PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers,MPI_REPLACE);
1528: PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers,MPI_REPLACE);
1530: /* Now count how many leafs we have on each processor. */
1531: leafCounter=0;
1532: for (i=0; i<pEnd-pStart; i++) {
1533: if (newOwners[i] >= 0) {
1534: if (newOwners[i] != rank) {
1535: leafCounter++;
1536: }
1537: } else {
1538: if (isLeaf[i]) {
1539: leafCounter++;
1540: }
1541: }
1542: }
1544: /* Now set up the new sf by creating the leaf arrays */
1545: PetscMalloc1(leafCounter, &leafsNew);
1546: PetscMalloc1(leafCounter, &leafLocationsNew);
1548: leafCounter = 0;
1549: counter = 0;
1550: for (i=0; i<pEnd-pStart; i++) {
1551: if (newOwners[i] >= 0) {
1552: if (newOwners[i] != rank) {
1553: leafsNew[leafCounter] = i;
1554: leafLocationsNew[leafCounter].rank = newOwners[i];
1555: leafLocationsNew[leafCounter].index = newNumbers[i];
1556: leafCounter++;
1557: }
1558: } else {
1559: if (isLeaf[i]) {
1560: leafsNew[leafCounter] = i;
1561: leafLocationsNew[leafCounter].rank = iremote[counter].rank;
1562: leafLocationsNew[leafCounter].index = iremote[counter].index;
1563: leafCounter++;
1564: }
1565: }
1566: if (isLeaf[i]) {
1567: counter++;
1568: }
1569: }
1571: PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
1572: PetscFree(newOwners);
1573: PetscFree(newNumbers);
1574: PetscFree(isLeaf);
1575: return 0;
1576: }
1578: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1579: {
1580: PetscInt *distribution, min, max, sum;
1581: PetscMPIInt rank, size;
1583: MPI_Comm_size(comm, &size);
1584: MPI_Comm_rank(comm, &rank);
1585: PetscCalloc1(size, &distribution);
1586: for (PetscInt i=0; i<n; i++) {
1587: if (part) distribution[part[i]] += vtxwgt[skip*i];
1588: else distribution[rank] += vtxwgt[skip*i];
1589: }
1590: MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
1591: min = distribution[0];
1592: max = distribution[0];
1593: sum = distribution[0];
1594: for (PetscInt i=1; i<size; i++) {
1595: if (distribution[i]<min) min=distribution[i];
1596: if (distribution[i]>max) max=distribution[i];
1597: sum += distribution[i];
1598: }
1599: PetscViewerASCIIPrintf(viewer, "Min: %D, Avg: %D, Max: %D, Balance: %f\n", min, sum/size, max, (max*1.*size)/sum);
1600: PetscFree(distribution);
1601: return 0;
1602: }
1604: #endif
1606: /*@
1607: DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the PointSF of the DM inplace.
1609: Input parameters:
1610: + dm - The DMPlex object.
1611: . entityDepth - depth of the entity to balance (0 -> balance vertices).
1612: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
1613: - parallel - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.
1615: Output parameters:
1616: . success - whether the graph partitioning was successful or not. If not, try useInitialGuess=True and parallel=True.
1618: Level: intermediate
1620: @*/
1622: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1623: {
1624: #if defined(PETSC_HAVE_PARMETIS)
1625: PetscSF sf;
1626: PetscInt ierr, i, j, idx, jdx;
1627: PetscInt eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1628: const PetscInt *degrees, *ilocal;
1629: const PetscSFNode *iremote;
1630: PetscBool *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1631: PetscInt numExclusivelyOwned, numNonExclusivelyOwned;
1632: PetscMPIInt rank, size;
1633: PetscInt *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1634: const PetscInt *cumSumVertices;
1635: PetscInt offset, counter;
1636: PetscInt lenadjncy;
1637: PetscInt *xadj, *adjncy, *vtxwgt;
1638: PetscInt lenxadj;
1639: PetscInt *adjwgt = NULL;
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, Apre;
1647: const char *prefix = NULL;
1648: PetscViewer viewer;
1649: PetscViewerFormat format;
1650: PetscLayout layout;
1652: if (success) *success = PETSC_FALSE;
1653: PetscObjectGetComm((PetscObject) dm, &comm);
1654: MPI_Comm_rank(comm, &rank);
1655: MPI_Comm_size(comm, &size);
1656: if (size==1) return 0;
1658: PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1660: PetscOptionsGetViewer(comm,((PetscObject)dm)->options, prefix,"-dm_rebalance_partition_view",&viewer,&format,NULL);
1661: if (viewer) {
1662: PetscViewerPushFormat(viewer,format);
1663: }
1665: /* Figure out all points in the plex that we are interested in balancing. */
1666: DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
1667: DMPlexGetChart(dm, &pStart, &pEnd);
1668: PetscMalloc1(pEnd-pStart, &toBalance);
1670: for (i=0; i<pEnd-pStart; i++) {
1671: toBalance[i] = (PetscBool)(i-pStart>=eBegin && i-pStart<eEnd);
1672: }
1674: /* There are three types of points:
1675: * exclusivelyOwned: points that are owned by this process and only seen by this process
1676: * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1677: * leaf: a point that is seen by this process but owned by a different process
1678: */
1679: DMGetPointSF(dm, &sf);
1680: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1681: PetscMalloc1(pEnd-pStart, &isLeaf);
1682: PetscMalloc1(pEnd-pStart, &isNonExclusivelyOwned);
1683: PetscMalloc1(pEnd-pStart, &isExclusivelyOwned);
1684: for (i=0; i<pEnd-pStart; i++) {
1685: isNonExclusivelyOwned[i] = PETSC_FALSE;
1686: isExclusivelyOwned[i] = PETSC_FALSE;
1687: isLeaf[i] = PETSC_FALSE;
1688: }
1690: /* start by marking all the leafs */
1691: for (i=0; i<nleafs; i++) {
1692: isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
1693: }
1695: /* for an owned point, we can figure out whether another processor sees it or
1696: * not by calculating its degree */
1697: PetscSFComputeDegreeBegin(sf, °rees);
1698: PetscSFComputeDegreeEnd(sf, °rees);
1700: numExclusivelyOwned = 0;
1701: numNonExclusivelyOwned = 0;
1702: for (i=0; i<pEnd-pStart; i++) {
1703: if (toBalance[i]) {
1704: if (degrees[i] > 0) {
1705: isNonExclusivelyOwned[i] = PETSC_TRUE;
1706: numNonExclusivelyOwned += 1;
1707: } else {
1708: if (!isLeaf[i]) {
1709: isExclusivelyOwned[i] = PETSC_TRUE;
1710: numExclusivelyOwned += 1;
1711: }
1712: }
1713: }
1714: }
1716: /* We are going to build a graph with one vertex per core representing the
1717: * exclusively owned points and then one vertex per nonExclusively owned
1718: * point. */
1720: PetscLayoutCreate(comm, &layout);
1721: PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
1722: PetscLayoutSetUp(layout);
1723: PetscLayoutGetRanges(layout, &cumSumVertices);
1725: PetscMalloc1(pEnd-pStart, &globalNumbersOfLocalOwnedVertices);
1726: for (i=0; i<pEnd-pStart; i++) {globalNumbersOfLocalOwnedVertices[i] = pStart - 1;}
1727: offset = cumSumVertices[rank];
1728: counter = 0;
1729: for (i=0; i<pEnd-pStart; i++) {
1730: if (toBalance[i]) {
1731: if (degrees[i] > 0) {
1732: globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1733: counter++;
1734: }
1735: }
1736: }
1738: /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
1739: PetscMalloc1(pEnd-pStart, &leafGlobalNumbers);
1740: PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers,MPI_REPLACE);
1741: PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers,MPI_REPLACE);
1743: /* Now start building the data structure for ParMETIS */
1745: MatCreate(comm, &Apre);
1746: MatSetType(Apre, MATPREALLOCATOR);
1747: MatSetSizes(Apre, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
1748: MatSetUp(Apre);
1750: for (i=0; i<pEnd-pStart; i++) {
1751: if (toBalance[i]) {
1752: idx = cumSumVertices[rank];
1753: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1754: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1755: else continue;
1756: MatSetValue(Apre, idx, jdx, 1, INSERT_VALUES);
1757: MatSetValue(Apre, jdx, idx, 1, INSERT_VALUES);
1758: }
1759: }
1761: MatAssemblyBegin(Apre, MAT_FINAL_ASSEMBLY);
1762: MatAssemblyEnd(Apre, MAT_FINAL_ASSEMBLY);
1764: MatCreate(comm, &A);
1765: MatSetType(A, MATMPIAIJ);
1766: MatSetSizes(A, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
1767: MatPreallocatorPreallocate(Apre, PETSC_FALSE, A);
1768: MatDestroy(&Apre);
1770: for (i=0; i<pEnd-pStart; i++) {
1771: if (toBalance[i]) {
1772: idx = cumSumVertices[rank];
1773: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1774: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1775: else continue;
1776: MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
1777: MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
1778: }
1779: }
1781: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1782: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1783: PetscFree(leafGlobalNumbers);
1784: PetscFree(globalNumbersOfLocalOwnedVertices);
1786: nparts = size;
1787: wgtflag = 2;
1788: numflag = 0;
1789: ncon = 2;
1790: real_t *tpwgts;
1791: PetscMalloc1(ncon * nparts, &tpwgts);
1792: for (i=0; i<ncon*nparts; i++) {
1793: tpwgts[i] = 1./(nparts);
1794: }
1796: PetscMalloc1(ncon, &ubvec);
1797: ubvec[0] = 1.01;
1798: ubvec[1] = 1.01;
1799: lenadjncy = 0;
1800: for (i=0; i<1+numNonExclusivelyOwned; i++) {
1801: PetscInt temp=0;
1802: MatGetRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
1803: lenadjncy += temp;
1804: MatRestoreRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
1805: }
1806: PetscMalloc1(lenadjncy, &adjncy);
1807: lenxadj = 2 + numNonExclusivelyOwned;
1808: PetscMalloc1(lenxadj, &xadj);
1809: xadj[0] = 0;
1810: counter = 0;
1811: for (i=0; i<1+numNonExclusivelyOwned; i++) {
1812: PetscInt temp=0;
1813: const PetscInt *cols;
1814: MatGetRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
1815: PetscArraycpy(&adjncy[counter], cols, temp);
1816: counter += temp;
1817: xadj[i+1] = counter;
1818: MatRestoreRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
1819: }
1821: PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);
1822: PetscMalloc1(ncon*(1 + numNonExclusivelyOwned), &vtxwgt);
1823: vtxwgt[0] = numExclusivelyOwned;
1824: if (ncon>1) vtxwgt[1] = 1;
1825: for (i=0; i<numNonExclusivelyOwned; i++) {
1826: vtxwgt[ncon*(i+1)] = 1;
1827: if (ncon>1) vtxwgt[ncon*(i+1)+1] = 0;
1828: }
1830: if (viewer) {
1831: PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %D on interface of mesh distribution.\n", entityDepth);
1832: PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %D\n", cumSumVertices[size]);
1833: }
1834: if (parallel) {
1835: PetscMalloc1(4, &options);
1836: options[0] = 1;
1837: options[1] = 0; /* Verbosity */
1838: options[2] = 0; /* Seed */
1839: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1840: if (viewer) PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n");
1841: if (useInitialGuess) {
1842: if (viewer) PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n");
1843: PetscStackPush("ParMETIS_V3_RefineKway");
1844: ParMETIS_V3_RefineKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1846: PetscStackPop;
1847: } else {
1848: PetscStackPush("ParMETIS_V3_PartKway");
1849: ParMETIS_V3_PartKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1850: PetscStackPop;
1852: }
1853: PetscFree(options);
1854: } else {
1855: if (viewer) PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n");
1856: Mat As;
1857: PetscInt numRows;
1858: PetscInt *partGlobal;
1859: MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &As);
1861: PetscInt *numExclusivelyOwnedAll;
1862: PetscMalloc1(size, &numExclusivelyOwnedAll);
1863: numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1864: MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,numExclusivelyOwnedAll,1,MPIU_INT,comm);
1866: MatGetSize(As, &numRows, NULL);
1867: PetscMalloc1(numRows, &partGlobal);
1868: if (rank == 0) {
1869: PetscInt *adjncy_g, *xadj_g, *vtxwgt_g;
1870: lenadjncy = 0;
1872: for (i=0; i<numRows; i++) {
1873: PetscInt temp=0;
1874: MatGetRow(As, i, &temp, NULL, NULL);
1875: lenadjncy += temp;
1876: MatRestoreRow(As, i, &temp, NULL, NULL);
1877: }
1878: PetscMalloc1(lenadjncy, &adjncy_g);
1879: lenxadj = 1 + numRows;
1880: PetscMalloc1(lenxadj, &xadj_g);
1881: xadj_g[0] = 0;
1882: counter = 0;
1883: for (i=0; i<numRows; i++) {
1884: PetscInt temp=0;
1885: const PetscInt *cols;
1886: MatGetRow(As, i, &temp, &cols, NULL);
1887: PetscArraycpy(&adjncy_g[counter], cols, temp);
1888: counter += temp;
1889: xadj_g[i+1] = counter;
1890: MatRestoreRow(As, i, &temp, &cols, NULL);
1891: }
1892: PetscMalloc1(2*numRows, &vtxwgt_g);
1893: for (i=0; i<size; i++) {
1894: vtxwgt_g[ncon*cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1895: if (ncon>1) vtxwgt_g[ncon*cumSumVertices[i]+1] = 1;
1896: for (j=cumSumVertices[i]+1; j<cumSumVertices[i+1]; j++) {
1897: vtxwgt_g[ncon*j] = 1;
1898: if (ncon>1) vtxwgt_g[2*j+1] = 0;
1899: }
1900: }
1901: PetscMalloc1(64, &options);
1902: METIS_SetDefaultOptions(options); /* initialize all defaults */
1904: options[METIS_OPTION_CONTIG] = 1;
1905: PetscStackPush("METIS_PartGraphKway");
1906: METIS_PartGraphKway(&numRows, &ncon, xadj_g, adjncy_g, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1907: PetscStackPop;
1909: PetscFree(options);
1910: PetscFree(xadj_g);
1911: PetscFree(adjncy_g);
1912: PetscFree(vtxwgt_g);
1913: }
1914: PetscFree(numExclusivelyOwnedAll);
1916: /* Now scatter the parts array. */
1917: {
1918: PetscMPIInt *counts, *mpiCumSumVertices;
1919: PetscMalloc1(size, &counts);
1920: PetscMalloc1(size+1, &mpiCumSumVertices);
1921: for (i=0; i<size; i++) {
1922: PetscMPIIntCast(cumSumVertices[i+1] - cumSumVertices[i], &(counts[i]));
1923: }
1924: for (i=0; i<=size; i++) {
1925: PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
1926: }
1927: MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
1928: PetscFree(counts);
1929: PetscFree(mpiCumSumVertices);
1930: }
1932: PetscFree(partGlobal);
1933: MatDestroy(&As);
1934: }
1936: MatDestroy(&A);
1937: PetscFree(ubvec);
1938: PetscFree(tpwgts);
1940: /* Now rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1942: PetscMalloc1(size, &firstVertices);
1943: PetscMalloc1(size, &renumbering);
1944: firstVertices[rank] = part[0];
1945: MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);
1946: for (i=0; i<size; i++) {
1947: renumbering[firstVertices[i]] = i;
1948: }
1949: for (i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) {
1950: part[i] = renumbering[part[i]];
1951: }
1952: /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1953: failed = (PetscInt)(part[0] != rank);
1954: MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
1956: PetscFree(firstVertices);
1957: PetscFree(renumbering);
1959: if (failedGlobal > 0) {
1960: PetscLayoutDestroy(&layout);
1961: PetscFree(xadj);
1962: PetscFree(adjncy);
1963: PetscFree(vtxwgt);
1964: PetscFree(toBalance);
1965: PetscFree(isLeaf);
1966: PetscFree(isNonExclusivelyOwned);
1967: PetscFree(isExclusivelyOwned);
1968: PetscFree(part);
1969: if (viewer) {
1970: PetscViewerPopFormat(viewer);
1971: PetscViewerDestroy(&viewer);
1972: }
1973: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1974: return 0;
1975: }
1977: /*Let's check how well we did distributing points*/
1978: if (viewer) {
1979: PetscViewerASCIIPrintf(viewer, "Comparing number of owned entities of depth %D on each process before rebalancing, after rebalancing, and after consistency checks.\n", entityDepth);
1980: PetscViewerASCIIPrintf(viewer, "Initial. ");
1981: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
1982: PetscViewerASCIIPrintf(viewer, "Rebalanced. ");
1983: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1984: }
1986: /* Now check that every vertex is owned by a process that it is actually connected to. */
1987: for (i=1; i<=numNonExclusivelyOwned; i++) {
1988: PetscInt loc = 0;
1989: PetscFindInt(cumSumVertices[part[i]], xadj[i+1]-xadj[i], &adjncy[xadj[i]], &loc);
1990: /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1991: if (loc<0) {
1992: part[i] = rank;
1993: }
1994: }
1996: /* Let's see how significant the influences of the previous fixing up step was.*/
1997: if (viewer) {
1998: PetscViewerASCIIPrintf(viewer, "After. ");
1999: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
2000: }
2002: PetscLayoutDestroy(&layout);
2003: PetscFree(xadj);
2004: PetscFree(adjncy);
2005: PetscFree(vtxwgt);
2007: /* Almost done, now rewrite the SF to reflect the new ownership. */
2008: {
2009: PetscInt *pointsToRewrite;
2010: PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
2011: counter = 0;
2012: for (i=0; i<pEnd-pStart; i++) {
2013: if (toBalance[i]) {
2014: if (isNonExclusivelyOwned[i]) {
2015: pointsToRewrite[counter] = i + pStart;
2016: counter++;
2017: }
2018: }
2019: }
2020: DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part+1, degrees);
2021: PetscFree(pointsToRewrite);
2022: }
2024: PetscFree(toBalance);
2025: PetscFree(isLeaf);
2026: PetscFree(isNonExclusivelyOwned);
2027: PetscFree(isExclusivelyOwned);
2028: PetscFree(part);
2029: if (viewer) {
2030: PetscViewerPopFormat(viewer);
2031: PetscViewerDestroy(&viewer);
2032: }
2033: if (success) *success = PETSC_TRUE;
2034: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
2035: return 0;
2036: #else
2037: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
2038: #endif
2039: }