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(MatColoring mc, PetscOptionItems *PetscOptionsObject)
 19: {
 20:   MC_JP *jp = (MC_JP *)mc->data;

 22:   PetscOptionsHeadBegin(PetscOptionsObject, "JP options");
 23:   PetscOptionsBool("-mat_coloring_jp_local", "Do an initial coloring of local columns", "", jp->local, &jp->local, NULL);
 24:   PetscOptionsHeadEnd();
 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++) dwts[i] = maxweights[i];
116:       if (oG) {
117:         PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0);
118:         PetscSFBcastBegin(sf, MPIU_REAL, dwts, owts, MPI_REPLACE);
119:         PetscSFBcastEnd(sf, MPIU_REAL, dwts, owts, MPI_REPLACE);
120:         PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0);
121:       }
122:     }
123:   }
124:   return 0;
125: }

127: static PetscErrorCode MCJPInitialLocalColor_Private(MatColoring mc, PetscInt *lperm, ISColoringValue *colors)
128: {
129:   PetscInt        j, i, s, e, n, bidx, cidx, idx, dist, distance = mc->dist;
130:   Mat             G = mc->mat, dG, oG;
131:   PetscInt       *seen;
132:   PetscInt       *idxbuf;
133:   PetscBool      *boundary;
134:   PetscInt       *distbuf;
135:   PetscInt       *colormask;
136:   PetscInt        ncols;
137:   const PetscInt *cols;
138:   PetscBool       isSeq, isMPI;
139:   Mat_MPIAIJ     *aij;
140:   Mat_SeqAIJ     *daij, *oaij;
141:   PetscInt       *di, *dj, dn;
142:   PetscInt       *oi;

144:   PetscLogEventBegin(MATCOLORING_Local, mc, 0, 0, 0);
145:   MatGetOwnershipRange(G, &s, &e);
146:   n = e - s;
147:   PetscObjectBaseTypeCompare((PetscObject)G, MATSEQAIJ, &isSeq);
148:   PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPI);

151:   /* get the inner matrix structure */
152:   oG = NULL;
153:   oi = NULL;
154:   if (isMPI) {
155:     aij  = (Mat_MPIAIJ *)G->data;
156:     dG   = aij->A;
157:     oG   = aij->B;
158:     daij = (Mat_SeqAIJ *)dG->data;
159:     oaij = (Mat_SeqAIJ *)oG->data;
160:     di   = daij->i;
161:     dj   = daij->j;
162:     oi   = oaij->i;
163:     MatGetSize(oG, &dn, NULL);
164:   } else {
165:     dG = G;
166:     MatGetSize(dG, NULL, &dn);
167:     daij = (Mat_SeqAIJ *)dG->data;
168:     di   = daij->i;
169:     dj   = daij->j;
170:   }
171:   PetscMalloc5(n, &colormask, n, &seen, n, &idxbuf, n, &distbuf, n, &boundary);
172:   for (i = 0; i < dn; i++) {
173:     seen[i]      = -1;
174:     colormask[i] = -1;
175:     boundary[i]  = PETSC_FALSE;
176:   }
177:   /* pass one -- figure out which ones are off-boundary in the distance-n sense */
178:   if (oG) {
179:     for (i = 0; i < dn; i++) {
180:       bidx = -1;
181:       /* nonempty off-diagonal, so this one is on the boundary */
182:       if (oi[i] != oi[i + 1]) {
183:         boundary[i] = PETSC_TRUE;
184:         continue;
185:       }
186:       ncols = di[i + 1] - di[i];
187:       cols  = &(dj[di[i]]);
188:       for (j = 0; j < ncols; j++) {
189:         bidx++;
190:         seen[cols[j]] = i;
191:         distbuf[bidx] = 1;
192:         idxbuf[bidx]  = cols[j];
193:       }
194:       while (bidx >= 0) {
195:         idx  = idxbuf[bidx];
196:         dist = distbuf[bidx];
197:         bidx--;
198:         if (dist < distance) {
199:           if (oi[idx + 1] != oi[idx]) {
200:             boundary[i] = PETSC_TRUE;
201:             break;
202:           }
203:           ncols = di[idx + 1] - di[idx];
204:           cols  = &(dj[di[idx]]);
205:           for (j = 0; j < ncols; j++) {
206:             if (seen[cols[j]] != i) {
207:               bidx++;
208:               seen[cols[j]] = i;
209:               idxbuf[bidx]  = cols[j];
210:               distbuf[bidx] = dist + 1;
211:             }
212:           }
213:         }
214:       }
215:     }
216:     for (i = 0; i < dn; i++) seen[i] = -1;
217:   }
218:   /* pass two -- color it by looking at nearby vertices and building a mask */
219:   for (i = 0; i < dn; i++) {
220:     cidx = lperm[i];
221:     if (!boundary[cidx]) {
222:       bidx  = -1;
223:       ncols = di[cidx + 1] - di[cidx];
224:       cols  = &(dj[di[cidx]]);
225:       for (j = 0; j < ncols; j++) {
226:         bidx++;
227:         seen[cols[j]] = cidx;
228:         distbuf[bidx] = 1;
229:         idxbuf[bidx]  = cols[j];
230:       }
231:       while (bidx >= 0) {
232:         idx  = idxbuf[bidx];
233:         dist = distbuf[bidx];
234:         bidx--;
235:         /* mask this color */
236:         if (colors[idx] < IS_COLORING_MAX) colormask[colors[idx]] = cidx;
237:         if (dist < distance) {
238:           ncols = di[idx + 1] - di[idx];
239:           cols  = &(dj[di[idx]]);
240:           for (j = 0; j < ncols; j++) {
241:             if (seen[cols[j]] != cidx) {
242:               bidx++;
243:               seen[cols[j]] = cidx;
244:               idxbuf[bidx]  = cols[j];
245:               distbuf[bidx] = dist + 1;
246:             }
247:           }
248:         }
249:       }
250:       /* find the lowest untaken color */
251:       for (j = 0; j < n; j++) {
252:         if (colormask[j] != cidx || j >= mc->maxcolors) {
253:           colors[cidx] = j;
254:           break;
255:         }
256:       }
257:     }
258:   }
259:   PetscFree5(colormask, seen, idxbuf, distbuf, boundary);
260:   PetscLogEventEnd(MATCOLORING_Local, mc, 0, 0, 0);
261:   return 0;
262: }

