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:    MatCoarsenApply_MIS_private - 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 MatCoarsenApply_MIS_private(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, 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++) lid_gid[kk] = gid;
 58:     VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
 59:     PetscMalloc1(num_fine_ghosts, &cpcol_gid);
 60:     PetscMalloc1(num_fine_ghosts, &cpcol_state);
 61:     PetscSFCreate(PetscObjectComm((PetscObject)Gmat), &sf);
 62:     MatGetLayouts(Gmat, &layout, NULL);
 63:     PetscSFSetGraphLayout(sf, layout, num_fine_ghosts, NULL, PETSC_COPY_VALUES, mpimat->garray);
 64:     PetscSFBcastBegin(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE);
 65:     PetscSFBcastEnd(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE);
 66:     for (kk = 0; kk < num_fine_ghosts; kk++) cpcol_state[kk] = MIS_NOT_DONE;
 67:   } else num_fine_ghosts = 0;

 69:   PetscMalloc1(nloc, &lid_cprowID);
 70:   PetscMalloc1(nloc, &lid_removed); /* explicit array needed */
 71:   if (strict_aggs) PetscMalloc1(nloc, &lid_parent_gid);
 72:   PetscMalloc1(nloc, &lid_state);

 74:   /* has ghost nodes for !strict and uses local indexing (yuck) */
 75:   PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts + nloc, &agg_lists);
 76:   if (a_locals_llist) *a_locals_llist = agg_lists;

 78:   /* need an inverse map - locals */
 79:   for (kk = 0; kk < nloc; kk++) {
 80:     lid_cprowID[kk] = -1;
 81:     lid_removed[kk] = PETSC_FALSE;
 82:     if (strict_aggs) lid_parent_gid[kk] = -1.0;
 83:     lid_state[kk] = MIS_NOT_DONE;
 84:   }
 85:   /* set index into cmpressed row 'lid_cprowID' */
 86:   if (matB) {
 87:     for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
 88:       lid              = matB->compressedrow.rindex[ix];
 89:       lid_cprowID[lid] = ix;
 90:     }
 91:   }
 92:   /* MIS */
 93:   nremoved = nDone = 0;
 94:   ISGetIndices(perm, &perm_ix);
 95:   while (nDone < nloc || PETSC_TRUE) { /* asynchronous not implemented */
 96:     /* check all vertices */
 97:     for (kk = 0; kk < nloc; kk++) {
 98:       lid   = perm_ix[kk];
 99:       state = lid_state[lid];
100:       if (lid_removed[lid]) continue;
101:       if (state == MIS_NOT_DONE) {
102:         /* parallel test, delete if selected ghost */
103:         isOK = PETSC_TRUE;
104:         if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
105:           ii  = matB->compressedrow.i;
106:           n   = ii[ix + 1] - ii[ix];
107:           idx = matB->j + ii[ix];
108:           for (j = 0; j < n; j++) {
109:             cpid   = idx[j]; /* compressed row ID in B mat */
110:             gid    = cpcol_gid[cpid];
111:             statej = cpcol_state[cpid];
112:             if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
113:               isOK = PETSC_FALSE;                        /* can not delete */
114:               break;
115:             }
116:           }
117:         }           /* parallel test */
118:         if (isOK) { /* select or remove this vertex */
119:           nDone++;
120:           /* check for singleton */
121:           ii = matA->i;
122:           n  = ii[lid + 1] - ii[lid];
123:           if (n < 2) {
124:             /* if I have any ghost adj then not a sing */
125:             ix = lid_cprowID[lid];
126:             if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) {
127:               nremoved++;
128:               lid_removed[lid] = PETSC_TRUE;
129:               /* should select this because it is technically in the MIS but lets not */
130:               continue; /* one local adj (me) and no ghost - singleton */
131:             }
132:           }
133:           /* SELECTED state encoded with global index */
134:           lid_state[lid] = lid + my0; /* needed???? */
135:           nselected++;
136:           if (strict_aggs) {
137:             PetscCDAppendID(agg_lists, lid, lid + my0);
138:           } else {
139:             PetscCDAppendID(agg_lists, lid, lid);
140:           }
141:           /* delete local adj */
142:           idx = matA->j + ii[lid];
143:           for (j = 0; j < n; j++) {
144:             lidj   = idx[j];
145:             statej = lid_state[lidj];
146:             if (statej == MIS_NOT_DONE) {
147:               nDone++;
148:               if (strict_aggs) {
149:                 PetscCDAppendID(agg_lists, lid, lidj + my0);
150:               } else {
151:                 PetscCDAppendID(agg_lists, lid, lidj);
152:               }
153:               lid_state[lidj] = MIS_DELETED; /* delete this */
154:             }
155:           }
156:           /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
157:           if (!strict_aggs) {
158:             if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
159:               ii  = matB->compressedrow.i;
160:               n   = ii[ix + 1] - ii[ix];
161:               idx = matB->j + ii[ix];
162:               for (j = 0; j < n; j++) {
163:                 cpid   = idx[j]; /* compressed row ID in B mat */
164:                 statej = cpcol_state[cpid];
165:                 if (statej == MIS_NOT_DONE) PetscCDAppendID(agg_lists, lid, nloc + cpid);
166:               }
167:             }
168:           }
169:         } /* selected */
170:       }   /* not done vertex */
171:     }     /* vertex loop */

