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: }