1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMPlexCreateNeighborCSR" 5 PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency) 6 { 7 const PetscInt maxFaceCases = 30; 8 PetscInt numFaceCases = 0; 9 PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */ 10 PetscInt *off, *adj; 11 PetscInt *neighborCells = NULL; 12 PetscInt dim, cellDim, depth = 0, faceDepth, cStart, cEnd, c, numCells, cell; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 /* For parallel partitioning, I think you have to communicate supports */ 17 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 18 cellDim = dim - cellHeight; 19 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 20 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 21 if (cEnd - cStart == 0) { 22 if (numVertices) *numVertices = 0; 23 if (offsets) *offsets = NULL; 24 if (adjacency) *adjacency = NULL; 25 PetscFunctionReturn(0); 26 } 27 numCells = cEnd - cStart; 28 faceDepth = depth - cellHeight; 29 if (dim == depth) { 30 PetscInt f, fStart, fEnd; 31 32 ierr = PetscCalloc1(numCells+1, &off);CHKERRQ(ierr); 33 /* Count neighboring cells */ 34 ierr = DMPlexGetHeightStratum(dm, cellHeight+1, &fStart, &fEnd);CHKERRQ(ierr); 35 for (f = fStart; f < fEnd; ++f) { 36 const PetscInt *support; 37 PetscInt supportSize; 38 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 39 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 40 if (supportSize == 2) { 41 ++off[support[0]-cStart+1]; 42 ++off[support[1]-cStart+1]; 43 } 44 } 45 /* Prefix sum */ 46 for (c = 1; c <= numCells; ++c) off[c] += off[c-1]; 47 if (adjacency) { 48 PetscInt *tmp; 49 50 ierr = PetscMalloc1(off[numCells], &adj);CHKERRQ(ierr); 51 ierr = PetscMalloc1((numCells+1), &tmp);CHKERRQ(ierr); 52 ierr = PetscMemcpy(tmp, off, (numCells+1) * sizeof(PetscInt));CHKERRQ(ierr); 53 /* Get neighboring cells */ 54 for (f = fStart; f < fEnd; ++f) { 55 const PetscInt *support; 56 PetscInt supportSize; 57 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 58 ierr = DMPlexGetSupport(dm, f, &support);CHKERRQ(ierr); 59 if (supportSize == 2) { 60 adj[tmp[support[0]-cStart]++] = support[1]; 61 adj[tmp[support[1]-cStart]++] = support[0]; 62 } 63 } 64 #if defined(PETSC_USE_DEBUG) 65 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); 66 #endif 67 ierr = PetscFree(tmp);CHKERRQ(ierr); 68 } 69 if (numVertices) *numVertices = numCells; 70 if (offsets) *offsets = off; 71 if (adjacency) *adjacency = adj; 72 PetscFunctionReturn(0); 73 } 74 /* Setup face recognition */ 75 if (faceDepth == 1) { 76 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 */ 77 78 for (c = cStart; c < cEnd; ++c) { 79 PetscInt corners; 80 81 ierr = DMPlexGetConeSize(dm, c, &corners);CHKERRQ(ierr); 82 if (!cornersSeen[corners]) { 83 PetscInt nFV; 84 85 if (numFaceCases >= maxFaceCases) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases"); 86 cornersSeen[corners] = 1; 87 88 ierr = DMPlexGetNumFaceVertices(dm, cellDim, corners, &nFV);CHKERRQ(ierr); 89 90 numFaceVertices[numFaceCases++] = nFV; 91 } 92 } 93 } 94 ierr = PetscCalloc1(numCells+1, &off);CHKERRQ(ierr); 95 /* Count neighboring cells */ 96 for (cell = cStart; cell < cEnd; ++cell) { 97 PetscInt numNeighbors = PETSC_DETERMINE, n; 98 99 ierr = DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, &numNeighbors, &neighborCells);CHKERRQ(ierr); 100 /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */ 101 for (n = 0; n < numNeighbors; ++n) { 102 PetscInt cellPair[2]; 103 PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE; 104 PetscInt meetSize = 0; 105 const PetscInt *meet = NULL; 106 107 cellPair[0] = cell; cellPair[1] = neighborCells[n]; 108 if (cellPair[0] == cellPair[1]) continue; 109 if (!found) { 110 ierr = DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr); 111 if (meetSize) { 112 PetscInt f; 113 114 for (f = 0; f < numFaceCases; ++f) { 115 if (numFaceVertices[f] == meetSize) { 116 found = PETSC_TRUE; 117 break; 118 } 119 } 120 } 121 ierr = DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr); 122 } 123 if (found) ++off[cell-cStart+1]; 124 } 125 } 126 /* Prefix sum */ 127 for (cell = 1; cell <= numCells; ++cell) off[cell] += off[cell-1]; 128 129 if (adjacency) { 130 ierr = PetscMalloc1(off[numCells], &adj);CHKERRQ(ierr); 131 /* Get neighboring cells */ 132 for (cell = cStart; cell < cEnd; ++cell) { 133 PetscInt numNeighbors = PETSC_DETERMINE, n; 134 PetscInt cellOffset = 0; 135 136 ierr = DMPlexGetAdjacency_Internal(dm, cell, PETSC_TRUE, PETSC_FALSE, &numNeighbors, &neighborCells);CHKERRQ(ierr); 137 /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */ 138 for (n = 0; n < numNeighbors; ++n) { 139 PetscInt cellPair[2]; 140 PetscBool found = faceDepth > 1 ? PETSC_TRUE : PETSC_FALSE; 141 PetscInt meetSize = 0; 142 const PetscInt *meet = NULL; 143 144 cellPair[0] = cell; cellPair[1] = neighborCells[n]; 145 if (cellPair[0] == cellPair[1]) continue; 146 if (!found) { 147 ierr = DMPlexGetMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr); 148 if (meetSize) { 149 PetscInt f; 150 151 for (f = 0; f < numFaceCases; ++f) { 152 if (numFaceVertices[f] == meetSize) { 153 found = PETSC_TRUE; 154 break; 155 } 156 } 157 } 158 ierr = DMPlexRestoreMeet(dm, 2, cellPair, &meetSize, &meet);CHKERRQ(ierr); 159 } 160 if (found) { 161 adj[off[cell-cStart]+cellOffset] = neighborCells[n]; 162 ++cellOffset; 163 } 164 } 165 } 166 } 167 ierr = PetscFree(neighborCells);CHKERRQ(ierr); 168 if (numVertices) *numVertices = numCells; 169 if (offsets) *offsets = off; 170 if (adjacency) *adjacency = adj; 171 PetscFunctionReturn(0); 172 } 173 174 #if defined(PETSC_HAVE_CHACO) 175 #if defined(PETSC_HAVE_UNISTD_H) 176 #include <unistd.h> 177 #endif 178 /* Chaco does not have an include file */ 179 PETSC_EXTERN int interface(int nvtxs, int *start, int *adjacency, int *vwgts, 180 float *ewgts, float *x, float *y, float *z, char *outassignname, 181 char *outfilename, short *assignment, int architecture, int ndims_tot, 182 int mesh_dims[3], double *goal, int global_method, int local_method, 183 int rqi_flag, int vmax, int ndims, double eigtol, long seed); 184 185 extern int FREE_GRAPH; 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "DMPlexPartition_Chaco" 189 PetscErrorCode DMPlexPartition_Chaco(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition) 190 { 191 enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3}; 192 MPI_Comm comm; 193 int nvtxs = numVertices; /* number of vertices in full graph */ 194 int *vwgts = NULL; /* weights for all vertices */ 195 float *ewgts = NULL; /* weights for all edges */ 196 float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */ 197 char *outassignname = NULL; /* name of assignment output file */ 198 char *outfilename = NULL; /* output file name */ 199 int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */ 200 int ndims_tot = 0; /* total number of cube dimensions to divide */ 201 int mesh_dims[3]; /* dimensions of mesh of processors */ 202 double *goal = NULL; /* desired set sizes for each set */ 203 int global_method = 1; /* global partitioning algorithm */ 204 int local_method = 1; /* local partitioning algorithm */ 205 int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */ 206 int vmax = 200; /* how many vertices to coarsen down to? */ 207 int ndims = 1; /* number of eigenvectors (2^d sets) */ 208 double eigtol = 0.001; /* tolerance on eigenvectors */ 209 long seed = 123636512; /* for random graph mutations */ 210 short int *assignment; /* Output partition */ 211 int fd_stdout, fd_pipe[2]; 212 PetscInt *points; 213 PetscMPIInt commSize; 214 int i, v, p; 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 219 ierr = MPI_Comm_size(comm, &commSize);CHKERRQ(ierr); 220 if (!numVertices) { 221 ierr = PetscSectionCreate(comm, partSection);CHKERRQ(ierr); 222 ierr = PetscSectionSetChart(*partSection, 0, commSize);CHKERRQ(ierr); 223 ierr = PetscSectionSetUp(*partSection);CHKERRQ(ierr); 224 ierr = ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, partition);CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 FREE_GRAPH = 0; /* Do not let Chaco free my memory */ 228 for (i = 0; i < start[numVertices]; ++i) ++adjacency[i]; 229 230 if (global_method == INERTIAL_METHOD) { 231 /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */ 232 SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported"); 233 } 234 mesh_dims[0] = commSize; 235 mesh_dims[1] = 1; 236 mesh_dims[2] = 1; 237 ierr = PetscMalloc1(nvtxs, &assignment);CHKERRQ(ierr); 238 /* Chaco outputs to stdout. We redirect this to a buffer. */ 239 /* TODO: check error codes for UNIX calls */ 240 #if defined(PETSC_HAVE_UNISTD_H) 241 { 242 int piperet; 243 piperet = pipe(fd_pipe); 244 if (piperet) SETERRQ(comm,PETSC_ERR_SYS,"Could not create pipe"); 245 fd_stdout = dup(1); 246 close(1); 247 dup2(fd_pipe[1], 1); 248 } 249 #endif 250 ierr = interface(nvtxs, (int*) start, (int*) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename, 251 assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag, 252 vmax, ndims, eigtol, seed); 253 #if defined(PETSC_HAVE_UNISTD_H) 254 { 255 char msgLog[10000]; 256 int count; 257 258 fflush(stdout); 259 count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char)); 260 if (count < 0) count = 0; 261 msgLog[count] = 0; 262 close(1); 263 dup2(fd_stdout, 1); 264 close(fd_stdout); 265 close(fd_pipe[0]); 266 close(fd_pipe[1]); 267 if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog); 268 } 269 #endif 270 /* Convert to PetscSection+IS */ 271 ierr = PetscSectionCreate(comm, partSection);CHKERRQ(ierr); 272 ierr = PetscSectionSetChart(*partSection, 0, commSize);CHKERRQ(ierr); 273 for (v = 0; v < nvtxs; ++v) { 274 ierr = PetscSectionAddDof(*partSection, assignment[v], 1);CHKERRQ(ierr); 275 } 276 ierr = PetscSectionSetUp(*partSection);CHKERRQ(ierr); 277 ierr = PetscMalloc1(nvtxs, &points);CHKERRQ(ierr); 278 for (p = 0, i = 0; p < commSize; ++p) { 279 for (v = 0; v < nvtxs; ++v) { 280 if (assignment[v] == p) points[i++] = v; 281 } 282 } 283 if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs); 284 ierr = ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);CHKERRQ(ierr); 285 if (global_method == INERTIAL_METHOD) { 286 /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */ 287 } 288 ierr = PetscFree(assignment);CHKERRQ(ierr); 289 for (i = 0; i < start[numVertices]; ++i) --adjacency[i]; 290 PetscFunctionReturn(0); 291 } 292 #endif 293 294 #if defined(PETSC_HAVE_PARMETIS) 295 #include <parmetis.h> 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "DMPlexPartition_ParMetis" 299 PetscErrorCode DMPlexPartition_ParMetis(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition) 300 { 301 MPI_Comm comm; 302 PetscInt nvtxs = numVertices; /* The number of vertices in full graph */ 303 PetscInt *vtxdist; /* Distribution of vertices across processes */ 304 PetscInt *xadj = start; /* Start of edge list for each vertex */ 305 PetscInt *adjncy = adjacency; /* Edge lists for all vertices */ 306 PetscInt *vwgt = NULL; /* Vertex weights */ 307 PetscInt *adjwgt = NULL; /* Edge weights */ 308 PetscInt wgtflag = 0; /* Indicates which weights are present */ 309 PetscInt numflag = 0; /* Indicates initial offset (0 or 1) */ 310 PetscInt ncon = 1; /* The number of weights per vertex */ 311 PetscInt nparts; /* The number of partitions */ 312 PetscReal *tpwgts; /* The fraction of vertex weights assigned to each partition */ 313 PetscReal *ubvec; /* The balance intolerance for vertex weights */ 314 PetscInt options[5]; /* Options */ 315 /* Outputs */ 316 PetscInt edgeCut; /* The number of edges cut by the partition */ 317 PetscInt *assignment, *points; 318 PetscMPIInt commSize, rank, p, v, i; 319 PetscErrorCode ierr; 320 321 PetscFunctionBegin; 322 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 323 ierr = MPI_Comm_size(comm, &commSize);CHKERRQ(ierr); 324 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 325 nparts = commSize; 326 options[0] = 0; /* Use all defaults */ 327 /* Calculate vertex distribution */ 328 ierr = PetscMalloc4(nparts+1,&vtxdist,nparts*ncon,&tpwgts,ncon,&ubvec,nvtxs,&assignment);CHKERRQ(ierr); 329 vtxdist[0] = 0; 330 ierr = MPI_Allgather(&nvtxs, 1, MPIU_INT, &vtxdist[1], 1, MPIU_INT, comm);CHKERRQ(ierr); 331 for (p = 2; p <= nparts; ++p) { 332 vtxdist[p] += vtxdist[p-1]; 333 } 334 /* Calculate weights */ 335 for (p = 0; p < nparts; ++p) { 336 tpwgts[p] = 1.0/nparts; 337 } 338 ubvec[0] = 1.05; 339 340 if (nparts == 1) { 341 ierr = PetscMemzero(assignment, nvtxs * sizeof(PetscInt)); 342 } else { 343 if (vtxdist[1] == vtxdist[nparts]) { 344 if (!rank) { 345 PetscStackPush("METIS_PartGraphKway"); 346 ierr = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, NULL, adjwgt, &nparts, tpwgts, ubvec, NULL, &edgeCut, assignment); 347 PetscStackPop; 348 if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in METIS_PartGraphKway()"); 349 } 350 } else { 351 PetscStackPush("ParMETIS_V3_PartKway"); 352 ierr = ParMETIS_V3_PartKway(vtxdist, xadj, adjncy, vwgt, adjwgt, &wgtflag, &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgeCut, assignment, &comm); 353 PetscStackPop; 354 if (ierr != METIS_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in ParMETIS_V3_PartKway()"); 355 } 356 } 357 /* Convert to PetscSection+IS */ 358 ierr = PetscSectionCreate(comm, partSection);CHKERRQ(ierr); 359 ierr = PetscSectionSetChart(*partSection, 0, commSize);CHKERRQ(ierr); 360 for (v = 0; v < nvtxs; ++v) { 361 ierr = PetscSectionAddDof(*partSection, assignment[v], 1);CHKERRQ(ierr); 362 } 363 ierr = PetscSectionSetUp(*partSection);CHKERRQ(ierr); 364 ierr = PetscMalloc1(nvtxs, &points);CHKERRQ(ierr); 365 for (p = 0, i = 0; p < commSize; ++p) { 366 for (v = 0; v < nvtxs; ++v) { 367 if (assignment[v] == p) points[i++] = v; 368 } 369 } 370 if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs); 371 ierr = ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);CHKERRQ(ierr); 372 ierr = PetscFree4(vtxdist,tpwgts,ubvec,assignment);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375 #endif 376 377 #undef __FUNCT__ 378 #define __FUNCT__ "DMPlexEnlargePartition" 379 /* Expand the partition by BFS on the adjacency graph */ 380 PetscErrorCode DMPlexEnlargePartition(DM dm, const PetscInt start[], const PetscInt adjacency[], PetscSection origPartSection, IS origPartition, PetscSection *partSection, IS *partition) 381 { 382 PetscHashI h; 383 const PetscInt *points; 384 PetscInt **tmpPoints, *newPoints, totPoints = 0; 385 PetscInt pStart, pEnd, part, q; 386 PetscBool useCone; 387 PetscErrorCode ierr; 388 389 PetscFunctionBegin; 390 PetscHashICreate(h); 391 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), partSection);CHKERRQ(ierr); 392 ierr = PetscSectionGetChart(origPartSection, &pStart, &pEnd);CHKERRQ(ierr); 393 ierr = PetscSectionSetChart(*partSection, pStart, pEnd);CHKERRQ(ierr); 394 ierr = ISGetIndices(origPartition, &points);CHKERRQ(ierr); 395 ierr = PetscMalloc1((pEnd - pStart), &tmpPoints);CHKERRQ(ierr); 396 ierr = DMPlexGetAdjacencyUseCone(dm, &useCone);CHKERRQ(ierr); 397 ierr = DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);CHKERRQ(ierr); 398 for (part = pStart; part < pEnd; ++part) { 399 PetscInt *adj = NULL; 400 PetscInt numPoints, nP, numNewPoints, off, p, n = 0; 401 402 PetscHashIClear(h); 403 ierr = PetscSectionGetDof(origPartSection, part, &numPoints);CHKERRQ(ierr); 404 ierr = PetscSectionGetOffset(origPartSection, part, &off);CHKERRQ(ierr); 405 /* Add all existing points to h */ 406 for (p = 0; p < numPoints; ++p) { 407 const PetscInt point = points[off+p]; 408 PetscHashIAdd(h, point, 1); 409 } 410 PetscHashISize(h, nP); 411 if (nP != numPoints) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid partition has %d points, but only %d were unique", numPoints, nP); 412 /* Add all points in next BFS level */ 413 for (p = 0; p < numPoints; ++p) { 414 const PetscInt point = points[off+p]; 415 PetscInt adjSize = PETSC_DETERMINE, a; 416 417 ierr = DMPlexGetAdjacency(dm, point, &adjSize, &adj);CHKERRQ(ierr); 418 for (a = 0; a < adjSize; ++a) PetscHashIAdd(h, adj[a], 1); 419 } 420 PetscHashISize(h, numNewPoints); 421 ierr = PetscSectionSetDof(*partSection, part, numNewPoints);CHKERRQ(ierr); 422 ierr = PetscMalloc1(numNewPoints, &tmpPoints[part]);CHKERRQ(ierr); 423 ierr = PetscHashIGetKeys(h, &n, tmpPoints[part]);CHKERRQ(ierr); 424 ierr = PetscFree(adj);CHKERRQ(ierr); 425 totPoints += numNewPoints; 426 } 427 ierr = DMPlexSetAdjacencyUseCone(dm, useCone);CHKERRQ(ierr); 428 ierr = ISRestoreIndices(origPartition, &points);CHKERRQ(ierr); 429 PetscHashIDestroy(h); 430 ierr = PetscSectionSetUp(*partSection);CHKERRQ(ierr); 431 ierr = PetscMalloc1(totPoints, &newPoints);CHKERRQ(ierr); 432 for (part = pStart, q = 0; part < pEnd; ++part) { 433 PetscInt numPoints, p; 434 435 ierr = PetscSectionGetDof(*partSection, part, &numPoints);CHKERRQ(ierr); 436 for (p = 0; p < numPoints; ++p, ++q) newPoints[q] = tmpPoints[part][p]; 437 ierr = PetscFree(tmpPoints[part]);CHKERRQ(ierr); 438 } 439 ierr = PetscFree(tmpPoints);CHKERRQ(ierr); 440 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), totPoints, newPoints, PETSC_OWN_POINTER, partition);CHKERRQ(ierr); 441 PetscFunctionReturn(0); 442 } 443 444 #undef __FUNCT__ 445 #define __FUNCT__ "DMPlexCreatePartition" 446 /* 447 DMPlexCreatePartition - Create a non-overlapping partition of the points at the given height 448 449 Collective on DM 450 451 Input Parameters: 452 + dm - The DM 453 . height - The height for points in the partition 454 - enlarge - Expand each partition with neighbors 455 456 Output Parameters: 457 + partSection - The PetscSection giving the division of points by partition 458 . partition - The list of points by partition 459 . origPartSection - If enlarge is true, the PetscSection giving the division of points before enlarging by partition, otherwise NULL 460 - origPartition - If enlarge is true, the list of points before enlarging by partition, otherwise NULL 461 462 Level: developer 463 464 .seealso DMPlexDistribute() 465 */ 466 PetscErrorCode DMPlexCreatePartition(DM dm, const char name[], PetscInt height, PetscBool enlarge, PetscSection *partSection, IS *partition, PetscSection *origPartSection, IS *origPartition) 467 { 468 char partname[1024]; 469 PetscBool isChaco = PETSC_FALSE, isMetis = PETSC_FALSE, flg; 470 PetscMPIInt size; 471 PetscErrorCode ierr; 472 473 PetscFunctionBegin; 474 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);CHKERRQ(ierr); 475 476 *origPartSection = NULL; 477 *origPartition = NULL; 478 if (size == 1) { 479 PetscInt *points; 480 PetscInt cStart, cEnd, c; 481 482 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 483 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), partSection);CHKERRQ(ierr); 484 ierr = PetscSectionSetChart(*partSection, 0, size);CHKERRQ(ierr); 485 ierr = PetscSectionSetDof(*partSection, 0, cEnd-cStart);CHKERRQ(ierr); 486 ierr = PetscSectionSetUp(*partSection);CHKERRQ(ierr); 487 ierr = PetscMalloc1((cEnd - cStart), &points);CHKERRQ(ierr); 488 for (c = cStart; c < cEnd; ++c) points[c] = c; 489 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), cEnd-cStart, points, PETSC_OWN_POINTER, partition);CHKERRQ(ierr); 490 PetscFunctionReturn(0); 491 } 492 ierr = PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_plex_partitioner", partname, 1024, &flg);CHKERRQ(ierr); 493 if (flg) name = partname; 494 if (name) { 495 ierr = PetscStrcmp(name, "chaco", &isChaco);CHKERRQ(ierr); 496 ierr = PetscStrcmp(name, "metis", &isMetis);CHKERRQ(ierr); 497 } 498 if (height == 0) { 499 PetscInt numVertices; 500 PetscInt *start = NULL; 501 PetscInt *adjacency = NULL; 502 503 ierr = DMPlexCreateNeighborCSR(dm, 0, &numVertices, &start, &adjacency);CHKERRQ(ierr); 504 if (!name || isChaco) { 505 #if defined(PETSC_HAVE_CHACO) 506 ierr = DMPlexPartition_Chaco(dm, numVertices, start, adjacency, partSection, partition);CHKERRQ(ierr); 507 #else 508 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Mesh partitioning needs external package support.\nPlease reconfigure with --download-chaco."); 509 #endif 510 } else if (isMetis) { 511 #if defined(PETSC_HAVE_PARMETIS) 512 ierr = DMPlexPartition_ParMetis(dm, numVertices, start, adjacency, partSection, partition);CHKERRQ(ierr); 513 #endif 514 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Unknown mesh partitioning package %s", name); 515 if (enlarge) { 516 *origPartSection = *partSection; 517 *origPartition = *partition; 518 519 ierr = DMPlexEnlargePartition(dm, start, adjacency, *origPartSection, *origPartition, partSection, partition);CHKERRQ(ierr); 520 } 521 ierr = PetscFree(start);CHKERRQ(ierr); 522 ierr = PetscFree(adjacency);CHKERRQ(ierr); 523 # if 0 524 } else if (height == 1) { 525 /* Build the dual graph for faces and partition the hypergraph */ 526 PetscInt numEdges; 527 528 buildFaceCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase()); 529 GraphPartitioner().partition(numEdges, start, adjacency, partition, manager); 530 destroyCSR(numEdges, start, adjacency); 531 #endif 532 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid partition height %D", height); 533 PetscFunctionReturn(0); 534 } 535 536 #undef __FUNCT__ 537 #define __FUNCT__ "DMPlexCreatePartitionClosure" 538 PetscErrorCode DMPlexCreatePartitionClosure(DM dm, PetscSection pointSection, IS pointPartition, PetscSection *section, IS *partition) 539 { 540 /* const PetscInt height = 0; */ 541 const PetscInt *partArray; 542 PetscInt *allPoints, *packPoints; 543 PetscInt rStart, rEnd, rank, pStart, pEnd, newSize; 544 PetscErrorCode ierr; 545 PetscBT bt; 546 PetscSegBuffer segpack,segpart; 547 548 PetscFunctionBegin; 549 ierr = PetscSectionGetChart(pointSection, &rStart, &rEnd);CHKERRQ(ierr); 550 ierr = ISGetIndices(pointPartition, &partArray);CHKERRQ(ierr); 551 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);CHKERRQ(ierr); 552 ierr = PetscSectionSetChart(*section, rStart, rEnd);CHKERRQ(ierr); 553 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 554 ierr = PetscBTCreate(pEnd-pStart,&bt);CHKERRQ(ierr); 555 ierr = PetscSegBufferCreate(sizeof(PetscInt),1000,&segpack);CHKERRQ(ierr); 556 ierr = PetscSegBufferCreate(sizeof(PetscInt),1000,&segpart);CHKERRQ(ierr); 557 for (rank = rStart; rank < rEnd; ++rank) { 558 PetscInt partSize = 0, numPoints, offset, p, *PETSC_RESTRICT placePoints; 559 560 ierr = PetscSectionGetDof(pointSection, rank, &numPoints);CHKERRQ(ierr); 561 ierr = PetscSectionGetOffset(pointSection, rank, &offset);CHKERRQ(ierr); 562 for (p = 0; p < numPoints; ++p) { 563 PetscInt point = partArray[offset+p], closureSize, c; 564 PetscInt *closure = NULL; 565 566 /* TODO Include support for height > 0 case */ 567 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 568 for (c=0; c<closureSize; c++) { 569 PetscInt cpoint = closure[c*2]; 570 if (!PetscBTLookupSet(bt,cpoint-pStart)) { 571 PetscInt *PETSC_RESTRICT pt; 572 partSize++; 573 ierr = PetscSegBufferGetInts(segpart,1,&pt);CHKERRQ(ierr); 574 *pt = cpoint; 575 } 576 } 577 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 578 } 579 ierr = PetscSectionSetDof(*section, rank, partSize);CHKERRQ(ierr); 580 ierr = PetscSegBufferGetInts(segpack,partSize,&placePoints);CHKERRQ(ierr); 581 ierr = PetscSegBufferExtractTo(segpart,placePoints);CHKERRQ(ierr); 582 ierr = PetscSortInt(partSize,placePoints);CHKERRQ(ierr); 583 for (p=0; p<partSize; p++) {ierr = PetscBTClear(bt,placePoints[p]-pStart);CHKERRQ(ierr);} 584 } 585 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 586 ierr = PetscSegBufferDestroy(&segpart);CHKERRQ(ierr); 587 588 ierr = PetscSectionSetUp(*section);CHKERRQ(ierr); 589 ierr = PetscSectionGetStorageSize(*section, &newSize);CHKERRQ(ierr); 590 ierr = PetscMalloc1(newSize, &allPoints);CHKERRQ(ierr); 591 592 ierr = PetscSegBufferExtractInPlace(segpack,&packPoints);CHKERRQ(ierr); 593 for (rank = rStart; rank < rEnd; ++rank) { 594 PetscInt numPoints, offset; 595 596 ierr = PetscSectionGetDof(*section, rank, &numPoints);CHKERRQ(ierr); 597 ierr = PetscSectionGetOffset(*section, rank, &offset);CHKERRQ(ierr); 598 ierr = PetscMemcpy(&allPoints[offset], packPoints, numPoints * sizeof(PetscInt));CHKERRQ(ierr); 599 packPoints += numPoints; 600 } 601 602 ierr = PetscSegBufferDestroy(&segpack);CHKERRQ(ierr); 603 ierr = ISRestoreIndices(pointPartition, &partArray);CHKERRQ(ierr); 604 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), newSize, allPoints, PETSC_OWN_POINTER, partition);CHKERRQ(ierr); 605 PetscFunctionReturn(0); 606 } 607