264: static PetscErrorCode MCJPMinColor_Private(MatColoring mc, ISColoringValue maxcolor, const ISColoringValue *colors, ISColoringValue *mincolors)
265: {
266:   MC_JP          *jp = (MC_JP *)mc->data;
267:   Mat             G  = mc->mat, dG, oG;
268:   PetscBool       isSeq, isMPI;
269:   Mat_MPIAIJ     *aij;
270:   Mat_SeqAIJ     *daij, *oaij;
271:   PetscInt       *di, *oi, *dj, *oj;
272:   PetscSF         sf = jp->sf;
273:   PetscLayout     layout;
274:   PetscInt        maskrounds, maskbase, maskradix;
275:   PetscInt        dn, on;
276:   PetscInt        i, j, l, k;
277:   PetscInt       *dmask = jp->dmask, *omask = jp->omask, *cmask = jp->cmask, curmask;
278:   PetscInt        ncols;
279:   const PetscInt *cols;

281:   maskradix  = sizeof(PetscInt) * 8;
282:   maskrounds = 1 + maxcolor / (maskradix);
283:   maskbase   = 0;
284:   PetscObjectBaseTypeCompare((PetscObject)G, MATSEQAIJ, &isSeq);
285:   PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPI);

288:   /* get the inner matrix structure */
289:   oG = NULL;
290:   oi = NULL;
291:   oj = NULL;
292:   if (isMPI) {
293:     aij  = (Mat_MPIAIJ *)G->data;
294:     dG   = aij->A;
295:     oG   = aij->B;
296:     daij = (Mat_SeqAIJ *)dG->data;
297:     oaij = (Mat_SeqAIJ *)oG->data;
298:     di   = daij->i;
299:     dj   = daij->j;
300:     oi   = oaij->i;
301:     oj   = oaij->j;
302:     MatGetSize(oG, &dn, &on);
303:     if (!sf) {
304:       PetscSFCreate(PetscObjectComm((PetscObject)mc), &sf);
305:       MatGetLayouts(G, &layout, NULL);
306:       PetscSFSetGraphLayout(sf, layout, on, NULL, PETSC_COPY_VALUES, aij->garray);
307:       jp->sf = sf;
308:     }
309:   } else {
310:     dG = G;
311:     MatGetSize(dG, NULL, &dn);
312:     daij = (Mat_SeqAIJ *)dG->data;
313:     di   = daij->i;
314:     dj   = daij->j;
315:   }
316:   for (i = 0; i < dn; i++) mincolors[i] = IS_COLORING_MAX;
317:   /* set up the distance-zero mask */
318:   if (!dmask) {
319:     PetscMalloc1(dn, &dmask);
320:     PetscMalloc1(dn, &cmask);
321:     jp->cmask = cmask;
322:     jp->dmask = dmask;
323:     if (oG) {
324:       PetscMalloc1(on, &omask);
325:       jp->omask = omask;
326:     }
327:   }
328:   /* the number of colors may be more than the number of bits in a PetscInt; take multiple rounds */
329:   for (k = 0; k < maskrounds; k++) {
330:     for (i = 0; i < dn; i++) {
331:       cmask[i] = 0;
332:       if (colors[i] < maskbase + maskradix && colors[i] >= maskbase) cmask[i] = 1 << (colors[i] - maskbase);
333:       dmask[i] = cmask[i];
334:     }
335:     if (oG) {
336:       PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0);
337:       PetscSFBcastBegin(sf, MPIU_INT, dmask, omask, MPI_REPLACE);
338:       PetscSFBcastEnd(sf, MPIU_INT, dmask, omask, MPI_REPLACE);
339:       PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0);
340:     }
341:     /* fill in the mask out to the distance of the coloring */
342:     for (l = 0; l < mc->dist; l++) {
343:       /* fill in the on-and-off diagonal mask */
344:       for (i = 0; i < dn; i++) {
345:         ncols = di[i + 1] - di[i];
346:         cols  = &(dj[di[i]]);
347:         for (j = 0; j < ncols; j++) cmask[i] = cmask[i] | dmask[cols[j]];
348:         if (oG) {
349:           ncols = oi[i + 1] - oi[i];
350:           cols  = &(oj[oi[i]]);
351:           for (j = 0; j < ncols; j++) cmask[i] = cmask[i] | omask[cols[j]];
352:         }
353:       }
354:       for (i = 0; i < dn; i++) dmask[i] = cmask[i];
355:       if (l < mc->dist - 1) {
356:         if (oG) {
357:           PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0);
358:           PetscSFBcastBegin(sf, MPIU_INT, dmask, omask, MPI_REPLACE);
359:           PetscSFBcastEnd(sf, MPIU_INT, dmask, omask, MPI_REPLACE);
360:           PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0);
361:         }
362:       }
363:     }
364:     /* read through the mask to see if we've discovered an acceptable color for any vertices in this round */
365:     for (i = 0; i < dn; i++) {
366:       if (mincolors[i] == IS_COLORING_MAX) {
367:         curmask = dmask[i];
368:         for (j = 0; j < maskradix; j++) {
369:           if (curmask % 2 == 0) {
370:             mincolors[i] = j + maskbase;
371:             break;
372:           }
373:           curmask = curmask >> 1;
374:         }
375:       }
376:     }
377:     /* do the next maskradix colors */
378:     maskbase += maskradix;
379:   }
380:   for (i = 0; i < dn; i++) {
381:     if (mincolors[i] == IS_COLORING_MAX) mincolors[i] = maxcolor + 1;
382:   }
383:   return 0;
384: }

