Actual source code: hierarchical.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  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 = 0;
 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:   PetscErrorCode                ierr;

 56:   PetscObjectGetComm((PetscObject)part,&comm);
 57:   MPI_Comm_size(comm,&size);
 58:   MPI_Comm_rank(comm,&rank);
 59:   PetscObjectTypeCompare((PetscObject)mat,MATMPIADJ,&flg);
 60:   if (flg) {
 61:     adj = mat;
 62:     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:    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:   if(!hpart->nfineparts) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," must set number of small subdomains for each big subdomain \n");
 74:   if(!hpart->ncoarseparts && !part->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE," did not either set number of coarse parts or total number of parts \n");

 76:   /* Partitioning the domain into one single subdomain is a trivial case, and we should just return  */
 77:   if (part->n==1) {
 78:     PetscCalloc1(bs*adj->rmap->n,&parts_indices);
 79:     ISCreateGeneral(comm,bs*adj->rmap->n,parts_indices,PETSC_OWN_POINTER,partitioning);
 80:     hpart->ncoarseparts = 1;
 81:     hpart->nfineparts = 1;
 82:     PetscStrallocpy("NONE",&hpart->coarseparttype);
 83:     PetscStrallocpy("NONE",&hpart->fineparttype);
 84:     MatDestroy(&adj);
 85:     return(0);
 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:   PetscCalloc1(hpart->ncoarseparts+1, &offsets);
 97:   PetscCalloc1(hpart->ncoarseparts, &part_weights);


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++)
112:     offsets[i] += offsets[i-1];

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

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

146:   PetscCalloc1(mat_localsize, &fineparts_indices_tmp);

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) {
152:     ISCreateGeneral(comm,mat_localsize,part->vertex_weights,PETSC_COPY_VALUES,&vweights);
153:   }

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

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

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

213:     MatDestroy(&sadj);

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

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

227:   if (part->vertex_weights) {
228:     ISDestroy(&vweights);
229:   }

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


250: PetscErrorCode MatPartitioningHierarchical_ReassembleFineparts(Mat adj, IS fineparts, ISLocalToGlobalMapping mapping, IS *sfineparts)
251: {
252:   PetscInt            *local_indices, *global_indices,*owners,*sfineparts_indices,localsize,i;
253:   const PetscInt      *ranges,*fineparts_indices;
254:   PetscMPIInt         rank;
255:   MPI_Comm            comm;
256:   PetscLayout         rmap;
257:   PetscSFNode        *remote;
258:   PetscSF             sf;
259:   PetscErrorCode      ierr;

263:   PetscObjectGetComm((PetscObject)adj,&comm);
264:   MPI_Comm_rank(comm,&rank);
265:   MatGetLayouts(adj,&rmap,NULL);
266:   ISGetLocalSize(fineparts,&localsize);
267:   PetscCalloc2(localsize,&global_indices,localsize,&local_indices);
268:   for(i=0; i<localsize; i++){
269:         local_indices[i] = i;
270:   }
271:   /* map local indices back to global so that we can permulate data globally */
272:   ISLocalToGlobalMappingApply(mapping,localsize,local_indices,global_indices);
273:   PetscCalloc1(localsize,&owners);
274:   /* find owners for global indices */
275:   for(i=0; i<localsize; i++){
276:         PetscLayoutFindOwner(rmap,global_indices[i],&owners[i]);
277:   }
278:   PetscLayoutGetRanges(rmap,&ranges);
279:   PetscCalloc1(ranges[rank+1]-ranges[rank],&sfineparts_indices);

281:   for (i=0; i<(ranges[rank+1]-ranges[rank]); i++) {
282:     sfineparts_indices[i] = -1;
283:   }

285:   ISGetIndices(fineparts,&fineparts_indices);
286:   PetscSFCreate(comm,&sf);
287:   PetscCalloc1(localsize,&remote);
288:   for(i=0; i<localsize; i++){
289:         remote[i].rank  = owners[i];
290:         remote[i].index = global_indices[i]-ranges[owners[i]];
291:   }
292:   PetscSFSetType(sf,PETSCSFBASIC);
293:   /* not sure how to add prefix to sf */
294:   PetscSFSetFromOptions(sf);
295:   PetscSFSetGraph(sf,localsize,localsize,NULL,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
296:   PetscSFReduceBegin(sf,MPIU_INT,fineparts_indices,sfineparts_indices,MPIU_REPLACE);
297:   PetscSFReduceEnd(sf,MPIU_INT,fineparts_indices,sfineparts_indices,MPIU_REPLACE);
298:   PetscSFDestroy(&sf);
299:   ISRestoreIndices(fineparts,&fineparts_indices);
300:   ISCreateGeneral(comm,ranges[rank+1]-ranges[rank],sfineparts_indices,PETSC_OWN_POINTER,sfineparts);
301:   PetscFree2(global_indices,local_indices);
302:   PetscFree(owners);
303:   return(0);
304: }


307: PetscErrorCode MatPartitioningHierarchical_AssembleSubdomain(Mat adj,IS vweights, IS destination,IS *svweights,Mat *sadj,ISLocalToGlobalMapping *mapping)
308: {
309:   IS              irows,icols;
310:   PetscInt        irows_ln;
311:   PetscMPIInt     rank;
312:   const PetscInt *irows_indices;
313:   MPI_Comm        comm;
314:   PetscErrorCode  ierr;

317:   PetscObjectGetComm((PetscObject)adj,&comm);
318:   MPI_Comm_rank(comm,&rank);
319:   /* figure out where data comes from  */
320:   ISBuildTwoSided(destination,NULL,&irows);
321:   ISDuplicate(irows,&icols);
322:   ISGetLocalSize(irows,&irows_ln);
323:   ISGetIndices(irows,&irows_indices);
324:   ISLocalToGlobalMappingCreate(comm,1,irows_ln,irows_indices,PETSC_COPY_VALUES,mapping);
325:   ISRestoreIndices(irows,&irows_indices);
326:   MatCreateSubMatrices(adj,1,&irows,&icols,MAT_INITIAL_MATRIX,&sadj);
327:   if (vweights && svweights) {
328:     ISCreateSubIS(vweights,irows,svweights);
329:   }
330:   ISDestroy(&irows);
331:   ISDestroy(&icols);
332:   return(0);
333: }


336: PetscErrorCode MatPartitioningHierarchical_DetermineDestination(MatPartitioning part, IS partitioning, PetscInt pstart, PetscInt pend, IS *destination)
337: {
338:   MPI_Comm            comm;
339:   PetscMPIInt         rank,size,target;
340:   PetscInt            plocalsize,*dest_indices,i;
341:   const PetscInt     *part_indices;
342:   PetscErrorCode      ierr;

345:   PetscObjectGetComm((PetscObject)part,&comm);
346:   MPI_Comm_rank(comm,&rank);
347:   MPI_Comm_size(comm,&size);
348:   if((pend-pstart)>size) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"range [%D, %D] should be smaller than or equal to size %D",pstart,pend,size);
349:   if(pstart>pend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP," pstart %D should be smaller than pend %D",pstart,pend);
350:   ISGetLocalSize(partitioning,&plocalsize);
351:   PetscCalloc1(plocalsize,&dest_indices);
352:   ISGetIndices(partitioning,&part_indices);
353:   for(i=0; i<plocalsize; i++){
354:         /* compute target */
355:     target = part_indices[i]-pstart;
356:     /* mark out of range entity as -1 */
357:     if(part_indices[i]<pstart || part_indices[i]>=pend) target = -1;
358:         dest_indices[i] = target;
359:   }
360:   ISCreateGeneral(comm,plocalsize,dest_indices,PETSC_OWN_POINTER,destination);
361:   return(0);
362: }


