Actual source code: misk.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 >= 0)

 11: /* ********************************************************************** */
 12: /* edge for priority queue */
 13: typedef struct edge_tag {
 14:   PetscReal weight;
 15:   PetscInt  lid0, gid1, cpid1;
 16: } Edge;

 18: static PetscErrorCode PetscCoarsenDataView_private(PetscCoarsenData *agg_lists, PetscViewer viewer)
 19: {
 20:   PetscCDIntNd *pos, *pos2;

 22:   for (PetscInt kk = 0; kk < agg_lists->size; kk++) {
 23:     PetscCDGetHeadPos(agg_lists, kk, &pos);
 24:     if ((pos2 = pos)) PetscViewerASCIISynchronizedPrintf(viewer, "selected local %d: ", (int)kk);
 25:     while (pos) {
 26:       PetscInt gid1;
 27:       PetscCDIntNdGetID(pos, &gid1);
 28:       PetscCDGetNextPos(agg_lists, kk, &pos);
 29:       PetscViewerASCIISynchronizedPrintf(viewer, " %d ", (int)gid1);
 30:     }
 31:     if (pos2) PetscViewerASCIISynchronizedPrintf(viewer, "\n");
 32:   }
 33:   return 0;
 34: }

 36: /* -------------------------------------------------------------------------- */
 37: /*
 38:   MatCoarsenApply_MISK_private - parallel heavy edge matching

 40:   Input Parameter:
 41:    . perm - permutation
 42:    . Gmat - global matrix of graph (data not defined)

 44:   Output Parameter:
 45:    . a_locals_llist - array of list of local nodes rooted at local node
 46: */
 47: static PetscErrorCode MatCoarsenApply_MISK_private(IS perm, const PetscInt misk, Mat Gmat, PetscCoarsenData **a_locals_llist)
 48: {
 49:   PetscBool   isMPI;
 50:   MPI_Comm    comm;
 51:   PetscMPIInt rank, size;
 52:   Mat         cMat, Prols[5], Rtot;
 53:   PetscScalar one = 1;

 59:   PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPI);
 60:   PetscObjectGetComm((PetscObject)Gmat, &comm);
 61:   MPI_Comm_rank(comm, &rank);
 62:   MPI_Comm_size(comm, &size);
 63:   PetscInfo(Gmat, "misk %d\n", (int)misk);
 64:   /* make a copy of the graph, this gets destroyed in iterates */
 65:   if (misk > 1) MatDuplicate(Gmat, MAT_COPY_VALUES, &cMat);
 66:   else cMat = Gmat;
 67:   for (PetscInt iterIdx = 0; iterIdx < misk; iterIdx++) {
 68:     Mat_SeqAIJ       *matA, *matB = NULL;
 69:     Mat_MPIAIJ       *mpimat = NULL;
 70:     const PetscInt   *perm_ix;
 71:     const PetscInt    nloc = cMat->rmap->n;
 72:     PetscCoarsenData *agg_lists;
 73:     PetscInt         *cpcol_gid = NULL, *cpcol_state, *lid_cprowID, *lid_state, *lid_parent_gid = NULL;
 74:     PetscInt          num_fine_ghosts, kk, n, ix, j, *idx, *ai, Iend, my0, nremoved, gid, lid, cpid, lidj, sgid, t1, t2, slid, nDone, nselected = 0, state;
 75:     PetscBool        *lid_removed, isOK;
 76:     PetscLayout       layout;
 77:     PetscSF           sf;

 79:     if (isMPI) {
 80:       mpimat = (Mat_MPIAIJ *)cMat->data;
 81:       matA   = (Mat_SeqAIJ *)mpimat->A->data;
 82:       matB   = (Mat_SeqAIJ *)mpimat->B->data;
 83:       /* force compressed storage of B */
 84:       MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, cMat->rmap->n, -1.0);
 85:     } else {
 86:       PetscBool isAIJ;
 87:       PetscObjectBaseTypeCompare((PetscObject)cMat, MATSEQAIJ, &isAIJ);
 89:       matA = (Mat_SeqAIJ *)cMat->data;
 90:     }
 91:     MatGetOwnershipRange(cMat, &my0, &Iend);
 92:     if (mpimat) {
 93:       PetscInt *lid_gid;
 94:       PetscMalloc1(nloc, &lid_gid); /* explicit array needed */
 95:       for (kk = 0, gid = my0; kk < nloc; kk++, gid++) lid_gid[kk] = gid;
 96:       VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
 97:       PetscMalloc1(num_fine_ghosts, &cpcol_gid);
 98:       PetscMalloc1(num_fine_ghosts, &cpcol_state);
 99:       PetscSFCreate(PetscObjectComm((PetscObject)cMat), &sf);
100:       MatGetLayouts(cMat, &layout, NULL);
101:       PetscSFSetGraphLayout(sf, layout, num_fine_ghosts, NULL, PETSC_COPY_VALUES, mpimat->garray);
102:       PetscSFBcastBegin(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE);
103:       PetscSFBcastEnd(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE);
104:       for (kk = 0; kk < num_fine_ghosts; kk++) cpcol_state[kk] = MIS_NOT_DONE;
105:       PetscFree(lid_gid);
106:     } else num_fine_ghosts = 0;

108:     PetscMalloc1(nloc, &lid_cprowID);
109:     PetscMalloc1(nloc, &lid_removed); /* explicit array needed */
110:     PetscMalloc1(nloc, &lid_parent_gid);
111:     PetscMalloc1(nloc, &lid_state);

113:     /* the data structure */
114:     PetscCDCreate(nloc, &agg_lists);
115:     /* need an inverse map - locals */
116:     for (kk = 0; kk < nloc; kk++) {
117:       lid_cprowID[kk]    = -1;
118:       lid_removed[kk]    = PETSC_FALSE;
119:       lid_parent_gid[kk] = -1.0;
120:       lid_state[kk]      = MIS_NOT_DONE;
121:     }
122:     /* set index into cmpressed row 'lid_cprowID' */
123:     if (matB) {
124:       for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
125:         lid              = matB->compressedrow.rindex[ix];
126:         lid_cprowID[lid] = ix;
127:       }
128:     }
129:     /* MIS */
130:     nremoved = nDone = 0;
131:     if (!iterIdx) ISGetIndices(perm, &perm_ix); // use permutation on first MIS
132:     else perm_ix = NULL;
133:     while (nDone < nloc || PETSC_TRUE) { /* asynchronous not implemented */
134:       /* check all vertices */
135:       for (kk = 0; kk < nloc; kk++) {
136:         lid   = perm_ix ? perm_ix[kk] : kk;
137:         state = lid_state[lid];
138:         if (lid_removed[lid]) continue;
139:         if (state == MIS_NOT_DONE) {
140:           /* parallel test, delete if selected ghost */
141:           isOK = PETSC_TRUE;
142:           /* parallel test */
143:           if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
144:             ai  = matB->compressedrow.i;
145:             n   = ai[ix + 1] - ai[ix];
146:             idx = matB->j + ai[ix];
147:             for (j = 0; j < n; j++) {
148:               cpid = idx[j]; /* compressed row ID in B mat */
149:               gid  = cpcol_gid[cpid];
150:               if (cpcol_state[cpid] == MIS_NOT_DONE && gid >= Iend) { /* or pe>rank */
151:                 isOK = PETSC_FALSE;                                   /* can not delete */
152:                 break;
153:               }
154:             }
155:           }
156:           if (isOK) { /* select or remove this vertex if it is a true singleton like a BC */
157:             nDone++;
158:             /* check for singleton */
159:             ai = matA->i;
160:             n  = ai[lid + 1] - ai[lid];
161:             if (n < 2) {
162:               /* if I have any ghost adj then not a singleton */
163:               ix = lid_cprowID[lid];
164:               if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) {
165:                 nremoved++;
166:                 lid_removed[lid] = PETSC_TRUE;
167:                 /* should select this because it is technically in the MIS but lets not */
168:                 continue; /* one local adj (me) and no ghost - singleton */
169:               }
170:             }
171:             /* SELECTED state encoded with global index */
172:             lid_state[lid] = nselected; // >= 0  is selected, cache for ordering coarse grid
173:             nselected++;
174:             PetscCDAppendID(agg_lists, lid, lid + my0);
175:             /* delete local adj */
176:             idx = matA->j + ai[lid];
177:             for (j = 0; j < n; j++) {
178:               lidj = idx[j];
179:               if (lid_state[lidj] == MIS_NOT_DONE) {
180:                 nDone++;
181:                 PetscCDAppendID(agg_lists, lid, lidj + my0);
182:                 lid_state[lidj] = MIS_DELETED; /* delete this */
183:               }
184:             }
185:           } /* selected */
186:         }   /* not done vertex */
187:       }     /* vertex loop */

189:       /* update ghost states and count todos */
190:       if (mpimat) {
191:         /* scatter states, check for done */
192:         PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE);
193:         PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE);
194:         ai = matB->compressedrow.i;
195:         for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
196:           lid   = matB->compressedrow.rindex[ix]; /* local boundary node */
197:           state = lid_state[lid];
198:           if (state == MIS_NOT_DONE) {
199:             /* look at ghosts */
200:             n   = ai[ix + 1] - ai[ix];
201:             idx = matB->j + ai[ix];
202:             for (j = 0; j < n; j++) {
203:               cpid = idx[j];                            /* compressed row ID in B mat */
204:               if (MIS_IS_SELECTED(cpcol_state[cpid])) { /* lid is now deleted by ghost */
205:                 nDone++;
206:                 lid_state[lid]      = MIS_DELETED; /* delete this */
207:                 sgid                = cpcol_gid[cpid];
208:                 lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */
209:                 break;
210:               }
211:             }
212:           }
213:         }
214:         /* all done? */
215:         t1 = nloc - nDone;
216:         MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm); /* synchronous version */
217:         if (!t2) break;
218:       } else break; /* no mpi - all done */
219:     }               /* outer parallel MIS loop */
220:     if (!iterIdx) ISRestoreIndices(perm, &perm_ix);
221:     PetscInfo(Gmat, "\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices.  %" PetscInt_FMT " selected.\n", nremoved, nloc, nselected);

223:     /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
224:     if (matB) {
225:       PetscInt *cpcol_sel_gid, *icpcol_gid;
226:       /* need to copy this to free buffer -- should do this globally */
227:       PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid);
228:       PetscMalloc1(num_fine_ghosts, &icpcol_gid);
229:       for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];
230:       /* get proc of deleted ghost */
231:       PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE);
232:       PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE);
233:       for (cpid = 0; cpid < num_fine_ghosts; cpid++) {
234:         sgid = cpcol_sel_gid[cpid];
235:         gid  = icpcol_gid[cpid];
236:         if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
237:           slid = sgid - my0;
238:           PetscCDAppendID(agg_lists, slid, gid);
239:         }
240:       }
241:       // done - cleanup
242:       PetscFree(icpcol_gid);
243:       PetscFree(cpcol_sel_gid);
244:       PetscSFDestroy(&sf);
245:       PetscFree(cpcol_gid);
246:       PetscFree(cpcol_state);
247:     }
248:     PetscFree(lid_cprowID);
249:     PetscFree(lid_removed);
250:     PetscFree(lid_parent_gid);
251:     PetscFree(lid_state);

253:     /* MIS done - make projection matrix - P */
254:     MatType jtype;
255:     MatGetType(Gmat, &jtype);
256:     MatCreate(comm, &Prols[iterIdx]);
257:     MatSetType(Prols[iterIdx], jtype);
258:     MatSetSizes(Prols[iterIdx], nloc, nselected, PETSC_DETERMINE, PETSC_DETERMINE);
259:     MatSeqAIJSetPreallocation(Prols[iterIdx], 1, NULL);
260:     MatMPIAIJSetPreallocation(Prols[iterIdx], 1, NULL, 1, NULL);
261:     {
262:       PetscCDIntNd *pos, *pos2;
263:       PetscInt      colIndex, Iend, fgid;
264:       MatGetOwnershipRangeColumn(Prols[iterIdx], &colIndex, &Iend);
265:       // TODO - order with permutation in lid_selected (reversed)
266:       for (PetscInt lid = 0; lid < agg_lists->size; lid++) {
267:         PetscCDGetHeadPos(agg_lists, lid, &pos);
268:         pos2 = pos;
269:         while (pos) {
270:           PetscCDIntNdGetID(pos, &fgid);
271:           PetscCDGetNextPos(agg_lists, lid, &pos);
272:           MatSetValues(Prols[iterIdx], 1, &fgid, 1, &colIndex, &one, INSERT_VALUES);
273:         }
274:         if (pos2) colIndex++;
275:       }
277:     }
278:     MatAssemblyBegin(Prols[iterIdx], MAT_FINAL_ASSEMBLY);
279:     MatAssemblyEnd(Prols[iterIdx], MAT_FINAL_ASSEMBLY);
280:     /* project to make new graph for next MIS, skip if last */
281:     if (iterIdx < misk - 1) {
282:       Mat new_mat;
283:       MatPtAP(cMat, Prols[iterIdx], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &new_mat);
284:       MatDestroy(&cMat);
285:       cMat = new_mat; // next iter
286:     } else if (cMat != Gmat) MatDestroy(&cMat);
287:     // cleanup
288:     PetscCDDestroy(agg_lists);
289:   } /* MIS-k iteration */
290:   /* make total prolongator Rtot = P_0 * P_1 * ... */
291:   Rtot = Prols[misk - 1]; // compose P then transpose to get R
292:   for (PetscInt iterIdx = misk - 1; iterIdx > 0; iterIdx--) {
293:     Mat P;
294:     MatMatMult(Prols[iterIdx - 1], Rtot, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &P);
295:     MatDestroy(&Prols[iterIdx - 1]);
296:     MatDestroy(&Rtot);
297:     Rtot = P;
298:   }
299:   MatTranspose(Rtot, MAT_INPLACE_MATRIX, &Rtot); // R now
300:   MatViewFromOptions(Rtot, NULL, "-misk_aggregation_view");
301:   /* make aggregates with Rtot - could use Rtot directly in theory but have to go through the aggrate list data structure */
302:   {
303:     PetscInt          Istart, Iend, ncols, NN, MM, jj = 0, max_osz = 0;
304:     const PetscInt    nloc = Gmat->rmap->n;
305:     PetscCoarsenData *agg_lists;
306:     Mat               mat;
307:     PetscCDCreate(nloc, &agg_lists);
308:     *a_locals_llist = agg_lists; // return
309:     MatGetOwnershipRange(Rtot, &Istart, &Iend);
310:     for (int grow = Istart, lid = 0; grow < Iend; grow++, lid++) {
311:       const PetscInt *idx;
312:       MatGetRow(Rtot, grow, &ncols, &idx, NULL);
313:       for (int jj = 0; jj < ncols; jj++) {
314:         PetscInt gcol = idx[jj];
315:         PetscCDAppendID(agg_lists, lid, gcol); // local row, global column
316:       }
317:       MatRestoreRow(Rtot, grow, &ncols, &idx, NULL);
318:     }
319:     MatDestroy(&Rtot);

321:     /* make fake matrix, get largest */
322:     for (int lid = 0; lid < nloc; lid++) {
323:       PetscCDSizeAt(agg_lists, lid, &jj);
324:       if (jj > max_osz) max_osz = jj;
325:     }
326:     MatGetSize(Gmat, &MM, &NN);
327:     if (max_osz > MM - nloc) max_osz = MM - nloc;
328:     MatGetOwnershipRange(Gmat, &Istart, NULL);
329:     MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, max_osz, NULL, &mat);
330:     for (PetscInt lid = 0, gidi = Istart; lid < nloc; lid++, gidi++) {
331:       PetscCDIntNd *pos;
332:       PetscCDGetHeadPos(agg_lists, lid, &pos);
333:       while (pos) {
334:         PetscInt gidj;
335:         PetscCDIntNdGetID(pos, &gidj);
336:         PetscCDGetNextPos(agg_lists, lid, &pos);
337:         if (gidj < Istart || gidj >= Istart + nloc) MatSetValues(mat, 1, &gidi, 1, &gidj, &one, ADD_VALUES);
338:       }
339:     }
340:     MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
341:     MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
342:     PetscCDSetMat(agg_lists, mat);
343:   }

345:   return 0;
346: }

348: /*
349:    Distance k MIS. k is in 'subctx'
350: */
351: static PetscErrorCode MatCoarsenApply_MISK(MatCoarsen coarse)
352: {
353:   Mat      mat = coarse->graph;
354:   PetscInt k;

356:   MatCoarsenMISKGetDistance(coarse, &k);
358:   if (!coarse->perm) {
359:     IS       perm;
360:     PetscInt n, m;

362:     MatGetLocalSize(mat, &m, &n);
363:     ISCreateStride(PetscObjectComm((PetscObject)mat), m, 0, 1, &perm);
364:     MatCoarsenApply_MISK_private(perm, (PetscInt)k, mat, &coarse->agg_lists);
365:     ISDestroy(&perm);
366:   } else {
367:     MatCoarsenApply_MISK_private(coarse->perm, (PetscInt)k, mat, &coarse->agg_lists);
368:   }
369:   return 0;
370: }

372: static PetscErrorCode MatCoarsenView_MISK(MatCoarsen coarse, PetscViewer viewer)
373: {
374:   PetscMPIInt rank;
375:   PetscBool   iascii;

377:   MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank);
378:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
379:   if (iascii) {
380:     PetscViewerASCIIPushSynchronized(viewer);
381:     PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] MISK aggregator\n", rank);
382:     if (!rank) PetscCoarsenDataView_private(coarse->agg_lists, viewer);
383:     PetscViewerFlush(viewer);
384:     PetscViewerASCIIPopSynchronized(viewer);
385:   }
386:   return 0;
387: }

389: static PetscErrorCode MatCoarsenSetFromOptions_MISK(MatCoarsen coarse, PetscOptionItems *PetscOptionsObject)
390: {
391:   PetscInt  k = 1;
392:   PetscBool flg;
393:   PetscOptionsHeadBegin(PetscOptionsObject, "MatCoarsen-MISk options");
394:   PetscOptionsInt("-mat_coarsen_misk_distance", "k distance for MIS", "", k, &k, &flg);
395:   if (flg) coarse->subctx = (void *)(size_t)k;
396:   PetscOptionsHeadEnd();
397:   return 0;
398: }

400: /*MC
401:    MATCOARSENMISK - A coarsener that uses MISK, a simple greedy coarsener

403:    Level: beginner

405:    Options Database Key:
406: .   -mat_coarsen_misk_distance <k> - distance for MIS

408: .seealso: `MatCoarsen`, `MatCoarsenMISKSetDistance()`, `MatCoarsenApply()`, `MatCoarsenSetType()`, `MatCoarsenType`, `MatCoarsenCreate()`
409: M*/

411: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MISK(MatCoarsen coarse)
412: {
413:   coarse->ops->apply          = MatCoarsenApply_MISK;
414:   coarse->ops->view           = MatCoarsenView_MISK;
415:   coarse->subctx              = (void *)(size_t)1;
416:   coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MISK;
417:   return 0;
418: }

420: /*@
421:    MatCoarsenMISKSetDistance - the distance to be used by MISK

423:    Collective

425:    Input Parameters:
426: +   coarsen - the coarsen
427: -   k - the distance

429:    Options Database Key:
430: .   -mat_coarsen_misk_distance <k> - distance for MIS

432:    Level: advanced

434: .seealso: `MATCOARSENMISK`, `MatCoarsen`, `MatCoarseSetFromOptions()`, `MatCoarsenSetType()`, `MatCoarsenRegister()`, `MatCoarsenCreate()`,
435:           `MatCoarsenDestroy()`, `MatCoarsenSetAdjacency()`, `MatCoarsenMISKGetDistance()`
436:           `MatCoarsenGetData()`
437: @*/
438: PetscErrorCode MatCoarsenMISKSetDistance(MatCoarsen crs, PetscInt k)
439: {
440:   crs->subctx = (void *)(size_t)k;
441:   return 0;
442: }

444: /*@
445:    MatCoarsenMISKGetDistance - gets the distance to be used by MISK

447:    Collective

449:    Input Parameter:
450: .   coarsen - the coarsen

452:    Output Parameter:
453: .   k - the distance

455:    Level: advanced

457: .seealso: `MATCOARSENMISK`, `MatCoarsen`, `MatCoarseSetFromOptions()`, `MatCoarsenSetType()`, `MatCoarsenRegister()`, `MatCoarsenCreate()`,
458:           `MatCoarsenDestroy()`, `MatCoarsenSetAdjacency()`, `MatCoarsenMISKGetDistance()`
459:           `MatCoarsenGetData()`
460: @*/
461: PetscErrorCode MatCoarsenMISKGetDistance(MatCoarsen crs, PetscInt *k)
462: {
463:   *k = (PetscInt)(size_t)crs->subctx;
464:   return 0;
465: }