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