Actual source code: mis.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/
3: #include <../src/mat/impls/aij/seq/aij.h>
4: #include <../src/mat/impls/aij/mpi/mpiaij.h>
7: /* typedef enum { NOT_DONE=-2, DELETED=-1, REMOVED=-3 } NState; */
8: /* use int instead of enum to facilitate passing them via Scatters */
9: typedef PetscInt NState;
10: static const NState NOT_DONE=-2;
11: static const NState DELETED =-1;
12: static const NState REMOVED =-3;
13: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)
15: /* -------------------------------------------------------------------------- */
16: /*
17: maxIndSetAgg - parallel maximal independent set (MIS) with data locality info. MatAIJ specific!!!
19: Input Parameter:
20: . perm - serial permutation of rows of local to process in MIS
21: . Gmat - glabal matrix of graph (data not defined)
22: . strict_aggs - flag for whether to keep strict (non overlapping) aggregates in 'llist';
23: . verbose -
24: Output Parameter:
25: . a_selected - IS of selected vertices, includes 'ghost' nodes at end with natural local indices
26: . a_locals_llist - array of list of nodes rooted at selected nodes
27: */
30: PetscErrorCode maxIndSetAgg(IS perm,Mat Gmat,PetscBool strict_aggs,PetscInt verbose,PetscCoarsenData **a_locals_llist)
31: {
32: PetscErrorCode ierr;
33: PetscBool isMPI;
34: Mat_SeqAIJ *matA, *matB = 0;
35: MPI_Comm comm;
36: Vec locState, ghostState;
37: PetscInt num_fine_ghosts,kk,n,ix,j,*idx,*ii,iter,Iend,my0,nremoved;
38: Mat_MPIAIJ *mpimat = 0;
39: PetscScalar *cpcol_gid,*cpcol_state;
40: PetscMPIInt mype,npe;
41: const PetscInt *perm_ix;
42: PetscInt nDone, nselected = 0;
43: const PetscInt nloc = Gmat->rmap->n;
44: PetscInt *lid_cprowID, *lid_gid;
45: PetscBool *lid_removed;
46: PetscScalar *lid_parent_gid = NULL; /* only used for strict aggs */
47: PetscScalar *lid_state;
48: PetscCoarsenData *agg_lists;
51: PetscObjectGetComm((PetscObject)Gmat,&comm);
52: MPI_Comm_rank(comm, &mype);
53: MPI_Comm_size(comm, &npe);
55: /* get submatrices */
56: PetscObjectTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPI);
57: if (isMPI) {
58: mpimat = (Mat_MPIAIJ*)Gmat->data;
59: matA = (Mat_SeqAIJ*)mpimat->A->data;
60: matB = (Mat_SeqAIJ*)mpimat->B->data;
61: /* force compressed storage of B */
62: matB->compressedrow.check = PETSC_TRUE;
63: MatCheckCompressedRow(mpimat->B,&matB->compressedrow,matB->i,Gmat->rmap->n,-1.0);
64: } else {
65: PetscBool isAIJ;
66: PetscObjectTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isAIJ);
67: matA = (Mat_SeqAIJ*)Gmat->data;
68: }
69: /* get vector */
70: MatGetVecs(Gmat, &locState, 0);
72: MatGetOwnershipRange(Gmat,&my0,&Iend);
74: if (mpimat) {
75: PetscInt gid;
76: for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
77: PetscScalar v = (PetscScalar)(gid);
78: VecSetValues(locState, 1, &gid, &v, INSERT_VALUES); /* set with GID */
79: }
80: VecAssemblyBegin(locState);
81: VecAssemblyEnd(locState);
82: VecScatterBegin(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
83: VecScatterEnd(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
84: VecGetArray(mpimat->lvec, &cpcol_gid); /* get proc ID in 'cpcol_gid' */
85: VecDuplicate(mpimat->lvec, &ghostState); /* need 2nd compressed col. of off proc data */
86: VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
87: VecSet(ghostState, (PetscScalar)((PetscReal)NOT_DONE)); /* set with UNKNOWN state */
88: } else num_fine_ghosts = 0;
90: PetscMalloc(nloc*sizeof(PetscInt), &lid_cprowID);
91: PetscMalloc((nloc+1)*sizeof(PetscInt), &lid_gid); /* explicit array needed */
92: PetscMalloc(nloc*sizeof(PetscBool), &lid_removed); /* explicit array needed */
93: if (strict_aggs) {
94: PetscMalloc((nloc+1)*sizeof(PetscScalar), &lid_parent_gid);
95: }
96: PetscMalloc((nloc+1)*sizeof(PetscScalar), &lid_state);
98: /* has ghost nodes for !strict and uses local indexing (yuck) */
99: PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts+nloc, &agg_lists);
100: if (a_locals_llist) *a_locals_llist = agg_lists;
102: /* need an inverse map - locals */
103: for (kk=0; kk<nloc; kk++) {
104: lid_cprowID[kk] = -1; lid_removed[kk] = PETSC_FALSE;
105: if (strict_aggs) {
106: lid_parent_gid[kk] = -1.0;
107: }
108: lid_gid[kk] = kk + my0;
109: lid_state[kk] = (PetscScalar)((PetscReal)NOT_DONE);
110: }
111: /* set index into cmpressed row 'lid_cprowID' */
112: if (matB) {
113: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
114: PetscInt lid = matB->compressedrow.rindex[ix];
115: lid_cprowID[lid] = ix;
116: }
117: }
118: /* MIS */
119: iter = nremoved = nDone = 0;
120: ISGetIndices(perm, &perm_ix);
121: while (nDone < nloc || PETSC_TRUE) { /* asyncronous not implemented */
122: iter++;
123: if (mpimat) {
124: VecGetArray(ghostState, &cpcol_state);
125: }
126: /* check all vertices */
127: for (kk=0; kk<nloc; kk++) {
128: PetscInt lid = perm_ix[kk];
129: NState state = (NState)PetscRealPart(lid_state[lid]);
130: if (lid_removed[lid]) continue;
131: if (state == NOT_DONE) {
132: /* parallel test, delete if selected ghost */
133: PetscBool isOK = PETSC_TRUE;
134: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
135: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
136: idx = matB->j + ii[ix];
137: for (j=0; j<n; j++) {
138: PetscInt cpid = idx[j]; /* compressed row ID in B mat */
139: PetscInt gid = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
140: NState statej = (NState)PetscRealPart(cpcol_state[cpid]);
141: if (statej == NOT_DONE && gid >= Iend) { /* should be (pe>mype), use gid as pe proxy */
142: isOK = PETSC_FALSE; /* can not delete */
143: break;
144: }
145: }
146: } /* parallel test */
147: if (isOK) { /* select or remove this vertex */
148: nDone++;
149: /* check for singleton */
150: ii = matA->i; n = ii[lid+1] - ii[lid];
151: if (n < 2) {
152: /* if I have any ghost adj then not a sing */
153: ix = lid_cprowID[lid];
154: if (ix==-1 || (matB->compressedrow.i[ix+1]-matB->compressedrow.i[ix])==0) {
155: nremoved++;
156: lid_removed[lid] = PETSC_TRUE;
157: /* should select this because it is technically in the MIS but lets not */
158: /* lid_state[lid] = (PetscScalar)(lid+my0); */
159: continue; /* one local adj (me) and no ghost - singleton */
160: }
161: }
162: /* SELECTED state encoded with global index */
163: lid_state[lid] = (PetscScalar)(lid+my0); /* needed???? */
164: nselected++;
165: if (strict_aggs) {
166: PetscCDAppendID(agg_lists, lid, lid+my0);
167: } else {
168: PetscCDAppendID(agg_lists, lid, lid);
169: }
170: /* delete local adj */
171: idx = matA->j + ii[lid];
172: for (j=0; j<n; j++) {
173: PetscInt lidj = idx[j];
174: NState statej = (NState)PetscRealPart(lid_state[lidj]);
175: if (statej == NOT_DONE) {
176: nDone++;
177: /* id_llist[lidj] = id_llist[lid]; id_llist[lid] = lidj; */ /* insert 'lidj' into head of llist */
178: if (strict_aggs) {
179: PetscCDAppendID(agg_lists, lid, lidj+my0);
180: } else {
181: PetscCDAppendID(agg_lists, lid, lidj);
182: }
183: lid_state[lidj] = (PetscScalar)(PetscReal)DELETED; /* delete this */
184: }
185: }
187: /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
188: if (!strict_aggs) {
189: if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
190: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
191: idx = matB->j + ii[ix];
192: for (j=0; j<n; j++) {
193: PetscInt cpid = idx[j]; /* compressed row ID in B mat */
194: NState statej = (NState)PetscRealPart(cpcol_state[cpid]);
195: if (statej == NOT_DONE) {
196: /* cpcol_state[cpid] = (PetscScalar)DELETED; this should happen later ... */
197: /* id_llist[lidj] = id_llist[lid]; id_llist[lid] = lidj; */ /* insert 'lidj' into head of llist */
198: PetscCDAppendID(agg_lists, lid, nloc+cpid);
199: }
200: }
201: }
202: }
203: } /* selected */
204: } /* not done vertex */
205: } /* vertex loop */
207: /* update ghost states and count todos */
208: if (mpimat) {
209: VecRestoreArray(ghostState, &cpcol_state);
210: /* put lid state in 'locState' */
211: VecSetValues(locState, nloc, lid_gid, lid_state, INSERT_VALUES);
212: VecAssemblyBegin(locState);
213: VecAssemblyEnd(locState);
214: /* scatter states, check for done */
215: VecScatterBegin(mpimat->Mvctx,locState,ghostState,INSERT_VALUES,SCATTER_FORWARD);
216: VecScatterEnd(mpimat->Mvctx,locState,ghostState,INSERT_VALUES,SCATTER_FORWARD);
217: /* delete locals from selected ghosts */
218: VecGetArray(ghostState, &cpcol_state);
219: ii = matB->compressedrow.i;
220: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
221: PetscInt lid = matB->compressedrow.rindex[ix]; /* local boundary node */
222: NState state = (NState)PetscRealPart(lid_state[lid]);
223: if (state == NOT_DONE) {
224: /* look at ghosts */
225: n = ii[ix+1] - ii[ix];
226: idx = matB->j + ii[ix];
227: for (j=0; j<n; j++) {
228: PetscInt cpid = idx[j]; /* compressed row ID in B mat */
229: NState statej = (NState)PetscRealPart(cpcol_state[cpid]);
230: if (IS_SELECTED(statej)) { /* lid is now deleted, do it */
231: nDone++;
232: lid_state[lid] = (PetscScalar)(PetscReal)DELETED; /* delete this */
233: if (!strict_aggs) {
234: PetscInt lidj = nloc + cpid;
235: /* id_llist[lid] = id_llist[lidj]; id_llist[lidj] = lid; */ /* insert 'lid' into head of ghost llist */
236: PetscCDAppendID(agg_lists, lidj, lid);
237: } else {
238: PetscInt sgid = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
239: lid_parent_gid[lid] = (PetscScalar)sgid; /* keep track of proc that I belong to */
240: }
241: break;
242: }
243: }
244: }
245: }
246: VecRestoreArray(ghostState, &cpcol_state);
248: /* all done? */
249: {
250: PetscInt t1, t2;
251: t1 = nloc - nDone;
252: MPI_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm); /* synchronous version */
253: if (t2 == 0) break;
254: }
255: } else break; /* all done */
256: } /* outer parallel MIS loop */
257: ISRestoreIndices(perm,&perm_ix);
259: if (verbose) {
260: if (verbose == 1) {
261: PetscPrintf(comm,"\t[%d]%s removed %d of %d vertices. %d selected.\n",mype,__FUNCT__,nremoved,nloc,nselected);
262: } else {
263: MPI_Allreduce(&nremoved, &n, 1, MPIU_INT, MPI_SUM, comm);
264: MatGetSize(Gmat, &kk, &j);
265: MPI_Allreduce(&nselected, &j, 1, MPIU_INT, MPI_SUM, comm);
266: PetscPrintf(comm,"\t[%d]%s removed %d of %d vertices. (%d local) %d selected.\n",mype,__FUNCT__,n,kk,nremoved,j);
267: }
268: }
270: /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
271: if (strict_aggs && matB) {
272: PetscScalar *cpcol_sel_gid;
273: PetscInt cpid,*icpcol_gid;
275: /* need to copy this to free buffer -- should do this globaly */
276: PetscMalloc(num_fine_ghosts*sizeof(PetscInt), &icpcol_gid);
277: for (cpid=0; cpid<num_fine_ghosts; cpid++) icpcol_gid[cpid] = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
279: /* get proc of deleted ghost */
280: VecSetValues(locState, nloc, lid_gid, lid_parent_gid, INSERT_VALUES);
281: VecAssemblyBegin(locState);
282: VecAssemblyEnd(locState);
283: VecScatterBegin(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
284: VecScatterEnd(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
285: VecGetArray(mpimat->lvec, &cpcol_sel_gid); /* has pe that owns ghost */
286: for (cpid=0; cpid<num_fine_ghosts; cpid++) {
287: PetscInt sgid = (PetscInt)PetscRealPart(cpcol_sel_gid[cpid]);
288: PetscInt gid = icpcol_gid[cpid];
289: if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
290: PetscInt slid = sgid - my0;
291: /* id_llist[lidj] = id_llist[lid]; id_llist[lid] = lidj; */ /* insert 'lidj' into head of llist */
292: PetscCDAppendID(agg_lists, slid, gid);
293: }
294: }
295: VecRestoreArray(mpimat->lvec, &cpcol_sel_gid);
296: PetscFree(icpcol_gid);
297: } else if (matB) {
298: VecRestoreArray(mpimat->lvec, &cpcol_gid);
299: }
301: /* cache IS of removed nodes, use 'lid_gid' */
302: /* for (kk=n=0,ix=my0;kk<nloc;kk++,ix++) { */
303: /* if (lid_removed[kk]) lid_gid[n++] = ix; */
304: /* } */
305: /* PetscCDSetRemovedIS(agg_lists, comm, n, lid_gid); */
307: PetscFree(lid_cprowID);
308: PetscFree(lid_gid);
309: PetscFree(lid_removed);
310: if (strict_aggs) {
311: PetscFree(lid_parent_gid);
312: }
313: PetscFree(lid_state);
315: if (mpimat) {
316: VecDestroy(&ghostState);
317: }
318: VecDestroy(&locState);
319: return(0);
320: }
322: typedef struct {
323: int dummy;
324: } MatCoarsen_MIS;
325: /*
326: MIS coarsen, simple greedy.
327: */
330: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
331: {
332: /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
334: Mat mat = coarse->graph;
338: if (!coarse->perm) {
339: IS perm;
340: PetscInt n,m;
341: MPI_Comm comm;
342: PetscObjectGetComm((PetscObject)mat,&comm);
343: MatGetLocalSize(mat, &m, &n);
344: ISCreateStride(comm, m, 0, 1, &perm);
345: maxIndSetAgg(perm, mat, coarse->strict_aggs, coarse->verbose, &coarse->agg_lists);
346: ISDestroy(&perm);
347: } else {
348: maxIndSetAgg(coarse->perm, mat, coarse->strict_aggs, coarse->verbose, &coarse->agg_lists);
349: }
350: return(0);
351: }
355: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)
356: {
357: /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
359: int rank;
360: PetscBool iascii;
364: MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
365: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
366: if (iascii) {
367: PetscViewerASCIISynchronizedPrintf(viewer," [%d] MIS aggregator\n",rank);
368: PetscViewerFlush(viewer);
369: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
370: }
371: return(0);
372: }
376: PetscErrorCode MatCoarsenDestroy_MIS(MatCoarsen coarse)
377: {
378: MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->subctx;
383: PetscFree(MIS);
384: return(0);
385: }
387: /*MC
388: MATCOARSENMIS - Creates a coarsen context via the external package MIS.
390: Collective on MPI_Comm
392: Input Parameter:
393: . coarse - the coarsen context
395: Options Database Keys:
396: + -mat_coarsen_MIS_xxx -
398: Level: beginner
400: .keywords: Coarsen, create, context
402: .seealso: MatCoarsenSetType(), MatCoarsenType
404: M*/
408: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
409: {
411: MatCoarsen_MIS *MIS;
414: PetscNewLog(coarse, MatCoarsen_MIS, &MIS);
415: coarse->subctx = (void*)MIS;
417: coarse->ops->apply = MatCoarsenApply_MIS;
418: coarse->ops->view = MatCoarsenView_MIS;
419: coarse->ops->destroy = MatCoarsenDestroy_MIS;
420: /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MIS; */
421: return(0);
422: }