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