Actual source code: mis.c

petsc-master 2018-06-22
Report Typos and Errors
  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: }