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