386: static PetscErrorCode MatColoringApply_JP(MatColoring mc, ISColoring *iscoloring)
387: {
388:   MC_JP           *jp = (MC_JP *)mc->data;
389:   PetscInt         i, nadded, nadded_total, nadded_total_old, ntotal, n;
390:   PetscInt         maxcolor_local = 0, maxcolor_global = 0, *lperm;
391:   PetscMPIInt      rank;
392:   PetscReal       *weights, *maxweights;
393:   ISColoringValue *color, *mincolor;

395:   MPI_Comm_rank(PetscObjectComm((PetscObject)mc), &rank);
396:   PetscLogEventBegin(MATCOLORING_Weights, mc, 0, 0, 0);
397:   MatColoringCreateWeights(mc, &weights, &lperm);
398:   PetscLogEventEnd(MATCOLORING_Weights, mc, 0, 0, 0);
399:   MatGetSize(mc->mat, NULL, &ntotal);
400:   MatGetLocalSize(mc->mat, NULL, &n);
401:   PetscMalloc1(n, &maxweights);
402:   PetscMalloc1(n, &color);
403:   PetscMalloc1(n, &mincolor);
404:   for (i = 0; i < n; i++) {
405:     color[i]    = IS_COLORING_MAX;
406:     mincolor[i] = 0;
407:   }
408:   nadded           = 0;
409:   nadded_total     = 0;
410:   nadded_total_old = 0;
411:   /* compute purely local vertices */
412:   if (jp->local) {
413:     MCJPInitialLocalColor_Private(mc, lperm, color);
414:     for (i = 0; i < n; i++) {
415:       if (color[i] < IS_COLORING_MAX) {
416:         nadded++;
417:         weights[i] = -1;
418:         if (color[i] > maxcolor_local) maxcolor_local = color[i];
419:       }
420:     }
421:     MPIU_Allreduce(&nadded, &nadded_total, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc));
422:     MPIU_Allreduce(&maxcolor_local, &maxcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc));
423:   }
424:   while (nadded_total < ntotal) {
425:     MCJPMinColor_Private(mc, (ISColoringValue)maxcolor_global, color, mincolor);
426:     MCJPGreatestWeight_Private(mc, weights, maxweights);
427:     for (i = 0; i < n; i++) {
428:       /* choose locally maximal vertices; weights less than zero are omitted from the graph */
429:       if (weights[i] >= maxweights[i] && weights[i] >= 0.) {
430:         /* assign the minimum possible color */
431:         if (mc->maxcolors > mincolor[i]) {
432:           color[i] = mincolor[i];
433:         } else {
434:           color[i] = mc->maxcolors;
435:         }
436:         if (color[i] > maxcolor_local) maxcolor_local = color[i];
437:         weights[i] = -1.;
438:         nadded++;
439:       }
440:     }
441:     MPIU_Allreduce(&maxcolor_local, &maxcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc));
442:     MPIU_Allreduce(&nadded, &nadded_total, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc));
444:     nadded_total_old = nadded_total;
445:   }
446:   PetscLogEventBegin(MATCOLORING_ISCreate, mc, 0, 0, 0);
447:   ISColoringCreate(PetscObjectComm((PetscObject)mc), maxcolor_global + 1, n, color, PETSC_OWN_POINTER, iscoloring);
448:   PetscLogEventEnd(MATCOLORING_ISCreate, mc, 0, 0, 0);
449:   PetscFree(jp->dwts);
450:   PetscFree(jp->dmask);
451:   PetscFree(jp->cmask);
452:   PetscFree(jp->owts);
453:   PetscFree(jp->omask);
454:   PetscFree(weights);
455:   PetscFree(lperm);
456:   PetscFree(maxweights);
457:   PetscFree(mincolor);
458:   PetscSFDestroy(&jp->sf);
459:   return 0;
460: }

462: /*MC
463:   MATCOLORINGJP - Parallel Jones-Plassmann coloring

465:    Level: beginner

467:    Options Database Key:
468: .  -mat_coloring_jp_local - perform a local coloring before applying the parallel algorithm

470:    Notes:
471:     This method uses a parallel Luby-style coloring with weights to choose an independent set of processor
472:    boundary vertices at each stage that may be assigned colors independently.

474:    Supports both distance one and distance two colorings.

476:    References:
477: .  * - M. Jones and P. Plassmann, "A parallel graph coloring heuristic," SIAM Journal on Scientific Computing, vol. 14, no. 3,
478:    pp. 654-669, 1993.

480: .seealso: `MatColoring`, `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`
481: M*/
482: PETSC_EXTERN PetscErrorCode MatColoringCreate_JP(MatColoring mc)
483: {
484:   MC_JP *jp;

486:   PetscNew(&jp);
487:   jp->sf                  = NULL;
488:   jp->dmask               = NULL;
489:   jp->omask               = NULL;
490:   jp->cmask               = NULL;
491:   jp->dwts                = NULL;
492:   jp->owts                = NULL;
493:   jp->local               = PETSC_TRUE;
494:   mc->data                = jp;
495:   mc->ops->apply          = MatColoringApply_JP;
496:   mc->ops->view           = NULL;
497:   mc->ops->destroy        = MatColoringDestroy_JP;
498:   mc->ops->setfromoptions = MatColoringSetFromOptions_JP;
499:   return 0;
500: }