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