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