Actual source code: plexpartition.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/partitionerimpl.h>
3: #include <petsc/private/hashseti.h>
5: PETSC_STATIC_INLINE PetscInt DMPlex_GlobalID(PetscInt point) { return point >= 0 ? point : -(point+1); }
7: static PetscErrorCode DMPlexCreatePartitionerGraph_Native(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
8: {
9: PetscInt dim, depth, p, pStart, pEnd, a, adjSize, idx, size;
10: PetscInt *adj = NULL, *vOffsets = NULL, *graph = NULL;
11: IS cellNumbering;
12: const PetscInt *cellNum;
13: PetscBool useCone, useClosure;
14: PetscSection section;
15: PetscSegBuffer adjBuffer;
16: PetscSF sfPoint;
17: PetscInt *adjCells = NULL, *remoteCells = NULL;
18: const PetscInt *local;
19: PetscInt nroots, nleaves, l;
20: PetscMPIInt rank;
24: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
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 */
38: if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
39: return(0);
40: }
41: DMGetPointSF(dm, &sfPoint);
42: DMPlexGetHeightStratum(dm, height, &pStart, &pEnd);
43: /* Build adjacency graph via a section/segbuffer */
44: PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);
45: PetscSectionSetChart(section, pStart, pEnd);
46: PetscSegBufferCreate(sizeof(PetscInt),1000,&adjBuffer);
47: /* Always use FVM adjacency to create partitioner graph */
48: DMGetBasicAdjacency(dm, &useCone, &useClosure);
49: DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
50: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, sfPoint, &cellNumbering);
51: if (globalNumbering) {
52: PetscObjectReference((PetscObject)cellNumbering);
53: *globalNumbering = cellNumbering;
54: }
55: ISGetIndices(cellNumbering, &cellNum);
56: /* For all boundary faces (including faces adjacent to a ghost cell), record the local cell in adjCells
57: Broadcast adjCells to remoteCells (to get cells from roots) and Reduce adjCells to remoteCells (to get cells from leaves)
58: */
59: PetscSFGetGraph(sfPoint, &nroots, &nleaves, &local, NULL);
60: if (nroots >= 0) {
61: PetscInt fStart, fEnd, f;
63: PetscCalloc2(nroots, &adjCells, nroots, &remoteCells);
64: DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
65: for (l = 0; l < nroots; ++l) adjCells[l] = -3;
66: for (f = fStart; f < fEnd; ++f) {
67: const PetscInt *support;
68: PetscInt supportSize;
70: DMPlexGetSupport(dm, f, &support);
71: DMPlexGetSupportSize(dm, f, &supportSize);
72: if (supportSize == 1) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
73: else if (supportSize == 2) {
74: PetscFindInt(support[0], nleaves, local, &p);
75: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[1]]);
76: PetscFindInt(support[1], nleaves, local, &p);
77: if (p >= 0) adjCells[f] = DMPlex_GlobalID(cellNum[support[0]]);
78: }
79: /* Handle non-conforming meshes */
80: if (supportSize > 2) {
81: PetscInt numChildren, i;
82: const PetscInt *children;
84: DMPlexGetTreeChildren(dm, f, &numChildren, &children);
85: for (i = 0; i < numChildren; ++i) {
86: const PetscInt child = children[i];
87: if (fStart <= child && child < fEnd) {
88: DMPlexGetSupport(dm, child, &support);
89: DMPlexGetSupportSize(dm, child, &supportSize);
90: if (supportSize == 1) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
91: else if (supportSize == 2) {
92: PetscFindInt(support[0], nleaves, local, &p);
93: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[1]]);
94: PetscFindInt(support[1], nleaves, local, &p);
95: if (p >= 0) adjCells[child] = DMPlex_GlobalID(cellNum[support[0]]);
96: }
97: }
98: }
99: }
100: }
101: for (l = 0; l < nroots; ++l) remoteCells[l] = -1;
102: PetscSFBcastBegin(dm->sf, MPIU_INT, adjCells, remoteCells);
103: PetscSFBcastEnd(dm->sf, MPIU_INT, adjCells, remoteCells);
104: PetscSFReduceBegin(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
105: PetscSFReduceEnd(dm->sf, MPIU_INT, adjCells, remoteCells, MPI_MAX);
106: }
107: /* Combine local and global adjacencies */
108: for (*numVertices = 0, p = pStart; p < pEnd; p++) {
109: /* Skip non-owned cells in parallel (ParMetis expects no overlap) */
110: if (nroots > 0) {if (cellNum[p] < 0) continue;}
111: /* Add remote cells */
112: if (remoteCells) {
113: const PetscInt gp = DMPlex_GlobalID(cellNum[p]);
114: PetscInt coneSize, numChildren, c, i;
115: const PetscInt *cone, *children;
117: DMPlexGetCone(dm, p, &cone);
118: DMPlexGetConeSize(dm, p, &coneSize);
119: for (c = 0; c < coneSize; ++c) {
120: const PetscInt point = cone[c];
121: if (remoteCells[point] >= 0 && remoteCells[point] != gp) {
122: PetscInt *PETSC_RESTRICT pBuf;
123: PetscSectionAddDof(section, p, 1);
124: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
125: *pBuf = remoteCells[point];
126: }
127: /* Handle non-conforming meshes */
128: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
129: for (i = 0; i < numChildren; ++i) {
130: const PetscInt child = children[i];
131: if (remoteCells[child] >= 0 && remoteCells[child] != gp) {
132: PetscInt *PETSC_RESTRICT pBuf;
133: PetscSectionAddDof(section, p, 1);
134: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
135: *pBuf = remoteCells[child];
136: }
137: }
138: }
139: }
140: /* Add local cells */
141: adjSize = PETSC_DETERMINE;
142: DMPlexGetAdjacency(dm, p, &adjSize, &adj);
143: for (a = 0; a < adjSize; ++a) {
144: const PetscInt point = adj[a];
145: if (point != p && pStart <= point && point < pEnd) {
146: PetscInt *PETSC_RESTRICT pBuf;
147: PetscSectionAddDof(section, p, 1);
148: PetscSegBufferGetInts(adjBuffer, 1, &pBuf);
149: *pBuf = DMPlex_GlobalID(cellNum[point]);
150: }
151: }
152: (*numVertices)++;
153: }
154: PetscFree(adj);
155: PetscFree2(adjCells, remoteCells);
156: DMSetBasicAdjacency(dm, useCone, useClosure);
158: /* Derive CSR graph from section/segbuffer */
159: PetscSectionSetUp(section);
160: PetscSectionGetStorageSize(section, &size);
161: PetscMalloc1(*numVertices+1, &vOffsets);
162: for (idx = 0, p = pStart; p < pEnd; p++) {
163: if (nroots > 0) {if (cellNum[p] < 0) continue;}
164: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
165: }
166: vOffsets[*numVertices] = size;
167: PetscSegBufferExtractAlloc(adjBuffer, &graph);
169: if (nroots >= 0) {
170: /* Filter out duplicate edges using section/segbuffer */
171: PetscSectionReset(section);
172: PetscSectionSetChart(section, 0, *numVertices);
173: for (p = 0; p < *numVertices; p++) {
174: PetscInt start = vOffsets[p], end = vOffsets[p+1];
175: PetscInt numEdges = end-start, *PETSC_RESTRICT edges;
176: PetscSortRemoveDupsInt(&numEdges, &graph[start]);
177: PetscSectionSetDof(section, p, numEdges);
178: PetscSegBufferGetInts(adjBuffer, numEdges, &edges);
179: PetscArraycpy(edges, &graph[start], numEdges);
180: }
181: PetscFree(vOffsets);
182: PetscFree(graph);
183: /* Derive CSR graph from section/segbuffer */
184: PetscSectionSetUp(section);
185: PetscSectionGetStorageSize(section, &size);
186: PetscMalloc1(*numVertices+1, &vOffsets);
187: for (idx = 0, p = 0; p < *numVertices; p++) {
188: PetscSectionGetOffset(section, p, &(vOffsets[idx++]));
189: }
190: vOffsets[*numVertices] = size;
191: PetscSegBufferExtractAlloc(adjBuffer, &graph);
192: } else {
193: /* Sort adjacencies (not strictly necessary) */
194: for (p = 0; p < *numVertices; p++) {
195: PetscInt start = vOffsets[p], end = vOffsets[p+1];
196: PetscSortInt(end-start, &graph[start]);
197: }
198: }
200: if (offsets) *offsets = vOffsets;
201: if (adjacency) *adjacency = graph;
203: /* Cleanup */
204: PetscSegBufferDestroy(&adjBuffer);
205: PetscSectionDestroy(§ion);
206: ISRestoreIndices(cellNumbering, &cellNum);
207: ISDestroy(&cellNumbering);
208: return(0);
209: }
211: static PetscErrorCode DMPlexCreatePartitionerGraph_ViaMat(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
212: {
213: Mat conn, CSR;
214: IS fis, cis, cis_own;
215: PetscSF sfPoint;
216: const PetscInt *rows, *cols, *ii, *jj;
217: PetscInt *idxs,*idxs2;
218: PetscInt dim, depth, floc, cloc, i, M, N, c, lm, m, cStart, cEnd, fStart, fEnd;
219: PetscMPIInt rank;
220: PetscBool flg;
224: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
225: DMGetDimension(dm, &dim);
226: DMPlexGetDepth(dm, &depth);
227: if (dim != depth) {
228: /* We do not handle the uninterpolated case here */
229: DMPlexCreateNeighborCSR(dm, height, numVertices, offsets, adjacency);
230: /* DMPlexCreateNeighborCSR does not make a numbering */
231: if (globalNumbering) {DMPlexCreateCellNumbering_Internal(dm, PETSC_TRUE, globalNumbering);}
232: /* Different behavior for empty graphs */
233: if (!*numVertices) {
234: PetscMalloc1(1, offsets);
235: (*offsets)[0] = 0;
236: }
237: /* Broken in parallel */
238: if (rank && *numVertices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parallel partitioning of uninterpolated meshes not supported");
239: return(0);
240: }
241: /* Interpolated and parallel case */
243: /* numbering */
244: DMGetPointSF(dm, &sfPoint);
245: DMPlexGetHeightStratum(dm, height, &cStart, &cEnd);
246: DMPlexGetHeightStratum(dm, height+1, &fStart, &fEnd);
247: DMPlexCreateNumbering_Plex(dm, cStart, cEnd, 0, &N, sfPoint, &cis);
248: DMPlexCreateNumbering_Plex(dm, fStart, fEnd, 0, &M, sfPoint, &fis);
249: if (globalNumbering) {
250: ISDuplicate(cis, globalNumbering);
251: }
253: /* get positive global ids and local sizes for facets and cells */
254: ISGetLocalSize(fis, &m);
255: ISGetIndices(fis, &rows);
256: PetscMalloc1(m, &idxs);
257: for (i = 0, floc = 0; i < m; i++) {
258: const PetscInt p = rows[i];
260: if (p < 0) {
261: idxs[i] = -(p+1);
262: } else {
263: idxs[i] = p;
264: floc += 1;
265: }
266: }
267: ISRestoreIndices(fis, &rows);
268: ISDestroy(&fis);
269: ISCreateGeneral(PETSC_COMM_SELF, m, idxs, PETSC_OWN_POINTER, &fis);
271: ISGetLocalSize(cis, &m);
272: ISGetIndices(cis, &cols);
273: PetscMalloc1(m, &idxs);
274: PetscMalloc1(m, &idxs2);
275: for (i = 0, cloc = 0; i < m; i++) {
276: const PetscInt p = cols[i];
278: if (p < 0) {
279: idxs[i] = -(p+1);
280: } else {
281: idxs[i] = p;
282: idxs2[cloc++] = p;
283: }
284: }
285: ISRestoreIndices(cis, &cols);
286: ISDestroy(&cis);
287: ISCreateGeneral(PetscObjectComm((PetscObject)dm), m, idxs, PETSC_OWN_POINTER, &cis);
288: ISCreateGeneral(PetscObjectComm((PetscObject)dm), cloc, idxs2, PETSC_OWN_POINTER, &cis_own);
290: /* Create matrix to hold F-C connectivity (MatMatTranspose Mult not supported for MPIAIJ) */
291: MatCreate(PetscObjectComm((PetscObject)dm), &conn);
292: MatSetSizes(conn, floc, cloc, M, N);
293: MatSetType(conn, MATMPIAIJ);
294: DMPlexGetMaxSizes(dm, NULL, &lm);
295: MPI_Allreduce(&lm, &m, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) dm));
296: MatMPIAIJSetPreallocation(conn, m, NULL, m, NULL);
298: /* Assemble matrix */
299: ISGetIndices(fis, &rows);
300: ISGetIndices(cis, &cols);
301: for (c = cStart; c < cEnd; c++) {
302: const PetscInt *cone;
303: PetscInt coneSize, row, col, f;
305: col = cols[c-cStart];
306: DMPlexGetCone(dm, c, &cone);
307: DMPlexGetConeSize(dm, c, &coneSize);
308: for (f = 0; f < coneSize; f++) {
309: const PetscScalar v = 1.0;
310: const PetscInt *children;
311: PetscInt numChildren, ch;
313: row = rows[cone[f]-fStart];
314: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
316: /* non-conforming meshes */
317: DMPlexGetTreeChildren(dm, cone[f], &numChildren, &children);
318: for (ch = 0; ch < numChildren; ch++) {
319: const PetscInt child = children[ch];
321: if (child < fStart || child >= fEnd) continue;
322: row = rows[child-fStart];
323: MatSetValues(conn, 1, &row, 1, &col, &v, INSERT_VALUES);
324: }
325: }
326: }
327: ISRestoreIndices(fis, &rows);
328: ISRestoreIndices(cis, &cols);
329: ISDestroy(&fis);
330: ISDestroy(&cis);
331: MatAssemblyBegin(conn, MAT_FINAL_ASSEMBLY);
332: MatAssemblyEnd(conn, MAT_FINAL_ASSEMBLY);
334: /* Get parallel CSR by doing conn^T * conn */
335: MatTransposeMatMult(conn, conn, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CSR);
336: MatDestroy(&conn);
338: /* extract local part of the CSR */
339: MatMPIAIJGetLocalMat(CSR, MAT_INITIAL_MATRIX, &conn);
340: MatDestroy(&CSR);
341: MatGetRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
342: if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
344: /* get back requested output */
345: if (numVertices) *numVertices = m;
346: if (offsets) {
347: PetscCalloc1(m+1, &idxs);
348: for (i = 1; i < m+1; i++) idxs[i] = ii[i] - i; /* ParMetis does not like self-connectivity */
349: *offsets = idxs;
350: }
351: if (adjacency) {
352: PetscMalloc1(ii[m] - m, &idxs);
353: ISGetIndices(cis_own, &rows);
354: for (i = 0, c = 0; i < m; i++) {
355: PetscInt j, g = rows[i];
357: for (j = ii[i]; j < ii[i+1]; j++) {
358: if (jj[j] == g) continue; /* again, self-connectivity */
359: idxs[c++] = jj[j];
360: }
361: }
362: if (c != ii[m] - m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %D != %D",c,ii[m]-m);
363: ISRestoreIndices(cis_own, &rows);
364: *adjacency = idxs;
365: }
367: /* cleanup */
368: ISDestroy(&cis_own);
369: MatRestoreRowIJ(conn, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg);
370: if (!flg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "No IJ format");
371: MatDestroy(&conn);
372: return(0);
373: }
375: /*@C
376: DMPlexCreatePartitionerGraph - Create a CSR graph of point connections for the partitioner
378: Input Parameters:
379: + dm - The mesh DM dm
380: - height - Height of the strata from which to construct the graph
382: Output Parameter:
383: + numVertices - Number of vertices in the graph
384: . offsets - Point offsets in the graph
385: . adjacency - Point connectivity in the graph
386: - 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.
388: The user can control the definition of adjacency for the mesh using DMSetAdjacency(). They should choose the combination appropriate for the function
389: representation on the mesh. If requested, globalNumbering needs to be destroyed by the caller; offsets and adjacency need to be freed with PetscFree().
391: Level: developer
393: .seealso: PetscPartitionerGetType(), PetscPartitionerCreate(), DMSetAdjacency()
394: @*/
395: PetscErrorCode DMPlexCreatePartitionerGraph(DM dm, PetscInt height, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency, IS *globalNumbering)
396: {
398: PetscBool usemat = PETSC_FALSE;
401: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_csr_via_mat", &usemat, NULL);
402: if (usemat) {
403: DMPlexCreatePartitionerGraph_ViaMat(dm, height, numVertices, offsets, adjacency, globalNumbering);
404: } else {
405: DMPlexCreatePartitionerGraph_Native(dm, height, numVertices, offsets, adjacency, globalNumbering);
406: }
407: return(0);
408: }
410: /*@C
411: DMPlexCreateNeighborCSR - Create a mesh graph (cell-cell adjacency) in parallel CSR format.
413: Collective on DM
415: Input Arguments:
416: + dm - The DMPlex
417: - cellHeight - The height of mesh points to treat as cells (default should be 0)
419: Output Arguments:
420: + numVertices - The number of local vertices in the graph, or cells in the mesh.
421: . offsets - The offset to the adjacency list for each cell
422: - adjacency - The adjacency list for all cells
424: Note: This is suitable for input to a mesh partitioner like ParMetis.
426: Level: advanced
428: .seealso: DMPlexCreate()
429: @*/
430: PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
431: {
432: const PetscInt maxFaceCases = 30;
433: PetscInt numFaceCases = 0;
434: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
435: PetscInt *off, *adj;
436: PetscInt *neighborCells = NULL;
437: PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell;
441: /* For parallel partitioning, I think you have to communicate supports */
442: DMGetDimension(dm, &dim);
443: cellDim = dim - cellHeight;
444: DMPlexGetDepth(dm, &depth);
445: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
446: if (cEnd - cStart == 0) {
447: if (numVertices) *numVertices = 0;
448: if (offsets) *offsets = NULL;
449: if (adjacency) *adjacency = NULL;
450: return(0);
451: }
452: numCells = cEnd - cStart;
453: faceDepth = depth - cellHeight;
454: if (dim == depth) {
455: PetscInt f, fStart, fEnd;
457: PetscCalloc1(numCells+1, &off);
458: /* Count neighboring cells */
459: DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);
460: for (f = fStart; f < fEnd; ++f) {
461: const PetscInt *support;
462: PetscInt supportSize;
463: DMPlexGetSupportSize(dm, f, &supportSize);
464: DMPlexGetSupport(dm, f, &support);
465: if (supportSize == 2) {
466: PetscInt numChildren;
468: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
469: if (!numChildren) {
470: ++off[support[0]-cStart+1];
471: ++off[support[1]-cStart+1];
472: }
473: }
474: }
475: /* Prefix sum */
476: for (c = 1; c <= numCells; ++c) off[c] += off[c-1];
477: if (adjacency) {
478: PetscInt *tmp;
480: PetscMalloc1(off[numCells], &adj);
481: PetscMalloc1(numCells+1, &tmp);
482: PetscArraycpy(tmp, off, numCells+1);
483: /* Get neighboring cells */
484: for (f = fStart; f < fEnd; ++f) {
485: const PetscInt *support;
486: PetscInt supportSize;
487: DMPlexGetSupportSize(dm, f, &supportSize);
488: DMPlexGetSupport(dm, f, &support);
489: if (supportSize == 2) {
490: PetscInt numChildren;
492: DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
493: if (!numChildren) {
494: adj[tmp[support[0]-cStart]++] = support[1];
495: adj[tmp[support[1]-cStart]++] = support[0];
496: }
497: }
498: }
499: if (PetscDefined(USE_DEBUG)) {
500: for (c = 0; c < cEnd-cStart; ++c) if (tmp[c] != off[c+1]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Offset %d != %d for cell %d", tmp[c], off[c], c+cStart);
501: }
502: PetscFree(tmp);
503: }
504: if (numVertices) *numVertices = numCells;
505: if (offsets) *offsets = off;
506: if (adjacency) *adjacency = adj;
507: return(0);
508: }
509: /* Setup face recognition */
510: if (faceDepth == 1) {
511: 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 */
513: for (c = cStart; c < cEnd; ++c) {
514: PetscInt corners;
516: DMPlexGetConeSize(dm, c, &corners);
517: if (!cornersSeen[corners]) {
518: PetscInt nFV;
520: if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
521: cornersSeen[corners] = 1;
523: DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);
525: numFaceVertices[numFaceCases++] = nFV;
526: }
527: }
528: }
529: PetscCalloc1(numCells+1, &off);
530: /* Count neighboring cells */
531: for (cell = cStart; cell < cEnd; ++cell) {
532: PetscInt numNeighbors = PETSC_DETERMINE, n;
534: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
535: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
536: for (n = 0; n < numNeighbors; ++n) {
537: PetscInt cellPair[2];
538: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
539: PetscInt meetSize = 0;
540: const PetscInt *meet = NULL;
542: cellPair[0] = cell; cellPair[1] = neighborCells[n];
543: if (cellPair[0] == cellPair[1]) continue;
544: if (!found) {
545: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
546: if (meetSize) {
547: PetscInt f;
549: for (f = 0; f < numFaceCases; ++f) {
550: if (numFaceVertices[f] == meetSize) {
551: found = PETSC_TRUE;
552: break;
553: }
554: }
555: }
556: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
557: }
558: if (found) ++off[cell-cStart+1];
559: }
560: }
561: /* Prefix sum */
562: for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1];
564: if (adjacency) {
565: PetscMalloc1(off[numCells], &adj);
566: /* Get neighboring cells */
567: for (cell = cStart; cell < cEnd; ++cell) {
568: PetscInt numNeighbors = PETSC_DETERMINE, n;
569: PetscInt cellOffset = 0;
571: DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &numNeighbors, &neighborCells);
572: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
573: for (n = 0; n < numNeighbors; ++n) {
574: PetscInt cellPair[2];
575: PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE;
576: PetscInt meetSize = 0;
577: const PetscInt *meet = NULL;
579: cellPair[0] = cell; cellPair[1] = neighborCells[n];
580: if (cellPair[0] == cellPair[1]) continue;
581: if (!found) {
582: DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);
583: if (meetSize) {
584: PetscInt f;
586: for (f = 0; f < numFaceCases; ++f) {
587: if (numFaceVertices[f] == meetSize) {
588: found = PETSC_TRUE;
589: break;
590: }
591: }
592: }
593: DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);
594: }
595: if (found) {
596: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
597: ++cellOffset;
598: }
599: }
600: }
601: }
602: PetscFree(neighborCells);
603: if (numVertices) *numVertices = numCells;
604: if (offsets) *offsets = off;
605: if (adjacency) *adjacency = adj;
606: return(0);
607: }
609: /*@
610: PetscPartitionerDMPlexPartition - Create a non-overlapping partition of the cells in the mesh
612: Collective on PetscPartitioner
614: Input Parameters:
615: + part - The PetscPartitioner
616: . targetSection - The PetscSection describing the absolute weight of each partition (can be NULL)
617: - dm - The mesh DM
619: Output Parameters:
620: + partSection - The PetscSection giving the division of points by partition
621: - partition - The list of points by partition
623: Notes:
624: If the DM has a local section associated, each point to be partitioned will be weighted by the total number of dofs identified
625: by the section in the transitive closure of the point.
627: Level: developer
629: .seealso DMPlexDistribute(), PetscPartitionerCreate(), PetscSectionCreate(), PetscSectionSetChart(), PetscPartitionerPartition()
630: @*/
631: PetscErrorCode PetscPartitionerDMPlexPartition(PetscPartitioner part, DM dm, PetscSection targetSection, PetscSection partSection, IS *partition)
632: {
633: PetscMPIInt size;
634: PetscBool isplex;
636: PetscSection vertSection = NULL;
644: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isplex);
645: if (!isplex) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not for type %s",((PetscObject)dm)->type_name);
646: MPI_Comm_size(PetscObjectComm((PetscObject) part), &size);
647: if (size == 1) {
648: PetscInt *points;
649: PetscInt cStart, cEnd, c;
651: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
652: PetscSectionReset(partSection);
653: PetscSectionSetChart(partSection, 0, size);
654: PetscSectionSetDof(partSection, 0, cEnd-cStart);
655: PetscSectionSetUp(partSection);
656: PetscMalloc1(cEnd-cStart, &points);
657: for (c = cStart; c < cEnd; ++c) points[c] = c;
658: ISCreateGeneral(PetscObjectComm((PetscObject) part), cEnd-cStart, points, PETSC_OWN_POINTER, partition);
659: return(0);
660: }
661: if (part->height == 0) {
662: PetscInt numVertices = 0;
663: PetscInt *start = NULL;
664: PetscInt *adjacency = NULL;
665: IS globalNumbering;
667: if (!part->noGraph || part->viewGraph) {
668: DMPlexCreatePartitionerGraph(dm, part->height, &numVertices, &start, &adjacency, &globalNumbering);
669: } else { /* only compute the number of owned local vertices */
670: const PetscInt *idxs;
671: PetscInt p, pStart, pEnd;
673: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
674: DMPlexCreateNumbering_Plex(dm, pStart, pEnd, 0, NULL, dm->sf, &globalNumbering);
675: ISGetIndices(globalNumbering, &idxs);
676: for (p = 0; p < pEnd - pStart; p++) numVertices += idxs[p] < 0 ? 0 : 1;
677: ISRestoreIndices(globalNumbering, &idxs);
678: }
679: if (part->usevwgt) {
680: PetscSection section = dm->localSection, clSection = NULL;
681: IS clPoints = NULL;
682: const PetscInt *gid,*clIdx;
683: PetscInt v, p, pStart, pEnd;
685: /* 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) */
686: /* We do this only if the local section has been set */
687: if (section) {
688: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, NULL);
689: if (!clSection) {
690: DMPlexCreateClosureIndex(dm,NULL);
691: }
692: PetscSectionGetClosureIndex(section, (PetscObject)dm, &clSection, &clPoints);
693: ISGetIndices(clPoints,&clIdx);
694: }
695: DMPlexGetHeightStratum(dm, part->height, &pStart, &pEnd);
696: PetscSectionCreate(PETSC_COMM_SELF, &vertSection);
697: PetscSectionSetChart(vertSection, 0, numVertices);
698: if (globalNumbering) {
699: ISGetIndices(globalNumbering,&gid);
700: } else gid = NULL;
701: for (p = pStart, v = 0; p < pEnd; ++p) {
702: PetscInt dof = 1;
704: /* skip cells in the overlap */
705: if (gid && gid[p-pStart] < 0) continue;
707: if (section) {
708: PetscInt cl, clSize, clOff;
710: dof = 0;
711: PetscSectionGetDof(clSection, p, &clSize);
712: PetscSectionGetOffset(clSection, p, &clOff);
713: for (cl = 0; cl < clSize; cl+=2) {
714: PetscInt clDof, clPoint = clIdx[clOff + cl]; /* odd indices are reserved for orientations */
716: PetscSectionGetDof(section, clPoint, &clDof);
717: dof += clDof;
718: }
719: }
720: if (!dof) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of dofs for point %D in the local section should be positive",p);
721: PetscSectionSetDof(vertSection, v, dof);
722: v++;
723: }
724: if (globalNumbering) {
725: ISRestoreIndices(globalNumbering,&gid);
726: }
727: if (clPoints) {
728: ISRestoreIndices(clPoints,&clIdx);
729: }
730: PetscSectionSetUp(vertSection);
731: }
732: PetscPartitionerPartition(part, size, numVertices, start, adjacency, vertSection, targetSection, partSection, partition);
733: PetscFree(start);
734: PetscFree(adjacency);
735: if (globalNumbering) { /* partition is wrt global unique numbering: change this to be wrt local numbering */
736: const PetscInt *globalNum;
737: const PetscInt *partIdx;
738: PetscInt *map, cStart, cEnd;
739: PetscInt *adjusted, i, localSize, offset;
740: IS newPartition;
742: ISGetLocalSize(*partition,&localSize);
743: PetscMalloc1(localSize,&adjusted);
744: ISGetIndices(globalNumbering,&globalNum);
745: ISGetIndices(*partition,&partIdx);
746: PetscMalloc1(localSize,&map);
747: DMPlexGetHeightStratum(dm, part->height, &cStart, &cEnd);
748: for (i = cStart, offset = 0; i < cEnd; i++) {
749: if (globalNum[i - cStart] >= 0) map[offset++] = i;
750: }
751: for (i = 0; i < localSize; i++) {
752: adjusted[i] = map[partIdx[i]];
753: }
754: PetscFree(map);
755: ISRestoreIndices(*partition,&partIdx);
756: ISRestoreIndices(globalNumbering,&globalNum);
757: ISCreateGeneral(PETSC_COMM_SELF,localSize,adjusted,PETSC_OWN_POINTER,&newPartition);
758: ISDestroy(&globalNumbering);
759: ISDestroy(partition);
760: *partition = newPartition;
761: }
762: } else SETERRQ1(PetscObjectComm((PetscObject) part), PETSC_ERR_ARG_OUTOFRANGE, "Invalid height %D for points to partition", part->height);
763: PetscSectionDestroy(&vertSection);
764: return(0);
765: }
767: /*@
768: DMPlexGetPartitioner - Get the mesh partitioner
770: Not collective
772: Input Parameter:
773: . dm - The DM
775: Output Parameter:
776: . part - The PetscPartitioner
778: Level: developer
780: Note: This gets a borrowed reference, so the user should not destroy this PetscPartitioner.
782: .seealso DMPlexDistribute(), DMPlexSetPartitioner(), PetscPartitionerDMPlexPartition(), PetscPartitionerCreate()
783: @*/
784: PetscErrorCode DMPlexGetPartitioner(DM dm, PetscPartitioner *part)
785: {
786: DM_Plex *mesh = (DM_Plex *) dm->data;
791: *part = mesh->partitioner;
792: return(0);
793: }
795: /*@
796: DMPlexSetPartitioner - Set the mesh partitioner
798: logically collective on DM
800: Input Parameters:
801: + dm - The DM
802: - part - The partitioner
804: Level: developer
806: Note: Any existing PetscPartitioner will be destroyed.
808: .seealso DMPlexDistribute(), DMPlexGetPartitioner(), PetscPartitionerCreate()
809: @*/
810: PetscErrorCode DMPlexSetPartitioner(DM dm, PetscPartitioner part)
811: {
812: DM_Plex *mesh = (DM_Plex *) dm->data;
818: PetscObjectReference((PetscObject)part);
819: PetscPartitionerDestroy(&mesh->partitioner);
820: mesh->partitioner = part;
821: return(0);
822: }
824: static PetscErrorCode DMPlexAddClosure_Private(DM dm, PetscHSetI ht, PetscInt point)
825: {
826: const PetscInt *cone;
827: PetscInt coneSize, c;
828: PetscBool missing;
832: PetscHSetIQueryAdd(ht, point, &missing);
833: if (missing) {
834: DMPlexGetCone(dm, point, &cone);
835: DMPlexGetConeSize(dm, point, &coneSize);
836: for (c = 0; c < coneSize; c++) {
837: DMPlexAddClosure_Private(dm, ht, cone[c]);
838: }
839: }
840: return(0);
841: }
843: PETSC_UNUSED static PetscErrorCode DMPlexAddClosure_Tree(DM dm, PetscHSetI ht, PetscInt point, PetscBool up, PetscBool down)
844: {
848: if (up) {
849: PetscInt parent;
851: DMPlexGetTreeParent(dm,point,&parent,NULL);
852: if (parent != point) {
853: PetscInt closureSize, *closure = NULL, i;
855: DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
856: for (i = 0; i < closureSize; i++) {
857: PetscInt cpoint = closure[2*i];
859: PetscHSetIAdd(ht, cpoint);
860: DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_TRUE,PETSC_FALSE);
861: }
862: DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);
863: }
864: }
865: if (down) {
866: PetscInt numChildren;
867: const PetscInt *children;
869: DMPlexGetTreeChildren(dm,point,&numChildren,&children);
870: if (numChildren) {
871: PetscInt i;
873: for (i = 0; i < numChildren; i++) {
874: PetscInt cpoint = children[i];
876: PetscHSetIAdd(ht, cpoint);
877: DMPlexAddClosure_Tree(dm,ht,cpoint,PETSC_FALSE,PETSC_TRUE);
878: }
879: }
880: }
881: return(0);
882: }
884: static PetscErrorCode DMPlexAddClosureTree_Up_Private(DM dm, PetscHSetI ht, PetscInt point)
885: {
886: PetscInt parent;
890: DMPlexGetTreeParent(dm, point, &parent,NULL);
891: if (point != parent) {
892: const PetscInt *cone;
893: PetscInt coneSize, c;
895: DMPlexAddClosureTree_Up_Private(dm, ht, parent);
896: DMPlexAddClosure_Private(dm, ht, parent);
897: DMPlexGetCone(dm, parent, &cone);
898: DMPlexGetConeSize(dm, parent, &coneSize);
899: for (c = 0; c < coneSize; c++) {
900: const PetscInt cp = cone[c];
902: DMPlexAddClosureTree_Up_Private(dm, ht, cp);
903: }
904: }
905: return(0);
906: }
908: static PetscErrorCode DMPlexAddClosureTree_Down_Private(DM dm, PetscHSetI ht, PetscInt point)
909: {
910: PetscInt i, numChildren;
911: const PetscInt *children;
915: DMPlexGetTreeChildren(dm, point, &numChildren, &children);
916: for (i = 0; i < numChildren; i++) {
917: PetscHSetIAdd(ht, children[i]);
918: }
919: return(0);
920: }
922: static PetscErrorCode DMPlexAddClosureTree_Private(DM dm, PetscHSetI ht, PetscInt point)
923: {
924: const PetscInt *cone;
925: PetscInt coneSize, c;
929: PetscHSetIAdd(ht, point);
930: DMPlexAddClosureTree_Up_Private(dm, ht, point);
931: DMPlexAddClosureTree_Down_Private(dm, ht, point);
932: DMPlexGetCone(dm, point, &cone);
933: DMPlexGetConeSize(dm, point, &coneSize);
934: for (c = 0; c < coneSize; c++) {
935: DMPlexAddClosureTree_Private(dm, ht, cone[c]);
936: }
937: return(0);
938: }
940: PetscErrorCode DMPlexClosurePoints_Private(DM dm, PetscInt numPoints, const PetscInt points[], IS *closureIS)
941: {
942: DM_Plex *mesh = (DM_Plex *)dm->data;
943: const PetscBool hasTree = (mesh->parentSection || mesh->childSection) ? PETSC_TRUE : PETSC_FALSE;
944: PetscInt nelems, *elems, off = 0, p;
945: PetscHSetI ht;
946: PetscErrorCode ierr;
949: PetscHSetICreate(&ht);
950: PetscHSetIResize(ht, numPoints*16);
951: if (!hasTree) {
952: for (p = 0; p < numPoints; ++p) {
953: DMPlexAddClosure_Private(dm, ht, points[p]);
954: }
955: } else {
956: #if 1
957: for (p = 0; p < numPoints; ++p) {
958: DMPlexAddClosureTree_Private(dm, ht, points[p]);
959: }
960: #else
961: PetscInt *closure = NULL, closureSize, c;
962: for (p = 0; p < numPoints; ++p) {
963: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
964: for (c = 0; c < closureSize*2; c += 2) {
965: PetscHSetIAdd(ht, closure[c]);
966: if (hasTree) {DMPlexAddClosure_Tree(dm, ht, closure[c], PETSC_TRUE, PETSC_TRUE);}
967: }
968: }
969: if (closure) {DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);}
970: #endif
971: }
972: PetscHSetIGetSize(ht, &nelems);
973: PetscMalloc1(nelems, &elems);
974: PetscHSetIGetElems(ht, &off, elems);
975: PetscHSetIDestroy(&ht);
976: PetscSortInt(nelems, elems);
977: ISCreateGeneral(PETSC_COMM_SELF, nelems, elems, PETSC_OWN_POINTER, closureIS);
978: return(0);
979: }
981: /*@
982: DMPlexPartitionLabelClosure - Add the closure of all points to the partition label
984: Input Parameters:
985: + dm - The DM
986: - label - DMLabel assinging ranks to remote roots
988: Level: developer
990: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
991: @*/
992: PetscErrorCode DMPlexPartitionLabelClosure(DM dm, DMLabel label)
993: {
994: IS rankIS, pointIS, closureIS;
995: const PetscInt *ranks, *points;
996: PetscInt numRanks, numPoints, r;
997: PetscErrorCode ierr;
1000: DMLabelGetValueIS(label, &rankIS);
1001: ISGetLocalSize(rankIS, &numRanks);
1002: ISGetIndices(rankIS, &ranks);
1003: for (r = 0; r < numRanks; ++r) {
1004: const PetscInt rank = ranks[r];
1005: DMLabelGetStratumIS(label, rank, &pointIS);
1006: ISGetLocalSize(pointIS, &numPoints);
1007: ISGetIndices(pointIS, &points);
1008: DMPlexClosurePoints_Private(dm, numPoints, points, &closureIS);
1009: ISRestoreIndices(pointIS, &points);
1010: ISDestroy(&pointIS);
1011: DMLabelSetStratumIS(label, rank, closureIS);
1012: ISDestroy(&closureIS);
1013: }
1014: ISRestoreIndices(rankIS, &ranks);
1015: ISDestroy(&rankIS);
1016: return(0);
1017: }
1019: /*@
1020: DMPlexPartitionLabelAdjacency - Add one level of adjacent points to the partition label
1022: Input Parameters:
1023: + dm - The DM
1024: - label - DMLabel assinging ranks to remote roots
1026: Level: developer
1028: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1029: @*/
1030: PetscErrorCode DMPlexPartitionLabelAdjacency(DM dm, DMLabel label)
1031: {
1032: IS rankIS, pointIS;
1033: const PetscInt *ranks, *points;
1034: PetscInt numRanks, numPoints, r, p, a, adjSize;
1035: PetscInt *adj = NULL;
1036: PetscErrorCode ierr;
1039: DMLabelGetValueIS(label, &rankIS);
1040: ISGetLocalSize(rankIS, &numRanks);
1041: ISGetIndices(rankIS, &ranks);
1042: for (r = 0; r < numRanks; ++r) {
1043: const PetscInt rank = ranks[r];
1045: DMLabelGetStratumIS(label, rank, &pointIS);
1046: ISGetLocalSize(pointIS, &numPoints);
1047: ISGetIndices(pointIS, &points);
1048: for (p = 0; p < numPoints; ++p) {
1049: adjSize = PETSC_DETERMINE;
1050: DMPlexGetAdjacency(dm, points[p], &adjSize, &adj);
1051: for (a = 0; a < adjSize; ++a) {DMLabelSetValue(label, adj[a], rank);}
1052: }
1053: ISRestoreIndices(pointIS, &points);
1054: ISDestroy(&pointIS);
1055: }
1056: ISRestoreIndices(rankIS, &ranks);
1057: ISDestroy(&rankIS);
1058: PetscFree(adj);
1059: return(0);
1060: }
1062: /*@
1063: DMPlexPartitionLabelPropagate - Propagate points in a partition label over the point SF
1065: Input Parameters:
1066: + dm - The DM
1067: - label - DMLabel assinging ranks to remote roots
1069: Level: developer
1071: Note: This is required when generating multi-level overlaps to capture
1072: overlap points from non-neighbouring partitions.
1074: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1075: @*/
1076: PetscErrorCode DMPlexPartitionLabelPropagate(DM dm, DMLabel label)
1077: {
1078: MPI_Comm comm;
1079: PetscMPIInt rank;
1080: PetscSF sfPoint;
1081: DMLabel lblRoots, lblLeaves;
1082: IS rankIS, pointIS;
1083: const PetscInt *ranks;
1084: PetscInt numRanks, r;
1085: PetscErrorCode ierr;
1088: PetscObjectGetComm((PetscObject) dm, &comm);
1089: MPI_Comm_rank(comm, &rank);
1090: DMGetPointSF(dm, &sfPoint);
1091: /* Pull point contributions from remote leaves into local roots */
1092: DMLabelGather(label, sfPoint, &lblLeaves);
1093: DMLabelGetValueIS(lblLeaves, &rankIS);
1094: ISGetLocalSize(rankIS, &numRanks);
1095: ISGetIndices(rankIS, &ranks);
1096: for (r = 0; r < numRanks; ++r) {
1097: const PetscInt remoteRank = ranks[r];
1098: if (remoteRank == rank) continue;
1099: DMLabelGetStratumIS(lblLeaves, remoteRank, &pointIS);
1100: DMLabelInsertIS(label, pointIS, remoteRank);
1101: ISDestroy(&pointIS);
1102: }
1103: ISRestoreIndices(rankIS, &ranks);
1104: ISDestroy(&rankIS);
1105: DMLabelDestroy(&lblLeaves);
1106: /* Push point contributions from roots into remote leaves */
1107: DMLabelDistribute(label, sfPoint, &lblRoots);
1108: DMLabelGetValueIS(lblRoots, &rankIS);
1109: ISGetLocalSize(rankIS, &numRanks);
1110: ISGetIndices(rankIS, &ranks);
1111: for (r = 0; r < numRanks; ++r) {
1112: const PetscInt remoteRank = ranks[r];
1113: if (remoteRank == rank) continue;
1114: DMLabelGetStratumIS(lblRoots, remoteRank, &pointIS);
1115: DMLabelInsertIS(label, pointIS, remoteRank);
1116: ISDestroy(&pointIS);
1117: }
1118: ISRestoreIndices(rankIS, &ranks);
1119: ISDestroy(&rankIS);
1120: DMLabelDestroy(&lblRoots);
1121: return(0);
1122: }
1124: /*@
1125: DMPlexPartitionLabelInvert - Create a partition label of remote roots from a local root label
1127: Input Parameters:
1128: + dm - The DM
1129: . rootLabel - DMLabel assinging ranks to local roots
1130: - processSF - A star forest mapping into the local index on each remote rank
1132: Output Parameter:
1133: . leafLabel - DMLabel assinging ranks to remote roots
1135: Note: The rootLabel defines a send pattern by mapping local points to remote target ranks. The
1136: resulting leafLabel is a receiver mapping of remote roots to their parent rank.
1138: Level: developer
1140: .seealso: DMPlexPartitionLabelCreateSF(), DMPlexDistribute(), DMPlexCreateOverlap()
1141: @*/
1142: PetscErrorCode DMPlexPartitionLabelInvert(DM dm, DMLabel rootLabel, PetscSF processSF, DMLabel leafLabel)
1143: {
1144: MPI_Comm comm;
1145: PetscMPIInt rank, size, r;
1146: PetscInt p, n, numNeighbors, numPoints, dof, off, rootSize, l, nleaves, leafSize;
1147: PetscSF sfPoint;
1148: PetscSection rootSection;
1149: PetscSFNode *rootPoints, *leafPoints;
1150: const PetscSFNode *remote;
1151: const PetscInt *local, *neighbors;
1152: IS valueIS;
1153: PetscBool mpiOverflow = PETSC_FALSE;
1154: PetscErrorCode ierr;
1157: PetscLogEventBegin(DMPLEX_PartLabelInvert,dm,0,0,0);
1158: PetscObjectGetComm((PetscObject) dm, &comm);
1159: MPI_Comm_rank(comm, &rank);
1160: MPI_Comm_size(comm, &size);
1161: DMGetPointSF(dm, &sfPoint);
1163: /* Convert to (point, rank) and use actual owners */
1164: PetscSectionCreate(comm, &rootSection);
1165: PetscSectionSetChart(rootSection, 0, size);
1166: DMLabelGetValueIS(rootLabel, &valueIS);
1167: ISGetLocalSize(valueIS, &numNeighbors);
1168: ISGetIndices(valueIS, &neighbors);
1169: for (n = 0; n < numNeighbors; ++n) {
1170: DMLabelGetStratumSize(rootLabel, neighbors[n], &numPoints);
1171: PetscSectionAddDof(rootSection, neighbors[n], numPoints);
1172: }
1173: PetscSectionSetUp(rootSection);
1174: PetscSectionGetStorageSize(rootSection, &rootSize);
1175: PetscMalloc1(rootSize, &rootPoints);
1176: PetscSFGetGraph(sfPoint, NULL, &nleaves, &local, &remote);
1177: for (n = 0; n < numNeighbors; ++n) {
1178: IS pointIS;
1179: const PetscInt *points;
1181: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1182: DMLabelGetStratumIS(rootLabel, neighbors[n], &pointIS);
1183: ISGetLocalSize(pointIS, &numPoints);
1184: ISGetIndices(pointIS, &points);
1185: for (p = 0; p < numPoints; ++p) {
1186: if (local) {PetscFindInt(points[p], nleaves, local, &l);}
1187: else {l = -1;}
1188: if (l >= 0) {rootPoints[off+p] = remote[l];}
1189: else {rootPoints[off+p].index = points[p]; rootPoints[off+p].rank = rank;}
1190: }
1191: ISRestoreIndices(pointIS, &points);
1192: ISDestroy(&pointIS);
1193: }
1195: /* Try to communicate overlap using All-to-All */
1196: if (!processSF) {
1197: PetscInt64 counter = 0;
1198: PetscBool locOverflow = PETSC_FALSE;
1199: PetscMPIInt *scounts, *sdispls, *rcounts, *rdispls;
1201: PetscCalloc4(size, &scounts, size, &sdispls, size, &rcounts, size, &rdispls);
1202: for (n = 0; n < numNeighbors; ++n) {
1203: PetscSectionGetDof(rootSection, neighbors[n], &dof);
1204: PetscSectionGetOffset(rootSection, neighbors[n], &off);
1205: #if defined(PETSC_USE_64BIT_INDICES)
1206: if (dof > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
1207: if (off > PETSC_MPI_INT_MAX) {locOverflow = PETSC_TRUE; break;}
1208: #endif
1209: scounts[neighbors[n]] = (PetscMPIInt) dof;
1210: sdispls[neighbors[n]] = (PetscMPIInt) off;
1211: }
1212: MPI_Alltoall(scounts, 1, MPI_INT, rcounts, 1, MPI_INT, comm);
1213: for (r = 0; r < size; ++r) { rdispls[r] = (int)counter; counter += rcounts[r]; }
1214: if (counter > PETSC_MPI_INT_MAX) locOverflow = PETSC_TRUE;
1215: MPI_Allreduce(&locOverflow, &mpiOverflow, 1, MPIU_BOOL, MPI_LOR, comm);
1216: if (!mpiOverflow) {
1217: PetscInfo(dm,"Using Alltoallv for mesh distribution\n");
1218: leafSize = (PetscInt) counter;
1219: PetscMalloc1(leafSize, &leafPoints);
1220: MPI_Alltoallv(rootPoints, scounts, sdispls, MPIU_2INT, leafPoints, rcounts, rdispls, MPIU_2INT, comm);
1221: }
1222: PetscFree4(scounts, sdispls, rcounts, rdispls);
1223: }
1225: /* Communicate overlap using process star forest */
1226: if (processSF || mpiOverflow) {
1227: PetscSF procSF;
1228: PetscSection leafSection;
1230: if (processSF) {
1231: PetscInfo(dm,"Using processSF for mesh distribution\n");
1232: PetscObjectReference((PetscObject)processSF);
1233: procSF = processSF;
1234: } else {
1235: PetscInfo(dm,"Using processSF for mesh distribution (MPI overflow)\n");
1236: PetscSFCreate(comm,&procSF);
1237: PetscSFSetGraphWithPattern(procSF,NULL,PETSCSF_PATTERN_ALLTOALL);
1238: }
1240: PetscSectionCreate(PetscObjectComm((PetscObject)dm), &leafSection);
1241: DMPlexDistributeData(dm, procSF, rootSection, MPIU_2INT, rootPoints, leafSection, (void**) &leafPoints);
1242: PetscSectionGetStorageSize(leafSection, &leafSize);
1243: PetscSectionDestroy(&leafSection);
1244: PetscSFDestroy(&procSF);
1245: }
1247: for (p = 0; p < leafSize; p++) {
1248: DMLabelSetValue(leafLabel, leafPoints[p].index, leafPoints[p].rank);
1249: }
1251: ISRestoreIndices(valueIS, &neighbors);
1252: ISDestroy(&valueIS);
1253: PetscSectionDestroy(&rootSection);
1254: PetscFree(rootPoints);
1255: PetscFree(leafPoints);
1256: PetscLogEventEnd(DMPLEX_PartLabelInvert,dm,0,0,0);
1257: return(0);
1258: }
1260: /*@
1261: DMPlexPartitionLabelCreateSF - Create a star forest from a label that assigns ranks to points
1263: Input Parameters:
1264: + dm - The DM
1265: - label - DMLabel assinging ranks to remote roots
1267: Output Parameter:
1268: . sf - The star forest communication context encapsulating the defined mapping
1270: Note: The incoming label is a receiver mapping of remote points to their parent rank.
1272: Level: developer
1274: .seealso: DMPlexDistribute(), DMPlexCreateOverlap()
1275: @*/
1276: PetscErrorCode DMPlexPartitionLabelCreateSF(DM dm, DMLabel label, PetscSF *sf)
1277: {
1278: PetscMPIInt rank;
1279: PetscInt n, numRemote, p, numPoints, pStart, pEnd, idx = 0, nNeighbors;
1280: PetscSFNode *remotePoints;
1281: IS remoteRootIS, neighborsIS;
1282: const PetscInt *remoteRoots, *neighbors;
1286: PetscLogEventBegin(DMPLEX_PartLabelCreateSF,dm,0,0,0);
1287: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1289: DMLabelGetValueIS(label, &neighborsIS);
1290: #if 0
1291: {
1292: IS is;
1293: ISDuplicate(neighborsIS, &is);
1294: ISSort(is);
1295: ISDestroy(&neighborsIS);
1296: neighborsIS = is;
1297: }
1298: #endif
1299: ISGetLocalSize(neighborsIS, &nNeighbors);
1300: ISGetIndices(neighborsIS, &neighbors);
1301: for (numRemote = 0, n = 0; n < nNeighbors; ++n) {
1302: DMLabelGetStratumSize(label, neighbors[n], &numPoints);
1303: numRemote += numPoints;
1304: }
1305: PetscMalloc1(numRemote, &remotePoints);
1306: /* Put owned points first */
1307: DMLabelGetStratumSize(label, rank, &numPoints);
1308: if (numPoints > 0) {
1309: DMLabelGetStratumIS(label, rank, &remoteRootIS);
1310: ISGetIndices(remoteRootIS, &remoteRoots);
1311: for (p = 0; p < numPoints; p++) {
1312: remotePoints[idx].index = remoteRoots[p];
1313: remotePoints[idx].rank = rank;
1314: idx++;
1315: }
1316: ISRestoreIndices(remoteRootIS, &remoteRoots);
1317: ISDestroy(&remoteRootIS);
1318: }
1319: /* Now add remote points */
1320: for (n = 0; n < nNeighbors; ++n) {
1321: const PetscInt nn = neighbors[n];
1323: DMLabelGetStratumSize(label, nn, &numPoints);
1324: if (nn == rank || numPoints <= 0) continue;
1325: DMLabelGetStratumIS(label, nn, &remoteRootIS);
1326: ISGetIndices(remoteRootIS, &remoteRoots);
1327: for (p = 0; p < numPoints; p++) {
1328: remotePoints[idx].index = remoteRoots[p];
1329: remotePoints[idx].rank = nn;
1330: idx++;
1331: }
1332: ISRestoreIndices(remoteRootIS, &remoteRoots);
1333: ISDestroy(&remoteRootIS);
1334: }
1335: PetscSFCreate(PetscObjectComm((PetscObject) dm), sf);
1336: DMPlexGetChart(dm, &pStart, &pEnd);
1337: PetscSFSetGraph(*sf, pEnd-pStart, numRemote, NULL, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
1338: PetscSFSetUp(*sf);
1339: ISDestroy(&neighborsIS);
1340: PetscLogEventEnd(DMPLEX_PartLabelCreateSF,dm,0,0,0);
1341: return(0);
1342: }
1344: #if defined(PETSC_HAVE_PARMETIS)
1345: #include <parmetis.h>
1346: #endif
1348: /* The two functions below are used by DMPlexRebalanceSharedPoints which errors
1349: * when PETSc is built without ParMETIS. To avoid -Wunused-function, we take
1350: * them out in that case. */
1351: #if defined(PETSC_HAVE_PARMETIS)
1352: /*@C
1354: DMPlexRewriteSF - Rewrites the ownership of the SF of a DM (in place).
1356: Input parameters:
1357: + dm - The DMPlex object.
1358: . n - The number of points.
1359: . pointsToRewrite - The points in the SF whose ownership will change.
1360: . targetOwners - New owner for each element in pointsToRewrite.
1361: - degrees - Degrees of the points in the SF as obtained by PetscSFComputeDegreeBegin/PetscSFComputeDegreeEnd.
1363: Level: developer
1365: @*/
1366: static PetscErrorCode DMPlexRewriteSF(DM dm, PetscInt n, PetscInt *pointsToRewrite, PetscInt *targetOwners, const PetscInt *degrees)
1367: {
1368: PetscInt ierr, pStart, pEnd, i, j, counter, leafCounter, sumDegrees, nroots, nleafs;
1369: PetscInt *cumSumDegrees, *newOwners, *newNumbers, *rankOnLeafs, *locationsOfLeafs, *remoteLocalPointOfLeafs, *points, *leafsNew;
1370: PetscSFNode *leafLocationsNew;
1371: const PetscSFNode *iremote;
1372: const PetscInt *ilocal;
1373: PetscBool *isLeaf;
1374: PetscSF sf;
1375: MPI_Comm comm;
1376: PetscMPIInt rank, size;
1379: PetscObjectGetComm((PetscObject) dm, &comm);
1380: MPI_Comm_rank(comm, &rank);
1381: MPI_Comm_size(comm, &size);
1382: DMPlexGetChart(dm, &pStart, &pEnd);
1384: DMGetPointSF(dm, &sf);
1385: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1386: PetscMalloc1(pEnd-pStart, &isLeaf);
1387: for (i=0; i<pEnd-pStart; i++) {
1388: isLeaf[i] = PETSC_FALSE;
1389: }
1390: for (i=0; i<nleafs; i++) {
1391: isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
1392: }
1394: PetscMalloc1(pEnd-pStart+1, &cumSumDegrees);
1395: cumSumDegrees[0] = 0;
1396: for (i=1; i<=pEnd-pStart; i++) {
1397: cumSumDegrees[i] = cumSumDegrees[i-1] + degrees[i-1];
1398: }
1399: sumDegrees = cumSumDegrees[pEnd-pStart];
1400: /* get the location of my leafs (we have sumDegrees many leafs pointing at our roots) */
1402: PetscMalloc1(sumDegrees, &locationsOfLeafs);
1403: PetscMalloc1(pEnd-pStart, &rankOnLeafs);
1404: for (i=0; i<pEnd-pStart; i++) {
1405: rankOnLeafs[i] = rank;
1406: }
1407: PetscSFGatherBegin(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1408: PetscSFGatherEnd(sf, MPIU_INT, rankOnLeafs, locationsOfLeafs);
1409: PetscFree(rankOnLeafs);
1411: /* get the remote local points of my leaves */
1412: PetscMalloc1(sumDegrees, &remoteLocalPointOfLeafs);
1413: PetscMalloc1(pEnd-pStart, &points);
1414: for (i=0; i<pEnd-pStart; i++) {
1415: points[i] = pStart+i;
1416: }
1417: PetscSFGatherBegin(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1418: PetscSFGatherEnd(sf, MPIU_INT, points, remoteLocalPointOfLeafs);
1419: PetscFree(points);
1420: /* Figure out the new owners of the vertices that are up for grabs and their numbers on the new owners */
1421: PetscMalloc1(pEnd-pStart, &newOwners);
1422: PetscMalloc1(pEnd-pStart, &newNumbers);
1423: for (i=0; i<pEnd-pStart; i++) {
1424: newOwners[i] = -1;
1425: newNumbers[i] = -1;
1426: }
1427: {
1428: PetscInt oldNumber, newNumber, oldOwner, newOwner;
1429: for (i=0; i<n; i++) {
1430: oldNumber = pointsToRewrite[i];
1431: newNumber = -1;
1432: oldOwner = rank;
1433: newOwner = targetOwners[i];
1434: if (oldOwner == newOwner) {
1435: newNumber = oldNumber;
1436: } else {
1437: for (j=0; j<degrees[oldNumber]; j++) {
1438: if (locationsOfLeafs[cumSumDegrees[oldNumber]+j] == newOwner) {
1439: newNumber = remoteLocalPointOfLeafs[cumSumDegrees[oldNumber]+j];
1440: break;
1441: }
1442: }
1443: }
1444: if (newNumber == -1) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Couldn't find the new owner of vertex.");
1446: newOwners[oldNumber] = newOwner;
1447: newNumbers[oldNumber] = newNumber;
1448: }
1449: }
1450: PetscFree(cumSumDegrees);
1451: PetscFree(locationsOfLeafs);
1452: PetscFree(remoteLocalPointOfLeafs);
1454: PetscSFBcastBegin(sf, MPIU_INT, newOwners, newOwners);
1455: PetscSFBcastEnd(sf, MPIU_INT, newOwners, newOwners);
1456: PetscSFBcastBegin(sf, MPIU_INT, newNumbers, newNumbers);
1457: PetscSFBcastEnd(sf, MPIU_INT, newNumbers, newNumbers);
1459: /* Now count how many leafs we have on each processor. */
1460: leafCounter=0;
1461: for (i=0; i<pEnd-pStart; i++) {
1462: if (newOwners[i] >= 0) {
1463: if (newOwners[i] != rank) {
1464: leafCounter++;
1465: }
1466: } else {
1467: if (isLeaf[i]) {
1468: leafCounter++;
1469: }
1470: }
1471: }
1473: /* Now set up the new sf by creating the leaf arrays */
1474: PetscMalloc1(leafCounter, &leafsNew);
1475: PetscMalloc1(leafCounter, &leafLocationsNew);
1477: leafCounter = 0;
1478: counter = 0;
1479: for (i=0; i<pEnd-pStart; i++) {
1480: if (newOwners[i] >= 0) {
1481: if (newOwners[i] != rank) {
1482: leafsNew[leafCounter] = i;
1483: leafLocationsNew[leafCounter].rank = newOwners[i];
1484: leafLocationsNew[leafCounter].index = newNumbers[i];
1485: leafCounter++;
1486: }
1487: } else {
1488: if (isLeaf[i]) {
1489: leafsNew[leafCounter] = i;
1490: leafLocationsNew[leafCounter].rank = iremote[counter].rank;
1491: leafLocationsNew[leafCounter].index = iremote[counter].index;
1492: leafCounter++;
1493: }
1494: }
1495: if (isLeaf[i]) {
1496: counter++;
1497: }
1498: }
1500: PetscSFSetGraph(sf, nroots, leafCounter, leafsNew, PETSC_OWN_POINTER, leafLocationsNew, PETSC_OWN_POINTER);
1501: PetscFree(newOwners);
1502: PetscFree(newNumbers);
1503: PetscFree(isLeaf);
1504: return(0);
1505: }
1507: static PetscErrorCode DMPlexViewDistribution(MPI_Comm comm, PetscInt n, PetscInt skip, PetscInt *vtxwgt, PetscInt *part, PetscViewer viewer)
1508: {
1509: PetscInt *distribution, min, max, sum, i, ierr;
1510: PetscMPIInt rank, size;
1512: MPI_Comm_size(comm, &size);
1513: MPI_Comm_rank(comm, &rank);
1514: PetscCalloc1(size, &distribution);
1515: for (i=0; i<n; i++) {
1516: if (part) distribution[part[i]] += vtxwgt[skip*i];
1517: else distribution[rank] += vtxwgt[skip*i];
1518: }
1519: MPI_Allreduce(MPI_IN_PLACE, distribution, size, MPIU_INT, MPI_SUM, comm);
1520: min = distribution[0];
1521: max = distribution[0];
1522: sum = distribution[0];
1523: for (i=1; i<size; i++) {
1524: if (distribution[i]<min) min=distribution[i];
1525: if (distribution[i]>max) max=distribution[i];
1526: sum += distribution[i];
1527: }
1528: PetscViewerASCIIPrintf(viewer, "Min: %D, Avg: %D, Max: %D, Balance: %f\n", min, sum/size, max, (max*1.*size)/sum);
1529: PetscFree(distribution);
1530: return(0);
1531: }
1533: #endif
1535: /*@
1536: DMPlexRebalanceSharedPoints - Redistribute points in the plex that are shared in order to achieve better balancing. This routine updates the PointSF of the DM inplace.
1538: Input parameters:
1539: + dm - The DMPlex object.
1540: . entityDepth - depth of the entity to balance (0 -> balance vertices).
1541: . useInitialGuess - whether to use the current distribution as initial guess (only used by ParMETIS).
1542: - parallel - whether to use ParMETIS and do the partition in parallel or whether to gather the graph onto a single process and use METIS.
1544: Output parameters:
1545: . success - whether the graph partitioning was successful or not. If not, try useInitialGuess=True and parallel=True.
1547: Level: intermediate
1549: @*/
1551: PetscErrorCode DMPlexRebalanceSharedPoints(DM dm, PetscInt entityDepth, PetscBool useInitialGuess, PetscBool parallel, PetscBool *success)
1552: {
1553: #if defined(PETSC_HAVE_PARMETIS)
1554: PetscSF sf;
1555: PetscInt ierr, i, j, idx, jdx;
1556: PetscInt eBegin, eEnd, nroots, nleafs, pStart, pEnd;
1557: const PetscInt *degrees, *ilocal;
1558: const PetscSFNode *iremote;
1559: PetscBool *toBalance, *isLeaf, *isExclusivelyOwned, *isNonExclusivelyOwned;
1560: PetscInt numExclusivelyOwned, numNonExclusivelyOwned;
1561: PetscMPIInt rank, size;
1562: PetscInt *globalNumbersOfLocalOwnedVertices, *leafGlobalNumbers;
1563: const PetscInt *cumSumVertices;
1564: PetscInt offset, counter;
1565: PetscInt lenadjncy;
1566: PetscInt *xadj, *adjncy, *vtxwgt;
1567: PetscInt lenxadj;
1568: PetscInt *adjwgt = NULL;
1569: PetscInt *part, *options;
1570: PetscInt nparts, wgtflag, numflag, ncon, edgecut;
1571: real_t *ubvec;
1572: PetscInt *firstVertices, *renumbering;
1573: PetscInt failed, failedGlobal;
1574: MPI_Comm comm;
1575: Mat A, Apre;
1576: const char *prefix = NULL;
1577: PetscViewer viewer;
1578: PetscViewerFormat format;
1579: PetscLayout layout;
1582: if (success) *success = PETSC_FALSE;
1583: PetscObjectGetComm((PetscObject) dm, &comm);
1584: MPI_Comm_rank(comm, &rank);
1585: MPI_Comm_size(comm, &size);
1586: if (size==1) return(0);
1588: PetscLogEventBegin(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1590: PetscOptionsGetViewer(comm,((PetscObject)dm)->options, prefix,"-dm_rebalance_partition_view",&viewer,&format,NULL);
1591: if (viewer) {
1592: PetscViewerPushFormat(viewer,format);
1593: }
1595: /* Figure out all points in the plex that we are interested in balancing. */
1596: DMPlexGetDepthStratum(dm, entityDepth, &eBegin, &eEnd);
1597: DMPlexGetChart(dm, &pStart, &pEnd);
1598: PetscMalloc1(pEnd-pStart, &toBalance);
1600: for (i=0; i<pEnd-pStart; i++) {
1601: toBalance[i] = (PetscBool)(i-pStart>=eBegin && i-pStart<eEnd);
1602: }
1604: /* There are three types of points:
1605: * exclusivelyOwned: points that are owned by this process and only seen by this process
1606: * nonExclusivelyOwned: points that are owned by this process but seen by at least another process
1607: * leaf: a point that is seen by this process but owned by a different process
1608: */
1609: DMGetPointSF(dm, &sf);
1610: PetscSFGetGraph(sf, &nroots, &nleafs, &ilocal, &iremote);
1611: PetscMalloc1(pEnd-pStart, &isLeaf);
1612: PetscMalloc1(pEnd-pStart, &isNonExclusivelyOwned);
1613: PetscMalloc1(pEnd-pStart, &isExclusivelyOwned);
1614: for (i=0; i<pEnd-pStart; i++) {
1615: isNonExclusivelyOwned[i] = PETSC_FALSE;
1616: isExclusivelyOwned[i] = PETSC_FALSE;
1617: isLeaf[i] = PETSC_FALSE;
1618: }
1620: /* start by marking all the leafs */
1621: for (i=0; i<nleafs; i++) {
1622: isLeaf[ilocal[i]-pStart] = PETSC_TRUE;
1623: }
1625: /* for an owned point, we can figure out whether another processor sees it or
1626: * not by calculating its degree */
1627: PetscSFComputeDegreeBegin(sf, °rees);
1628: PetscSFComputeDegreeEnd(sf, °rees);
1630: numExclusivelyOwned = 0;
1631: numNonExclusivelyOwned = 0;
1632: for (i=0; i<pEnd-pStart; i++) {
1633: if (toBalance[i]) {
1634: if (degrees[i] > 0) {
1635: isNonExclusivelyOwned[i] = PETSC_TRUE;
1636: numNonExclusivelyOwned += 1;
1637: } else {
1638: if (!isLeaf[i]) {
1639: isExclusivelyOwned[i] = PETSC_TRUE;
1640: numExclusivelyOwned += 1;
1641: }
1642: }
1643: }
1644: }
1646: /* We are going to build a graph with one vertex per core representing the
1647: * exclusively owned points and then one vertex per nonExclusively owned
1648: * point. */
1650: PetscLayoutCreate(comm, &layout);
1651: PetscLayoutSetLocalSize(layout, 1 + numNonExclusivelyOwned);
1652: PetscLayoutSetUp(layout);
1653: PetscLayoutGetRanges(layout, &cumSumVertices);
1655: PetscMalloc1(pEnd-pStart, &globalNumbersOfLocalOwnedVertices);
1656: for (i=0; i<pEnd-pStart; i++) {globalNumbersOfLocalOwnedVertices[i] = pStart - 1;}
1657: offset = cumSumVertices[rank];
1658: counter = 0;
1659: for (i=0; i<pEnd-pStart; i++) {
1660: if (toBalance[i]) {
1661: if (degrees[i] > 0) {
1662: globalNumbersOfLocalOwnedVertices[i] = counter + 1 + offset;
1663: counter++;
1664: }
1665: }
1666: }
1668: /* send the global numbers of vertices I own to the leafs so that they know to connect to it */
1669: PetscMalloc1(pEnd-pStart, &leafGlobalNumbers);
1670: PetscSFBcastBegin(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);
1671: PetscSFBcastEnd(sf, MPIU_INT, globalNumbersOfLocalOwnedVertices, leafGlobalNumbers);
1673: /* Now start building the data structure for ParMETIS */
1675: MatCreate(comm, &Apre);
1676: MatSetType(Apre, MATPREALLOCATOR);
1677: MatSetSizes(Apre, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
1678: MatSetUp(Apre);
1680: for (i=0; i<pEnd-pStart; i++) {
1681: if (toBalance[i]) {
1682: idx = cumSumVertices[rank];
1683: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1684: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1685: else continue;
1686: MatSetValue(Apre, idx, jdx, 1, INSERT_VALUES);
1687: MatSetValue(Apre, jdx, idx, 1, INSERT_VALUES);
1688: }
1689: }
1691: MatAssemblyBegin(Apre, MAT_FINAL_ASSEMBLY);
1692: MatAssemblyEnd(Apre, MAT_FINAL_ASSEMBLY);
1694: MatCreate(comm, &A);
1695: MatSetType(A, MATMPIAIJ);
1696: MatSetSizes(A, 1+numNonExclusivelyOwned, 1+numNonExclusivelyOwned, cumSumVertices[size], cumSumVertices[size]);
1697: MatPreallocatorPreallocate(Apre, PETSC_FALSE, A);
1698: MatDestroy(&Apre);
1700: for (i=0; i<pEnd-pStart; i++) {
1701: if (toBalance[i]) {
1702: idx = cumSumVertices[rank];
1703: if (isNonExclusivelyOwned[i]) jdx = globalNumbersOfLocalOwnedVertices[i];
1704: else if (isLeaf[i]) jdx = leafGlobalNumbers[i];
1705: else continue;
1706: MatSetValue(A, idx, jdx, 1, INSERT_VALUES);
1707: MatSetValue(A, jdx, idx, 1, INSERT_VALUES);
1708: }
1709: }
1711: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1712: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1713: PetscFree(leafGlobalNumbers);
1714: PetscFree(globalNumbersOfLocalOwnedVertices);
1716: nparts = size;
1717: wgtflag = 2;
1718: numflag = 0;
1719: ncon = 2;
1720: real_t *tpwgts;
1721: PetscMalloc1(ncon * nparts, &tpwgts);
1722: for (i=0; i<ncon*nparts; i++) {
1723: tpwgts[i] = 1./(nparts);
1724: }
1726: PetscMalloc1(ncon, &ubvec);
1727: ubvec[0] = 1.01;
1728: ubvec[1] = 1.01;
1729: lenadjncy = 0;
1730: for (i=0; i<1+numNonExclusivelyOwned; i++) {
1731: PetscInt temp=0;
1732: MatGetRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
1733: lenadjncy += temp;
1734: MatRestoreRow(A, cumSumVertices[rank] + i, &temp, NULL, NULL);
1735: }
1736: PetscMalloc1(lenadjncy, &adjncy);
1737: lenxadj = 2 + numNonExclusivelyOwned;
1738: PetscMalloc1(lenxadj, &xadj);
1739: xadj[0] = 0;
1740: counter = 0;
1741: for (i=0; i<1+numNonExclusivelyOwned; i++) {
1742: PetscInt temp=0;
1743: const PetscInt *cols;
1744: MatGetRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
1745: PetscArraycpy(&adjncy[counter], cols, temp);
1746: counter += temp;
1747: xadj[i+1] = counter;
1748: MatRestoreRow(A, cumSumVertices[rank] + i, &temp, &cols, NULL);
1749: }
1751: PetscMalloc1(cumSumVertices[rank+1]-cumSumVertices[rank], &part);
1752: PetscMalloc1(ncon*(1 + numNonExclusivelyOwned), &vtxwgt);
1753: vtxwgt[0] = numExclusivelyOwned;
1754: if (ncon>1) vtxwgt[1] = 1;
1755: for (i=0; i<numNonExclusivelyOwned; i++) {
1756: vtxwgt[ncon*(i+1)] = 1;
1757: if (ncon>1) vtxwgt[ncon*(i+1)+1] = 0;
1758: }
1760: if (viewer) {
1761: PetscViewerASCIIPrintf(viewer, "Attempt rebalancing of shared points of depth %D on interface of mesh distribution.\n", entityDepth);
1762: PetscViewerASCIIPrintf(viewer, "Size of generated auxiliary graph: %D\n", cumSumVertices[size]);
1763: }
1764: if (parallel) {
1765: PetscMalloc1(4, &options);
1766: options[0] = 1;
1767: options[1] = 0; /* Verbosity */
1768: options[2] = 0; /* Seed */
1769: options[3] = PARMETIS_PSR_COUPLED; /* Seed */
1770: if (viewer) { PetscViewerASCIIPrintf(viewer, "Using ParMETIS to partition graph.\n"); }
1771: if (useInitialGuess) {
1772: if (viewer) { PetscViewerASCIIPrintf(viewer, "Using current distribution of points as initial guess.\n"); }
1773: PetscStackPush("ParMETIS_V3_RefineKway");
1774: ParMETIS_V3_RefineKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1775: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_RefineKway()");
1776: PetscStackPop;
1777: } else {
1778: PetscStackPush("ParMETIS_V3_PartKway");
1779: ParMETIS_V3_PartKway((PetscInt*)cumSumVertices, xadj, adjncy, vtxwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &comm);
1780: PetscStackPop;
1781: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()");
1782: }
1783: PetscFree(options);
1784: } else {
1785: if (viewer) { PetscViewerASCIIPrintf(viewer, "Using METIS to partition graph.\n"); }
1786: Mat As;
1787: PetscInt numRows;
1788: PetscInt *partGlobal;
1789: MatCreateRedundantMatrix(A, size, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &As);
1791: PetscInt *numExclusivelyOwnedAll;
1792: PetscMalloc1(size, &numExclusivelyOwnedAll);
1793: numExclusivelyOwnedAll[rank] = numExclusivelyOwned;
1794: MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,numExclusivelyOwnedAll,1,MPIU_INT,comm);
1796: MatGetSize(As, &numRows, NULL);
1797: PetscMalloc1(numRows, &partGlobal);
1798: if (!rank) {
1799: PetscInt *adjncy_g, *xadj_g, *vtxwgt_g;
1800: lenadjncy = 0;
1802: for (i=0; i<numRows; i++) {
1803: PetscInt temp=0;
1804: MatGetRow(As, i, &temp, NULL, NULL);
1805: lenadjncy += temp;
1806: MatRestoreRow(As, i, &temp, NULL, NULL);
1807: }
1808: PetscMalloc1(lenadjncy, &adjncy_g);
1809: lenxadj = 1 + numRows;
1810: PetscMalloc1(lenxadj, &xadj_g);
1811: xadj_g[0] = 0;
1812: counter = 0;
1813: for (i=0; i<numRows; i++) {
1814: PetscInt temp=0;
1815: const PetscInt *cols;
1816: MatGetRow(As, i, &temp, &cols, NULL);
1817: PetscArraycpy(&adjncy_g[counter], cols, temp);
1818: counter += temp;
1819: xadj_g[i+1] = counter;
1820: MatRestoreRow(As, i, &temp, &cols, NULL);
1821: }
1822: PetscMalloc1(2*numRows, &vtxwgt_g);
1823: for (i=0; i<size; i++){
1824: vtxwgt_g[ncon*cumSumVertices[i]] = numExclusivelyOwnedAll[i];
1825: if (ncon>1) vtxwgt_g[ncon*cumSumVertices[i]+1] = 1;
1826: for (j=cumSumVertices[i]+1; j<cumSumVertices[i+1]; j++) {
1827: vtxwgt_g[ncon*j] = 1;
1828: if (ncon>1) vtxwgt_g[2*j+1] = 0;
1829: }
1830: }
1831: PetscMalloc1(64, &options);
1832: METIS_SetDefaultOptions(options); /* initialize all defaults */
1833: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_SetDefaultOptions()");
1834: options[METIS_OPTION_CONTIG] = 1;
1835: PetscStackPush("METIS_PartGraphKway");
1836: METIS_PartGraphKway(&numRows, &ncon, xadj_g, adjncy_g, vtxwgt_g, NULL, NULL, &nparts, tpwgts, ubvec, options, &edgecut, partGlobal);
1837: PetscStackPop;
1838: if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()");
1839: PetscFree(options);
1840: PetscFree(xadj_g);
1841: PetscFree(adjncy_g);
1842: PetscFree(vtxwgt_g);
1843: }
1844: PetscFree(numExclusivelyOwnedAll);
1846: /* Now scatter the parts array. */
1847: {
1848: PetscMPIInt *counts, *mpiCumSumVertices;
1849: PetscMalloc1(size, &counts);
1850: PetscMalloc1(size+1, &mpiCumSumVertices);
1851: for (i=0; i<size; i++) {
1852: PetscMPIIntCast(cumSumVertices[i+1] - cumSumVertices[i], &(counts[i]));
1853: }
1854: for (i=0; i<=size; i++) {
1855: PetscMPIIntCast(cumSumVertices[i], &(mpiCumSumVertices[i]));
1856: }
1857: MPI_Scatterv(partGlobal, counts, mpiCumSumVertices, MPIU_INT, part, counts[rank], MPIU_INT, 0, comm);
1858: PetscFree(counts);
1859: PetscFree(mpiCumSumVertices);
1860: }
1862: PetscFree(partGlobal);
1863: MatDestroy(&As);
1864: }
1866: MatDestroy(&A);
1867: PetscFree(ubvec);
1868: PetscFree(tpwgts);
1870: /* Now rename the result so that the vertex resembling the exclusively owned points stays on the same rank */
1872: PetscMalloc1(size, &firstVertices);
1873: PetscMalloc1(size, &renumbering);
1874: firstVertices[rank] = part[0];
1875: MPI_Allgather(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,firstVertices,1,MPIU_INT,comm);
1876: for (i=0; i<size; i++) {
1877: renumbering[firstVertices[i]] = i;
1878: }
1879: for (i=0; i<cumSumVertices[rank+1]-cumSumVertices[rank]; i++) {
1880: part[i] = renumbering[part[i]];
1881: }
1882: /* Check if the renumbering worked (this can fail when ParMETIS gives fewer partitions than there are processes) */
1883: failed = (PetscInt)(part[0] != rank);
1884: MPI_Allreduce(&failed, &failedGlobal, 1, MPIU_INT, MPI_SUM, comm);
1886: PetscFree(firstVertices);
1887: PetscFree(renumbering);
1889: if (failedGlobal > 0) {
1890: PetscLayoutDestroy(&layout);
1891: PetscFree(xadj);
1892: PetscFree(adjncy);
1893: PetscFree(vtxwgt);
1894: PetscFree(toBalance);
1895: PetscFree(isLeaf);
1896: PetscFree(isNonExclusivelyOwned);
1897: PetscFree(isExclusivelyOwned);
1898: PetscFree(part);
1899: if (viewer) {
1900: PetscViewerPopFormat(viewer);
1901: PetscViewerDestroy(&viewer);
1902: }
1903: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1904: return(0);
1905: }
1907: /*Let's check how well we did distributing points*/
1908: if (viewer) {
1909: PetscViewerASCIIPrintf(viewer, "Comparing number of owned entities of depth %D on each process before rebalancing, after rebalancing, and after consistency checks.\n", entityDepth);
1910: PetscViewerASCIIPrintf(viewer, "Initial. ");
1911: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, NULL, viewer);
1912: PetscViewerASCIIPrintf(viewer, "Rebalanced. ");
1913: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1914: }
1916: /* Now check that every vertex is owned by a process that it is actually connected to. */
1917: for (i=1; i<=numNonExclusivelyOwned; i++) {
1918: PetscInt loc = 0;
1919: PetscFindInt(cumSumVertices[part[i]], xadj[i+1]-xadj[i], &adjncy[xadj[i]], &loc);
1920: /* If not, then just set the owner to the original owner (hopefully a rare event, it means that a vertex has been isolated) */
1921: if (loc<0) {
1922: part[i] = rank;
1923: }
1924: }
1926: /* Let's see how significant the influences of the previous fixing up step was.*/
1927: if (viewer) {
1928: PetscViewerASCIIPrintf(viewer, "After. ");
1929: DMPlexViewDistribution(comm, cumSumVertices[rank+1]-cumSumVertices[rank], ncon, vtxwgt, part, viewer);
1930: }
1932: PetscLayoutDestroy(&layout);
1933: PetscFree(xadj);
1934: PetscFree(adjncy);
1935: PetscFree(vtxwgt);
1937: /* Almost done, now rewrite the SF to reflect the new ownership. */
1938: {
1939: PetscInt *pointsToRewrite;
1940: PetscMalloc1(numNonExclusivelyOwned, &pointsToRewrite);
1941: counter = 0;
1942: for (i=0; i<pEnd-pStart; i++) {
1943: if (toBalance[i]) {
1944: if (isNonExclusivelyOwned[i]) {
1945: pointsToRewrite[counter] = i + pStart;
1946: counter++;
1947: }
1948: }
1949: }
1950: DMPlexRewriteSF(dm, numNonExclusivelyOwned, pointsToRewrite, part+1, degrees);
1951: PetscFree(pointsToRewrite);
1952: }
1954: PetscFree(toBalance);
1955: PetscFree(isLeaf);
1956: PetscFree(isNonExclusivelyOwned);
1957: PetscFree(isExclusivelyOwned);
1958: PetscFree(part);
1959: if (viewer) {
1960: PetscViewerPopFormat(viewer);
1961: PetscViewerDestroy(&viewer);
1962: }
1963: if (success) *success = PETSC_TRUE;
1964: PetscLogEventEnd(DMPLEX_RebalanceSharedPoints, dm, 0, 0, 0);
1965: return(0);
1966: #else
1967: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-parmetis.");
1968: #endif
1969: }