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: 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;
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;
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;
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 */
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));
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];
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;
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: }
253: /*
254: MIS coarsen, simple greedy.
255: */
256: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
257: {
258: Mat mat = coarse->graph;
260: PetscFunctionBegin;
261: if (!coarse->perm) {
262: IS perm;
263: PetscInt n, m;
264: MPI_Comm comm;
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: }
277: static PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse, PetscViewer viewer)
278: {
279: PetscMPIInt rank;
280: PetscBool iascii;
281: PetscViewerFormat format;
283: PetscFunctionBegin;
284: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
285: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
286: PetscCall(PetscViewerGetFormat(viewer, &format));
287: if (iascii && 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: }
314: /*MC
315: MATCOARSENMIS - Creates a coarsening object that uses a maximal independent set (MIS) algorithm
317: Collective
319: Input Parameter:
320: . coarse - the coarsen context
322: Level: beginner
324: .seealso: `MatCoarsen`, `MatCoarsenApply()`, `MatCoarsenGetData()`, `MatCoarsenSetType()`, `MatCoarsenType`
325: M*/
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: }