Actual source code: hierarchical.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

210:     PetscCall(MatDestroy(&sadj));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

510:    Collective

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

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

521:    Level: beginner

523:    Note:
524:    See {cite}`kong2016highly` and {cite}`kongstognergastonpetersonpermannslaughtermartineau2018`.

526: .seealso: `MatPartitioningSetType()`, `MatPartitioningType`, `MATPARTITIONINGMETIS`, `MATPARTITIONINGPARMETIS`,
527: M*/

529: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Hierarchical(MatPartitioning part)
530: {
531:   MatPartitioning_Hierarchical *hpart;

533:   PetscFunctionBegin;
534:   PetscCall(PetscNew(&hpart));
535:   part->data = (void *)hpart;

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

546:   part->ops->apply          = MatPartitioningApply_Hierarchical;
547:   part->ops->view           = MatPartitioningView_Hierarchical;
548:   part->ops->destroy        = MatPartitioningDestroy_Hierarchical;
549:   part->ops->setfromoptions = MatPartitioningSetFromOptions_Hierarchical;
550:   part->ops->improve        = MatPartitioningImprove_Hierarchical;
551:   PetscFunctionReturn(PETSC_SUCCESS);
552: }