1 /* 2 * overlapsplit.c: increase the overlap of a 'big' subdomain across several processor cores 3 * 4 * Author: Fande Kong <fdkong.jd@gmail.com> 5 */ 6 7 #include <petscsf.h> 8 #include <petsc/private/matimpl.h> 9 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatIncreaseOverlapSplit_Single" 13 14 /* 15 * Increase overlap for the sub-matrix across sub communicator 16 * sub-matrix could be a graph or numerical matrix 17 * */ 18 PetscErrorCode MatIncreaseOverlapSplit_Single(Mat mat,IS *is,PetscInt ov) 19 { 20 PetscInt i,nindx,*indices_sc,*indices_ov,localsize,*localsizes_sc,localsize_tmp; 21 PetscInt *indices_ov_rd,nroots,nleaves,*localoffsets,*indices_recv,*sources_sc,*sources_sc_rd; 22 const PetscInt *indices; 23 PetscMPIInt srank,ssize,issamecomm,k,grank; 24 IS is_sc,allis_sc,partitioning; 25 MPI_Comm gcomm,dcomm,scomm; 26 PetscSF sf; 27 PetscSFNode *remote; 28 Mat *smat; 29 MatPartitioning part; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 /* get a sub communicator before call individual MatIncreaseOverlap 34 * since the sub communicator may be changed. 35 * */ 36 ierr = PetscObjectGetComm((PetscObject)(*is),&dcomm);CHKERRQ(ierr); 37 /*make a copy before the original one is deleted*/ 38 ierr = PetscCommDuplicate(dcomm,&scomm,NULL);CHKERRQ(ierr); 39 /*get a global communicator, where mat should be a global matrix */ 40 ierr = PetscObjectGetComm((PetscObject)mat,&gcomm);CHKERRQ(ierr); 41 /*increase overlap on each individual subdomain*/ 42 ierr = (*mat->ops->increaseoverlap)(mat,1,is,ov);CHKERRQ(ierr); 43 /*compare communicators */ 44 ierr = MPI_Comm_compare(gcomm,scomm,&issamecomm);CHKERRQ(ierr); 45 /* if the sub-communicator is the same as the global communicator, 46 * user does not want to use a sub-communicator 47 * */ 48 if(issamecomm == MPI_IDENT || issamecomm == MPI_CONGRUENT) PetscFunctionReturn(0); 49 /* if the sub-communicator is petsc_comm_self, 50 * user also does not care the sub-communicator 51 * */ 52 ierr = MPI_Comm_compare(scomm,PETSC_COMM_SELF,&issamecomm);CHKERRQ(ierr); 53 if(issamecomm == MPI_IDENT || issamecomm == MPI_CONGRUENT){PetscFunctionReturn(0);} 54 /*local rank, size in a sub-communicator */ 55 ierr = MPI_Comm_rank(scomm,&srank);CHKERRQ(ierr); 56 ierr = MPI_Comm_size(scomm,&ssize);CHKERRQ(ierr); 57 ierr = MPI_Comm_rank(gcomm,&grank);CHKERRQ(ierr); 58 /*create a new IS based on sub-communicator 59 * since the old IS is often based on petsc_comm_self 60 * */ 61 ierr = ISGetLocalSize(*is,&nindx);CHKERRQ(ierr); 62 ierr = PetscCalloc1(nindx,&indices_sc);CHKERRQ(ierr); 63 ierr = ISGetIndices(*is,&indices);CHKERRQ(ierr); 64 ierr = PetscMemcpy(indices_sc,indices,sizeof(PetscInt)*nindx);CHKERRQ(ierr); 65 ierr = ISRestoreIndices(*is,&indices);CHKERRQ(ierr); 66 /*we do not need any more*/ 67 ierr = ISDestroy(is);CHKERRQ(ierr); 68 /*create a index set based on the sub communicator */ 69 ierr = ISCreateGeneral(scomm,nindx,indices_sc,PETSC_OWN_POINTER,&is_sc);CHKERRQ(ierr); 70 /*gather all indices within the sub communicator*/ 71 ierr = ISAllGather(is_sc,&allis_sc);CHKERRQ(ierr); 72 ierr = ISDestroy(&is_sc);CHKERRQ(ierr); 73 /* gather local sizes */ 74 ierr = PetscMalloc1(ssize,&localsizes_sc);CHKERRQ(ierr); 75 /*get individual local sizes for all index sets*/ 76 ierr = MPI_Gather(&nindx,1,MPIU_INT,localsizes_sc,1,MPIU_INT,0,scomm);CHKERRQ(ierr); 77 /*only root does these computations */ 78 if(!srank){ 79 /*get local size for the big index set*/ 80 ierr = ISGetLocalSize(allis_sc,&localsize);CHKERRQ(ierr); 81 ierr = PetscCalloc2(localsize,&indices_ov,localsize,&sources_sc);CHKERRQ(ierr); 82 ierr = PetscCalloc2(localsize,&indices_ov_rd,localsize,&sources_sc_rd);CHKERRQ(ierr); 83 ierr = ISGetIndices(allis_sc,&indices);CHKERRQ(ierr); 84 ierr = PetscMemcpy(indices_ov,indices,sizeof(PetscInt)*localsize);CHKERRQ(ierr); 85 ierr = ISRestoreIndices(allis_sc,&indices);CHKERRQ(ierr); 86 /*we do not need it any more */ 87 ierr = ISDestroy(&allis_sc);CHKERRQ(ierr); 88 /*assign corresponding sources */ 89 localsize_tmp = 0; 90 for(k=0; k<ssize; k++){ 91 for(i=0; i<localsizes_sc[k]; i++){ 92 sources_sc[localsize_tmp++] = k; 93 } 94 } 95 /*record where indices come from */ 96 ierr = PetscSortIntWithArray(localsize,indices_ov,sources_sc);CHKERRQ(ierr); 97 /*count local sizes for reduced indices */ 98 ierr = PetscMemzero(localsizes_sc,sizeof(PetscInt)*ssize);CHKERRQ(ierr); 99 /*initialize the first entity*/ 100 if(localsize){ 101 indices_ov_rd[0] = indices_ov[0]; 102 sources_sc_rd[0] = sources_sc[0]; 103 localsizes_sc[sources_sc[0]]++; 104 } 105 localsize_tmp = 1; 106 /*remove duplicate integers */ 107 for(i=1; i<localsize; i++){ 108 if(indices_ov[i] != indices_ov[i-1]){ 109 indices_ov_rd[localsize_tmp] = indices_ov[i]; 110 sources_sc_rd[localsize_tmp++] = sources_sc[i]; 111 localsizes_sc[sources_sc[i]]++; 112 } 113 } 114 ierr = PetscFree2(indices_ov,sources_sc);CHKERRQ(ierr); 115 ierr = PetscCalloc1(ssize+1,&localoffsets);CHKERRQ(ierr); 116 for(k=0; k<ssize; k++){ 117 localoffsets[k+1] = localoffsets[k] + localsizes_sc[k]; 118 } 119 /*construct a star forest to send data back */ 120 nleaves = localoffsets[ssize]; 121 ierr = PetscMemzero(localoffsets,(ssize+1)*sizeof(PetscInt));CHKERRQ(ierr); 122 nroots = localsizes_sc[srank]; 123 ierr = PetscCalloc1(nleaves,&remote);CHKERRQ(ierr); 124 for(i=0; i<nleaves; i++){ 125 remote[i].rank = sources_sc_rd[i]; 126 remote[i].index = localoffsets[sources_sc_rd[i]]++; 127 } 128 ierr = PetscFree(localoffsets);CHKERRQ(ierr); 129 }else{ 130 ierr = ISDestroy(&allis_sc);CHKERRQ(ierr); 131 /*Allocate a 'zero' pointer */ 132 ierr = PetscCalloc1(0,&remote);CHKERRQ(ierr); 133 nleaves = 0; 134 indices_ov_rd = 0; 135 sources_sc_rd = 0; 136 } 137 /*scatter sizes to everybody */ 138 ierr = MPI_Scatter(localsizes_sc,1, MPIU_INT,&nroots,1, MPIU_INT,0,scomm);CHKERRQ(ierr); 139 /*free memory */ 140 ierr = PetscFree(localsizes_sc);CHKERRQ(ierr); 141 ierr = PetscCalloc1(nroots,&indices_recv);CHKERRQ(ierr); 142 /*ierr = MPI_Comm_dup(scomm,&dcomm);CHKERRQ(ierr);*/ 143 /*set data back to every body */ 144 ierr = PetscSFCreate(scomm,&sf);CHKERRQ(ierr); 145 ierr = PetscSFSetType(sf,PETSCSFBASIC);CHKERRQ(ierr); 146 ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr); 147 ierr = PetscSFSetGraph(sf,nroots,nleaves,PETSC_NULL,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);CHKERRQ(ierr); 148 ierr = PetscSFReduceBegin(sf,MPIU_INT,indices_ov_rd,indices_recv,MPIU_REPLACE);CHKERRQ(ierr); 149 ierr = PetscSFReduceEnd(sf,MPIU_INT,indices_ov_rd,indices_recv,MPIU_REPLACE);CHKERRQ(ierr); 150 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 151 /* free memory */ 152 ierr = PetscFree2(indices_ov_rd,sources_sc_rd);CHKERRQ(ierr); 153 /*create a index set*/ 154 ierr = ISCreateGeneral(scomm,nroots,indices_recv,PETSC_OWN_POINTER,&is_sc);CHKERRQ(ierr); 155 /*construct a parallel submatrix */ 156 ierr = MatGetSubMatricesMPI(mat,1,&is_sc,&is_sc,MAT_INITIAL_MATRIX,&smat);CHKERRQ(ierr); 157 /* we do not need them any more */ 158 ierr = ISDestroy(&allis_sc);CHKERRQ(ierr); 159 /*create a partitioner to repartition the sub-matrix*/ 160 ierr = MatPartitioningCreate(scomm,&part);CHKERRQ(ierr); 161 ierr = MatPartitioningSetAdjacency(part,smat[0]);CHKERRQ(ierr); 162 #if PETSC_HAVE_PARMETIS 163 /* if there exists a ParMETIS installation, we try to use ParMETIS 164 * because a repartition routine possibly work better 165 * */ 166 ierr = MatPartitioningSetType(part,MATPARTITIONINGPARMETIS);CHKERRQ(ierr); 167 /*try to use reparition function, instead of partition function */ 168 ierr = MatPartitioningParmetisSetRepartition(part);CHKERRQ(ierr); 169 #else 170 /*we at least provide a default partitioner to rebalance the computation */ 171 ierr = MatPartitioningSetType(part,MATPARTITIONINGAVERAGE);CHKERRQ(ierr); 172 #endif 173 /*user can pick up any partitioner by using an option*/ 174 ierr = MatPartitioningSetFromOptions(part);CHKERRQ(ierr); 175 /* apply partition */ 176 ierr = MatPartitioningApply(part,&partitioning);CHKERRQ(ierr); 177 ierr = MatPartitioningDestroy(&part);CHKERRQ(ierr); 178 ierr = MatDestroy(&(smat[0]));CHKERRQ(ierr); 179 ierr = PetscFree(smat);CHKERRQ(ierr); 180 /* get local rows including overlap */ 181 ierr = ISBuildTwoSided(partitioning,is_sc,is);CHKERRQ(ierr); 182 /* destroy */ 183 ierr = ISDestroy(&is_sc);CHKERRQ(ierr); 184 ierr = ISDestroy(&partitioning);CHKERRQ(ierr); 185 ierr = PetscCommDestroy(&scomm);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 190