Actual source code: jp.c
2: #include <../src/mat/impls/aij/mpi/mpiaij.h>
3: #include <petscsf.h>
5: typedef struct {
6: PetscSF sf;
7: PetscReal *dwts, *owts;
8: PetscInt *dmask, *omask, *cmask;
9: PetscBool local;
10: } MC_JP;
12: static PetscErrorCode MatColoringDestroy_JP(MatColoring mc)
13: {
14: PetscFunctionBegin;
15: PetscCall(PetscFree(mc->data));
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: static PetscErrorCode MatColoringSetFromOptions_JP(MatColoring mc, PetscOptionItems *PetscOptionsObject)
20: {
21: MC_JP *jp = (MC_JP *)mc->data;
23: PetscFunctionBegin;
24: PetscOptionsHeadBegin(PetscOptionsObject, "JP options");
25: PetscCall(PetscOptionsBool("-mat_coloring_jp_local", "Do an initial coloring of local columns", "", jp->local, &jp->local, NULL));
26: PetscOptionsHeadEnd();
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode MCJPGreatestWeight_Private(MatColoring mc, const PetscReal *weights, PetscReal *maxweights)
31: {
32: MC_JP *jp = (MC_JP *)mc->data;
33: Mat G = mc->mat, dG, oG;
34: PetscBool isSeq, isMPI;
35: Mat_MPIAIJ *aij;
36: Mat_SeqAIJ *daij, *oaij;
37: PetscInt *di, *oi, *dj, *oj;
38: PetscSF sf = jp->sf;
39: PetscLayout layout;
40: PetscInt dn, on;
41: PetscInt i, j, l;
42: PetscReal *dwts = jp->dwts, *owts = jp->owts;
43: PetscInt ncols;
44: const PetscInt *cols;
46: PetscFunctionBegin;
47: PetscCall(PetscObjectTypeCompare((PetscObject)G, MATSEQAIJ, &isSeq));
48: PetscCall(PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPI));
49: PetscCheck(isSeq || isMPI, PetscObjectComm((PetscObject)G), PETSC_ERR_ARG_WRONGSTATE, "MatColoringDegrees requires an MPI/SEQAIJ Matrix");
51: /* get the inner matrix structure */
52: oG = NULL;
53: oi = NULL;
54: oj = NULL;
55: if (isMPI) {
56: aij = (Mat_MPIAIJ *)G->data;
57: dG = aij->A;
58: oG = aij->B;
59: daij = (Mat_SeqAIJ *)dG->data;
60: oaij = (Mat_SeqAIJ *)oG->data;
61: di = daij->i;
62: dj = daij->j;
63: oi = oaij->i;
64: oj = oaij->j;
65: PetscCall(MatGetSize(oG, &dn, &on));
66: if (!sf) {
67: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)mc), &sf));
68: PetscCall(MatGetLayouts(G, &layout, NULL));
69: PetscCall(PetscSFSetGraphLayout(sf, layout, on, NULL, PETSC_COPY_VALUES, aij->garray));
70: jp->sf = sf;
71: }
72: } else {
73: dG = G;
74: PetscCall(MatGetSize(dG, NULL, &dn));
75: daij = (Mat_SeqAIJ *)dG->data;
76: di = daij->i;
77: dj = daij->j;
78: }
79: /* set up the distance-zero weights */
80: if (!dwts) {
81: PetscCall(PetscMalloc1(dn, &dwts));
82: jp->dwts = dwts;
83: if (oG) {
84: PetscCall(PetscMalloc1(on, &owts));
85: jp->owts = owts;
86: }
87: }
88: for (i = 0; i < dn; i++) {
89: maxweights[i] = weights[i];
90: dwts[i] = maxweights[i];
91: }
92: /* get the off-diagonal weights */
93: if (oG) {
94: PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
95: PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, dwts, owts, MPI_REPLACE));
96: PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, dwts, owts, MPI_REPLACE));
97: PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
98: }
99: /* check for the maximum out to the distance of the coloring */
100: for (l = 0; l < mc->dist; l++) {
101: /* check for on-diagonal greater weights */
102: for (i = 0; i < dn; i++) {
103: ncols = di[i + 1] - di[i];
104: cols = &(dj[di[i]]);
105: for (j = 0; j < ncols; j++) {
106: if (dwts[cols[j]] > maxweights[i]) maxweights[i] = dwts[cols[j]];
107: }
108: /* check for off-diagonal greater weights */
109: if (oG) {
110: ncols = oi[i + 1] - oi[i];
111: cols = &(oj[oi[i]]);
112: for (j = 0; j < ncols; j++) {
113: if (owts[cols[j]] > maxweights[i]) maxweights[i] = owts[cols[j]];
114: }
115: }
116: }
117: if (l < mc->dist - 1) {
118: for (i = 0; i < dn; i++) dwts[i] = maxweights[i];
119: if (oG) {
120: PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
121: PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, dwts, owts, MPI_REPLACE));
122: PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, dwts, owts, MPI_REPLACE));
123: PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
124: }
125: }
126: }
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode MCJPInitialLocalColor_Private(MatColoring mc, PetscInt *lperm, ISColoringValue *colors)
131: {
132: PetscInt j, i, s, e, n, bidx, cidx, idx, dist, distance = mc->dist;
133: Mat G = mc->mat, dG, oG;
134: PetscInt *seen;
135: PetscInt *idxbuf;
136: PetscBool *boundary;
137: PetscInt *distbuf;
138: PetscInt *colormask;
139: PetscInt ncols;
140: const PetscInt *cols;
141: PetscBool isSeq, isMPI;
142: Mat_MPIAIJ *aij;
143: Mat_SeqAIJ *daij, *oaij;
144: PetscInt *di, *dj, dn;
145: PetscInt *oi;
147: PetscFunctionBegin;
148: PetscCall(PetscLogEventBegin(MATCOLORING_Local, mc, 0, 0, 0));
149: PetscCall(MatGetOwnershipRange(G, &s, &e));
150: n = e - s;
151: PetscCall(PetscObjectBaseTypeCompare((PetscObject)G, MATSEQAIJ, &isSeq));
152: PetscCall(PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPI));
153: PetscCheck(isSeq || isMPI, PetscObjectComm((PetscObject)G), PETSC_ERR_ARG_WRONGSTATE, "MatColoringDegrees requires an MPI/SEQAIJ Matrix");
155: /* get the inner matrix structure */
156: oG = NULL;
157: oi = NULL;
158: if (isMPI) {
159: aij = (Mat_MPIAIJ *)G->data;
160: dG = aij->A;
161: oG = aij->B;
162: daij = (Mat_SeqAIJ *)dG->data;
163: oaij = (Mat_SeqAIJ *)oG->data;
164: di = daij->i;
165: dj = daij->j;
166: oi = oaij->i;
167: PetscCall(MatGetSize(oG, &dn, NULL));
168: } else {
169: dG = G;
170: PetscCall(MatGetSize(dG, NULL, &dn));
171: daij = (Mat_SeqAIJ *)dG->data;
172: di = daij->i;
173: dj = daij->j;
174: }
175: PetscCall(PetscMalloc5(n, &colormask, n, &seen, n, &idxbuf, n, &distbuf, n, &boundary));
176: for (i = 0; i < dn; i++) {
177: seen[i] = -1;
178: colormask[i] = -1;
179: boundary[i] = PETSC_FALSE;
180: }
181: /* pass one -- figure out which ones are off-boundary in the distance-n sense */
182: if (oG) {
183: for (i = 0; i < dn; i++) {
184: bidx = -1;
185: /* nonempty off-diagonal, so this one is on the boundary */
186: if (oi[i] != oi[i + 1]) {
187: boundary[i] = PETSC_TRUE;
188: continue;
189: }
190: ncols = di[i + 1] - di[i];
191: cols = &(dj[di[i]]);
192: for (j = 0; j < ncols; j++) {
193: bidx++;
194: seen[cols[j]] = i;
195: distbuf[bidx] = 1;
196: idxbuf[bidx] = cols[j];
197: }
198: while (bidx >= 0) {
199: idx = idxbuf[bidx];
200: dist = distbuf[bidx];
201: bidx--;
202: if (dist < distance) {
203: if (oi[idx + 1] != oi[idx]) {
204: boundary[i] = PETSC_TRUE;
205: break;
206: }
207: ncols = di[idx + 1] - di[idx];
208: cols = &(dj[di[idx]]);
209: for (j = 0; j < ncols; j++) {
210: if (seen[cols[j]] != i) {
211: bidx++;
212: seen[cols[j]] = i;
213: idxbuf[bidx] = cols[j];
214: distbuf[bidx] = dist + 1;
215: }
216: }
217: }
218: }
219: }
220: for (i = 0; i < dn; i++) seen[i] = -1;
221: }
222: /* pass two -- color it by looking at nearby vertices and building a mask */
223: for (i = 0; i < dn; i++) {
224: cidx = lperm[i];
225: if (!boundary[cidx]) {
226: bidx = -1;
227: ncols = di[cidx + 1] - di[cidx];
228: cols = &(dj[di[cidx]]);
229: for (j = 0; j < ncols; j++) {
230: bidx++;
231: seen[cols[j]] = cidx;
232: distbuf[bidx] = 1;
233: idxbuf[bidx] = cols[j];
234: }
235: while (bidx >= 0) {
236: idx = idxbuf[bidx];
237: dist = distbuf[bidx];
238: bidx--;
239: /* mask this color */
240: if (colors[idx] < IS_COLORING_MAX) colormask[colors[idx]] = cidx;
241: if (dist < distance) {
242: ncols = di[idx + 1] - di[idx];
243: cols = &(dj[di[idx]]);
244: for (j = 0; j < ncols; j++) {
245: if (seen[cols[j]] != cidx) {
246: bidx++;
247: seen[cols[j]] = cidx;
248: idxbuf[bidx] = cols[j];
249: distbuf[bidx] = dist + 1;
250: }
251: }
252: }
253: }
254: /* find the lowest untaken color */
255: for (j = 0; j < n; j++) {
256: if (colormask[j] != cidx || j >= mc->maxcolors) {
257: colors[cidx] = j;
258: break;
259: }
260: }
261: }
262: }
263: PetscCall(PetscFree5(colormask, seen, idxbuf, distbuf, boundary));
264: PetscCall(PetscLogEventEnd(MATCOLORING_Local, mc, 0, 0, 0));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: static PetscErrorCode MCJPMinColor_Private(MatColoring mc, ISColoringValue maxcolor, const ISColoringValue *colors, ISColoringValue *mincolors)
269: {
270: MC_JP *jp = (MC_JP *)mc->data;
271: Mat G = mc->mat, dG, oG;
272: PetscBool isSeq, isMPI;
273: Mat_MPIAIJ *aij;
274: Mat_SeqAIJ *daij, *oaij;
275: PetscInt *di, *oi, *dj, *oj;
276: PetscSF sf = jp->sf;
277: PetscLayout layout;
278: PetscInt maskrounds, maskbase, maskradix;
279: PetscInt dn, on;
280: PetscInt i, j, l, k;
281: PetscInt *dmask = jp->dmask, *omask = jp->omask, *cmask = jp->cmask, curmask;
282: PetscInt ncols;
283: const PetscInt *cols;
285: PetscFunctionBegin;
286: maskradix = sizeof(PetscInt) * 8;
287: maskrounds = 1 + maxcolor / (maskradix);
288: maskbase = 0;
289: PetscCall(PetscObjectBaseTypeCompare((PetscObject)G, MATSEQAIJ, &isSeq));
290: PetscCall(PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPI));
291: PetscCheck(isSeq || isMPI, PetscObjectComm((PetscObject)G), PETSC_ERR_ARG_WRONGSTATE, "MatColoringDegrees requires an MPI/SEQAIJ Matrix");
293: /* get the inner matrix structure */
294: oG = NULL;
295: oi = NULL;
296: oj = NULL;
297: if (isMPI) {
298: aij = (Mat_MPIAIJ *)G->data;
299: dG = aij->A;
300: oG = aij->B;
301: daij = (Mat_SeqAIJ *)dG->data;
302: oaij = (Mat_SeqAIJ *)oG->data;
303: di = daij->i;
304: dj = daij->j;
305: oi = oaij->i;
306: oj = oaij->j;
307: PetscCall(MatGetSize(oG, &dn, &on));
308: if (!sf) {
309: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)mc), &sf));
310: PetscCall(MatGetLayouts(G, &layout, NULL));
311: PetscCall(PetscSFSetGraphLayout(sf, layout, on, NULL, PETSC_COPY_VALUES, aij->garray));
312: jp->sf = sf;
313: }
314: } else {
315: dG = G;
316: PetscCall(MatGetSize(dG, NULL, &dn));
317: daij = (Mat_SeqAIJ *)dG->data;
318: di = daij->i;
319: dj = daij->j;
320: }
321: for (i = 0; i < dn; i++) mincolors[i] = IS_COLORING_MAX;
322: /* set up the distance-zero mask */
323: if (!dmask) {
324: PetscCall(PetscMalloc1(dn, &dmask));
325: PetscCall(PetscMalloc1(dn, &cmask));
326: jp->cmask = cmask;
327: jp->dmask = dmask;
328: if (oG) {
329: PetscCall(PetscMalloc1(on, &omask));
330: jp->omask = omask;
331: }
332: }
333: /* the number of colors may be more than the number of bits in a PetscInt; take multiple rounds */
334: for (k = 0; k < maskrounds; k++) {
335: for (i = 0; i < dn; i++) {
336: cmask[i] = 0;
337: if (colors[i] < maskbase + maskradix && colors[i] >= maskbase) cmask[i] = 1 << (colors[i] - maskbase);
338: dmask[i] = cmask[i];
339: }
340: if (oG) {
341: PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
342: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dmask, omask, MPI_REPLACE));
343: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dmask, omask, MPI_REPLACE));
344: PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
345: }
346: /* fill in the mask out to the distance of the coloring */
347: for (l = 0; l < mc->dist; l++) {
348: /* fill in the on-and-off diagonal mask */
349: for (i = 0; i < dn; i++) {
350: ncols = di[i + 1] - di[i];
351: cols = &(dj[di[i]]);
352: for (j = 0; j < ncols; j++) cmask[i] = cmask[i] | dmask[cols[j]];
353: if (oG) {
354: ncols = oi[i + 1] - oi[i];
355: cols = &(oj[oi[i]]);
356: for (j = 0; j < ncols; j++) cmask[i] = cmask[i] | omask[cols[j]];
357: }
358: }
359: for (i = 0; i < dn; i++) dmask[i] = cmask[i];
360: if (l < mc->dist - 1) {
361: if (oG) {
362: PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
363: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dmask, omask, MPI_REPLACE));
364: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dmask, omask, MPI_REPLACE));
365: PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
366: }
367: }
368: }
369: /* read through the mask to see if we've discovered an acceptable color for any vertices in this round */
370: for (i = 0; i < dn; i++) {
371: if (mincolors[i] == IS_COLORING_MAX) {
372: curmask = dmask[i];
373: for (j = 0; j < maskradix; j++) {
374: if (curmask % 2 == 0) {
375: mincolors[i] = j + maskbase;
376: break;
377: }
378: curmask = curmask >> 1;
379: }
380: }
381: }
382: /* do the next maskradix colors */
383: maskbase += maskradix;
384: }
385: for (i = 0; i < dn; i++) {
386: if (mincolors[i] == IS_COLORING_MAX) mincolors[i] = maxcolor + 1;
387: }
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode MatColoringApply_JP(MatColoring mc, ISColoring *iscoloring)
392: {
393: MC_JP *jp = (MC_JP *)mc->data;
394: PetscInt i, nadded, nadded_total, nadded_total_old, ntotal, n;
395: PetscInt maxcolor_local = 0, maxcolor_global = 0, *lperm;
396: PetscMPIInt rank;
397: PetscReal *weights, *maxweights;
398: ISColoringValue *color, *mincolor;
400: PetscFunctionBegin;
401: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mc), &rank));
402: PetscCall(PetscLogEventBegin(MATCOLORING_Weights, mc, 0, 0, 0));
403: PetscCall(MatColoringCreateWeights(mc, &weights, &lperm));
404: PetscCall(PetscLogEventEnd(MATCOLORING_Weights, mc, 0, 0, 0));
405: PetscCall(MatGetSize(mc->mat, NULL, &ntotal));
406: PetscCall(MatGetLocalSize(mc->mat, NULL, &n));
407: PetscCall(PetscMalloc1(n, &maxweights));
408: PetscCall(PetscMalloc1(n, &color));
409: PetscCall(PetscMalloc1(n, &mincolor));
410: for (i = 0; i < n; i++) {
411: color[i] = IS_COLORING_MAX;
412: mincolor[i] = 0;
413: }
414: nadded = 0;
415: nadded_total = 0;
416: nadded_total_old = 0;
417: /* compute purely local vertices */
418: if (jp->local) {
419: PetscCall(MCJPInitialLocalColor_Private(mc, lperm, color));
420: for (i = 0; i < n; i++) {
421: if (color[i] < IS_COLORING_MAX) {
422: nadded++;
423: weights[i] = -1;
424: if (color[i] > maxcolor_local) maxcolor_local = color[i];
425: }
426: }
427: PetscCall(MPIU_Allreduce(&nadded, &nadded_total, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
428: PetscCall(MPIU_Allreduce(&maxcolor_local, &maxcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
429: }
430: while (nadded_total < ntotal) {
431: PetscCall(MCJPMinColor_Private(mc, (ISColoringValue)maxcolor_global, color, mincolor));
432: PetscCall(MCJPGreatestWeight_Private(mc, weights, maxweights));
433: for (i = 0; i < n; i++) {
434: /* choose locally maximal vertices; weights less than zero are omitted from the graph */
435: if (weights[i] >= maxweights[i] && weights[i] >= 0.) {
436: /* assign the minimum possible color */
437: if (mc->maxcolors > mincolor[i]) {
438: color[i] = mincolor[i];
439: } else {
440: color[i] = mc->maxcolors;
441: }
442: if (color[i] > maxcolor_local) maxcolor_local = color[i];
443: weights[i] = -1.;
444: nadded++;
445: }
446: }
447: PetscCall(MPIU_Allreduce(&maxcolor_local, &maxcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
448: PetscCall(MPIU_Allreduce(&nadded, &nadded_total, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
449: PetscCheck(nadded_total != nadded_total_old, PetscObjectComm((PetscObject)mc), PETSC_ERR_NOT_CONVERGED, "JP didn't make progress");
450: nadded_total_old = nadded_total;
451: }
452: PetscCall(PetscLogEventBegin(MATCOLORING_ISCreate, mc, 0, 0, 0));
453: PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)mc), maxcolor_global + 1, n, color, PETSC_OWN_POINTER, iscoloring));
454: PetscCall(PetscLogEventEnd(MATCOLORING_ISCreate, mc, 0, 0, 0));
455: PetscCall(PetscFree(jp->dwts));
456: PetscCall(PetscFree(jp->dmask));
457: PetscCall(PetscFree(jp->cmask));
458: PetscCall(PetscFree(jp->owts));
459: PetscCall(PetscFree(jp->omask));
460: PetscCall(PetscFree(weights));
461: PetscCall(PetscFree(lperm));
462: PetscCall(PetscFree(maxweights));
463: PetscCall(PetscFree(mincolor));
464: PetscCall(PetscSFDestroy(&jp->sf));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*MC
469: MATCOLORINGJP - Parallel Jones-Plassmann coloring
471: Level: beginner
473: Options Database Key:
474: . -mat_coloring_jp_local - perform a local coloring before applying the parallel algorithm
476: Notes:
477: This method uses a parallel Luby-style coloring with weights to choose an independent set of processor
478: boundary vertices at each stage that may be assigned colors independently.
480: Supports both distance one and distance two colorings.
482: References:
483: . * - M. Jones and P. Plassmann, "A parallel graph coloring heuristic," SIAM Journal on Scientific Computing, vol. 14, no. 3,
484: pp. 654-669, 1993.
486: .seealso: `MatColoring`, `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`
487: M*/
488: PETSC_EXTERN PetscErrorCode MatColoringCreate_JP(MatColoring mc)
489: {
490: MC_JP *jp;
492: PetscFunctionBegin;
493: PetscCall(PetscNew(&jp));
494: jp->sf = NULL;
495: jp->dmask = NULL;
496: jp->omask = NULL;
497: jp->cmask = NULL;
498: jp->dwts = NULL;
499: jp->owts = NULL;
500: jp->local = PETSC_TRUE;
501: mc->data = jp;
502: mc->ops->apply = MatColoringApply_JP;
503: mc->ops->view = NULL;
504: mc->ops->destroy = MatColoringDestroy_JP;
505: mc->ops->setfromoptions = MatColoringSetFromOptions_JP;
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }