Actual source code: mis.c

petsc-3.4.4 2014-03-13
  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: }