Actual source code: overlapsplit.c

  1: /*
  2:  * Increase the overlap of a 'big' subdomain across several processor cores
  3:  *
  4:  * Author: Fande Kong <fdkong.jd@gmail.com>
  5:  */

  7: #include <petscsf.h>
  8: #include <petsc/private/matimpl.h>

 10: /*
 11:  * Increase overlap for the sub-matrix across sub communicator
 12:  * sub-matrix could be a graph or numerical matrix
 13:  * */
 14: PetscErrorCode  MatIncreaseOverlapSplit_Single(Mat mat,IS *is,PetscInt ov)
 15: {
 16:   PetscInt         i,nindx,*indices_sc,*indices_ov,localsize,*localsizes_sc,localsize_tmp;
 17:   PetscInt         *indices_ov_rd,nroots,nleaves,*localoffsets,*indices_recv,*sources_sc,*sources_sc_rd;
 18:   const PetscInt   *indices;
 19:   PetscMPIInt      srank,ssize,issamecomm,k,grank;
 20:   IS               is_sc,allis_sc,partitioning;
 21:   MPI_Comm         gcomm,dcomm,scomm;
 22:   PetscSF          sf;
 23:   PetscSFNode      *remote;
 24:   Mat              *smat;
 25:   MatPartitioning  part;

 27:   /* get a sub communicator before call individual MatIncreaseOverlap
 28:    * since the sub communicator may be changed.
 29:    * */
 30:   PetscObjectGetComm((PetscObject)(*is),&dcomm);
 31:   /* make a copy before the original one is deleted */
 32:   PetscCommDuplicate(dcomm,&scomm,NULL);
 33:   /* get a global communicator, where mat should be a global matrix  */
 34:   PetscObjectGetComm((PetscObject)mat,&gcomm);
 35:   (*mat->ops->increaseoverlap)(mat,1,is,ov);
 36:   MPI_Comm_compare(gcomm,scomm,&issamecomm);
 37:   /* if the sub-communicator is the same as the global communicator,
 38:    * user does not want to use a sub-communicator
 39:    * */
 40:   if (issamecomm == MPI_IDENT || issamecomm == MPI_CONGRUENT) {
 41:         PetscCommDestroy(&scomm);
 42:         return 0;
 43:   }
 44:   /* if the sub-communicator is petsc_comm_self,
 45:    * user also does not care the sub-communicator
 46:    * */
 47:   MPI_Comm_compare(scomm,PETSC_COMM_SELF,&issamecomm);
 48:   if (issamecomm == MPI_IDENT || issamecomm == MPI_CONGRUENT) {
 49:     PetscCommDestroy(&scomm);
 50:     return 0;
 51:   }
 52:   MPI_Comm_rank(scomm,&srank);
 53:   MPI_Comm_size(scomm,&ssize);
 54:   MPI_Comm_rank(gcomm,&grank);
 55:   /* create a new IS based on sub-communicator
 56:    * since the old IS is often based on petsc_comm_self
 57:    * */
 58:   ISGetLocalSize(*is,&nindx);
 59:   PetscMalloc1(nindx,&indices_sc);
 60:   ISGetIndices(*is,&indices);
 61:   PetscArraycpy(indices_sc,indices,nindx);
 62:   ISRestoreIndices(*is,&indices);
 63:   /* we do not need any more */
 64:   ISDestroy(is);
 65:   /* create a index set based on the sub communicator  */
 66:   ISCreateGeneral(scomm,nindx,indices_sc,PETSC_OWN_POINTER,&is_sc);
 67:   /* gather all indices within  the sub communicator */
 68:   ISAllGather(is_sc,&allis_sc);
 69:   ISDestroy(&is_sc);
 70:   /* gather local sizes */
 71:   PetscMalloc1(ssize,&localsizes_sc);
 72:   /* get individual local sizes for all index sets */
 73:   MPI_Gather(&nindx,1,MPIU_INT,localsizes_sc,1,MPIU_INT,0,scomm);
 74:   /* only root does these computations */
 75:   if (!srank) {
 76:    /* get local size for the big index set */
 77:    ISGetLocalSize(allis_sc,&localsize);
 78:    PetscCalloc2(localsize,&indices_ov,localsize,&sources_sc);
 79:    PetscCalloc2(localsize,&indices_ov_rd,localsize,&sources_sc_rd);
 80:    ISGetIndices(allis_sc,&indices);
 81:    PetscArraycpy(indices_ov,indices,localsize);
 82:    ISRestoreIndices(allis_sc,&indices);
 83:    ISDestroy(&allis_sc);
 84:    /* assign corresponding sources */
 85:    localsize_tmp = 0;
 86:    for (k=0; k<ssize; k++) {
 87:      for (i=0; i<localsizes_sc[k]; i++) {
 88:        sources_sc[localsize_tmp++] = k;
 89:      }
 90:    }
 91:    /* record where indices come from */
 92:    PetscSortIntWithArray(localsize,indices_ov,sources_sc);
 93:    /* count local sizes for reduced indices */
 94:    PetscArrayzero(localsizes_sc,ssize);
 95:    /* initialize the first entity */
 96:    if (localsize) {
 97:      indices_ov_rd[0] = indices_ov[0];
 98:      sources_sc_rd[0] = sources_sc[0];
 99:      localsizes_sc[sources_sc[0]]++;
100:    }
101:    localsize_tmp = 1;
102:    /* remove duplicate integers */
103:    for (i=1; i<localsize; i++) {
104:      if (indices_ov[i] != indices_ov[i-1]) {
105:        indices_ov_rd[localsize_tmp]   = indices_ov[i];
106:        sources_sc_rd[localsize_tmp++] = sources_sc[i];
107:        localsizes_sc[sources_sc[i]]++;
108:      }
109:    }
110:    PetscFree2(indices_ov,sources_sc);
111:    PetscCalloc1(ssize+1,&localoffsets);
112:    for (k=0; k<ssize; k++) {
113:      localoffsets[k+1] = localoffsets[k] + localsizes_sc[k];
114:    }
115:    nleaves = localoffsets[ssize];
116:    PetscArrayzero(localoffsets,ssize+1);
117:    nroots  = localsizes_sc[srank];
118:    PetscMalloc1(nleaves,&remote);
119:    for (i=0; i<nleaves; i++) {
120:      remote[i].rank  = sources_sc_rd[i];
121:      remote[i].index = localoffsets[sources_sc_rd[i]]++;
122:    }
123:    PetscFree(localoffsets);
124:   } else {
125:    ISDestroy(&allis_sc);
126:    /* Allocate a 'zero' pointer to avoid using uninitialized variable  */
127:    PetscCalloc1(0,&remote);
128:    nleaves       = 0;
129:    indices_ov_rd = NULL;
130:    sources_sc_rd = NULL;
131:   }
132:   /* scatter sizes to everybody */
133:   MPI_Scatter(localsizes_sc,1, MPIU_INT,&nroots,1, MPIU_INT,0,scomm);
134:   PetscFree(localsizes_sc);
135:   PetscCalloc1(nroots,&indices_recv);
136:   /* set data back to every body */
137:   PetscSFCreate(scomm,&sf);
138:   PetscSFSetType(sf,PETSCSFBASIC);
139:   PetscSFSetFromOptions(sf);
140:   PetscSFSetGraph(sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
141:   PetscSFReduceBegin(sf,MPIU_INT,indices_ov_rd,indices_recv,MPI_REPLACE);
142:   PetscSFReduceEnd(sf,MPIU_INT,indices_ov_rd,indices_recv,MPI_REPLACE);
143:   PetscSFDestroy(&sf);
144:   PetscFree2(indices_ov_rd,sources_sc_rd);
145:   ISCreateGeneral(scomm,nroots,indices_recv,PETSC_OWN_POINTER,&is_sc);
146:   MatCreateSubMatricesMPI(mat,1,&is_sc,&is_sc,MAT_INITIAL_MATRIX,&smat);
147:   ISDestroy(&allis_sc);
148:   /* create a partitioner to repartition the sub-matrix */
149:   MatPartitioningCreate(scomm,&part);
150:   MatPartitioningSetAdjacency(part,smat[0]);
151: #if defined(PETSC_HAVE_PARMETIS)
152:   /* if there exists a ParMETIS installation, we try to use ParMETIS
153:    * because a repartition routine possibly work better
154:    * */
155:   MatPartitioningSetType(part,MATPARTITIONINGPARMETIS);
156:   /* try to use reparition function, instead of partition function */
157:   MatPartitioningParmetisSetRepartition(part);
158: #else
159:   /* we at least provide a default partitioner to rebalance the computation  */
160:   MatPartitioningSetType(part,MATPARTITIONINGAVERAGE);
161: #endif
162:   /* user can pick up any partitioner by using an option */
163:   MatPartitioningSetFromOptions(part);
164:   MatPartitioningApply(part,&partitioning);
165:   MatPartitioningDestroy(&part);
166:   MatDestroy(&(smat[0]));
167:   PetscFree(smat);
168:   /* get local rows including  overlap */
169:   ISBuildTwoSided(partitioning,is_sc,is);
170:   ISDestroy(&is_sc);
171:   ISDestroy(&partitioning);
172:   PetscCommDestroy(&scomm);
173:   return 0;
174: }