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(PetscOptionItems *PetscOptionsObject,MatColoring mc)
19: {
20: MC_JP *jp = (MC_JP*)mc->data;
22: PetscOptionsHead(PetscOptionsObject,"JP options");
23: PetscOptionsBool("-mat_coloring_jp_local","Do an initial coloring of local columns","",jp->local,&jp->local,NULL);
24: PetscOptionsTail();
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++) {
116: dwts[i] = maxweights[i];
117: }
118: if (oG) {
119: PetscLogEventBegin(MATCOLORING_Comm,mc,0,0,0);
120: PetscSFBcastBegin(sf,MPIU_REAL,dwts,owts,MPI_REPLACE);
121: PetscSFBcastEnd(sf,MPIU_REAL,dwts,owts,MPI_REPLACE);
122: PetscLogEventEnd(MATCOLORING_Comm,mc,0,0,0);
123: }
124: }
125: }
126: return 0;
127: }
129: static PetscErrorCode MCJPInitialLocalColor_Private(MatColoring mc,PetscInt *lperm,ISColoringValue *colors)
130: {
131: PetscInt j,i,s,e,n,bidx,cidx,idx,dist,distance=mc->dist;
132: Mat G=mc->mat,dG,oG;
133: PetscInt *seen;
134: PetscInt *idxbuf;
135: PetscBool *boundary;
136: PetscInt *distbuf;
137: PetscInt *colormask;
138: PetscInt ncols;
139: const PetscInt *cols;
140: PetscBool isSeq,isMPI;
141: Mat_MPIAIJ *aij;
142: Mat_SeqAIJ *daij,*oaij;
143: PetscInt *di,*dj,dn;
144: PetscInt *oi;
146: PetscLogEventBegin(MATCOLORING_Local,mc,0,0,0);
147: MatGetOwnershipRange(G,&s,&e);
148: n=e-s;
149: PetscObjectBaseTypeCompare((PetscObject)G,MATSEQAIJ,&isSeq);
150: PetscObjectTypeCompare((PetscObject)G,MATMPIAIJ,&isMPI);
153: /* get the inner matrix structure */
154: oG = NULL;
155: oi = NULL;
156: if (isMPI) {
157: aij = (Mat_MPIAIJ*)G->data;
158: dG = aij->A;
159: oG = aij->B;
160: daij = (Mat_SeqAIJ*)dG->data;
161: oaij = (Mat_SeqAIJ*)oG->data;
162: di = daij->i;
163: dj = daij->j;
164: oi = oaij->i;
165: MatGetSize(oG,&dn,NULL);
166: } else {
167: dG = G;
168: MatGetSize(dG,NULL,&dn);
169: daij = (Mat_SeqAIJ*)dG->data;
170: di = daij->i;
171: dj = daij->j;
172: }
173: PetscMalloc5(n,&colormask,n,&seen,n,&idxbuf,n,&distbuf,n,&boundary);
174: for (i=0;i<dn;i++) {
175: seen[i]=-1;
176: colormask[i] = -1;
177: boundary[i] = PETSC_FALSE;
178: }
179: /* pass one -- figure out which ones are off-boundary in the distance-n sense */
180: if (oG) {
181: for (i=0;i<dn;i++) {
182: bidx=-1;
183: /* nonempty off-diagonal, so this one is on the boundary */
184: if (oi[i]!=oi[i+1]) {
185: boundary[i] = PETSC_TRUE;
186: continue;
187: }
188: ncols = di[i+1]-di[i];
189: cols = &(dj[di[i]]);
190: for (j=0;j<ncols;j++) {
191: bidx++;
192: seen[cols[j]] = i;
193: distbuf[bidx] = 1;
194: idxbuf[bidx] = cols[j];
195: }
196: while (bidx >= 0) {
197: idx = idxbuf[bidx];
198: dist = distbuf[bidx];
199: bidx--;
200: if (dist < distance) {
201: if (oi[idx+1]!=oi[idx]) {
202: boundary[i] = PETSC_TRUE;
203: break;
204: }
205: ncols = di[idx+1]-di[idx];
206: cols = &(dj[di[idx]]);
207: for (j=0;j<ncols;j++) {
208: if (seen[cols[j]] != i) {
209: bidx++;
210: seen[cols[j]] = i;
211: idxbuf[bidx] = cols[j];
212: distbuf[bidx] = dist+1;
213: }
214: }
215: }
216: }
217: }
218: for (i=0;i<dn;i++) {
219: seen[i]=-1;
220: }
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) {
241: colormask[colors[idx]] = cidx;
242: }
243: if (dist < distance) {
244: ncols = di[idx+1]-di[idx];
245: cols = &(dj[di[idx]]);
246: for (j=0;j<ncols;j++) {
247: if (seen[cols[j]] != cidx) {
248: bidx++;
249: seen[cols[j]] = cidx;
250: idxbuf[bidx] = cols[j];
251: distbuf[bidx] = dist+1;
252: }
253: }
254: }
255: }
256: /* find the lowest untaken color */
257: for (j=0;j<n;j++) {
258: if (colormask[j] != cidx || j >= mc->maxcolors) {
259: colors[cidx] = j;
260: break;
261: }
262: }
263: }
264: }
265: PetscFree5(colormask,seen,idxbuf,distbuf,boundary);
266: PetscLogEventEnd(MATCOLORING_Local,mc,0,0,0);
267: return 0;
268: }
270: static PetscErrorCode MCJPMinColor_Private(MatColoring mc,ISColoringValue maxcolor,const ISColoringValue *colors,ISColoringValue *mincolors)
271: {
272: MC_JP *jp = (MC_JP*)mc->data;
273: Mat G=mc->mat,dG,oG;
274: PetscBool isSeq,isMPI;
275: Mat_MPIAIJ *aij;
276: Mat_SeqAIJ *daij,*oaij;
277: PetscInt *di,*oi,*dj,*oj;
278: PetscSF sf=jp->sf;
279: PetscLayout layout;
280: PetscInt maskrounds,maskbase,maskradix;
281: PetscInt dn,on;
282: PetscInt i,j,l,k;
283: PetscInt *dmask=jp->dmask,*omask=jp->omask,*cmask=jp->cmask,curmask;
284: PetscInt ncols;
285: const PetscInt *cols;
287: maskradix = sizeof(PetscInt)*8;
288: maskrounds = 1 + maxcolor / (maskradix);
289: maskbase = 0;
290: PetscObjectBaseTypeCompare((PetscObject)G,MATSEQAIJ,&isSeq);
291: PetscObjectTypeCompare((PetscObject)G,MATMPIAIJ,&isMPI);
294: /* get the inner matrix structure */
295: oG = NULL;
296: oi = NULL;
297: oj = NULL;
298: if (isMPI) {
299: aij = (Mat_MPIAIJ*)G->data;
300: dG = aij->A;
301: oG = aij->B;
302: daij = (Mat_SeqAIJ*)dG->data;
303: oaij = (Mat_SeqAIJ*)oG->data;
304: di = daij->i;
305: dj = daij->j;
306: oi = oaij->i;
307: oj = oaij->j;
308: MatGetSize(oG,&dn,&on);
309: if (!sf) {
310: PetscSFCreate(PetscObjectComm((PetscObject)mc),&sf);
311: MatGetLayouts(G,&layout,NULL);
312: PetscSFSetGraphLayout(sf,layout,on,NULL,PETSC_COPY_VALUES,aij->garray);
313: jp->sf = sf;
314: }
315: } else {
316: dG = G;
317: MatGetSize(dG,NULL,&dn);
318: daij = (Mat_SeqAIJ*)dG->data;
319: di = daij->i;
320: dj = daij->j;
321: }
322: for (i=0;i<dn;i++) {
323: mincolors[i] = IS_COLORING_MAX;
324: }
325: /* set up the distance-zero mask */
326: if (!dmask) {
327: PetscMalloc1(dn,&dmask);
328: PetscMalloc1(dn,&cmask);
329: jp->cmask = cmask;
330: jp->dmask = dmask;
331: if (oG) {
332: PetscMalloc1(on,&omask);
333: jp->omask = omask;
334: }
335: }
336: /* the number of colors may be more than the number of bits in a PetscInt; take multiple rounds */
337: for (k=0;k<maskrounds;k++) {
338: for (i=0;i<dn;i++) {
339: cmask[i] = 0;
340: if (colors[i] < maskbase+maskradix && colors[i] >= maskbase)
341: cmask[i] = 1 << (colors[i]-maskbase);
342: dmask[i] = cmask[i];
343: }
344: if (oG) {
345: PetscLogEventBegin(MATCOLORING_Comm,mc,0,0,0);
346: PetscSFBcastBegin(sf,MPIU_INT,dmask,omask,MPI_REPLACE);
347: PetscSFBcastEnd(sf,MPIU_INT,dmask,omask,MPI_REPLACE);
348: PetscLogEventEnd(MATCOLORING_Comm,mc,0,0,0);
349: }
350: /* fill in the mask out to the distance of the coloring */
351: for (l=0;l<mc->dist;l++) {
352: /* fill in the on-and-off diagonal mask */
353: for (i=0;i<dn;i++) {
354: ncols = di[i+1]-di[i];
355: cols = &(dj[di[i]]);
356: for (j=0;j<ncols;j++) {
357: cmask[i] = cmask[i] | dmask[cols[j]];
358: }
359: if (oG) {
360: ncols = oi[i+1]-oi[i];
361: cols = &(oj[oi[i]]);
362: for (j=0;j<ncols;j++) {
363: cmask[i] = cmask[i] | omask[cols[j]];
364: }
365: }
366: }
367: for (i=0;i<dn;i++) {
368: dmask[i]=cmask[i];
369: }
370: if (l < mc->dist-1) {
371: if (oG) {
372: PetscLogEventBegin(MATCOLORING_Comm,mc,0,0,0);
373: PetscSFBcastBegin(sf,MPIU_INT,dmask,omask,MPI_REPLACE);
374: PetscSFBcastEnd(sf,MPIU_INT,dmask,omask,MPI_REPLACE);
375: PetscLogEventEnd(MATCOLORING_Comm,mc,0,0,0);
376: }
377: }
378: }
379: /* read through the mask to see if we've discovered an acceptable color for any vertices in this round */
380: for (i=0;i<dn;i++) {
381: if (mincolors[i] == IS_COLORING_MAX) {
382: curmask = dmask[i];
383: for (j=0;j<maskradix;j++) {
384: if (curmask % 2 == 0) {
385: mincolors[i] = j+maskbase;
386: break;
387: }
388: curmask = curmask >> 1;
389: }
390: }
391: }
392: /* do the next maskradix colors */
393: maskbase += maskradix;
394: }
395: for (i=0;i<dn;i++) {
396: if (mincolors[i] == IS_COLORING_MAX) {
397: mincolors[i] = maxcolor+1;
398: }
399: }
400: return 0;
401: }
403: static PetscErrorCode MatColoringApply_JP(MatColoring mc,ISColoring *iscoloring)
404: {
405: MC_JP *jp = (MC_JP*)mc->data;
406: PetscInt i,nadded,nadded_total,nadded_total_old,ntotal,n,round;
407: PetscInt maxcolor_local=0,maxcolor_global = 0,*lperm;
408: PetscMPIInt rank;
409: PetscReal *weights,*maxweights;
410: ISColoringValue *color,*mincolor;
412: MPI_Comm_rank(PetscObjectComm((PetscObject)mc),&rank);
413: PetscLogEventBegin(MATCOLORING_Weights,mc,0,0,0);
414: MatColoringCreateWeights(mc,&weights,&lperm);
415: PetscLogEventEnd(MATCOLORING_Weights,mc,0,0,0);
416: MatGetSize(mc->mat,NULL,&ntotal);
417: MatGetLocalSize(mc->mat,NULL,&n);
418: PetscMalloc1(n,&maxweights);
419: PetscMalloc1(n,&color);
420: PetscMalloc1(n,&mincolor);
421: for (i=0;i<n;i++) {
422: color[i] = IS_COLORING_MAX;
423: mincolor[i] = 0;
424: }
425: nadded=0;
426: nadded_total=0;
427: nadded_total_old=0;
428: /* compute purely local vertices */
429: if (jp->local) {
430: MCJPInitialLocalColor_Private(mc,lperm,color);
431: for (i=0;i<n;i++) {
432: if (color[i] < IS_COLORING_MAX) {
433: nadded++;
434: weights[i] = -1;
435: if (color[i] > maxcolor_local) maxcolor_local = color[i];
436: }
437: }
438: MPIU_Allreduce(&nadded,&nadded_total,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)mc));
439: MPIU_Allreduce(&maxcolor_local,&maxcolor_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)mc));
440: }
441: round = 0;
442: while (nadded_total < ntotal) {
443: MCJPMinColor_Private(mc,(ISColoringValue)maxcolor_global,color,mincolor);
444: MCJPGreatestWeight_Private(mc,weights,maxweights);
445: for (i=0;i<n;i++) {
446: /* choose locally maximal vertices; weights less than zero are omitted from the graph */
447: if (weights[i] >= maxweights[i] && weights[i] >= 0.) {
448: /* assign the minimum possible color */
449: if (mc->maxcolors > mincolor[i]) {
450: color[i] = mincolor[i];
451: } else {
452: color[i] = mc->maxcolors;
453: }
454: if (color[i] > maxcolor_local) maxcolor_local = color[i];
455: weights[i] = -1.;
456: nadded++;
457: }
458: }
459: MPIU_Allreduce(&maxcolor_local,&maxcolor_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)mc));
460: MPIU_Allreduce(&nadded,&nadded_total,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)mc));
462: nadded_total_old = nadded_total;
463: round++;
464: }
465: PetscLogEventBegin(MATCOLORING_ISCreate,mc,0,0,0);
466: ISColoringCreate(PetscObjectComm((PetscObject)mc),maxcolor_global+1,n,color,PETSC_OWN_POINTER,iscoloring);
467: PetscLogEventEnd(MATCOLORING_ISCreate,mc,0,0,0);
468: PetscFree(jp->dwts);
469: PetscFree(jp->dmask);
470: PetscFree(jp->cmask);
471: PetscFree(jp->owts);
472: PetscFree(jp->omask);
473: PetscFree(weights);
474: PetscFree(lperm);
475: PetscFree(maxweights);
476: PetscFree(mincolor);
477: PetscSFDestroy(&jp->sf);
478: return 0;
479: }
481: /*MC
482: MATCOLORINGJP - Parallel Jones-Plassmann Coloring
484: Level: beginner
486: Notes:
487: This method uses a parallel Luby-style coloring with weights to choose an independent set of processor
488: boundary vertices at each stage that may be assigned colors independently.
490: Supports both distance one and distance two colorings.
492: References:
493: . * - M. Jones and P. Plassmann, "A parallel graph coloring heuristic," SIAM Journal on Scientific Computing, vol. 14, no. 3,
494: pp. 654-669, 1993.
496: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType()
497: M*/
498: PETSC_EXTERN PetscErrorCode MatColoringCreate_JP(MatColoring mc)
499: {
500: MC_JP *jp;
502: PetscNewLog(mc,&jp);
503: jp->sf = NULL;
504: jp->dmask = NULL;
505: jp->omask = NULL;
506: jp->cmask = NULL;
507: jp->dwts = NULL;
508: jp->owts = NULL;
509: jp->local = PETSC_TRUE;
510: mc->data = jp;
511: mc->ops->apply = MatColoringApply_JP;
512: mc->ops->view = NULL;
513: mc->ops->destroy = MatColoringDestroy_JP;
514: mc->ops->setfromoptions = MatColoringSetFromOptions_JP;
515: return 0;
516: }