Actual source code: hierarchical.c


  2: #include <../src/mat/impls/adj/mpi/mpiadj.h>
  3: #include <petscsf.h>
  4: #include <petsc/private/matimpl.h>

  6: /*
  7:   It is a hierarchical partitioning. The partitioner has two goals:
  8:   (1) Most of current partitioners fail at a large scale. The hierarchical partitioning
  9:   strategy is trying to produce large number of subdomains when number of processor cores is large.
 10:   (2) PCGASM needs one 'big' subdomain across multi-cores. The partitioner provides two
 11:   consistent partitions, coarse parts and fine parts. A coarse part is a 'big' subdomain consisting
 12:   of several small subdomains.
 13: */

 15: PetscErrorCode MatPartitioningHierarchical_DetermineDestination(MatPartitioning, IS, PetscInt, PetscInt, IS *);
 16: PetscErrorCode MatPartitioningHierarchical_AssembleSubdomain(Mat, IS, IS, IS *, Mat *, ISLocalToGlobalMapping *);
 17: PetscErrorCode MatPartitioningHierarchical_ReassembleFineparts(Mat, IS, ISLocalToGlobalMapping, IS *);

 19: typedef struct {
 20:   char           *fineparttype;   /* partitioner on fine level */
 21:   char           *coarseparttype; /* partitioner on coarse level */
 22:   PetscInt        nfineparts;     /* number of fine parts on each coarse subdomain */
 23:   PetscInt        ncoarseparts;   /* number of coarse parts */
 24:   IS              coarseparts;    /* partitioning on coarse level */
 25:   IS              fineparts;      /* partitioning on fine level */
 26:   MatPartitioning coarseMatPart;  /* MatPartititioning on coarse level (first level) */
 27:   MatPartitioning fineMatPart;    /* MatPartitioning on fine level (second level) */
 28:   MatPartitioning improver;       /* Improve the quality of a partition */
 29: } MatPartitioning_Hierarchical;

 31: /*
 32:    Uses a hierarchical partitioning strategy to partition the matrix in parallel.
 33:    Use this interface to make the partitioner consistent with others
 34: */
 35: static PetscErrorCode MatPartitioningApply_Hierarchical(MatPartitioning part, IS *partitioning)
 36: {
 37:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;
 38:   const PetscInt               *fineparts_indices, *coarseparts_indices;
 39:   PetscInt                     *fineparts_indices_tmp;
 40:   PetscInt                     *parts_indices, i, j, mat_localsize, *offsets;
 41:   Mat                           mat = part->adj, adj, sadj;
 42:   PetscReal                    *part_weights;
 43:   PetscBool                     flg;
 44:   PetscInt                      bs                    = 1;
 45:   PetscInt                     *coarse_vertex_weights = NULL;
 46:   PetscMPIInt                   size, rank;
 47:   MPI_Comm                      comm, scomm;
 48:   IS                            destination, fineparts_temp, vweights, svweights;
 49:   PetscInt                      nsvwegihts, *fp_vweights;
 50:   const PetscInt               *svweights_indices;
 51:   ISLocalToGlobalMapping        mapping;
 52:   const char                   *prefix;
 53:   PetscBool                     use_edge_weights;

 55:   PetscFunctionBegin;
 56:   PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
 57:   PetscCallMPI(MPI_Comm_size(comm, &size));
 58:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 59:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIADJ, &flg));
 60:   if (flg) {
 61:     adj = mat;
 62:     PetscCall(PetscObjectReference((PetscObject)adj));
 63:   } else {
 64:     /* bs indicates if the converted matrix is "reduced" from the original and hence the
 65:        resulting partition results need to be stretched to match the original matrix */
 66:     PetscCall(MatConvert(mat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
 67:     if (adj->rmap->n > 0) bs = mat->rmap->n / adj->rmap->n;
 68:   }
 69:   /* local size of mat */
 70:   mat_localsize = adj->rmap->n;
 71:   /* check parameters */
 72:   /* how many small subdomains we want from a given 'big' suddomain */
 73:   PetscCheck(hpart->nfineparts, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, " must set number of small subdomains for each big subdomain ");
 74:   PetscCheck(hpart->ncoarseparts || part->n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, " did not either set number of coarse parts or total number of parts ");

 76:   /* Partitioning the domain into one single subdomain is a trivial case, and we should just return  */
 77:   if (part->n == 1) {
 78:     PetscCall(PetscCalloc1(bs * adj->rmap->n, &parts_indices));
 79:     PetscCall(ISCreateGeneral(comm, bs * adj->rmap->n, parts_indices, PETSC_OWN_POINTER, partitioning));
 80:     hpart->ncoarseparts = 1;
 81:     hpart->nfineparts   = 1;
 82:     PetscCall(PetscStrallocpy("NONE", &hpart->coarseparttype));
 83:     PetscCall(PetscStrallocpy("NONE", &hpart->fineparttype));
 84:     PetscCall(MatDestroy(&adj));
 85:     PetscFunctionReturn(PETSC_SUCCESS);
 86:   }

 88:   if (part->n) {
 89:     hpart->ncoarseparts = part->n / hpart->nfineparts;

 91:     if (part->n % hpart->nfineparts != 0) hpart->ncoarseparts++;
 92:   } else {
 93:     part->n = hpart->ncoarseparts * hpart->nfineparts;
 94:   }

 96:   PetscCall(PetscMalloc1(hpart->ncoarseparts + 1, &offsets));
 97:   PetscCall(PetscMalloc1(hpart->ncoarseparts, &part_weights));

 99:   offsets[0] = 0;
100:   if (part->n % hpart->nfineparts != 0) offsets[1] = part->n % hpart->nfineparts;
101:   else offsets[1] = hpart->nfineparts;

103:   part_weights[0] = ((PetscReal)offsets[1]) / part->n;

105:   for (i = 2; i <= hpart->ncoarseparts; i++) {
106:     offsets[i]          = hpart->nfineparts;
107:     part_weights[i - 1] = ((PetscReal)offsets[i]) / part->n;
108:   }

110:   offsets[0] = 0;
111:   for (i = 1; i <= hpart->ncoarseparts; i++) offsets[i] += offsets[i - 1];

113:   /* If these exists a mat partitioner, we should delete it */
114:   PetscCall(MatPartitioningDestroy(&hpart->coarseMatPart));
115:   PetscCall(MatPartitioningCreate(comm, &hpart->coarseMatPart));
116:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)part, &prefix));
117:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)hpart->coarseMatPart, prefix));
118:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)hpart->coarseMatPart, "hierarch_coarse_"));
119:   /* if did not set partitioning type yet, use parmetis by default */
120:   if (!hpart->coarseparttype) {
121: #if defined(PETSC_HAVE_PARMETIS)
122:     PetscCall(MatPartitioningSetType(hpart->coarseMatPart, MATPARTITIONINGPARMETIS));
123:     PetscCall(PetscStrallocpy(MATPARTITIONINGPARMETIS, &hpart->coarseparttype));
124: #elif defined(PETSC_HAVE_PTSCOTCH)
125:     PetscCall(MatPartitioningSetType(hpart->coarseMatPart, MATPARTITIONINGPTSCOTCH));
126:     PetscCall(PetscStrallocpy(MATPARTITIONINGPTSCOTCH, &hpart->coarseparttype));
127: #else
128:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Requires PETSc be installed with ParMetis or run with -mat_partitioning_hierarchical_coarseparttype partitiontype");
129: #endif
130:   } else {
131:     PetscCall(MatPartitioningSetType(hpart->coarseMatPart, hpart->coarseparttype));
132:   }
133:   PetscCall(MatPartitioningSetAdjacency(hpart->coarseMatPart, adj));
134:   PetscCall(MatPartitioningSetNParts(hpart->coarseMatPart, hpart->ncoarseparts));
135:   /* copy over vertex weights */
136:   if (part->vertex_weights) {
137:     PetscCall(PetscMalloc1(mat_localsize, &coarse_vertex_weights));
138:     PetscCall(PetscArraycpy(coarse_vertex_weights, part->vertex_weights, mat_localsize));
139:     PetscCall(MatPartitioningSetVertexWeights(hpart->coarseMatPart, coarse_vertex_weights));
140:   }
141:   /* Copy use_edge_weights flag from part to coarse part */
142:   PetscCall(MatPartitioningGetUseEdgeWeights(part, &use_edge_weights));
143:   PetscCall(MatPartitioningSetUseEdgeWeights(hpart->coarseMatPart, use_edge_weights));

145:   PetscCall(MatPartitioningSetPartitionWeights(hpart->coarseMatPart, part_weights));
146:   PetscCall(MatPartitioningApply(hpart->coarseMatPart, &hpart->coarseparts));

148:   /* Wrap the original vertex weights into an index set so that we can extract the corresponding
149:    * vertex weights for each big subdomain using ISCreateSubIS().
150:    * */
151:   if (part->vertex_weights) PetscCall(ISCreateGeneral(comm, mat_localsize, part->vertex_weights, PETSC_COPY_VALUES, &vweights));

153:   PetscCall(PetscCalloc1(mat_localsize, &fineparts_indices_tmp));
154:   for (i = 0; i < hpart->ncoarseparts; i += size) {
155:     /* Determine where we want to send big subdomains */
156:     PetscCall(MatPartitioningHierarchical_DetermineDestination(part, hpart->coarseparts, i, i + size, &destination));
157:     /* Assemble a submatrix and its vertex weights for partitioning subdomains  */
158:     PetscCall(MatPartitioningHierarchical_AssembleSubdomain(adj, part->vertex_weights ? vweights : NULL, destination, part->vertex_weights ? &svweights : NULL, &sadj, &mapping));
159:     /* We have to create a new array to hold vertex weights since coarse partitioner needs to own the vertex-weights array */
160:     if (part->vertex_weights) {
161:       PetscCall(ISGetLocalSize(svweights, &nsvwegihts));
162:       PetscCall(PetscMalloc1(nsvwegihts, &fp_vweights));
163:       PetscCall(ISGetIndices(svweights, &svweights_indices));
164:       PetscCall(PetscArraycpy(fp_vweights, svweights_indices, nsvwegihts));
165:       PetscCall(ISRestoreIndices(svweights, &svweights_indices));
166:       PetscCall(ISDestroy(&svweights));
167:     }

169:     PetscCall(ISDestroy(&destination));
170:     PetscCall(PetscObjectGetComm((PetscObject)sadj, &scomm));

172:     /*
173:      * If the number of big subdomains is smaller than the number of processor cores, the higher ranks do not
174:      * need to do partitioning
175:      * */
176:     if ((i + rank) < hpart->ncoarseparts) {
177:       PetscCall(MatPartitioningDestroy(&hpart->fineMatPart));
178:       /* create a fine partitioner */
179:       PetscCall(MatPartitioningCreate(scomm, &hpart->fineMatPart));
180:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)hpart->fineMatPart, prefix));
181:       PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)hpart->fineMatPart, "hierarch_fine_"));
182:       /* if do not set partitioning type, use parmetis by default */
183:       if (!hpart->fineparttype) {
184: #if defined(PETSC_HAVE_PARMETIS)
185:         PetscCall(MatPartitioningSetType(hpart->fineMatPart, MATPARTITIONINGPARMETIS));
186:         PetscCall(PetscStrallocpy(MATPARTITIONINGPARMETIS, &hpart->fineparttype));
187: #elif defined(PETSC_HAVE_PTSCOTCH)
188:         PetscCall(MatPartitioningSetType(hpart->fineMatPart, MATPARTITIONINGPTSCOTCH));
189:         PetscCall(PetscStrallocpy(MATPARTITIONINGPTSCOTCH, &hpart->fineparttype));
190: #elif defined(PETSC_HAVE_CHACO)
191:         PetscCall(MatPartitioningSetType(hpart->fineMatPart, MATPARTITIONINGCHACO));
192:         PetscCall(PetscStrallocpy(MATPARTITIONINGCHACO, &hpart->fineparttype));
193: #elif defined(PETSC_HAVE_PARTY)
194:         PetscCall(MatPartitioningSetType(hpart->fineMatPart, MATPARTITIONINGPARTY));
195:         PetscCall(PetscStrallocpy(PETSC_HAVE_PARTY, &hpart->fineparttype));
196: #else
197:         SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Requires PETSc be installed with ParMetis or run with -mat_partitioning_hierarchical_coarseparttype partitiontype");
198: #endif
199:       } else {
200:         PetscCall(MatPartitioningSetType(hpart->fineMatPart, hpart->fineparttype));
201:       }
202:       PetscCall(MatPartitioningSetUseEdgeWeights(hpart->fineMatPart, use_edge_weights));
203:       PetscCall(MatPartitioningSetAdjacency(hpart->fineMatPart, sadj));
204:       PetscCall(MatPartitioningSetNParts(hpart->fineMatPart, offsets[rank + 1 + i] - offsets[rank + i]));
205:       if (part->vertex_weights) PetscCall(MatPartitioningSetVertexWeights(hpart->fineMatPart, fp_vweights));
206:       PetscCall(MatPartitioningApply(hpart->fineMatPart, &fineparts_temp));
207:     } else {
208:       PetscCall(ISCreateGeneral(scomm, 0, NULL, PETSC_OWN_POINTER, &fineparts_temp));
209:     }

211:     PetscCall(MatDestroy(&sadj));

213:     /* Send partition back to the original owners */
214:     PetscCall(MatPartitioningHierarchical_ReassembleFineparts(adj, fineparts_temp, mapping, &hpart->fineparts));
215:     PetscCall(ISGetIndices(hpart->fineparts, &fineparts_indices));
216:     for (j = 0; j < mat_localsize; j++)
217:       if (fineparts_indices[j] >= 0) fineparts_indices_tmp[j] = fineparts_indices[j];

219:     PetscCall(ISRestoreIndices(hpart->fineparts, &fineparts_indices));
220:     PetscCall(ISDestroy(&hpart->fineparts));
221:     PetscCall(ISDestroy(&fineparts_temp));
222:     PetscCall(ISLocalToGlobalMappingDestroy(&mapping));
223:   }

225:   if (part->vertex_weights) PetscCall(ISDestroy(&vweights));

227:   PetscCall(ISCreateGeneral(comm, mat_localsize, fineparts_indices_tmp, PETSC_OWN_POINTER, &hpart->fineparts));
228:   PetscCall(ISGetIndices(hpart->fineparts, &fineparts_indices));
229:   PetscCall(ISGetIndices(hpart->coarseparts, &coarseparts_indices));
230:   PetscCall(PetscMalloc1(bs * adj->rmap->n, &parts_indices));
231:   /* Modify the local indices to the global indices by combing the coarse partition and the fine partitions */
232:   for (i = 0; i < adj->rmap->n; i++) {
233:     for (j = 0; j < bs; j++) parts_indices[bs * i + j] = fineparts_indices[i] + offsets[coarseparts_indices[i]];
234:   }
235:   PetscCall(ISRestoreIndices(hpart->fineparts, &fineparts_indices));
236:   PetscCall(ISRestoreIndices(hpart->coarseparts, &coarseparts_indices));
237:   PetscCall(PetscFree(offsets));
238:   PetscCall(ISCreateGeneral(comm, bs * adj->rmap->n, parts_indices, PETSC_OWN_POINTER, partitioning));
239:   PetscCall(MatDestroy(&adj));
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: PetscErrorCode MatPartitioningHierarchical_ReassembleFineparts(Mat adj, IS fineparts, ISLocalToGlobalMapping mapping, IS *sfineparts)
244: {
245:   PetscInt       *local_indices, *global_indices, *sfineparts_indices, localsize, i;
246:   const PetscInt *ranges, *fineparts_indices;
247:   PetscMPIInt     rank, *owners;
248:   MPI_Comm        comm;
249:   PetscLayout     rmap;
250:   PetscSFNode    *remote;
251:   PetscSF         sf;

253:   PetscFunctionBegin;
255:   PetscCall(PetscObjectGetComm((PetscObject)adj, &comm));
256:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
257:   PetscCall(MatGetLayouts(adj, &rmap, NULL));
258:   PetscCall(ISGetLocalSize(fineparts, &localsize));
259:   PetscCall(PetscMalloc2(localsize, &global_indices, localsize, &local_indices));
260:   for (i = 0; i < localsize; i++) local_indices[i] = i;
261:   /* map local indices back to global so that we can permulate data globally */
262:   PetscCall(ISLocalToGlobalMappingApply(mapping, localsize, local_indices, global_indices));
263:   PetscCall(PetscCalloc1(localsize, &owners));
264:   /* find owners for global indices */
265:   for (i = 0; i < localsize; i++) PetscCall(PetscLayoutFindOwner(rmap, global_indices[i], &owners[i]));
266:   PetscCall(PetscLayoutGetRanges(rmap, &ranges));
267:   PetscCall(PetscMalloc1(ranges[rank + 1] - ranges[rank], &sfineparts_indices));

269:   for (i = 0; i < ranges[rank + 1] - ranges[rank]; i++) sfineparts_indices[i] = -1;

271:   PetscCall(ISGetIndices(fineparts, &fineparts_indices));
272:   PetscCall(PetscSFCreate(comm, &sf));
273:   PetscCall(PetscMalloc1(localsize, &remote));
274:   for (i = 0; i < localsize; i++) {
275:     remote[i].rank  = owners[i];
276:     remote[i].index = global_indices[i] - ranges[owners[i]];
277:   }
278:   PetscCall(PetscSFSetType(sf, PETSCSFBASIC));
279:   /* not sure how to add prefix to sf */
280:   PetscCall(PetscSFSetFromOptions(sf));
281:   PetscCall(PetscSFSetGraph(sf, localsize, localsize, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
282:   PetscCall(PetscSFReduceBegin(sf, MPIU_INT, fineparts_indices, sfineparts_indices, MPI_REPLACE));
283:   PetscCall(PetscSFReduceEnd(sf, MPIU_INT, fineparts_indices, sfineparts_indices, MPI_REPLACE));
284:   PetscCall(PetscSFDestroy(&sf));
285:   PetscCall(ISRestoreIndices(fineparts, &fineparts_indices));
286:   PetscCall(ISCreateGeneral(comm, ranges[rank + 1] - ranges[rank], sfineparts_indices, PETSC_OWN_POINTER, sfineparts));
287:   PetscCall(PetscFree2(global_indices, local_indices));
288:   PetscCall(PetscFree(owners));
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: PetscErrorCode MatPartitioningHierarchical_AssembleSubdomain(Mat adj, IS vweights, IS destination, IS *svweights, Mat *sadj, ISLocalToGlobalMapping *mapping)
293: {
294:   IS              irows, icols;
295:   PetscInt        irows_ln;
296:   PetscMPIInt     rank;
297:   const PetscInt *irows_indices;
298:   MPI_Comm        comm;

300:   PetscFunctionBegin;
301:   PetscCall(PetscObjectGetComm((PetscObject)adj, &comm));
302:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
303:   /* figure out where data comes from  */
304:   PetscCall(ISBuildTwoSided(destination, NULL, &irows));
305:   PetscCall(ISDuplicate(irows, &icols));
306:   PetscCall(ISGetLocalSize(irows, &irows_ln));
307:   PetscCall(ISGetIndices(irows, &irows_indices));
308:   PetscCall(ISLocalToGlobalMappingCreate(comm, 1, irows_ln, irows_indices, PETSC_COPY_VALUES, mapping));
309:   PetscCall(ISRestoreIndices(irows, &irows_indices));
310:   PetscCall(MatCreateSubMatrices(adj, 1, &irows, &icols, MAT_INITIAL_MATRIX, &sadj));
311:   if (vweights && svweights) PetscCall(ISCreateSubIS(vweights, irows, svweights));
312:   PetscCall(ISDestroy(&irows));
313:   PetscCall(ISDestroy(&icols));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: PetscErrorCode MatPartitioningHierarchical_DetermineDestination(MatPartitioning part, IS partitioning, PetscInt pstart, PetscInt pend, IS *destination)
318: {
319:   MPI_Comm        comm;
320:   PetscMPIInt     rank, size, target;
321:   PetscInt        plocalsize, *dest_indices, i;
322:   const PetscInt *part_indices;

324:   PetscFunctionBegin;
325:   PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
326:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
327:   PetscCallMPI(MPI_Comm_size(comm, &size));
328:   PetscCheck((pend - pstart) <= size, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "range [%" PetscInt_FMT ", %" PetscInt_FMT "] should be smaller than or equal to size %d", pstart, pend, size);
329:   PetscCheck(pstart <= pend, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, " pstart %" PetscInt_FMT " should be smaller than pend %" PetscInt_FMT, pstart, pend);
330:   PetscCall(ISGetLocalSize(partitioning, &plocalsize));
331:   PetscCall(PetscMalloc1(plocalsize, &dest_indices));
332:   PetscCall(ISGetIndices(partitioning, &part_indices));
333:   for (i = 0; i < plocalsize; i++) {
334:     /* compute target */
335:     target = part_indices[i] - pstart;
336:     /* mark out of range entity as -1 */
337:     if (part_indices[i] < pstart || part_indices[i] >= pend) target = -1;
338:     dest_indices[i] = target;
339:   }
340:   PetscCall(ISCreateGeneral(comm, plocalsize, dest_indices, PETSC_OWN_POINTER, destination));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }

344: PetscErrorCode MatPartitioningView_Hierarchical(MatPartitioning part, PetscViewer viewer)
345: {
346:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;
347:   PetscMPIInt                   rank;
348:   PetscBool                     iascii;
349:   PetscViewer                   sviewer;

351:   PetscFunctionBegin;
352:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)part), &rank));
353:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
354:   if (iascii) {
355:     PetscCall(PetscViewerASCIIPrintf(viewer, " Number of coarse parts: %" PetscInt_FMT "\n", hpart->ncoarseparts));
356:     PetscCall(PetscViewerASCIIPrintf(viewer, " Coarse partitioner: %s\n", hpart->coarseparttype));
357:     if (hpart->coarseMatPart) {
358:       PetscCall(PetscViewerASCIIPushTab(viewer));
359:       PetscCall(MatPartitioningView(hpart->coarseMatPart, viewer));
360:       PetscCall(PetscViewerASCIIPopTab(viewer));
361:     }
362:     PetscCall(PetscViewerASCIIPrintf(viewer, " Number of fine parts: %" PetscInt_FMT "\n", hpart->nfineparts));
363:     PetscCall(PetscViewerASCIIPrintf(viewer, " Fine partitioner: %s\n", hpart->fineparttype));
364:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
365:     if (rank == 0 && hpart->fineMatPart) {
366:       PetscCall(PetscViewerASCIIPushTab(viewer));
367:       PetscCall(MatPartitioningView(hpart->fineMatPart, sviewer));
368:       PetscCall(PetscViewerASCIIPopTab(viewer));
369:     }
370:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
371:   }
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: PetscErrorCode MatPartitioningHierarchicalGetFineparts(MatPartitioning part, IS *fineparts)
376: {
377:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;

379:   PetscFunctionBegin;
380:   *fineparts = hpart->fineparts;
381:   PetscCall(PetscObjectReference((PetscObject)hpart->fineparts));
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: PetscErrorCode MatPartitioningHierarchicalGetCoarseparts(MatPartitioning part, IS *coarseparts)
386: {
387:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;

389:   PetscFunctionBegin;
390:   *coarseparts = hpart->coarseparts;
391:   PetscCall(PetscObjectReference((PetscObject)hpart->coarseparts));
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: PetscErrorCode MatPartitioningHierarchicalSetNcoarseparts(MatPartitioning part, PetscInt ncoarseparts)
396: {
397:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;

399:   PetscFunctionBegin;
400:   hpart->ncoarseparts = ncoarseparts;
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: PetscErrorCode MatPartitioningHierarchicalSetNfineparts(MatPartitioning part, PetscInt nfineparts)
405: {
406:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;

408:   PetscFunctionBegin;
409:   hpart->nfineparts = nfineparts;
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: PetscErrorCode MatPartitioningSetFromOptions_Hierarchical(MatPartitioning part, PetscOptionItems *PetscOptionsObject)
414: {
415:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;
416:   char                          value[1024];
417:   PetscBool                     flag = PETSC_FALSE;

419:   PetscFunctionBegin;
420:   PetscOptionsHeadBegin(PetscOptionsObject, "Set hierarchical partitioning options");
421:   PetscCall(PetscOptionsString("-mat_partitioning_hierarchical_coarseparttype", "coarse part type", NULL, NULL, value, sizeof(value), &flag));
422:   if (flag) {
423:     PetscCall(PetscFree(hpart->coarseparttype));
424:     PetscCall(PetscStrallocpy(value, &hpart->coarseparttype));
425:   }
426:   PetscCall(PetscOptionsString("-mat_partitioning_hierarchical_fineparttype", "fine part type", NULL, NULL, value, sizeof(value), &flag));
427:   if (flag) {
428:     PetscCall(PetscFree(hpart->fineparttype));
429:     PetscCall(PetscStrallocpy(value, &hpart->fineparttype));
430:   }
431:   PetscCall(PetscOptionsInt("-mat_partitioning_hierarchical_ncoarseparts", "number of coarse parts", NULL, hpart->ncoarseparts, &hpart->ncoarseparts, &flag));
432:   PetscCall(PetscOptionsInt("-mat_partitioning_hierarchical_nfineparts", "number of fine parts", NULL, hpart->nfineparts, &hpart->nfineparts, &flag));
433:   PetscOptionsHeadEnd();
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

437: PetscErrorCode MatPartitioningDestroy_Hierarchical(MatPartitioning part)
438: {
439:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;

441:   PetscFunctionBegin;
442:   if (hpart->coarseparttype) PetscCall(PetscFree(hpart->coarseparttype));
443:   if (hpart->fineparttype) PetscCall(PetscFree(hpart->fineparttype));
444:   PetscCall(ISDestroy(&hpart->fineparts));
445:   PetscCall(ISDestroy(&hpart->coarseparts));
446:   PetscCall(MatPartitioningDestroy(&hpart->coarseMatPart));
447:   PetscCall(MatPartitioningDestroy(&hpart->fineMatPart));
448:   PetscCall(MatPartitioningDestroy(&hpart->improver));
449:   PetscCall(PetscFree(hpart));
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

453: /*
454:    Improves the quality  of a partition
455: */
456: static PetscErrorCode MatPartitioningImprove_Hierarchical(MatPartitioning part, IS *partitioning)
457: {
458:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical *)part->data;
459:   Mat                           mat   = part->adj, adj;
460:   PetscBool                     flg;
461:   const char                   *prefix;
462: #if defined(PETSC_HAVE_PARMETIS)
463:   PetscInt *vertex_weights;
464: #endif

466:   PetscFunctionBegin;
467:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIADJ, &flg));
468:   if (flg) {
469:     adj = mat;
470:     PetscCall(PetscObjectReference((PetscObject)adj));
471:   } else {
472:     /* bs indicates if the converted matrix is "reduced" from the original and hence the
473:        resulting partition results need to be stretched to match the original matrix */
474:     PetscCall(MatConvert(mat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
475:   }

477:   /* If there exists a mat partitioner, we should delete it */
478:   PetscCall(MatPartitioningDestroy(&hpart->improver));
479:   PetscCall(MatPartitioningCreate(PetscObjectComm((PetscObject)part), &hpart->improver));
480:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)part, &prefix));
481:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)hpart->improver, prefix));
482:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)hpart->improver, "hierarch_improver_"));
483:   /* Only parmetis supports to refine a partition */
484: #if defined(PETSC_HAVE_PARMETIS)
485:   PetscCall(MatPartitioningSetType(hpart->improver, MATPARTITIONINGPARMETIS));
486:   PetscCall(MatPartitioningSetAdjacency(hpart->improver, adj));
487:   PetscCall(MatPartitioningSetNParts(hpart->improver, part->n));
488:   /* copy over vertex weights */
489:   if (part->vertex_weights) {
490:     PetscCall(PetscMalloc1(adj->rmap->n, &vertex_weights));
491:     PetscCall(PetscArraycpy(vertex_weights, part->vertex_weights, adj->rmap->n));
492:     PetscCall(MatPartitioningSetVertexWeights(hpart->improver, vertex_weights));
493:   }
494:   PetscCall(MatPartitioningImprove(hpart->improver, partitioning));
495:   PetscCall(MatDestroy(&adj));
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: #else
498:   SETERRQ(PetscObjectComm((PetscObject)adj), PETSC_ERR_SUP, "Requires PETSc be installed with ParMetis");
499: #endif
500: }

502: /*MC
503:    MATPARTITIONINGHIERARCH - Creates a partitioning context via hierarchical partitioning strategy.
504:    The graph is partitioned into a number of subgraphs, and each subgraph is further split into a few smaller
505:    subgraphs. The idea can be applied in a recursive manner. It is useful when you want to partition the graph
506:    into a large number of subgraphs (often more than 10K) since partitions obtained with existing partitioners
507:    such as ParMETIS and PTScotch are far from ideal. The hierarchical partitioning also tries to avoid off-node
508:    communication as much as possible for multi-core processor. Another user case for the hierarchical partitioning
509:    is to improve `PCGASM` convergence by generating multi-rank connected subdomain.

511:    Collective

513:    Input Parameter:
514: .  part - the partitioning context

516:    Options Database Keys:
517: +     -mat_partitioning_hierarchical_coarseparttype - partitioner type at the first level and parmetis is used by default
518: .     -mat_partitioning_hierarchical_fineparttype - partitioner type at the second level and parmetis is used by default
519: .     -mat_partitioning_hierarchical_ncoarseparts - number of subgraphs is required at the first level, which is often the number of compute nodes
520: -     -mat_partitioning_hierarchical_nfineparts - number of smaller subgraphs for each subgraph, which is often the number of cores per compute node

522:    Level: beginner

524:    References:
525: +  * - Fande Kong, Xiao-Chuan Cai, A highly scalable multilevel Schwarz method with boundary geometry preserving coarse spaces for 3D elasticity
526:       problems on domains with complex geometry,   SIAM Journal on Scientific Computing 38 (2), C73-C95, 2016
527: -  * - Fande Kong, Roy H. Stogner, Derek Gaston, John W. Peterson, Cody J. Permann, Andrew E. Slaughter, and Richard C. Martineau,
528:       A general-purpose hierarchical mesh partitioning method with node balancing strategies for large-scale numerical simulations,
529:       arXiv preprint arXiv:1809.02666CoRR, 2018.

531: .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MATPARTITIONINGMETIS`, `MATPARTITIONINGPARMETIS`,
532: M*/

534: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Hierarchical(MatPartitioning part)
535: {
536:   MatPartitioning_Hierarchical *hpart;

538:   PetscFunctionBegin;
539:   PetscCall(PetscNew(&hpart));
540:   part->data = (void *)hpart;

542:   hpart->fineparttype   = NULL; /* fine level (second) partitioner */
543:   hpart->coarseparttype = NULL; /* coarse level (first) partitioner */
544:   hpart->nfineparts     = 1;    /* we do not further partition coarse partition any more by default */
545:   hpart->ncoarseparts   = 0;    /* number of coarse parts (first level) */
546:   hpart->coarseparts    = NULL;
547:   hpart->fineparts      = NULL;
548:   hpart->coarseMatPart  = NULL;
549:   hpart->fineMatPart    = NULL;

551:   part->ops->apply          = MatPartitioningApply_Hierarchical;
552:   part->ops->view           = MatPartitioningView_Hierarchical;
553:   part->ops->destroy        = MatPartitioningDestroy_Hierarchical;
554:   part->ops->setfromoptions = MatPartitioningSetFromOptions_Hierarchical;
555:   part->ops->improve        = MatPartitioningImprove_Hierarchical;
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }