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