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: maxIndSetAgg - 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 maxIndSetAgg(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,iter,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++) {
58: lid_gid[kk] = gid;
59: }
60: VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
61: PetscMalloc1(num_fine_ghosts,&cpcol_gid);
62: PetscMalloc1(num_fine_ghosts,&cpcol_state);
63: PetscSFCreate(PetscObjectComm((PetscObject)Gmat),&sf);
64: MatGetLayouts(Gmat,&layout,NULL);
65: PetscSFSetGraphLayout(sf,layout,num_fine_ghosts,NULL,PETSC_COPY_VALUES,mpimat->garray);
66: PetscSFBcastBegin(sf,MPIU_INT,lid_gid,cpcol_gid,MPI_REPLACE);
67: PetscSFBcastEnd(sf,MPIU_INT,lid_gid,cpcol_gid,MPI_REPLACE);
68: for (kk=0;kk<num_fine_ghosts;kk++) {
69: cpcol_state[kk]=MIS_NOT_DONE;
70: }
71: } else num_fine_ghosts = 0;
73: PetscMalloc1(nloc, &lid_cprowID);
74: PetscMalloc1(nloc, &lid_removed); /* explicit array needed */
75: if (strict_aggs) {
76: PetscMalloc1(nloc,&lid_parent_gid);
77: }
78: PetscMalloc1(nloc,&lid_state);
80: /* has ghost nodes for !strict and uses local indexing (yuck) */
81: PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts+nloc, &agg_lists);
82: if (a_locals_llist) *a_locals_llist = agg_lists;
84: /* need an inverse map - locals */
85: for (kk=0; kk<nloc; kk++) {
86: lid_cprowID[kk] = -1; lid_removed[kk] = PETSC_FALSE;
87: if (strict_aggs) {
88: lid_parent_gid[kk] = -1.0;
89: }
90: lid_state[kk] = MIS_NOT_DONE;
91: }
92: /* set index into cmpressed row 'lid_cprowID' */
93: if (matB) {
94: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
95: lid = matB->compressedrow.rindex[ix];
96: lid_cprowID[lid] = ix;
97: }
98: }
99: /* MIS */
100: iter = nremoved = nDone = 0;
101: ISGetIndices(perm, &perm_ix);
102: while (nDone < nloc || PETSC_TRUE) { /* asyncronous not implemented */
103: iter++;
104: /* check all vertices */
105: for (kk=0; kk<nloc; kk++) {
106: lid = perm_ix[kk];
107: state = lid_state[lid];
108: if (lid_removed[lid]) continue;
109: if (state == MIS_NOT_DONE) {
110: /* parallel test, delete if selected ghost */
111: isOK = PETSC_TRUE;
112: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
113: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
114: idx = matB->j + ii[ix];
115: for (j=0; j<n; j++) {
116: cpid = idx[j]; /* compressed row ID in B mat */
117: gid = cpcol_gid[cpid];
118: statej = cpcol_state[cpid];
119: if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
120: isOK = PETSC_FALSE; /* can not delete */
121: break;
122: }
123: }
124: } /* parallel test */
125: if (isOK) { /* select or remove this vertex */
126: nDone++;
127: /* check for singleton */
128: ii = matA->i; n = ii[lid+1] - ii[lid];
129: if (n < 2) {
130: /* if I have any ghost adj then not a sing */
131: ix = lid_cprowID[lid];
132: if (ix==-1 || !(matB->compressedrow.i[ix+1]-matB->compressedrow.i[ix])) {
133: nremoved++;
134: lid_removed[lid] = PETSC_TRUE;
135: /* should select this because it is technically in the MIS but lets not */
136: continue; /* one local adj (me) and no ghost - singleton */
137: }
138: }
139: /* SELECTED state encoded with global index */
140: lid_state[lid] = lid+my0; /* needed???? */
141: nselected++;
142: if (strict_aggs) {
143: PetscCDAppendID(agg_lists, lid, lid+my0);
144: } else {
145: PetscCDAppendID(agg_lists, lid, lid);
146: }
147: /* delete local adj */
148: idx = matA->j + ii[lid];
149: for (j=0; j<n; j++) {
150: lidj = idx[j];
151: statej = lid_state[lidj];
152: if (statej == MIS_NOT_DONE) {
153: nDone++;
154: if (strict_aggs) {
155: PetscCDAppendID(agg_lists, lid, lidj+my0);
156: } else {
157: PetscCDAppendID(agg_lists, lid, lidj);
158: }
159: lid_state[lidj] = MIS_DELETED; /* delete this */
160: }
161: }
162: /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
163: if (!strict_aggs) {
164: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
165: ii = matB->compressedrow.i; 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) {
171: PetscCDAppendID(agg_lists, lid, nloc+cpid);
172: }
173: }
174: }
175: }
176: } /* selected */
177: } /* not done vertex */
178: } /* vertex loop */
180: /* update ghost states and count todos */
181: if (mpimat) {
182: /* scatter states, check for done */
183: PetscSFBcastBegin(sf,MPIU_INT,lid_state,cpcol_state,MPI_REPLACE);
184: PetscSFBcastEnd(sf,MPIU_INT,lid_state,cpcol_state,MPI_REPLACE);
185: ii = matB->compressedrow.i;
186: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
187: lid = matB->compressedrow.rindex[ix]; /* local boundary node */
188: state = lid_state[lid];
189: if (state == MIS_NOT_DONE) {
190: /* look at ghosts */
191: n = ii[ix+1] - ii[ix];
192: idx = matB->j + ii[ix];
193: for (j=0; j<n; j++) {
194: cpid = idx[j]; /* compressed row ID in B mat */
195: statej = cpcol_state[cpid];
196: if (MIS_IS_SELECTED(statej)) { /* lid is now deleted, do it */
197: nDone++;
198: lid_state[lid] = MIS_DELETED; /* delete this */
199: if (!strict_aggs) {
200: lidj = nloc + cpid;
201: PetscCDAppendID(agg_lists, lidj, lid);
202: } else {
203: sgid = cpcol_gid[cpid];
204: lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */
205: }
206: break;
207: }
208: }
209: }
210: }
211: /* all done? */
212: t1 = nloc - nDone;
213: MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm); /* synchronous version */
214: if (!t2) break;
215: } else break; /* all done */
216: } /* outer parallel MIS loop */
217: ISRestoreIndices(perm,&perm_ix);
218: PetscInfo(Gmat,"\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices. %" PetscInt_FMT " selected.\n",nremoved,nloc,nselected);
220: /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
221: if (strict_aggs && matB) {
222: /* need to copy this to free buffer -- should do this globaly */
223: PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid);
224: PetscMalloc1(num_fine_ghosts, &icpcol_gid);
225: for (cpid=0; cpid<num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];
227: /* get proc of deleted ghost */
228: PetscSFBcastBegin(sf,MPIU_INT,lid_parent_gid,cpcol_sel_gid,MPI_REPLACE);
229: PetscSFBcastEnd(sf,MPIU_INT,lid_parent_gid,cpcol_sel_gid,MPI_REPLACE);
230: for (cpid=0; cpid<num_fine_ghosts; cpid++) {
231: sgid = cpcol_sel_gid[cpid];
232: gid = icpcol_gid[cpid];
233: if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
234: slid = sgid - my0;
235: PetscCDAppendID(agg_lists, slid, gid);
236: }
237: }
238: PetscFree(icpcol_gid);
239: PetscFree(cpcol_sel_gid);
240: }
241: if (mpimat) {
242: PetscSFDestroy(&sf);
243: PetscFree(cpcol_gid);
244: PetscFree(cpcol_state);
245: }
246: PetscFree(lid_cprowID);
247: PetscFree(lid_gid);
248: PetscFree(lid_removed);
249: if (strict_aggs) {
250: PetscFree(lid_parent_gid);
251: }
252: PetscFree(lid_state);
253: return 0;
254: }
256: /*
257: MIS coarsen, simple greedy.
258: */
259: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
260: {
261: Mat mat = coarse->graph;
263: if (!coarse->perm) {
264: IS perm;
265: PetscInt n,m;
266: MPI_Comm comm;
268: PetscObjectGetComm((PetscObject)mat,&comm);
269: MatGetLocalSize(mat, &m, &n);
270: ISCreateStride(comm, m, 0, 1, &perm);
271: maxIndSetAgg(perm, mat, coarse->strict_aggs, &coarse->agg_lists);
272: ISDestroy(&perm);
273: } else {
274: maxIndSetAgg(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists);
275: }
276: return 0;
277: }
279: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)
280: {
281: PetscMPIInt rank;
282: PetscBool iascii;
284: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
285: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
286: if (iascii) {
287: PetscViewerASCIIPushSynchronized(viewer);
288: PetscViewerASCIISynchronizedPrintf(viewer," [%d] MIS aggregator\n",rank);
289: PetscViewerFlush(viewer);
290: PetscViewerASCIIPopSynchronized(viewer);
291: }
292: return 0;
293: }
295: /*MC
296: MATCOARSENMIS - Creates a coarsen context via the external package MIS.
298: Collective
300: Input Parameter:
301: . coarse - the coarsen context
303: Options Database Keys:
304: . -mat_coarsen_MIS_xxx -
306: Level: beginner
308: .seealso: MatCoarsenSetType(), MatCoarsenType
310: M*/
312: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
313: {
314: coarse->ops->apply = MatCoarsenApply_MIS;
315: coarse->ops->view = MatCoarsenView_MIS;
316: /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MIS; */
317: return 0;
318: }