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