Actual source code: jp.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  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:
500:     This method uses a parallel Luby-style coloring with weights to choose an independent set of processor
501:    boundary vertices at each stage that may be assigned colors independently.

503:    Supports both distance one and distance two colorings.

505:    References:
506: .  1. - M. Jones and P. Plassmann, "A parallel graph coloring heuristic," SIAM Journal on Scientific Computing, vol. 14, no. 3,
507:    pp. 654-669, 1993.

509: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType()
510: M*/
511: PETSC_EXTERN PetscErrorCode MatColoringCreate_JP(MatColoring mc)
512: {
513:   MC_JP          *jp;

517:   PetscNewLog(mc,&jp);
518:   jp->sf                  = NULL;
519:   jp->dmask               = NULL;
520:   jp->omask               = NULL;
521:   jp->cmask               = NULL;
522:   jp->dwts                = NULL;
523:   jp->owts                = NULL;
524:   jp->local               = PETSC_TRUE;
525:   mc->data                = jp;
526:   mc->ops->apply          = MatColoringApply_JP;
527:   mc->ops->view           = NULL;
528:   mc->ops->destroy        = MatColoringDestroy_JP;
529:   mc->ops->setfromoptions = MatColoringSetFromOptions_JP;
530:   return(0);
531: }