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