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