173:     /* update ghost states and count todos */
174:     if (mpimat) {
175:       /* scatter states, check for done */
176:       PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE);
177:       PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE);
178:       ii = matB->compressedrow.i;
179:       for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
180:         lid   = matB->compressedrow.rindex[ix]; /* local boundary node */
181:         state = lid_state[lid];
182:         if (state == MIS_NOT_DONE) {
183:           /* look at ghosts */
184:           n   = ii[ix + 1] - ii[ix];
185:           idx = matB->j + ii[ix];
186:           for (j = 0; j < n; j++) {
187:             cpid   = idx[j]; /* compressed row ID in B mat */
188:             statej = cpcol_state[cpid];
189:             if (MIS_IS_SELECTED(statej)) { /* lid is now deleted, do it */
190:               nDone++;
191:               lid_state[lid] = MIS_DELETED; /* delete this */
192:               if (!strict_aggs) {
193:                 lidj = nloc + cpid;
194:                 PetscCDAppendID(agg_lists, lidj, lid);
195:               } else {
196:                 sgid                = cpcol_gid[cpid];
197:                 lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */
198:               }
199:               break;
200:             }
201:           }
202:         }
203:       }
204:       /* all done? */
205:       t1 = nloc - nDone;
206:       MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm); /* synchronous version */
207:       if (!t2) break;
208:     } else break; /* all done */
209:   }               /* outer parallel MIS loop */
210:   ISRestoreIndices(perm, &perm_ix);
211:   PetscInfo(Gmat, "\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices.  %" PetscInt_FMT " selected.\n", nremoved, nloc, nselected);

213:   /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
214:   if (strict_aggs && matB) {
215:     /* need to copy this to free buffer -- should do this globally */
216:     PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid);
217:     PetscMalloc1(num_fine_ghosts, &icpcol_gid);
218:     for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];

220:     /* get proc of deleted ghost */
221:     PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE);
222:     PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE);
223:     for (cpid = 0; cpid < num_fine_ghosts; cpid++) {
224:       sgid = cpcol_sel_gid[cpid];
225:       gid  = icpcol_gid[cpid];
226:       if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
227:         slid = sgid - my0;
228:         PetscCDAppendID(agg_lists, slid, gid);
229:       }
230:     }
231:     PetscFree(icpcol_gid);
232:     PetscFree(cpcol_sel_gid);
233:   }
234:   if (mpimat) {
235:     PetscSFDestroy(&sf);
236:     PetscFree(cpcol_gid);
237:     PetscFree(cpcol_state);
238:   }
239:   PetscFree(lid_cprowID);
240:   PetscFree(lid_gid);
241:   PetscFree(lid_removed);
242:   if (strict_aggs) PetscFree(lid_parent_gid);
243:   PetscFree(lid_state);
244:   return 0;
245: }

247: /*
248:    MIS coarsen, simple greedy.
249: */
250: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
251: {
252:   Mat mat = coarse->graph;

254:   if (!coarse->perm) {
255:     IS       perm;
256:     PetscInt n, m;
257:     MPI_Comm comm;

259:     PetscObjectGetComm((PetscObject)mat, &comm);
260:     MatGetLocalSize(mat, &m, &n);
261:     ISCreateStride(comm, m, 0, 1, &perm);
262:     MatCoarsenApply_MIS_private(perm, mat, coarse->strict_aggs, &coarse->agg_lists);
263:     ISDestroy(&perm);
264:   } else {
265:     MatCoarsenApply_MIS_private(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists);
266:   }
267:   return 0;
268: }

270: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse, PetscViewer viewer)
271: {
272:   PetscMPIInt rank;
273:   PetscBool   iascii;

275:   MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank);
276:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
277:   if (iascii) {
278:     PetscViewerASCIIPushSynchronized(viewer);
279:     PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] MIS aggregator\n", rank);
280:     if (!rank) {
281:       PetscCDIntNd *pos, *pos2;
282:       for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) {
283:         PetscCDGetHeadPos(coarse->agg_lists, kk, &pos);
284:         if ((pos2 = pos)) PetscViewerASCIISynchronizedPrintf(viewer, "selected %d: ", (int)kk);
285:         while (pos) {
286:           PetscInt gid1;
287:           PetscCDIntNdGetID(pos, &gid1);
288:           PetscCDGetNextPos(coarse->agg_lists, kk, &pos);
289:           PetscViewerASCIISynchronizedPrintf(viewer, " %d ", (int)gid1);
290:         }
291:         if (pos2) PetscViewerASCIISynchronizedPrintf(viewer, "\n");
292:       }
293:     }
294:     PetscViewerFlush(viewer);
295:     PetscViewerASCIIPopSynchronized(viewer);
296:   }
297:   return 0;
298: }

300: /*MC
301:    MATCOARSENMIS - Creates a coarsening with a maximal independent set (MIS) algorithm

303:    Collective

305:    Input Parameter:
306: .  coarse - the coarsen context

308:    Level: beginner

310: .seealso: `MatCoarsen`, `MatCoarsenApply()`, `MatCoarsenGetData()`, `MatCoarsenSetType()`, `MatCoarsenType`
311: M*/

313: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
314: {
315:   coarse->ops->apply = MatCoarsenApply_MIS;
316:   coarse->ops->view  = MatCoarsenView_MIS;
317:   return 0;
318: }