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: 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: }