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