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: }