Actual source code: pmetis.c
2: #include <../src/mat/impls/adj/mpi/mpiadj.h>
4: /*
5: Currently using ParMetis-4.0.2
6: */
8: #include <parmetis.h>
10: /*
11: The first 5 elements of this structure are the input control array to Metis
12: */
13: typedef struct {
14: PetscInt cuts; /* number of cuts made (output) */
15: PetscInt foldfactor;
16: PetscInt parallel; /* use parallel partitioner for coarse problem */
17: PetscInt indexing; /* 0 indicates C indexing, 1 Fortran */
18: PetscInt printout; /* indicates if one wishes Metis to print info */
19: PetscBool repartition;
20: } MatPartitioning_Parmetis;
22: #define PetscCallPARMETIS(n,func) do { \
26: } while (0)
28: #define PetscStackCallParmetis_(name,func,args) do { \
29: PetscStackPush(name); \
30: int status = func args; \
31: PetscStackPop; \
32: status,name; \
33: } while (0)
35: #define PetscStackCallParmetis(func,args) PetscStackCallParmetis_(PetscStringize(func),func,args)
37: static PetscErrorCode MatPartitioningApply_Parmetis_Private(MatPartitioning part, PetscBool useND, PetscBool isImprove, IS *partitioning)
38: {
39: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis*)part->data;
40: PetscInt *locals = NULL;
41: Mat mat = part->adj,amat,pmat;
42: PetscBool flg;
43: PetscInt bs = 1;
47: PetscObjectTypeCompare((PetscObject)mat,MATMPIADJ,&flg);
48: if (flg) {
49: amat = mat;
50: PetscObjectReference((PetscObject)amat);
51: } else {
52: /* bs indicates if the converted matrix is "reduced" from the original and hence the
53: resulting partition results need to be stretched to match the original matrix */
54: MatConvert(mat,MATMPIADJ,MAT_INITIAL_MATRIX,&amat);
55: if (amat->rmap->n > 0) bs = mat->rmap->n/amat->rmap->n;
56: }
57: MatMPIAdjCreateNonemptySubcommMat(amat,&pmat);
58: MPI_Barrier(PetscObjectComm((PetscObject)part));
60: if (pmat) {
61: MPI_Comm pcomm,comm;
62: Mat_MPIAdj *adj = (Mat_MPIAdj*)pmat->data;
63: PetscInt *vtxdist = pmat->rmap->range;
64: PetscInt *xadj = adj->i;
65: PetscInt *adjncy = adj->j;
66: PetscInt *NDorder = NULL;
67: PetscInt itmp = 0,wgtflag=0, numflag=0, ncon=1, nparts=part->n, options[24], i, j;
68: real_t *tpwgts,*ubvec,itr=0.1;
70: PetscObjectGetComm((PetscObject)pmat,&pcomm);
71: if (PetscDefined(USE_DEBUG)) {
72: /* check that matrix has no diagonal entries */
73: PetscInt rstart;
74: MatGetOwnershipRange(pmat,&rstart,NULL);
75: for (i=0; i<pmat->rmap->n; i++) {
76: for (j=xadj[i]; j<xadj[i+1]; j++) {
78: }
79: }
80: }
82: PetscMalloc1(pmat->rmap->n,&locals);
84: if (isImprove) {
85: PetscInt i;
86: const PetscInt *part_indices;
88: ISGetIndices(*partitioning,&part_indices);
89: for (i=0; i<pmat->rmap->n; i++) locals[i] = part_indices[i*bs];
90: ISRestoreIndices(*partitioning,&part_indices);
91: ISDestroy(partitioning);
92: }
94: if (adj->values && part->use_edge_weights && !part->vertex_weights) wgtflag = 1;
95: if (part->vertex_weights && !adj->values) wgtflag = 2;
96: if (part->vertex_weights && adj->values && part->use_edge_weights) wgtflag = 3;
98: if (PetscLogPrintInfo) {itmp = pmetis->printout; pmetis->printout = 127;}
99: PetscMalloc1(ncon*nparts,&tpwgts);
100: for (i=0; i<ncon; i++) {
101: for (j=0; j<nparts; j++) {
102: if (part->part_weights) {
103: tpwgts[i*nparts+j] = part->part_weights[i*nparts+j];
104: } else {
105: tpwgts[i*nparts+j] = 1./nparts;
106: }
107: }
108: }
109: PetscMalloc1(ncon,&ubvec);
110: for (i=0; i<ncon; i++) ubvec[i] = 1.05;
111: /* This sets the defaults */
112: options[0] = 0;
113: for (i=1; i<24; i++) options[i] = -1;
114: /* Duplicate the communicator to be sure that ParMETIS attribute caching does not interfere with PETSc. */
115: MPI_Comm_dup(pcomm,&comm);
116: if (useND) {
117: PetscInt *sizes, *seps, log2size, subd, *level;
118: PetscMPIInt size;
119: idx_t mtype = PARMETIS_MTYPE_GLOBAL, rtype = PARMETIS_SRTYPE_2PHASE, p_nseps = 1, s_nseps = 1;
120: real_t ubfrac = 1.05;
122: MPI_Comm_size(comm,&size);
123: PetscMalloc1(pmat->rmap->n,&NDorder);
124: PetscMalloc3(2*size,&sizes,4*size,&seps,size,&level);
125: PetscStackCallParmetis(ParMETIS_V32_NodeND,((idx_t*)vtxdist,(idx_t*)xadj,(idx_t*)adjncy,(idx_t*)part->vertex_weights,(idx_t*)&numflag,&mtype,&rtype,&p_nseps,&s_nseps,&ubfrac,NULL/* seed */,NULL/* dbglvl */,(idx_t*)NDorder,(idx_t*)(sizes),&comm));
126: log2size = PetscLog2Real(size);
127: subd = PetscPowInt(2,log2size);
128: MatPartitioningSizesToSep_Private(subd,sizes,seps,level);
129: for (i=0;i<pmat->rmap->n;i++) {
130: PetscInt loc;
132: PetscFindInt(NDorder[i],2*subd,seps,&loc);
133: if (loc < 0) {
134: loc = -(loc+1);
135: if (loc%2) { /* part of subdomain */
136: locals[i] = loc/2;
137: } else {
138: PetscFindInt(NDorder[i],2*(subd-1),seps+2*subd,&loc);
139: loc = loc < 0 ? -(loc+1)/2 : loc/2;
140: locals[i] = level[loc];
141: }
142: } else locals[i] = loc/2;
143: }
144: PetscFree3(sizes,seps,level);
145: } else {
146: if (pmetis->repartition) {
147: PetscStackCallParmetis(ParMETIS_V3_AdaptiveRepart,((idx_t*)vtxdist,(idx_t*)xadj,(idx_t*)adjncy,(idx_t*)part->vertex_weights,(idx_t*)part->vertex_weights,(idx_t*)adj->values,(idx_t*)&wgtflag,(idx_t*)&numflag,(idx_t*)&ncon,(idx_t*)&nparts,tpwgts,ubvec,&itr,(idx_t*)options,(idx_t*)&pmetis->cuts,(idx_t*)locals,&comm));
148: } else if (isImprove) {
149: PetscStackCallParmetis(ParMETIS_V3_RefineKway,((idx_t*)vtxdist,(idx_t*)xadj,(idx_t*)adjncy,(idx_t*)part->vertex_weights,(idx_t*)adj->values,(idx_t*)&wgtflag,(idx_t*)&numflag,(idx_t*)&ncon,(idx_t*)&nparts,tpwgts,ubvec,(idx_t*)options,(idx_t*)&pmetis->cuts,(idx_t*)locals,&comm));
150: } else {
151: PetscStackCallParmetis(ParMETIS_V3_PartKway,((idx_t*)vtxdist,(idx_t*)xadj,(idx_t*)adjncy,(idx_t*)part->vertex_weights,(idx_t*)adj->values,(idx_t*)&wgtflag,(idx_t*)&numflag,(idx_t*)&ncon,(idx_t*)&nparts,tpwgts,ubvec,(idx_t*)options,(idx_t*)&pmetis->cuts,(idx_t*)locals,&comm));
152: }
153: }
154: MPI_Comm_free(&comm);
156: PetscFree(tpwgts);
157: PetscFree(ubvec);
158: if (PetscLogPrintInfo) pmetis->printout = itmp;
160: if (bs > 1) {
161: PetscInt i,j,*newlocals;
162: PetscMalloc1(bs*pmat->rmap->n,&newlocals);
163: for (i=0; i<pmat->rmap->n; i++) {
164: for (j=0; j<bs; j++) {
165: newlocals[bs*i + j] = locals[i];
166: }
167: }
168: PetscFree(locals);
169: ISCreateGeneral(PetscObjectComm((PetscObject)part),bs*pmat->rmap->n,newlocals,PETSC_OWN_POINTER,partitioning);
170: } else {
171: ISCreateGeneral(PetscObjectComm((PetscObject)part),pmat->rmap->n,locals,PETSC_OWN_POINTER,partitioning);
172: }
173: if (useND) {
174: IS ndis;
176: if (bs > 1) {
177: ISCreateBlock(PetscObjectComm((PetscObject)part),bs,pmat->rmap->n,NDorder,PETSC_OWN_POINTER,&ndis);
178: } else {
179: ISCreateGeneral(PetscObjectComm((PetscObject)part),pmat->rmap->n,NDorder,PETSC_OWN_POINTER,&ndis);
180: }
181: ISSetPermutation(ndis);
182: PetscObjectCompose((PetscObject)(*partitioning),"_petsc_matpartitioning_ndorder",(PetscObject)ndis);
183: ISDestroy(&ndis);
184: }
185: } else {
186: ISCreateGeneral(PetscObjectComm((PetscObject)part),0,NULL,PETSC_COPY_VALUES,partitioning);
187: if (useND) {
188: IS ndis;
190: if (bs > 1) {
191: ISCreateBlock(PetscObjectComm((PetscObject)part),bs,0,NULL,PETSC_COPY_VALUES,&ndis);
192: } else {
193: ISCreateGeneral(PetscObjectComm((PetscObject)part),0,NULL,PETSC_COPY_VALUES,&ndis);
194: }
195: ISSetPermutation(ndis);
196: PetscObjectCompose((PetscObject)(*partitioning),"_petsc_matpartitioning_ndorder",(PetscObject)ndis);
197: ISDestroy(&ndis);
198: }
199: }
200: MatDestroy(&pmat);
201: MatDestroy(&amat);
202: return 0;
203: }
205: /*
206: Uses the ParMETIS parallel matrix partitioner to compute a nested dissection ordering of the matrix in parallel
207: */
208: static PetscErrorCode MatPartitioningApplyND_Parmetis(MatPartitioning part, IS *partitioning)
209: {
210: MatPartitioningApply_Parmetis_Private(part, PETSC_TRUE, PETSC_FALSE, partitioning);
211: return 0;
212: }
214: /*
215: Uses the ParMETIS parallel matrix partitioner to partition the matrix in parallel
216: */
217: static PetscErrorCode MatPartitioningApply_Parmetis(MatPartitioning part, IS *partitioning)
218: {
219: MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_FALSE, partitioning);
220: return 0;
221: }
223: /*
224: Uses the ParMETIS to improve the quality of a partition
225: */
226: static PetscErrorCode MatPartitioningImprove_Parmetis(MatPartitioning part, IS *partitioning)
227: {
228: MatPartitioningApply_Parmetis_Private(part, PETSC_FALSE, PETSC_TRUE, partitioning);
229: return 0;
230: }
232: PetscErrorCode MatPartitioningView_Parmetis(MatPartitioning part,PetscViewer viewer)
233: {
234: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis*)part->data;
235: PetscMPIInt rank;
236: PetscBool iascii;
238: MPI_Comm_rank(PetscObjectComm((PetscObject)part),&rank);
239: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
240: if (iascii) {
241: if (pmetis->parallel == 2) {
242: PetscViewerASCIIPrintf(viewer," Using parallel coarse grid partitioner\n");
243: } else {
244: PetscViewerASCIIPrintf(viewer," Using sequential coarse grid partitioner\n");
245: }
246: PetscViewerASCIIPrintf(viewer," Using %" PetscInt_FMT " fold factor\n",pmetis->foldfactor);
247: PetscViewerASCIIPushSynchronized(viewer);
248: PetscViewerASCIISynchronizedPrintf(viewer," [%d]Number of cuts found %" PetscInt_FMT "\n",rank,pmetis->cuts);
249: PetscViewerFlush(viewer);
250: PetscViewerASCIIPopSynchronized(viewer);
251: }
252: return 0;
253: }
255: /*@
256: MatPartitioningParmetisSetCoarseSequential - Use the sequential code to
257: do the partitioning of the coarse grid.
259: Logically Collective on MatPartitioning
261: Input Parameter:
262: . part - the partitioning context
264: Level: advanced
266: @*/
267: PetscErrorCode MatPartitioningParmetisSetCoarseSequential(MatPartitioning part)
268: {
269: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis*)part->data;
271: pmetis->parallel = 1;
272: return 0;
273: }
275: /*@
276: MatPartitioningParmetisSetRepartition - Repartition
277: current mesh to rebalance computation.
279: Logically Collective on MatPartitioning
281: Input Parameter:
282: . part - the partitioning context
284: Level: advanced
286: @*/
287: PetscErrorCode MatPartitioningParmetisSetRepartition(MatPartitioning part)
288: {
289: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis*)part->data;
291: pmetis->repartition = PETSC_TRUE;
292: return 0;
293: }
295: /*@
296: MatPartitioningParmetisGetEdgeCut - Returns the number of edge cuts in the vertex partition.
298: Input Parameter:
299: . part - the partitioning context
301: Output Parameter:
302: . cut - the edge cut
304: Level: advanced
306: @*/
307: PetscErrorCode MatPartitioningParmetisGetEdgeCut(MatPartitioning part, PetscInt *cut)
308: {
309: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis*) part->data;
311: *cut = pmetis->cuts;
312: return 0;
313: }
315: PetscErrorCode MatPartitioningSetFromOptions_Parmetis(PetscOptionItems *PetscOptionsObject,MatPartitioning part)
316: {
317: PetscBool flag = PETSC_FALSE;
319: PetscOptionsHead(PetscOptionsObject,"Set ParMeTiS partitioning options");
320: PetscOptionsBool("-mat_partitioning_parmetis_coarse_sequential","Use sequential coarse partitioner","MatPartitioningParmetisSetCoarseSequential",flag,&flag,NULL);
321: if (flag) {
322: MatPartitioningParmetisSetCoarseSequential(part);
323: }
324: PetscOptionsBool("-mat_partitioning_parmetis_repartition","","MatPartitioningParmetisSetRepartition",flag,&flag,NULL);
325: if (flag) {
326: MatPartitioningParmetisSetRepartition(part);
327: }
328: PetscOptionsTail();
329: return 0;
330: }
332: PetscErrorCode MatPartitioningDestroy_Parmetis(MatPartitioning part)
333: {
334: MatPartitioning_Parmetis *pmetis = (MatPartitioning_Parmetis*)part->data;
336: PetscFree(pmetis);
337: return 0;
338: }
340: /*MC
341: MATPARTITIONINGPARMETIS - Creates a partitioning context via the external package PARMETIS.
343: Collective
345: Input Parameter:
346: . part - the partitioning context
348: Options Database Keys:
349: . -mat_partitioning_parmetis_coarse_sequential - use sequential PARMETIS coarse partitioner
351: Level: beginner
353: Notes:
354: See https://www-users.cs.umn.edu/~karypis/metis/
356: .seealso: MatPartitioningSetType(), MatPartitioningType
358: M*/
360: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Parmetis(MatPartitioning part)
361: {
362: MatPartitioning_Parmetis *pmetis;
364: PetscNewLog(part,&pmetis);
365: part->data = (void*)pmetis;
367: pmetis->cuts = 0; /* output variable */
368: pmetis->foldfactor = 150; /*folding factor */
369: pmetis->parallel = 2; /* use parallel partitioner for coarse grid */
370: pmetis->indexing = 0; /* index numbering starts from 0 */
371: pmetis->printout = 0; /* print no output while running */
372: pmetis->repartition = PETSC_FALSE;
374: part->ops->apply = MatPartitioningApply_Parmetis;
375: part->ops->applynd = MatPartitioningApplyND_Parmetis;
376: part->ops->improve = MatPartitioningImprove_Parmetis;
377: part->ops->view = MatPartitioningView_Parmetis;
378: part->ops->destroy = MatPartitioningDestroy_Parmetis;
379: part->ops->setfromoptions = MatPartitioningSetFromOptions_Parmetis;
380: return 0;
381: }
383: /*@
384: MatMeshToVertexGraph - This routine does not exist because ParMETIS does not provide the functionality. Uses the ParMETIS package to
385: convert a Mat that represents a mesh to a Mat the represents the graph of the coupling
386: between vertices of the cells and is suitable for partitioning with the MatPartitioning object. Use this to partition
387: vertices of a mesh. More likely you should use MatMeshToCellGraph()
389: Collective on Mat
391: Input Parameters:
392: + mesh - the graph that represents the mesh
393: - ncommonnodes - mesh elements that share this number of common nodes are considered neighbors, use 2 for triangles and
394: quadrilaterials, 3 for tetrahedrals and 4 for hexahedrals
396: Output Parameter:
397: . dual - the dual graph
399: Notes:
400: Currently requires ParMetis to be installed and uses ParMETIS_V3_Mesh2Dual()
402: The columns of each row of the Mat mesh are the global vertex numbers of the vertices of that rows cell. The number of rows in mesh is
403: number of cells, the number of columns is the number of vertices.
405: Level: advanced
407: .seealso: MatMeshToCellGraph(), MatCreateMPIAdj(), MatPartitioningCreate()
409: @*/
410: PetscErrorCode MatMeshToVertexGraph(Mat mesh,PetscInt ncommonnodes,Mat *dual)
411: {
412: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"ParMETIS does not provide this functionality");
413: }
415: /*@
416: MatMeshToCellGraph - Uses the ParMETIS package to convert a Mat that represents a mesh to a Mat the represents the graph of the coupling
417: between cells (the "dual" graph) and is suitable for partitioning with the MatPartitioning object. Use this to partition
418: cells of a mesh.
420: Collective on Mat
422: Input Parameters:
423: + mesh - the graph that represents the mesh
424: - ncommonnodes - mesh elements that share this number of common nodes are considered neighbors, use 2 for triangles and
425: quadrilaterials, 3 for tetrahedrals and 4 for hexahedrals
427: Output Parameter:
428: . dual - the dual graph
430: Notes:
431: Currently requires ParMetis to be installed and uses ParMETIS_V3_Mesh2Dual()
433: $ Each row of the mesh object represents a single cell in the mesh. For triangles it has 3 entries, quadrilaterials 4 entries,
434: $ tetrahedrals 4 entries and hexahedrals 8 entries. You can mix triangles and quadrilaterals in the same mesh, but cannot
435: $ mix tetrahedrals and hexahedrals
436: $ The columns of each row of the Mat mesh are the global vertex numbers of the vertices of that row's cell.
437: $ The number of rows in mesh is number of cells, the number of columns is the number of vertices.
439: Level: advanced
441: .seealso: MatMeshToVertexGraph(), MatCreateMPIAdj(), MatPartitioningCreate()
443: @*/
444: PetscErrorCode MatMeshToCellGraph(Mat mesh,PetscInt ncommonnodes,Mat *dual)
445: {
446: PetscInt *newxadj,*newadjncy;
447: PetscInt numflag=0;
448: Mat_MPIAdj *adj = (Mat_MPIAdj*)mesh->data,*newadj;
449: PetscBool flg;
450: MPI_Comm comm;
452: PetscObjectTypeCompare((PetscObject)mesh,MATMPIADJ,&flg);
455: PetscObjectGetComm((PetscObject)mesh,&comm);
456: PetscStackCallParmetis(ParMETIS_V3_Mesh2Dual,((idx_t*)mesh->rmap->range,(idx_t*)adj->i,(idx_t*)adj->j,(idx_t*)&numflag,(idx_t*)&ncommonnodes,(idx_t**)&newxadj,(idx_t**)&newadjncy,&comm));
457: MatCreateMPIAdj(PetscObjectComm((PetscObject)mesh),mesh->rmap->n,mesh->rmap->N,newxadj,newadjncy,NULL,dual);
458: newadj = (Mat_MPIAdj*)(*dual)->data;
460: newadj->freeaijwithfree = PETSC_TRUE; /* signal the matrix should be freed with system free since space was allocated by ParMETIS */
461: return 0;
462: }