365: PetscErrorCode MatPartitioningView_Hierarchical(MatPartitioning part,PetscViewer viewer)
366: {
367:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
368:   PetscErrorCode           ierr;
369:   PetscMPIInt              rank;
370:   PetscBool                iascii;
371:   PetscViewer              sviewer;

374:   MPI_Comm_rank(PetscObjectComm((PetscObject)part),&rank);
375:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
376:   if(iascii){
377:     PetscViewerASCIIPrintf(viewer," Number of coarse parts: %D\n",hpart->ncoarseparts);
378:     PetscViewerASCIIPrintf(viewer," Coarse partitioner: %s\n",hpart->coarseparttype);
379:     if (hpart->coarseMatPart) {
380:       PetscViewerASCIIPushTab(viewer);
381:       MatPartitioningView(hpart->coarseMatPart,viewer);
382:       PetscViewerASCIIPopTab(viewer);
383:     }
384:     PetscViewerASCIIPrintf(viewer," Number of fine parts: %D\n",hpart->nfineparts);
385:     PetscViewerASCIIPrintf(viewer," Fine partitioner: %s\n",hpart->fineparttype);
386:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
387:     if (!rank && hpart->fineMatPart) {
388:       PetscViewerASCIIPushTab(viewer);
389:       MatPartitioningView(hpart->fineMatPart,sviewer);
390:       PetscViewerASCIIPopTab(viewer);
391:     }
392:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
393:   }
394:   return(0);
395: }


398: PetscErrorCode MatPartitioningHierarchicalGetFineparts(MatPartitioning part,IS *fineparts)
399: {
400:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
401:   PetscErrorCode                ierr;

404:   *fineparts = hpart->fineparts;
405:   PetscObjectReference((PetscObject)hpart->fineparts);
406:   return(0);
407: }

409: PetscErrorCode MatPartitioningHierarchicalGetCoarseparts(MatPartitioning part,IS *coarseparts)
410: {
411:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
412:   PetscErrorCode                ierr;

415:   *coarseparts = hpart->coarseparts;
416:   PetscObjectReference((PetscObject)hpart->coarseparts);
417:   return(0);
418: }

420: PetscErrorCode MatPartitioningHierarchicalSetNcoarseparts(MatPartitioning part, PetscInt ncoarseparts)
421: {
422:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;

425:   hpart->ncoarseparts = ncoarseparts;
426:   return(0);
427: }

429: PetscErrorCode MatPartitioningHierarchicalSetNfineparts(MatPartitioning part, PetscInt nfineparts)
430: {
431:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;

434:   hpart->nfineparts = nfineparts;
435:   return(0);
436: }

438: PetscErrorCode MatPartitioningSetFromOptions_Hierarchical(PetscOptionItems *PetscOptionsObject,MatPartitioning part)
439: {
440:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
442:   char           value[1024];
443:   PetscBool      flag = PETSC_FALSE;

446:   PetscOptionsHead(PetscOptionsObject,"Set hierarchical partitioning options");
447:   PetscOptionsString("-mat_partitioning_hierarchical_coarseparttype","coarse part type",NULL,NULL,value,1024,&flag);
448:   if(flag){
449:    PetscCalloc1(1024,&hpart->coarseparttype);
450:    PetscStrcpy(hpart->coarseparttype,value);
451:   }
452:   PetscOptionsString("-mat_partitioning_hierarchical_fineparttype","fine part type",NULL,NULL,value,1024,&flag);
453:   if(flag){
454:     PetscCalloc1(1024,&hpart->fineparttype);
455:     PetscStrcpy(hpart->fineparttype,value);
456:   }
457:   PetscOptionsInt("-mat_partitioning_hierarchical_ncoarseparts","number of coarse parts",NULL,hpart->ncoarseparts,&hpart->ncoarseparts,&flag);
458:   PetscOptionsInt("-mat_partitioning_hierarchical_nfineparts","number of fine parts",NULL,hpart->nfineparts,&hpart->nfineparts,&flag);
459:   PetscOptionsTail();
460:   return(0);
461: }


464: PetscErrorCode MatPartitioningDestroy_Hierarchical(MatPartitioning part)
465: {
466:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
467:   PetscErrorCode           ierr;

470:   if(hpart->coarseparttype) {PetscFree(hpart->coarseparttype);}
471:   if(hpart->fineparttype) {PetscFree(hpart->fineparttype);}
472:   ISDestroy(&hpart->fineparts);
473:   ISDestroy(&hpart->coarseparts);
474:   MatPartitioningDestroy(&hpart->coarseMatPart);
475:   MatPartitioningDestroy(&hpart->fineMatPart);
476:   MatPartitioningDestroy(&hpart->improver);
477:   PetscFree(hpart);
478:   return(0);
479: }

