Actual source code: hierarchical.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: #include <../src/mat/impls/adj/mpi/mpiadj.h>    /*I "petscmat.h" I*/
  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 part, IS partitioning, PetscInt pstart, PetscInt pend, IS *destination);
 16: PetscErrorCode MatPartitioningHierarchical_AssembleSubdomain(Mat adj,IS destination,Mat *sadj, ISLocalToGlobalMapping *mapping);
 17: PetscErrorCode MatPartitioningHierarchical_ReassembleFineparts(Mat adj, IS fineparts, ISLocalToGlobalMapping mapping, IS *sfineparts);

 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_Hierarchical;

 28: /*
 29:    Uses a hierarchical partitioning strategy to partition the matrix in parallel.
 30:    Use this interface to make the partitioner consistent with others
 31: */
 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                     *parts_indices,i,j,mat_localsize;
 39:   Mat                           mat    = part->adj,adj,sadj;
 40:   PetscBool                     flg;
 41:   PetscInt                      bs     = 1;
 42:   MatPartitioning               finePart, coarsePart;
 43:   PetscInt                     *coarse_vertex_weights = 0;
 44:   PetscMPIInt                   size,rank;
 45:   MPI_Comm                      comm,scomm;
 46:   IS                            destination,fineparts_temp;
 47:   ISLocalToGlobalMapping        mapping;
 48:   PetscErrorCode                ierr;

 51:   PetscObjectGetComm((PetscObject)part,&comm);
 52:   MPI_Comm_size(comm,&size);
 53:   MPI_Comm_rank(comm,&rank);
 54:   PetscObjectTypeCompare((PetscObject)mat,MATMPIADJ,&flg);
 55:   if (flg) {
 56:     adj = mat;
 57:     PetscObjectReference((PetscObject)adj);
 58:   }else {
 59:     /* bs indicates if the converted matrix is "reduced" from the original and hence the
 60:        resulting partition results need to be stretched to match the original matrix */
 61:    MatConvert(mat,MATMPIADJ,MAT_INITIAL_MATRIX,&adj);
 62:    if (adj->rmap->n > 0) bs = mat->rmap->n/adj->rmap->n;
 63:   }
 64:   /* local size of mat */
 65:   mat_localsize = adj->rmap->n;
 66:   /* check parameters */
 67:   /* how many small subdomains we want from a given 'big' suddomain */
 68:   if(!hpart->Nfineparts) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," must set number of small subdomains for each big subdomain \n");
 69:   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");
 70:   if(part->n && part->n%hpart->Nfineparts!=0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,
 71:                    " total number of parts %D can not be divided by number of fine parts %D\n",part->n,hpart->Nfineparts);
 72:   if(part->n){
 73:     hpart->Ncoarseparts = part->n/hpart->Nfineparts;
 74:   }else{
 75:         part->n = hpart->Ncoarseparts*hpart->Nfineparts;
 76:   }
 77:    /* we do not support this case currently, but this restriction should be
 78:      * removed in the further
 79:      * */
 80:   if(hpart->Ncoarseparts>size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP," we do not support number of coarse parts %D > size %D \n",hpart->Ncoarseparts,size);
 81:   MatPartitioningCreate(comm,&coarsePart);
 82:     /* if did not set partitioning type yet, use parmetis by default */
 83:   if(!hpart->coarseparttype){
 84:         MatPartitioningSetType(coarsePart,MATPARTITIONINGPARMETIS);
 85:   }else{
 86:         MatPartitioningSetType(coarsePart,hpart->coarseparttype);
 87:   }
 88:   MatPartitioningSetAdjacency(coarsePart,adj);
 89:   MatPartitioningSetNParts(coarsePart, hpart->Ncoarseparts);
 90:   /* copy over vertex weights */
 91:   if(part->vertex_weights){
 92:    PetscMalloc(sizeof(PetscInt)*mat_localsize,&coarse_vertex_weights);
 93:    PetscMemcpy(coarse_vertex_weights,part->vertex_weights,sizeof(PetscInt)*mat_localsize);
 94:    MatPartitioningSetVertexWeights(coarsePart,coarse_vertex_weights);
 95:   }
 96:   /* It looks nontrivial to support part weights,
 97:    * I will return back to implement it when have
 98:    * an idea.
 99:    *  */
100:   MatPartitioningApply(coarsePart,&hpart->coarseparts);
101:   MatPartitioningDestroy(&coarsePart);
102:   /* In the current implementation, destination should be the same as hpart->coarseparts,
103:    * and this interface is preserved to deal with the case hpart->coarseparts>size in the
104:    * future.
105:    * */
106:   MatPartitioningHierarchical_DetermineDestination(part,hpart->coarseparts,0,hpart->Ncoarseparts,&destination);
107:   /* assemble a submatrix for partitioning subdomains  */
108:   MatPartitioningHierarchical_AssembleSubdomain(adj,destination,&sadj,&mapping);
109:   ISDestroy(&destination);
110:   PetscObjectGetComm((PetscObject)sadj,&scomm);
111:   /* create a fine partitioner */
112:   MatPartitioningCreate(scomm,&finePart);
113:   /* if do not set partitioning type, use parmetis by default */
114:   if(!hpart->fineparttype){
115:     MatPartitioningSetType(finePart,MATPARTITIONINGPARMETIS);
116:   }else{
117:     MatPartitioningSetType(finePart,hpart->fineparttype);
118:   }
119:   MatPartitioningSetAdjacency(finePart,sadj);
120:   MatPartitioningSetNParts(finePart, hpart->Nfineparts);
121:   MatPartitioningApply(finePart,&fineparts_temp);
122:   MatDestroy(&sadj);
123:   MatPartitioningDestroy(&finePart);
124:   MatPartitioningHierarchical_ReassembleFineparts(adj,fineparts_temp,mapping,&hpart->fineparts);
125:   ISDestroy(&fineparts_temp);
126:   ISLocalToGlobalMappingDestroy(&mapping);

128:   ISGetIndices(hpart->fineparts,&fineparts_indices);
129:   ISGetIndices(hpart->coarseparts,&coarseparts_indices);
130:   PetscMalloc1(bs*adj->rmap->n,&parts_indices);
131:   for(i=0; i<adj->rmap->n; i++){
132:     for(j=0; j<bs; j++){
133:       parts_indices[bs*i+j] = fineparts_indices[i]+coarseparts_indices[i]*hpart->Nfineparts;
134:     }
135:   }
136:   ISCreateGeneral(comm,bs*adj->rmap->n,parts_indices,PETSC_OWN_POINTER,partitioning);
137:   MatDestroy(&adj);
138:   return(0);
139: }


144: PetscErrorCode MatPartitioningHierarchical_ReassembleFineparts(Mat adj, IS fineparts, ISLocalToGlobalMapping mapping, IS *sfineparts)
145: {
146:   PetscInt            *local_indices, *global_indices,*owners,*sfineparts_indices,localsize,i;
147:   const PetscInt      *ranges,*fineparts_indices;
148:   PetscMPIInt         rank;
149:   MPI_Comm            comm;
150:   PetscLayout         rmap;
151:   PetscSFNode        *remote;
152:   PetscSF             sf;
153:   PetscErrorCode      ierr;

156:   PetscObjectGetComm((PetscObject)adj,&comm);
157:   MPI_Comm_rank(comm,&rank);
158:   MatGetLayouts(adj,&rmap,PETSC_NULL);
159:   ISGetLocalSize(fineparts,&localsize);
160:   PetscCalloc2(localsize,&global_indices,localsize,&local_indices);
161:   for(i=0; i<localsize; i++){
162:         local_indices[i] = i;
163:   }
164:   /* map local indices back to global so that we can permulate data globally */
165:   ISLocalToGlobalMappingApply(mapping,localsize,local_indices,global_indices);
166:   PetscCalloc1(localsize,&owners);
167:   /* find owners for global indices */
168:   for(i=0; i<localsize; i++){
169:         PetscLayoutFindOwner(rmap,global_indices[i],&owners[i]);
170:   }
171:   PetscLayoutGetRanges(rmap,&ranges);
172:   PetscCalloc1(ranges[rank+1]-ranges[rank],&sfineparts_indices);
173:   ISGetIndices(fineparts,&fineparts_indices);
174:   PetscSFCreate(comm,&sf);
175:   PetscCalloc1(localsize,&remote);
176:   for(i=0; i<localsize; i++){
177:         remote[i].rank  = owners[i];
178:         remote[i].index = global_indices[i]-ranges[owners[i]];
179:   }
180:   PetscSFSetType(sf,PETSCSFBASIC);
181:   /* not sure how to add prefix to sf */
182:   PetscSFSetFromOptions(sf);
183:   PetscSFSetGraph(sf,localsize,localsize,PETSC_NULL,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
184:   PetscSFReduceBegin(sf,MPIU_INT,fineparts_indices,sfineparts_indices,MPIU_REPLACE);
185:   PetscSFReduceEnd(sf,MPIU_INT,fineparts_indices,sfineparts_indices,MPIU_REPLACE);
186:   PetscSFDestroy(&sf);
187:   ISRestoreIndices(fineparts,&fineparts_indices);
188:   ISCreateGeneral(comm,ranges[rank+1]-ranges[rank],sfineparts_indices,PETSC_OWN_POINTER,sfineparts);
189:   PetscFree2(global_indices,local_indices);
190:   PetscFree(owners);
191:   return(0);
192: }


197: PetscErrorCode MatPartitioningHierarchical_AssembleSubdomain(Mat adj,IS destination,Mat *sadj, ISLocalToGlobalMapping *mapping)
198: {
199:   IS              irows,icols;
200:   PetscInt        irows_ln;
201:   PetscMPIInt     rank;
202:   const PetscInt *irows_indices;
203:   MPI_Comm        comm;
204:   PetscErrorCode  ierr;

207:   PetscObjectGetComm((PetscObject)adj,&comm);
208:   MPI_Comm_rank(comm,&rank);
209:   /* figure out where data comes from  */
210:   ISBuildTwoSided(destination,NULL,&irows);
211:   ISDuplicate(irows,&icols);
212:   ISGetLocalSize(irows,&irows_ln);
213:   ISGetIndices(irows,&irows_indices);
214:   ISLocalToGlobalMappingCreate(comm,1,irows_ln,irows_indices,PETSC_COPY_VALUES,mapping);
215:   ISRestoreIndices(irows,&irows_indices);
216:   MatGetSubMatrices(adj,1,&irows,&icols,MAT_INITIAL_MATRIX,&sadj);
217:   ISDestroy(&irows);
218:   ISDestroy(&icols);
219:   return(0);
220: }


225: PetscErrorCode MatPartitioningHierarchical_DetermineDestination(MatPartitioning part, IS partitioning, PetscInt pstart, PetscInt pend, IS *destination)
226: {
227:   MPI_Comm            comm;
228:   PetscMPIInt         rank,size,target;
229:   PetscInt            plocalsize,*dest_indices,i;
230:   const PetscInt     *part_indices;
231:   PetscErrorCode      ierr;

234:   PetscObjectGetComm((PetscObject)part,&comm);
235:   MPI_Comm_rank(comm,&rank);
236:   MPI_Comm_size(comm,&size);
237:   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);
238:   if(pstart>pend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP," pstart %D should be smaller than pend %D",pstart,pend);
239:   ISGetLocalSize(partitioning,&plocalsize);
240:   PetscCalloc1(plocalsize,&dest_indices);
241:   ISGetIndices(partitioning,&part_indices);
242:   for(i=0; i<plocalsize; i++){
243:         /* compute target */
244:     target = part_indices[i]-pstart;
245:     /* mark out of range entity as -1 */
246:     if(part_indices[i]<pstart || part_indices[i]>pend) target = -1;
247:         dest_indices[i] = target;
248:   }
249:   ISCreateGeneral(comm,plocalsize,dest_indices,PETSC_OWN_POINTER,destination);
250:   return(0);
251: }


256: PetscErrorCode MatPartitioningView_Hierarchical(MatPartitioning part,PetscViewer viewer)
257: {
258:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
259:   PetscErrorCode           ierr;
260:   PetscMPIInt              rank;
261:   PetscBool                iascii;

264:   MPI_Comm_rank(PetscObjectComm((PetscObject)part),&rank);
265:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
266:   if(iascii){
267:          PetscViewerASCIIPrintf(viewer," Fine partitioner %s \n",hpart->fineparttype);
268:          PetscViewerASCIIPrintf(viewer," Coarse partitioner %s \n",hpart->coarseparttype);
269:          PetscViewerASCIIPrintf(viewer," Number of coarse parts %D \n",hpart->Ncoarseparts);
270:          PetscViewerASCIIPrintf(viewer," Number of fine parts %D \n",hpart->Nfineparts);
271:   }
272:   return(0);
273: }


278: PetscErrorCode MatPartitioningHierarchicalGetFineparts(MatPartitioning part,IS *fineparts)
279: {
280:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
281:   PetscErrorCode                ierr;

284:   *fineparts = hpart->fineparts;
285:   PetscObjectReference((PetscObject)hpart->fineparts);
286:   return(0);
287: }

291: PetscErrorCode MatPartitioningHierarchicalGetCoarseparts(MatPartitioning part,IS *coarseparts)
292: {
293:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
294:   PetscErrorCode                ierr;

297:   *coarseparts = hpart->coarseparts;
298:   PetscObjectReference((PetscObject)hpart->coarseparts);
299:   return(0);
300: }

304: PetscErrorCode MatPartitioningHierarchicalSetNcoarseparts(MatPartitioning part, PetscInt Ncoarseparts)
305: {
306:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;

309:   hpart->Ncoarseparts = Ncoarseparts;
310:   return(0);
311: }

315: PetscErrorCode MatPartitioningHierarchicalSetNfineparts(MatPartitioning part, PetscInt Nfineparts)
316: {
317:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;

320:   hpart->Nfineparts = Nfineparts;
321:   return(0);
322: }

326: PetscErrorCode MatPartitioningSetFromOptions_Hierarchical(PetscOptionItems *PetscOptionsObject,MatPartitioning part)
327: {
328:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
330:   char           value[1024];
331:   PetscBool      flag = PETSC_FALSE;

334:   PetscOptionsHead(PetscOptionsObject,"Set hierarchical partitioning options");
335:   PetscOptionsString("-mat_partitioning_hierarchical_coarseparttype","coarse part type",PETSC_NULL,PETSC_NULL,value,1024,&flag);
336:   if(flag){
337:    PetscCalloc1(1024,&hpart->coarseparttype);
338:    PetscStrcpy(hpart->coarseparttype,value);
339:   }
340:   PetscOptionsString("-mat_partitioning_hierarchical_fineparttype","fine part type",PETSC_NULL,PETSC_NULL,value,1024,&flag);
341:   if(flag){
342:     PetscCalloc1(1024,&hpart->fineparttype);
343:     PetscStrcpy(hpart->fineparttype,value);
344:   }
345:   PetscOptionsInt("-mat_partitioning_hierarchical_Ncoarseparts","number of coarse parts",PETSC_NULL,0,&hpart->Ncoarseparts,&flag);
346:   PetscOptionsInt("-mat_partitioning_hierarchical_Nfineparts","number of fine parts",PETSC_NULL,1,&hpart->Nfineparts,&flag);
347:   PetscOptionsTail();
348:   return(0);
349: }


354: PetscErrorCode MatPartitioningDestroy_Hierarchical(MatPartitioning part)
355: {
356:   MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
357:   PetscErrorCode           ierr;

360:   if(hpart->coarseparttype) {PetscFree(hpart->coarseparttype);}
361:   if(hpart->fineparttype) {PetscFree(hpart->fineparttype);}
362:   ISDestroy(&hpart->fineparts);
363:   ISDestroy(&hpart->coarseparts);
364:   PetscFree(hpart);
365:   return(0);
366: }


369: /*MC
370:    MATPARTITIONINGHIERARCHPART - Creates a partitioning context via hierarchical partitioning strategy.

372:    Collective on MPI_Comm

374:    Input Parameter:
375: .  part - the partitioning context

377:    Options Database Keys:

379:    Level: beginner

381: .keywords: Partitioning, create, context

383: .seealso: MatPartitioningSetType(), MatPartitioningType

385: M*/

389: PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Hierarchical(MatPartitioning part)
390: {
391:   PetscErrorCode                ierr;
392:   MatPartitioning_Hierarchical *hpart;

395:   PetscNewLog(part,&hpart);
396:   part->data = (void*)hpart;

398:   hpart->fineparttype       = 0; /* fine level partitioner */
399:   hpart->coarseparttype     = 0; /* coarse level partitioner */
400:   hpart->Nfineparts         = 1; /* we do not further partition coarse partition any more by default */
401:   hpart->Ncoarseparts       = 0; /* number of coarse parts (first level) */
402:   hpart->coarseparts        = 0;
403:   hpart->fineparts          = 0;

405:   part->ops->apply          = MatPartitioningApply_Hierarchical;
406:   part->ops->view           = MatPartitioningView_Hierarchical;
407:   part->ops->destroy        = MatPartitioningDestroy_Hierarchical;
408:   part->ops->setfromoptions = MatPartitioningSetFromOptions_Hierarchical;
409:   return(0);
410: }