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