Actual source code: mis.c
petsc-3.3-p7 2013-05-11
1:
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>
5: #include <assert.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( const IS perm,
31: const Mat Gmat,
32: const PetscBool strict_aggs,
33: const PetscInt verbose,
34: PetscCoarsenData **a_locals_llist
35: )
36: {
38: PetscBool isMPI;
39: Mat_SeqAIJ *matA, *matB = 0;
40: MPI_Comm wcomm = ((PetscObject)Gmat)->comm;
41: Vec locState, ghostState;
42: PetscInt num_fine_ghosts,kk,n,ix,j,*idx,*ii,iter,Iend,my0,nremoved;
43: Mat_MPIAIJ *mpimat = 0;
44: PetscScalar *cpcol_gid,*cpcol_state;
45: PetscMPIInt mype,npe;
46: const PetscInt *perm_ix;
47: PetscInt nDone, nselected = 0;
48: const PetscInt nloc = Gmat->rmap->n;
49: PetscInt *lid_cprowID, *lid_gid;
50: PetscBool *lid_removed;
51: PetscScalar *lid_parent_gid = PETSC_NULL; /* only used for strict aggs */
52: PetscScalar *lid_state;
53: PetscCoarsenData *agg_lists;
56: MPI_Comm_rank( wcomm, &mype );
57: MPI_Comm_size( wcomm, &npe );
59: /* get submatrices */
60: PetscObjectTypeCompare( (PetscObject)Gmat, MATMPIAIJ, &isMPI );
61: if (isMPI) {
62: mpimat = (Mat_MPIAIJ*)Gmat->data;
63: matA = (Mat_SeqAIJ*)mpimat->A->data;
64: matB = (Mat_SeqAIJ*)mpimat->B->data;
65: /* force compressed storage of B */
66: matB->compressedrow.check = PETSC_TRUE;
67: MatCheckCompressedRow(mpimat->B,&matB->compressedrow,matB->i,Gmat->rmap->n,-1.0);
68: assert( matB->compressedrow.use );
69: } else {
70: PetscBool isAIJ;
71: PetscObjectTypeCompare( (PetscObject)Gmat, MATSEQAIJ, &isAIJ );
72: assert(isAIJ);
73: matA = (Mat_SeqAIJ*)Gmat->data;
74: }
75: assert( matA && !matA->compressedrow.use );
76: assert( matB==0 || matB->compressedrow.use );
77: /* get vector */
78: MatGetVecs( Gmat, &locState, 0 );
80: MatGetOwnershipRange(Gmat,&my0,&Iend);
82: if( mpimat ) {
83: PetscInt gid;
84: for(kk=0,gid=my0;kk<nloc;kk++,gid++) {
85: PetscScalar v = (PetscScalar)(gid);
86: VecSetValues( locState, 1, &gid, &v, INSERT_VALUES ); /* set with GID */
87: }
88: VecAssemblyBegin( locState );
89: VecAssemblyEnd( locState );
90: VecScatterBegin(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
91: VecScatterEnd(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
92: VecGetArray( mpimat->lvec, &cpcol_gid ); /* get proc ID in 'cpcol_gid' */
93: VecDuplicate( mpimat->lvec, &ghostState ); /* need 2nd compressed col. of off proc data */
94: VecGetLocalSize( mpimat->lvec, &num_fine_ghosts );
95: VecSet( ghostState, (PetscScalar)((PetscReal)NOT_DONE) ); /* set with UNKNOWN state */
96: }
97: else num_fine_ghosts = 0;
98:
99: PetscMalloc( nloc*sizeof(PetscInt), &lid_cprowID );
100: PetscMalloc( (nloc+1)*sizeof(PetscInt), &lid_gid ); /* explicit array needed */
101: PetscMalloc( nloc*sizeof(PetscBool), &lid_removed ); /* explicit array needed */
102: if( strict_aggs ) {
103: PetscMalloc( (nloc+1)*sizeof(PetscScalar), &lid_parent_gid );
104: }
105: PetscMalloc( (nloc+1)*sizeof(PetscScalar), &lid_state );
107: /* has ghost nodes for !strict and uses local indexing (yuck) */
108: PetscCDCreate( strict_aggs ? nloc : num_fine_ghosts+nloc, &agg_lists );
109: if(a_locals_llist) *a_locals_llist = agg_lists;
111: /* need an inverse map - locals */
112: for(kk=0;kk<nloc;kk++) {
113: lid_cprowID[kk] = -1; lid_removed[kk] = PETSC_FALSE;
114: if( strict_aggs ) {
115: lid_parent_gid[kk] = -1.0;
116: }
117: lid_gid[kk] = kk + my0;
118: lid_state[kk] = (PetscScalar)((PetscReal)NOT_DONE);
119: }
120: /* set index into cmpressed row 'lid_cprowID' */
121: if( matB ) {
122: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
123: PetscInt lid = matB->compressedrow.rindex[ix];
124: lid_cprowID[lid] = ix;
125: }
126: }
127: /* MIS */
128: iter = nremoved = nDone = 0;
129: ISGetIndices( perm, &perm_ix );
130: while ( nDone < nloc || PETSC_TRUE ) { /* asyncronous not implemented */
131: iter++;
132: if( mpimat ) {
133: VecGetArray( ghostState, &cpcol_state );
134: }
135: /* check all vertices */
136: for(kk=0;kk<nloc;kk++){
137: PetscInt lid = perm_ix[kk];
138: NState state = (NState)PetscRealPart(lid_state[lid]);
139: if(lid_removed[lid]) continue;
140: if( state == NOT_DONE ) {
141: /* parallel test, delete if selected ghost */
142: PetscBool isOK = PETSC_TRUE;
143: if( (ix=lid_cprowID[lid]) != -1 ) { /* if I have any ghost neighbors */
144: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
145: idx = matB->j + ii[ix];
146: for( j=0 ; j<n ; j++ ) {
147: PetscInt cpid = idx[j]; /* compressed row ID in B mat */
148: PetscInt gid = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
149: NState statej = (NState)PetscRealPart(cpcol_state[cpid]);
150: if( statej == NOT_DONE && gid >= Iend ) { /* should be (pe>mype), use gid as pe proxy */
151: isOK = PETSC_FALSE; /* can not delete */
152: break;
153: }
154: else assert( !IS_SELECTED(statej) ); /* lid is now deleted, do it */
155: }
156: } /* parallel test */
157: if( isOK ){ /* select or remove this vertex */
158: nDone++;
159: /* check for singleton */
160: ii = matA->i; n = ii[lid+1] - ii[lid];
161: if( n < 2 ) {
162: /* if I have any ghost adj then not a sing */
163: ix = lid_cprowID[lid];
164: if( ix==-1 || (matB->compressedrow.i[ix+1]-matB->compressedrow.i[ix])==0 ){
165: nremoved++;
166: lid_removed[lid] = PETSC_TRUE;
167: /* should select this because it is technically in the MIS but lets not */
168: /* lid_state[lid] = (PetscScalar)(lid+my0); */
169: continue; /* one local adj (me) and no ghost - singleton */
170: }
171: }
172: /* SELECTED state encoded with global index */
173: lid_state[lid] = (PetscScalar)(lid+my0); /* needed???? */
174: nselected++;
175: if( strict_aggs ) {
176: PetscCDAppendID( agg_lists, lid, lid+my0 );
177: }
178: else {
179: PetscCDAppendID( agg_lists, lid, lid );
180: }
181: /* delete local adj */
182: idx = matA->j + ii[lid];
183: for (j=0; j<n; j++) {
184: PetscInt lidj = idx[j];
185: NState statej = (NState)PetscRealPart(lid_state[lidj]);
186: if( statej == NOT_DONE ){
187: nDone++;
188: /* id_llist[lidj] = id_llist[lid]; id_llist[lid] = lidj; */ /* insert 'lidj' into head of llist */
189: if( strict_aggs ) {
190: PetscCDAppendID( agg_lists, lid, lidj+my0 );
191: }
192: else {
193: PetscCDAppendID( agg_lists, lid, lidj );
194: }
195: lid_state[lidj] = (PetscScalar)(PetscReal)DELETED; /* delete this */
196: }
197: }
198:
199: /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
200: if( !strict_aggs ) {
201: if( (ix=lid_cprowID[lid]) != -1 ) { /* if I have any ghost neighbors */
202: ii = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
203: idx = matB->j + ii[ix];
204: for( j=0 ; j<n ; j++ ) {
205: PetscInt cpid = idx[j]; /* compressed row ID in B mat */
206: NState statej = (NState)PetscRealPart(cpcol_state[cpid]); assert( !IS_SELECTED(statej) );
207: if( statej == NOT_DONE ) {
208: /* cpcol_state[cpid] = (PetscScalar)DELETED; this should happen later ... */
209: /* id_llist[lidj] = id_llist[lid]; id_llist[lid] = lidj; */ /* insert 'lidj' into head of llist */
210: PetscCDAppendID( agg_lists, lid, nloc+cpid );
211: }
212: }
213: }
214: }
215: } /* selected */
216: } /* not done vertex */
217: } /* vertex loop */
218:
219: /* update ghost states and count todos */
220: if( mpimat ) {
221: VecRestoreArray( ghostState, &cpcol_state );
222: /* put lid state in 'locState' */
223: VecSetValues( locState, nloc, lid_gid, lid_state, INSERT_VALUES );
224: VecAssemblyBegin( locState );
225: VecAssemblyEnd( locState );
226: /* scatter states, check for done */
227: VecScatterBegin(mpimat->Mvctx,locState,ghostState,INSERT_VALUES,SCATTER_FORWARD);
228:
229: VecScatterEnd(mpimat->Mvctx,locState,ghostState,INSERT_VALUES,SCATTER_FORWARD);
230:
231: /* delete locals from selected ghosts */
232: VecGetArray( ghostState, &cpcol_state );
233: ii = matB->compressedrow.i;
234: for (ix=0; ix<matB->compressedrow.nrows; ix++) {
235: PetscInt lid = matB->compressedrow.rindex[ix]; /* local boundary node */
236: NState state = (NState)PetscRealPart(lid_state[lid]);
237: if( state == NOT_DONE ) {
238: /* look at ghosts */
239: n = ii[ix+1] - ii[ix];
240: idx = matB->j + ii[ix];
241: for( j=0 ; j<n ; j++ ) {
242: PetscInt cpid = idx[j]; /* compressed row ID in B mat */
243: NState statej = (NState)PetscRealPart(cpcol_state[cpid]);
244: if( IS_SELECTED(statej) ) { /* lid is now deleted, do it */
245: nDone++;
246: lid_state[lid] = (PetscScalar)(PetscReal)DELETED; /* delete this */
247: if( !strict_aggs ) {
248: PetscInt lidj = nloc + cpid;
249: /* id_llist[lid] = id_llist[lidj]; id_llist[lidj] = lid; */ /* insert 'lid' into head of ghost llist */
250: PetscCDAppendID( agg_lists, lidj, lid );
251: }
252: else {
253: PetscInt sgid = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
254: lid_parent_gid[lid] = (PetscScalar)sgid; /* keep track of proc that I belong to */
255: }
256: break;
257: }
258: }
259: }
260: }
261: VecRestoreArray( ghostState, &cpcol_state );
262:
263: /* all done? */
264: {
265: PetscInt t1, t2;
266: t1 = nloc - nDone; assert(t1>=0);
267: MPI_Allreduce( &t1, &t2, 1, MPIU_INT, MPI_SUM, wcomm ); /* synchronous version */
268: if( t2 == 0 ) break;
269: }
270: }
271: else break; /* all done */
272: } /* outer parallel MIS loop */
273: ISRestoreIndices(perm,&perm_ix);
275: if( verbose ) {
276: if( verbose == 1 ) {
277: PetscPrintf(wcomm,"\t[%d]%s removed %d of %d vertices. %d selected.\n",mype,__FUNCT__,nremoved,nloc,nselected);
278: }
279: else {
280: MPI_Allreduce( &nremoved, &n, 1, MPIU_INT, MPI_SUM, wcomm );
281: MatGetSize( Gmat, &kk, &j );
282: MPI_Allreduce( &nselected, &j, 1, MPIU_INT, MPI_SUM, wcomm );
283: PetscPrintf(wcomm,"\t[%d]%s removed %d of %d vertices. (%d local) %d selected.\n",mype,__FUNCT__,n,kk,nremoved,j);
284: }
285: }
287: /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
288: if( strict_aggs && matB ) {
289: PetscScalar *cpcol_sel_gid;
290: PetscInt cpid,*icpcol_gid;
292: /* need to copy this to free buffer -- should do this globaly */
293: PetscMalloc( num_fine_ghosts*sizeof(PetscInt), &icpcol_gid );
294: for(cpid=0; cpid<num_fine_ghosts; cpid++) icpcol_gid[cpid] = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
295:
296: /* get proc of deleted ghost */
297: VecSetValues(locState, nloc, lid_gid, lid_parent_gid, INSERT_VALUES);
298: VecAssemblyBegin( locState );
299: VecAssemblyEnd( locState );
300: VecScatterBegin(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
301: VecScatterEnd(mpimat->Mvctx,locState,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);
302: VecGetArray( mpimat->lvec, &cpcol_sel_gid ); /* has pe that owns ghost */
303: for(cpid=0; cpid<num_fine_ghosts; cpid++) {
304: PetscInt sgid = (PetscInt)PetscRealPart(cpcol_sel_gid[cpid]);
305: PetscInt gid = icpcol_gid[cpid];
306: if( sgid >= my0 && sgid < Iend ) { /* I own this deleted */
307: PetscInt slid = sgid - my0;
308: /* id_llist[lidj] = id_llist[lid]; id_llist[lid] = lidj; */ /* insert 'lidj' into head of llist */
309: PetscCDAppendID( agg_lists, slid, gid );
310: assert(IS_SELECTED( (NState)PetscRealPart(lid_state[slid]) ));
311: }
312: }
313: VecRestoreArray( mpimat->lvec, &cpcol_sel_gid );
314: PetscFree( icpcol_gid );
315: }
316: else if( matB ) {
317: VecRestoreArray( mpimat->lvec, &cpcol_gid );
318: }
320: /* cache IS of removed nodes, use 'lid_gid' */
321: for(kk=n=0,ix=my0;kk<nloc;kk++,ix++) {
322: if( lid_removed[kk] ) lid_gid[n++] = ix;
323: }
324: assert(n==nremoved);
325: PetscCDSetRemovedIS( agg_lists, wcomm, n, lid_gid );
327: PetscFree( lid_cprowID );
328: PetscFree( lid_gid );
329: PetscFree( lid_removed );
330: if( strict_aggs ) {
331: PetscFree( lid_parent_gid );
332: }
333: PetscFree( lid_state );
335: if(mpimat){
336: VecDestroy( &ghostState );
337: }
338: VecDestroy( &locState );
340: return(0);
341: }
343: typedef struct {
344: int dummy;
345: } MatCoarsen_MIS;
346: /*
347: MIS coarsen, simple greedy.
348: */
351: static PetscErrorCode MatCoarsenApply_MIS( MatCoarsen coarse )
352: {
353: /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
354: PetscErrorCode ierr;
355: Mat mat = coarse->graph;
356:
359: if(!coarse->perm) {
360: IS perm;
361: PetscInt n,m;
362: MPI_Comm wcomm = ((PetscObject)mat)->comm;
363: MatGetLocalSize( mat, &m, &n );
364: ISCreateStride( wcomm, m, 0, 1, &perm );
365: maxIndSetAgg( perm, mat, coarse->strict_aggs, coarse->verbose, &coarse->agg_lists );
366:
367: ISDestroy( &perm );
368: }
369: else {
370: maxIndSetAgg(coarse->perm, mat, coarse->strict_aggs, coarse->verbose, &coarse->agg_lists );
371:
372: }
373: return(0);
374: }
378: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)
379: {
380: /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS *)coarse->; */
382: int rank;
383: PetscBool iascii;
387: MPI_Comm_rank(((PetscObject)coarse)->comm,&rank);
388: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
389: if (iascii) {
390: PetscViewerASCIISynchronizedPrintf(viewer," [%d] MIS aggregator\n",rank);
391: PetscViewerFlush(viewer);
392: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
393: }
394: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for this MIS coarsener",
395: ((PetscObject)viewer)->type_name);
397: return(0);
398: }
402: PetscErrorCode MatCoarsenDestroy_MIS ( MatCoarsen coarse )
403: {
404: MatCoarsen_MIS *MIS = (MatCoarsen_MIS *)coarse->subctx;
409: PetscFree(MIS);
410:
411: return(0);
412: }
414: /*MC
415: MATCOARSENMIS - Creates a coarsen context via the external package MIS.
417: Collective on MPI_Comm
419: Input Parameter:
420: . coarse - the coarsen context
422: Options Database Keys:
423: + -mat_coarsen_MIS_xxx -
425: Level: beginner
427: .keywords: Coarsen, create, context
429: .seealso: MatCoarsenSetType(), MatCoarsenType
431: M*/
433: EXTERN_C_BEGIN
436: PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
437: {
439: MatCoarsen_MIS *MIS;
442: PetscNewLog( coarse, MatCoarsen_MIS, &MIS );
443: coarse->subctx = (void*)MIS;
445: coarse->ops->apply = MatCoarsenApply_MIS;
446: coarse->ops->view = MatCoarsenView_MIS;
447: coarse->ops->destroy = MatCoarsenDestroy_MIS;
448: /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MIS; */
449: return(0);
450: }
451: EXTERN_C_END