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: /*
13: MatCoarsenApply_MIS_private - parallel maximal independent set (MIS) with data locality info. MatAIJ specific!!!
15: Input Parameter:
16: . perm - serial permutation of rows of local to process in MIS
17: . Gmat - global matrix of graph (data not defined)
18: . strict_aggs - flag for whether to keep strict (non overlapping) aggregates in 'llist';
20: Output Parameter:
21: . a_selected - IS of selected vertices, includes 'ghost' nodes at end with natural local indices
22: . a_locals_llist - array of list of nodes rooted at selected nodes
23: */
24: PetscErrorCode MatCoarsenApply_MIS_private(IS perm, Mat Gmat, PetscBool strict_aggs, PetscCoarsenData **a_locals_llist)
25: {
26: Mat_SeqAIJ *matA, *matB = NULL;
27: Mat_MPIAIJ *mpimat = NULL;
28: MPI_Comm comm;
29: 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;
30: PetscInt *cpcol_gid, *cpcol_state, *lid_cprowID, *lid_gid, *cpcol_sel_gid, *icpcol_gid, *lid_state, *lid_parent_gid = NULL;
31: PetscBool *lid_removed;
32: PetscBool isMPI, isAIJ, isOK;
33: const PetscInt *perm_ix;
34: const PetscInt nloc = Gmat->rmap->n;
35: PetscCoarsenData *agg_lists;
36: PetscLayout layout;
37: PetscSF sf;
39: PetscObjectGetComm((PetscObject)Gmat, &comm);
41: /* get submatrices */
42: PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPI);
43: if (isMPI) {
44: mpimat = (Mat_MPIAIJ *)Gmat->data;
45: matA = (Mat_SeqAIJ *)mpimat->A->data;
46: matB = (Mat_SeqAIJ *)mpimat->B->data;
47: /* force compressed storage of B */
48: MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, Gmat->rmap->n, -1.0);
49: } else {
50: PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isAIJ);
52: matA = (Mat_SeqAIJ *)Gmat->data;
53: }
54: MatGetOwnershipRange(Gmat, &my0, &Iend);
55: PetscMalloc1(nloc, &lid_gid); /* explicit array needed */
56: if (mpimat) {
57: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) lid_gid[kk] = gid;
58: VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
59: PetscMalloc1(num_fine_ghosts, &cpcol_gid);
60: PetscMalloc1(num_fine_ghosts, &cpcol_state);
61: PetscSFCreate(PetscObjectComm((PetscObject)Gmat), &sf);
62: MatGetLayouts(Gmat, &layout, NULL);
63: PetscSFSetGraphLayout(sf, layout, num_fine_ghosts, NULL, PETSC_COPY_VALUES, mpimat->garray);
64: PetscSFBcastBegin(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE);
65: PetscSFBcastEnd(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE);
66: for (kk = 0; kk < num_fine_ghosts; kk++) cpcol_state[kk] = MIS_NOT_DONE;
67: } else num_fine_ghosts = 0;
69: PetscMalloc1(nloc, &lid_cprowID);
70: PetscMalloc1(nloc, &lid_removed); /* explicit array needed */
71: if (strict_aggs) PetscMalloc1(nloc, &lid_parent_gid);
72: PetscMalloc1(nloc, &lid_state);
74: /* has ghost nodes for !strict and uses local indexing (yuck) */
75: PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts + nloc, &agg_lists);
76: if (a_locals_llist) *a_locals_llist = agg_lists;
78: /* need an inverse map - locals */
79: for (kk = 0; kk < nloc; kk++) {
80: lid_cprowID[kk] = -1;
81: lid_removed[kk] = PETSC_FALSE;
82: if (strict_aggs) lid_parent_gid[kk] = -1.0;
83: lid_state[kk] = MIS_NOT_DONE;
84: }
85: /* set index into cmpressed row 'lid_cprowID' */
86: if (matB) {
87: for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
88: lid = matB->compressedrow.rindex[ix];
89: lid_cprowID[lid] = ix;
90: }
91: }
92: /* MIS */
93: nremoved = nDone = 0;
94: ISGetIndices(perm, &perm_ix);
95: while (nDone < nloc || PETSC_TRUE) { /* asynchronous not implemented */
96: /* check all vertices */
97: for (kk = 0; kk < nloc; kk++) {
98: lid = perm_ix[kk];
99: state = lid_state[lid];
100: if (lid_removed[lid]) continue;
101: if (state == MIS_NOT_DONE) {
102: /* parallel test, delete if selected ghost */
103: isOK = PETSC_TRUE;
104: if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
105: ii = matB->compressedrow.i;
106: n = ii[ix + 1] - ii[ix];
107: idx = matB->j + ii[ix];
108: for (j = 0; j < n; j++) {
109: cpid = idx[j]; /* compressed row ID in B mat */
110: gid = cpcol_gid[cpid];
111: statej = cpcol_state[cpid];
112: if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
113: isOK = PETSC_FALSE; /* can not delete */
114: break;
115: }
116: }
117: } /* parallel test */
118: if (isOK) { /* select or remove this vertex */
119: nDone++;
120: /* check for singleton */
121: ii = matA->i;
122: n = ii[lid + 1] - ii[lid];
123: if (n < 2) {
124: /* if I have any ghost adj then not a sing */
125: ix = lid_cprowID[lid];
126: if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) {
127: nremoved++;
128: lid_removed[lid] = PETSC_TRUE;
129: /* should select this because it is technically in the MIS but lets not */
130: continue; /* one local adj (me) and no ghost - singleton */
131: }
132: }
133: /* SELECTED state encoded with global index */
134: lid_state[lid] = lid + my0; /* needed???? */
135: nselected++;
136: if (strict_aggs) {
137: PetscCDAppendID(agg_lists, lid, lid + my0);
138: } else {
139: 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: PetscCDAppendID(agg_lists, lid, lidj + my0);
150: } else {
151: 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) 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 (mpimat) {
175: /* scatter states, check for done */
176: PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE);
177: 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: 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: 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: ISRestoreIndices(perm, &perm_ix);
211: PetscInfo(Gmat, "\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: PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid);
217: PetscMalloc1(num_fine_ghosts, &icpcol_gid);
218: for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];
220: /* get proc of deleted ghost */
221: PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE);
222: PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE);
223: for (cpid = 0; cpid < num_fine_ghosts; cpid++) {
224: sgid = cpcol_sel_gid[cpid];
225: gid = icpcol_gid[cpid];
226: if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
227: slid = sgid - my0;
228: PetscCDAppendID(agg_lists, slid, gid);
229: }
230: }
231: PetscFree(icpcol_gid);
232: PetscFree(cpcol_sel_gid);
233: }
234: if (mpimat) {
235: PetscSFDestroy(&sf);
236: PetscFree(cpcol_gid);
237: PetscFree(cpcol_state);
238: }
239: PetscFree(lid_cprowID);
240: PetscFree(lid_gid);
241: PetscFree(lid_removed);
242: if (strict_aggs) PetscFree(lid_parent_gid);
243: PetscFree(lid_state);
244: return 0;
245: }
247: /*
248: MIS coarsen, simple greedy.
249: */
250: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
251: {
252: Mat mat = coarse->graph;
254: if (!coarse->perm) {
255: IS perm;
256: PetscInt n, m;
257: MPI_Comm comm;
259: PetscObjectGetComm((PetscObject)mat, &comm);
260: MatGetLocalSize(mat, &m, &n);
261: ISCreateStride(comm, m, 0, 1, &perm);
262: MatCoarsenApply_MIS_private(perm, mat, coarse->strict_aggs, &coarse->agg_lists);
263: ISDestroy(&perm);
264: } else {
265: MatCoarsenApply_MIS_private(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists);
266: }
267: return 0;
268: }
270: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse, PetscViewer viewer)
271: {
272: PetscMPIInt rank;
273: PetscBool iascii;
275: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank);
276: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
277: if (iascii) {
278: PetscViewerASCIIPushSynchronized(viewer);
279: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] MIS aggregator\n", rank);
280: if (!rank) {
281: PetscCDIntNd *pos, *pos2;
282: for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) {
283: PetscCDGetHeadPos(coarse->agg_lists, kk, &pos);
284: if ((pos2 = pos)) PetscViewerASCIISynchronizedPrintf(viewer, "selected %d: ", (int)kk);
285: while (pos) {
286: PetscInt gid1;
287: PetscCDIntNdGetID(pos, &gid1);
288: PetscCDGetNextPos(coarse->agg_lists, kk, &pos);
289: PetscViewerASCIISynchronizedPrintf(viewer, " %d ", (int)gid1);
290: }
291: if (pos2) PetscViewerASCIISynchronizedPrintf(viewer, "\n");
292: }
293: }
294: PetscViewerFlush(viewer);
295: PetscViewerASCIIPopSynchronized(viewer);
296: }
297: return 0;
298: }
300: /*MC
301: MATCOARSENMIS - Creates a coarsening with a maximal independent set (MIS) algorithm
303: Collective
305: Input Parameter:
306: . coarse - the coarsen context
308: Level: beginner
310: .seealso: `MatCoarsen`, `MatCoarsenApply()`, `MatCoarsenGetData()`, `MatCoarsenSetType()`, `MatCoarsenType`
311: M*/
313: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
314: {
315: coarse->ops->apply = MatCoarsenApply_MIS;
316: coarse->ops->view = MatCoarsenView_MIS;
317: return 0;
318: }