481: /*
482:    Improves the quality  of a partition
483: */
484: static PetscErrorCode MatPartitioningImprove_Hierarchical(MatPartitioning part, IS *partitioning)
485: {
486:   PetscErrorCode               ierr;
487:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
488:   Mat                           mat = part->adj, adj;
489:   PetscBool                    flg;
490:   PetscInt                     *vertex_weights;
491:   const char                   *prefix;

494:   PetscObjectTypeCompare((PetscObject)mat,MATMPIADJ,&flg);
495:   if (flg) {
496:     adj = mat;
497:     PetscObjectReference((PetscObject)adj);
498:   }else {
499:     /* bs indicates if the converted matrix is "reduced" from the original and hence the
500:        resulting partition results need to be stretched to match the original matrix */
501:    MatConvert(mat,MATMPIADJ,MAT_INITIAL_MATRIX,&adj);
502:   }

504:   /* If there exists a mat partitioner, we should delete it */
505:   MatPartitioningDestroy(&hpart->improver);
506:   MatPartitioningCreate(PetscObjectComm((PetscObject)part),&hpart->improver);
507:   PetscObjectGetOptionsPrefix((PetscObject)part,&prefix);
508:   PetscObjectSetOptionsPrefix((PetscObject)hpart->improver,prefix);
509:   PetscObjectAppendOptionsPrefix((PetscObject)hpart->improver,"hierarch_improver_");
510:   /* Only parmetis supports to refine a partition */
511: #if defined(PETSC_HAVE_PARMETIS)
512:     MatPartitioningSetType(hpart->improver,MATPARTITIONINGPARMETIS);
513: #else
514:     SETERRQ(PetscObjectComm((PetscObject)adj),PETSC_ERR_SUP,"Requires PETSc be installed with ParMetis\n");
515: #endif

517:   MatPartitioningSetAdjacency(hpart->improver,adj);
518:   MatPartitioningSetNParts(hpart->improver, part->n);
519:   /* copy over vertex weights */
520:   if(part->vertex_weights){
521:     PetscMalloc1(adj->rmap->n,&vertex_weights);
522:     PetscMemcpy(vertex_weights,part->vertex_weights,sizeof(PetscInt)*adj->rmap->n);
523:     MatPartitioningSetVertexWeights(hpart->improver,vertex_weights);
524:   }
525:   MatPartitioningImprove(hpart->improver,partitioning);
526:   MatDestroy(&adj);
527:   return(0);
528: }


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

540:    Collective on MPI_Comm

542:    Input Parameter:
543: .  part - the partitioning context

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

551:    Level: beginner

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

560: .keywords: Partitioning, create, context

562: .seealso: MatPartitioningSetType(), MatPartitioningType

564: M*/

566: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Hierarchical(MatPartitioning part)
567: {
568:   PetscErrorCode                ierr;
569:   MatPartitioning_Hierarchical *hpart;

572:   PetscNewLog(part,&hpart);
573:   part->data = (void*)hpart;

575:   hpart->fineparttype       = 0; /* fine level (second) partitioner */
576:   hpart->coarseparttype     = 0; /* coarse level (first) partitioner */
577:   hpart->nfineparts         = 1; /* we do not further partition coarse partition any more by default */
578:   hpart->ncoarseparts       = 0; /* number of coarse parts (first level) */
579:   hpart->coarseparts        = 0;
580:   hpart->fineparts          = 0;
581:   hpart->coarseMatPart      = 0;
582:   hpart->fineMatPart        = 0;

584:   part->ops->apply          = MatPartitioningApply_Hierarchical;
585:   part->ops->view           = MatPartitioningView_Hierarchical;
586:   part->ops->destroy        = MatPartitioningDestroy_Hierarchical;
587:   part->ops->setfromoptions = MatPartitioningSetFromOptions_Hierarchical;
588:   part->ops->improve        = MatPartitioningImprove_Hierarchical;
589:   return(0);
590: }