1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2 #include <../src/mat/impls/aij/seq/aij.h>
3 #include <../src/mat/impls/aij/mpi/mpiaij.h>
4 #include <petscsf.h>
5
6 #define MIS_NOT_DONE -2
7 #define MIS_DELETED -1
8 #define MIS_REMOVED -3
9 #define MIS_IS_SELECTED(s) (s != MIS_DELETED && s != MIS_NOT_DONE && s != MIS_REMOVED)
10
11 /*
12 MatCoarsenApply_MIS_private - parallel maximal independent set (MIS) with data locality info. MatAIJ specific!!!
13
14 Input Parameter:
15 . perm - serial permutation of rows of local to process in MIS
16 . Gmat - global matrix of graph (data not defined)
17 . strict_aggs - flag for whether to keep strict (non overlapping) aggregates in 'llist';
18
19 Output Parameter:
20 . a_selected - IS of selected vertices, includes 'ghost' nodes at end with natural local indices
21 . a_locals_llist - array of list of nodes rooted at selected nodes
22 */
MatCoarsenApply_MIS_private(IS perm,Mat Gmat,PetscBool strict_aggs,PetscCoarsenData ** a_locals_llist)23 static PetscErrorCode MatCoarsenApply_MIS_private(IS perm, Mat Gmat, PetscBool strict_aggs, PetscCoarsenData **a_locals_llist)
24 {
25 Mat_SeqAIJ *matA, *matB = NULL;
26 Mat_MPIAIJ *mpimat = NULL;
27 MPI_Comm comm;
28 PetscInt num_fine_ghosts, kk, n, ix, j, *idx, *ii, Iend, my0, nremoved, gid, lid, cpid, lidj, sgid, t1, t2, slid, nDone, nselected = 0, state, statej;
29 PetscInt *cpcol_gid, *cpcol_state, *lid_cprowID, *lid_gid, *cpcol_sel_gid, *icpcol_gid, *lid_state, *lid_parent_gid = NULL, nrm_tot = 0;
30 PetscBool *lid_removed;
31 PetscBool isMPI, isAIJ, isOK;
32 const PetscInt *perm_ix;
33 const PetscInt nloc = Gmat->rmap->n;
34 PetscCoarsenData *agg_lists;
35 PetscLayout layout;
36 PetscSF sf;
37 IS info_is;
38
39 PetscFunctionBegin;
40 PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
41 PetscCall(ISCreate(comm, &info_is));
42 PetscCall(PetscInfo(info_is, "mis: nloc = %d\n", (int)nloc));
43 /* get submatrices */
44 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPI));
45 if (isMPI) {
46 mpimat = (Mat_MPIAIJ *)Gmat->data;
47 matA = (Mat_SeqAIJ *)mpimat->A->data;
48 matB = (Mat_SeqAIJ *)mpimat->B->data;
49 /* force compressed storage of B */
50 PetscCall(MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, Gmat->rmap->n, -1.0));
51 } else {
52 matA = (Mat_SeqAIJ *)Gmat->data;
53 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isAIJ));
54 PetscCheck(isAIJ, comm, PETSC_ERR_PLIB, "Require AIJ matrix.");
55 }
56 PetscCall(MatGetOwnershipRange(Gmat, &my0, &Iend));
57 PetscCall(PetscMalloc4(nloc, &lid_gid, nloc, &lid_cprowID, nloc, &lid_removed, nloc, &lid_state));
58 if (strict_aggs) PetscCall(PetscMalloc1(nloc, &lid_parent_gid));
59 if (isMPI) {
60 for (kk = 0, gid = my0; kk < nloc; kk++, gid++) lid_gid[kk] = gid;
61 PetscCall(VecGetLocalSize(mpimat->lvec, &num_fine_ghosts));
62 PetscCall(PetscMalloc2(num_fine_ghosts, &cpcol_gid, num_fine_ghosts, &cpcol_state));
63 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)Gmat), &sf));
64 PetscCall(MatGetLayouts(Gmat, &layout, NULL));
65 PetscCall(PetscSFSetGraphLayout(sf, layout, num_fine_ghosts, NULL, PETSC_COPY_VALUES, mpimat->garray));
66 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE));
67 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE));
68 for (kk = 0; kk < num_fine_ghosts; kk++) cpcol_state[kk] = MIS_NOT_DONE;
69 } else num_fine_ghosts = 0;
70
71 /* has ghost nodes for !strict and uses local indexing (yuck) */
72 PetscCall(PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts + nloc, &agg_lists));
73 if (a_locals_llist) *a_locals_llist = agg_lists;
74
75 /* need an inverse map - locals */
76 for (kk = 0; kk < nloc; kk++) {
77 lid_cprowID[kk] = -1;
78 lid_removed[kk] = PETSC_FALSE;
79 if (strict_aggs) lid_parent_gid[kk] = -1.0;
80 lid_state[kk] = MIS_NOT_DONE;
81 }
82 /* set index into cmpressed row 'lid_cprowID' */
83 if (matB) {
84 for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
85 lid = matB->compressedrow.rindex[ix];
86 if (lid >= 0) lid_cprowID[lid] = ix;
87 }
88 }
89 /* MIS */
90 nremoved = nDone = 0;
91
92 PetscCall(ISGetIndices(perm, &perm_ix));
93 while (nDone < nloc || PETSC_TRUE) { /* asynchronous not implemented */
94 /* check all vertices */
95 for (kk = 0; kk < nloc; kk++) {
96 lid = perm_ix[kk];
97 state = lid_state[lid];
98 if (lid_removed[lid]) continue;
99 if (state == MIS_NOT_DONE) {
100 /* parallel test, delete if selected ghost */
101 isOK = PETSC_TRUE;
102 if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
103 ii = matB->compressedrow.i;
104 n = ii[ix + 1] - ii[ix];
105 idx = matB->j + ii[ix];
106 for (j = 0; j < n; j++) {
107 cpid = idx[j]; /* compressed row ID in B mat */
108 gid = cpcol_gid[cpid];
109 statej = cpcol_state[cpid];
110 PetscCheck(!MIS_IS_SELECTED(statej), PETSC_COMM_SELF, PETSC_ERR_SUP, "selected ghost: %" PetscInt_FMT, gid);
111 if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
112 isOK = PETSC_FALSE; /* can not delete */
113 break;
114 }
115 }
116 } /* parallel test */
117 if (isOK) { /* select or remove this vertex */
118 nDone++;
119 /* check for singleton */
120 ii = matA->i;
121 n = ii[lid + 1] - ii[lid];
122 if (n < 2) {
123 /* if I have any ghost adj then not a sing */
124 ix = lid_cprowID[lid];
125 if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) {
126 nremoved++;
127 nrm_tot++;
128 lid_removed[lid] = PETSC_TRUE;
129 continue;
130 // lid_state[lidj] = MIS_REMOVED; /* add singleton to MIS (can cause low rank with elasticity on fine grid) */
131 }
132 }
133 /* SELECTED state encoded with global index */
134 lid_state[lid] = lid + my0;
135 nselected++;
136 if (strict_aggs) {
137 PetscCall(PetscCDAppendID(agg_lists, lid, lid + my0));
138 } else {
139 PetscCall(PetscCDAppendID(agg_lists, lid, lid));
140 }
141 /* delete local adj */
142 idx = matA->j + ii[lid];
143 for (j = 0; j < n; j++) {
144 lidj = idx[j];
145 statej = lid_state[lidj];
146 if (statej == MIS_NOT_DONE) {
147 nDone++;
148 if (strict_aggs) {
149 PetscCall(PetscCDAppendID(agg_lists, lid, lidj + my0));
150 } else {
151 PetscCall(PetscCDAppendID(agg_lists, lid, lidj));
152 }
153 lid_state[lidj] = MIS_DELETED; /* delete this */
154 }
155 }
156 /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
157 if (!strict_aggs) {
158 if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
159 ii = matB->compressedrow.i;
160 n = ii[ix + 1] - ii[ix];
161 idx = matB->j + ii[ix];
162 for (j = 0; j < n; j++) {
163 cpid = idx[j]; /* compressed row ID in B mat */
164 statej = cpcol_state[cpid];
165 if (statej == MIS_NOT_DONE) PetscCall(PetscCDAppendID(agg_lists, lid, nloc + cpid));
166 }
167 }
168 }
169 } /* selected */
170 } /* not done vertex */
171 } /* vertex loop */
172
173 /* update ghost states and count todos */
174 if (isMPI) {
175 /* scatter states, check for done */
176 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE));
177 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE));
178 ii = matB->compressedrow.i;
179 for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
180 lid = matB->compressedrow.rindex[ix]; /* local boundary node */
181 state = lid_state[lid];
182 if (state == MIS_NOT_DONE) {
183 /* look at ghosts */
184 n = ii[ix + 1] - ii[ix];
185 idx = matB->j + ii[ix];
186 for (j = 0; j < n; j++) {
187 cpid = idx[j]; /* compressed row ID in B mat */
188 statej = cpcol_state[cpid];
189 if (MIS_IS_SELECTED(statej)) { /* lid is now deleted, do it */
190 nDone++;
191 lid_state[lid] = MIS_DELETED; /* delete this */
192 if (!strict_aggs) {
193 lidj = nloc + cpid;
194 PetscCall(PetscCDAppendID(agg_lists, lidj, lid));
195 } else {
196 sgid = cpcol_gid[cpid];
197 lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */
198 }
199 break;
200 }
201 }
202 }
203 }
204 /* all done? */
205 t1 = nloc - nDone;
206 PetscCallMPI(MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm)); /* synchronous version */
207 if (!t2) break;
208 } else break; /* all done */
209 } /* outer parallel MIS loop */
210 PetscCall(ISRestoreIndices(perm, &perm_ix));
211 PetscCall(PetscInfo(info_is, "\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices. %" PetscInt_FMT " selected.\n", nremoved, nloc, nselected));
212
213 /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
214 if (strict_aggs && matB) {
215 /* need to copy this to free buffer -- should do this globally */
216 PetscCall(PetscMalloc2(num_fine_ghosts, &cpcol_sel_gid, num_fine_ghosts, &icpcol_gid));
217 for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];
218
219 /* get proc of deleted ghost */
220 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE));
221 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE));
222 for (cpid = 0; cpid < num_fine_ghosts; cpid++) {
223 sgid = cpcol_sel_gid[cpid];
224 gid = icpcol_gid[cpid];
225 if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
226 slid = sgid - my0;
227 PetscCall(PetscCDAppendID(agg_lists, slid, gid));
228 }
229 }
230 PetscCall(PetscFree2(cpcol_sel_gid, icpcol_gid));
231 }
232 if (isMPI) {
233 PetscCall(PetscSFDestroy(&sf));
234 PetscCall(PetscFree2(cpcol_gid, cpcol_state));
235 }
236 PetscCall(PetscFree4(lid_gid, lid_cprowID, lid_removed, lid_state));
237 if (strict_aggs) {
238 // check sizes -- all vertices must get in graph
239 PetscInt aa[2] = {0, nrm_tot}, bb[2], MM;
240
241 PetscCall(PetscFree(lid_parent_gid));
242 PetscCall(MatGetSize(Gmat, &MM, NULL));
243 // check sizes -- all vertices must get in graph
244 PetscCall(PetscCDCount(agg_lists, &aa[0]));
245 PetscCallMPI(MPIU_Allreduce(aa, bb, 2, MPIU_INT, MPI_SUM, comm));
246 if (MM != bb[0]) PetscCall(PetscInfo(info_is, "Warning: N = %" PetscInt_FMT ", sum of aggregates %" PetscInt_FMT ", %" PetscInt_FMT " removed total\n", MM, bb[0], bb[1]));
247 PetscCheck(MM >= bb[0], comm, PETSC_ERR_PLIB, "Sum of aggs is too large");
248 }
249 PetscCall(ISDestroy(&info_is));
250 PetscFunctionReturn(PETSC_SUCCESS);
251 }
252
253 /*
254 MIS coarsen, simple greedy.
255 */
MatCoarsenApply_MIS(MatCoarsen coarse)256 static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
257 {
258 Mat mat = coarse->graph;
259
260 PetscFunctionBegin;
261 if (!coarse->perm) {
262 IS perm;
263 PetscInt n, m;
264 MPI_Comm comm;
265
266 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
267 PetscCall(MatGetLocalSize(mat, &m, &n));
268 PetscCall(ISCreateStride(comm, m, 0, 1, &perm));
269 PetscCall(MatCoarsenApply_MIS_private(perm, mat, coarse->strict_aggs, &coarse->agg_lists));
270 PetscCall(ISDestroy(&perm));
271 } else {
272 PetscCall(MatCoarsenApply_MIS_private(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists));
273 }
274 PetscFunctionReturn(PETSC_SUCCESS);
275 }
276
MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)277 static PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse, PetscViewer viewer)
278 {
279 PetscMPIInt rank;
280 PetscBool isascii;
281 PetscViewerFormat format;
282
283 PetscFunctionBegin;
284 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
285 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
286 PetscCall(PetscViewerGetFormat(viewer, &format));
287 if (isascii && format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
288 if (coarse->agg_lists) {
289 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
290 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d] MIS aggregator\n", rank));
291 if (!rank) {
292 PetscCDIntNd *pos, *pos2;
293 for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) {
294 PetscCall(PetscCDGetHeadPos(coarse->agg_lists, kk, &pos));
295 if ((pos2 = pos)) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "selected %" PetscInt_FMT ": ", kk));
296 while (pos) {
297 PetscInt gid1;
298 PetscCall(PetscCDIntNdGetID(pos, &gid1));
299 PetscCall(PetscCDGetNextPos(coarse->agg_lists, kk, &pos));
300 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " ", gid1));
301 }
302 if (pos2) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
303 }
304 }
305 PetscCall(PetscViewerFlush(viewer));
306 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
307 } else {
308 PetscCall(PetscViewerASCIIPrintf(viewer, " MIS aggregator lists are not available\n"));
309 }
310 }
311 PetscFunctionReturn(PETSC_SUCCESS);
312 }
313
314 /*MC
315 MATCOARSENMIS - Creates a coarsening object that uses a maximal independent set (MIS) algorithm
316
317 Collective
318
319 Input Parameter:
320 . coarse - the coarsen context
321
322 Level: beginner
323
324 .seealso: `MatCoarsen`, `MatCoarsenApply()`, `MatCoarsenGetData()`, `MatCoarsenSetType()`, `MatCoarsenType`
325 M*/
MatCoarsenCreate_MIS(MatCoarsen coarse)326 PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
327 {
328 PetscFunctionBegin;
329 coarse->ops->apply = MatCoarsenApply_MIS;
330 coarse->ops->view = MatCoarsenView_MIS;
331 PetscFunctionReturn(PETSC_SUCCESS);
332 }
333