Actual source code: mis.c
petsc-3.11.4 2019-09-28
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 - glabal 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: PetscErrorCode ierr;
27: Mat_SeqAIJ *matA,*matB=NULL;
28: Mat_MPIAIJ *mpimat=NULL;
29: MPI_Comm comm;
30: 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;
31: PetscInt *cpcol_gid,*cpcol_state,*lid_cprowID,*lid_gid,*cpcol_sel_gid,*icpcol_gid,*lid_state,*lid_parent_gid=NULL;
32: PetscBool *lid_removed;
33: PetscBool isMPI,isAIJ,isOK;
34: const PetscInt *perm_ix;
35: const PetscInt nloc = Gmat->rmap->n;
36: PetscCoarsenData *agg_lists;
37: PetscLayout layout;
38: PetscSF sf;
41: PetscObjectGetComm((PetscObject)Gmat,&comm);
43: /* get submatrices */
44: 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: MatCheckCompressedRow(mpimat->B,matB->nonzerorowcnt,&matB->compressedrow,matB->i,Gmat->rmap->n,-1.0);
51: } else {
52: PetscObjectBaseTypeCompare((PetscObject)Gmat,MATSEQAIJ,&isAIJ);
53: if (!isAIJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Require AIJ matrix.");
54: matA = (Mat_SeqAIJ*)Gmat->data;
55: }
56: MatGetOwnershipRange(Gmat,&my0,&Iend);
57: PetscMalloc1(nloc,&lid_gid); /* explicit array needed */
58: if (mpimat) {
59: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
60: lid_gid[kk] = gid;
61: }
62: VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
63: PetscMalloc1(num_fine_ghosts,&cpcol_gid);
64: PetscMalloc1(num_fine_ghosts,&cpcol_state);
65: PetscSFCreate(PetscObjectComm((PetscObject)Gmat),&sf);
66: MatGetLayouts(Gmat,&layout,NULL);
67: PetscSFSetGraphLayout(sf,layout,num_fine_ghosts,NULL,PETSC_COPY_VALUES,mpimat->garray);
68: PetscSFBcastBegin(sf,MPIU_INT,lid_gid,cpcol_gid);
69: PetscSFBcastEnd(sf,MPIU_INT,lid_gid,cpcol_gid);
70: for (kk=0;kk<num_fine_ghosts;kk++) {
71: cpcol_state[kk]=MIS_NOT_DONE;
72: }
73: } else num_fine_ghosts = 0;
75: PetscMalloc1(nloc, &lid_cprowID);
76: PetscMalloc1(nloc, &lid_removed); /* explicit array needed */
77: if (strict_aggs) {
78: PetscMalloc1(nloc,&lid_parent_gid);
79: }
80: PetscMalloc1(nloc,&lid_state);
82: /* has ghost nodes for !strict and uses local indexing (yuck) */
83: PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts+nloc, &agg_lists);
84: if (a_locals_llist) *a_locals_llist = agg_lists;
86: /* need an inverse map - locals */
87: for (kk=0; kk<nloc; kk++) {
88: lid_cprowID[kk] = -1; lid_removed[kk] = PETSC_FALSE;
89: if (strict_aggs) {
90: lid_parent_gid[kk] = -1.0;
91: }
92: lid_state[kk] = MIS_NOT_DONE;
93: }
94: /* set index into cmpressed row 'lid_cprowID' */
95: if (matB) {
96: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
97: lid = matB->compressedrow.rindex[ix];
98: lid_cprowID[lid] = ix;
99: }
100: }
101: /* MIS */
102: iter = nremoved = nDone = 0;
103: ISGetIndices(perm, &perm_ix);
104: while (nDone < nloc || PETSC_TRUE) { /* asyncronous not implemented */
105: iter++;
106: /* check all vertices */
107: for (kk=0; kk<nloc; kk++) {
108: lid = perm_ix[kk];
109: state = lid_state[lid];
110: if (lid_removed[lid]) continue;
111: if (state == MIS_NOT_DONE) {
112: /* parallel test, delete if selected ghost */
113: isOK = PETSC_TRUE;
114: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
115: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
116: idx = matB->j + ii[ix];
117: for (j=0; j<n; j++) {
118: cpid = idx[j]; /* compressed row ID in B mat */
119: gid = cpcol_gid[cpid];
120: statej = cpcol_state[cpid];
121: if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
122: isOK = PETSC_FALSE; /* can not delete */
123: break;
124: }
125: }
126: } /* parallel test */
127: if (isOK) { /* select or remove this vertex */
128: nDone++;
129: /* check for singleton */
130: ii = matA->i; n = ii[lid+1] - ii[lid];
131: if (n < 2) {
132: /* if I have any ghost adj then not a sing */
133: ix = lid_cprowID[lid];
134: if (ix==-1 || !(matB->compressedrow.i[ix+1]-matB->compressedrow.i[ix])) {
135: nremoved++;
136: lid_removed[lid] = PETSC_TRUE;
137: /* should select this because it is technically in the MIS but lets not */
138: continue; /* one local adj (me) and no ghost - singleton */
139: }
140: }
141: /* SELECTED state encoded with global index */
142: lid_state[lid] = lid+my0; /* needed???? */
143: nselected++;
144: if (strict_aggs) {
145: PetscCDAppendID(agg_lists, lid, lid+my0);
146: } else {
147: PetscCDAppendID(agg_lists, lid, lid);
148: }
149: /* delete local adj */
150: idx = matA->j + ii[lid];
151: for (j=0; j<n; j++) {
152: lidj = idx[j];
153: statej = lid_state[lidj];
154: if (statej == MIS_NOT_DONE) {
155: nDone++;
156: if (strict_aggs) {
157: PetscCDAppendID(agg_lists, lid, lidj+my0);
158: } else {
159: PetscCDAppendID(agg_lists, lid, lidj);
160: }
161: lid_state[lidj] = MIS_DELETED; /* delete this */
162: }
163: }
164: /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
165: if (!strict_aggs) {
166: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
167: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
168: idx = matB->j + ii[ix];
169: for (j=0; j<n; j++) {
170: cpid = idx[j]; /* compressed row ID in B mat */
171: statej = cpcol_state[cpid];
172: if (statej == MIS_NOT_DONE) {
173: PetscCDAppendID(agg_lists, lid, nloc+cpid);
174: }
175: }
176: }
177: }
178: } /* selected */
179: } /* not done vertex */
180: } /* vertex loop */
182: /* update ghost states and count todos */
183: if (mpimat) {
184: /* scatter states, check for done */
185: PetscSFBcastBegin(sf,MPIU_INT,lid_state,cpcol_state);
186: PetscSFBcastEnd(sf,MPIU_INT,lid_state,cpcol_state);
187: ii = matB->compressedrow.i;
188: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
189: lid = matB->compressedrow.rindex[ix]; /* local boundary node */
190: state = lid_state[lid];
191: if (state == MIS_NOT_DONE) {
192: /* look at ghosts */
193: n = ii[ix+1] - ii[ix];
194: idx = matB->j + ii[ix];
195: for (j=0; j<n; j++) {
196: cpid = idx[j]; /* compressed row ID in B mat */
197: statej = cpcol_state[cpid];
198: if (MIS_IS_SELECTED(statej)) { /* lid is now deleted, do it */
199: nDone++;
200: lid_state[lid] = MIS_DELETED; /* delete this */
201: if (!strict_aggs) {
202: lidj = nloc + cpid;
203: PetscCDAppendID(agg_lists, lidj, lid);
204: } else {
205: sgid = cpcol_gid[cpid];
206: lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */
207: }
208: break;
209: }
210: }
211: }
212: }
213: /* all done? */
214: t1 = nloc - nDone;
215: MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm); /* synchronous version */
216: if (!t2) break;
217: } else break; /* all done */
218: } /* outer parallel MIS loop */
219: ISRestoreIndices(perm,&perm_ix);
220: PetscInfo3(Gmat,"\t removed %D of %D vertices. %D selected.\n",nremoved,nloc,nselected);
222: /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
223: if (strict_aggs && matB) {
224: /* need to copy this to free buffer -- should do this globaly */
225: PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid);
226: PetscMalloc1(num_fine_ghosts, &icpcol_gid);
227: for (cpid=0; cpid<num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];
229: /* get proc of deleted ghost */
230: PetscSFBcastBegin(sf,MPIU_INT,lid_parent_gid,cpcol_sel_gid);
231: PetscSFBcastEnd(sf,MPIU_INT,lid_parent_gid,cpcol_sel_gid);
232: for (cpid=0; cpid<num_fine_ghosts; cpid++) {
233: sgid = cpcol_sel_gid[cpid];
234: gid = icpcol_gid[cpid];
235: if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
236: slid = sgid - my0;
237: PetscCDAppendID(agg_lists, slid, gid);
238: }
239: }
240: PetscFree(icpcol_gid);
241: PetscFree(cpcol_sel_gid);
242: }
243: if (mpimat) {
244: PetscSFDestroy(&sf);
245: PetscFree(cpcol_gid);
246: PetscFree(cpcol_state);
247: }
248: PetscFree(lid_cprowID);
249: PetscFree(lid_gid);
250: PetscFree(lid_removed);
251: if (strict_aggs) {
252: PetscFree(lid_parent_gid);
253: }
254: PetscFree(lid_state);
255: return(0);
256: }
258: typedef struct {
259: int dummy;
260: } MatCoarsen_MIS;
261: /*
262: MIS coarsen, simple greedy.
263: */
264: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
265: {
266: /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
268: Mat mat = coarse->graph;
272: if (!coarse->perm) {
273: IS perm;
274: PetscInt n,m;
275: MPI_Comm comm;
276: PetscObjectGetComm((PetscObject)mat,&comm);
277: MatGetLocalSize(mat, &m, &n);
278: ISCreateStride(comm, m, 0, 1, &perm);
279: maxIndSetAgg(perm, mat, coarse->strict_aggs, &coarse->agg_lists);
280: ISDestroy(&perm);
281: } else {
282: maxIndSetAgg(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists);
283: }
284: return(0);
285: }
287: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)
288: {
289: /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
291: PetscMPIInt rank;
292: PetscBool iascii;
296: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
297: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
298: if (iascii) {
299: PetscViewerASCIIPushSynchronized(viewer);
300: PetscViewerASCIISynchronizedPrintf(viewer," [%d] MIS aggregator\n",rank);
301: PetscViewerFlush(viewer);
302: PetscViewerASCIIPopSynchronized(viewer);
303: }
304: return(0);
305: }
307: PetscErrorCode MatCoarsenDestroy_MIS(MatCoarsen coarse)
308: {
309: MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->subctx;
314: PetscFree(MIS);
315: return(0);
316: }
318: /*MC
319: MATCOARSENMIS - Creates a coarsen context via the external package MIS.
321: Collective on MPI_Comm
323: Input Parameter:
324: . coarse - the coarsen context
326: Options Database Keys:
327: + -mat_coarsen_MIS_xxx -
329: Level: beginner
331: .keywords: Coarsen, create, context
333: .seealso: MatCoarsenSetType(), MatCoarsenType
335: M*/
337: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
338: {
340: MatCoarsen_MIS *MIS;
343: PetscNewLog(coarse,&MIS);
344: coarse->subctx = (void*)MIS;
346: coarse->ops->apply = MatCoarsenApply_MIS;
347: coarse->ops->view = MatCoarsenView_MIS;
348: coarse->ops->destroy = MatCoarsenDestroy_MIS;
349: /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MIS; */
350: return(0);
351: }