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