Actual source code: greedy.c

  1: #include <petsc/private/matimpl.h>
  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:   PetscBool symmetric;
  8: } MC_Greedy;

 10: static PetscErrorCode MatColoringDestroy_Greedy(MatColoring mc)
 11: {
 12:   PetscFunctionBegin;
 13:   PetscCall(PetscFree(mc->data));
 14:   PetscFunctionReturn(PETSC_SUCCESS);
 15: }

 17: static PetscErrorCode GreedyColoringLocalDistanceOne_Private(MatColoring mc, PetscReal *wts, PetscInt *lperm, ISColoringValue *colors)
 18: {
 19:   PetscInt        i, j, k, s, e, n, no, nd, nd_global, n_global, idx, ncols, maxcolors, masksize, ccol, *mask;
 20:   Mat             m   = mc->mat;
 21:   Mat_MPIAIJ     *aij = (Mat_MPIAIJ *)m->data;
 22:   Mat             md = NULL, mo = NULL;
 23:   const PetscInt *md_i, *mo_i, *md_j, *mo_j;
 24:   PetscBool       isMPIAIJ, isSEQAIJ;
 25:   PetscInt        pcol;
 26:   const PetscInt *cidx;
 27:   PetscInt       *lcolors, *ocolors;
 28:   PetscReal      *owts = NULL;
 29:   PetscSF         sf;
 30:   PetscLayout     layout;

 32:   PetscFunctionBegin;
 33:   PetscCall(MatGetSize(m, &n_global, NULL));
 34:   PetscCall(MatGetOwnershipRange(m, &s, &e));
 35:   n         = e - s;
 36:   masksize  = 20;
 37:   nd_global = 0;
 38:   /* get the matrix communication structures */
 39:   PetscCall(PetscObjectTypeCompare((PetscObject)m, MATMPIAIJ, &isMPIAIJ));
 40:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)m, MATSEQAIJ, &isSEQAIJ));
 41:   PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_WRONG, "Matrix must be AIJ for greedy coloring");
 42:   if (isMPIAIJ) {
 43:     /* get the CSR data for on and off-diagonal portions of m */
 44:     Mat_SeqAIJ *dseq;
 45:     Mat_SeqAIJ *oseq;
 46:     md   = aij->A;
 47:     dseq = (Mat_SeqAIJ *)md->data;
 48:     mo   = aij->B;
 49:     oseq = (Mat_SeqAIJ *)mo->data;
 50:     md_i = dseq->i;
 51:     md_j = dseq->j;
 52:     mo_i = oseq->i;
 53:     mo_j = oseq->j;
 54:   } else {
 55:     /* get the CSR data for m */
 56:     Mat_SeqAIJ *dseq;
 57:     /* no off-processor nodes */
 58:     md   = m;
 59:     dseq = (Mat_SeqAIJ *)md->data;
 60:     mo   = NULL;
 61:     no   = 0;
 62:     md_i = dseq->i;
 63:     md_j = dseq->j;
 64:     mo_i = NULL;
 65:     mo_j = NULL;
 66:   }

 68:   PetscCall(MatColoringGetMaxColors(mc, &maxcolors));
 69:   if (mo) {
 70:     PetscCall(VecGetSize(aij->lvec, &no));
 71:     PetscCall(PetscMalloc2(no, &ocolors, no, &owts));
 72:     for (i = 0; i < no; i++) ocolors[i] = maxcolors;
 73:   }

 75:   PetscCall(PetscMalloc1(masksize, &mask));
 76:   PetscCall(PetscMalloc1(n, &lcolors));
 77:   for (i = 0; i < n; i++) lcolors[i] = maxcolors;
 78:   for (i = 0; i < masksize; i++) mask[i] = -1;
 79:   if (mo) {
 80:     /* transfer neighbor weights */
 81:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)m), &sf));
 82:     PetscCall(MatGetLayouts(m, &layout, NULL));
 83:     PetscCall(PetscSFSetGraphLayout(sf, layout, no, NULL, PETSC_COPY_VALUES, aij->garray));
 84:     PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
 85:     PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
 86:   }
 87:   while (nd_global < n_global) {
 88:     nd = n;
 89:     /* assign lowest possible color to each local vertex */
 90:     PetscCall(PetscLogEventBegin(MATCOLORING_Local, mc, 0, 0, 0));
 91:     for (i = 0; i < n; i++) {
 92:       idx = lperm[i];
 93:       if (lcolors[idx] == maxcolors) {
 94:         ncols = md_i[idx + 1] - md_i[idx];
 95:         cidx  = &(md_j[md_i[idx]]);
 96:         for (j = 0; j < ncols; j++) {
 97:           if (lcolors[cidx[j]] != maxcolors) {
 98:             ccol = lcolors[cidx[j]];
 99:             if (ccol >= masksize) {
100:               PetscInt *newmask;
101:               PetscCall(PetscMalloc1(masksize * 2, &newmask));
102:               for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
103:               for (k = 0; k < masksize; k++) newmask[k] = mask[k];
104:               PetscCall(PetscFree(mask));
105:               mask = newmask;
106:               masksize *= 2;
107:             }
108:             mask[ccol] = idx;
109:           }
110:         }
111:         if (mo) {
112:           ncols = mo_i[idx + 1] - mo_i[idx];
113:           cidx  = &(mo_j[mo_i[idx]]);
114:           for (j = 0; j < ncols; j++) {
115:             if (ocolors[cidx[j]] != maxcolors) {
116:               ccol = ocolors[cidx[j]];
117:               if (ccol >= masksize) {
118:                 PetscInt *newmask;
119:                 PetscCall(PetscMalloc1(masksize * 2, &newmask));
120:                 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
121:                 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
122:                 PetscCall(PetscFree(mask));
123:                 mask = newmask;
124:                 masksize *= 2;
125:               }
126:               mask[ccol] = idx;
127:             }
128:           }
129:         }
130:         for (j = 0; j < masksize; j++) {
131:           if (mask[j] != idx) break;
132:         }
133:         pcol = j;
134:         if (pcol > maxcolors) pcol = maxcolors;
135:         lcolors[idx] = pcol;
136:       }
137:     }
138:     PetscCall(PetscLogEventEnd(MATCOLORING_Local, mc, 0, 0, 0));
139:     if (mo) {
140:       /* transfer neighbor colors */
141:       PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
142:       PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lcolors, ocolors, MPI_REPLACE));
143:       PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lcolors, ocolors, MPI_REPLACE));
144:       /* check for conflicts -- this is merely checking if any adjacent off-processor rows have the same color and marking the ones that are lower weight locally for changing */
145:       for (i = 0; i < n; i++) {
146:         ncols = mo_i[i + 1] - mo_i[i];
147:         cidx  = &(mo_j[mo_i[i]]);
148:         for (j = 0; j < ncols; j++) {
149:           /* in the case of conflicts, the highest weight one stays and the others go */
150:           if ((ocolors[cidx[j]] == lcolors[i]) && (owts[cidx[j]] > wts[i]) && lcolors[i] < maxcolors) {
151:             lcolors[i] = maxcolors;
152:             nd--;
153:           }
154:         }
155:       }
156:       nd_global = 0;
157:     }
158:     PetscCallMPI(MPIU_Allreduce(&nd, &nd_global, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
159:   }
160:   for (i = 0; i < n; i++) colors[i] = (ISColoringValue)lcolors[i];
161:   PetscCall(PetscFree(mask));
162:   PetscCall(PetscFree(lcolors));
163:   if (mo) {
164:     PetscCall(PetscFree2(ocolors, owts));
165:     PetscCall(PetscSFDestroy(&sf));
166:   }
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }

170: static PetscErrorCode GreedyColoringLocalDistanceTwo_Private(MatColoring mc, PetscReal *wts, PetscInt *lperm, ISColoringValue *colors)
171: {
172:   MC_Greedy       *gr = (MC_Greedy *)mc->data;
173:   PetscInt         i, j, k, l, s, e, n, nd, nd_global, n_global, idx, ncols, maxcolors, mcol, mcol_global, nd1cols, *mask, masksize, *d1cols, *bad, *badnext, nbad, badsize, ccol, no, cbad;
174:   Mat              m   = mc->mat, mt;
175:   Mat_MPIAIJ      *aij = (Mat_MPIAIJ *)m->data;
176:   Mat              md = NULL, mo = NULL;
177:   const PetscInt  *md_i, *mo_i, *md_j, *mo_j;
178:   const PetscInt  *rmd_i, *rmo_i, *rmd_j, *rmo_j;
179:   PetscBool        isMPIAIJ, isSEQAIJ;
180:   PetscInt         pcol, *dcolors, *ocolors;
181:   ISColoringValue *badidx;
182:   const PetscInt  *cidx;
183:   PetscReal       *owts, *colorweights;
184:   PetscInt        *oconf, *conf;
185:   PetscSF          sf;
186:   PetscLayout      layout;

188:   PetscFunctionBegin;
189:   PetscCall(MatGetSize(m, &n_global, NULL));
190:   PetscCall(MatGetOwnershipRange(m, &s, &e));
191:   n         = e - s;
192:   nd_global = 0;
193:   /* get the matrix communication structures */
194:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)m, MATMPIAIJ, &isMPIAIJ));
195:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)m, MATSEQAIJ, &isSEQAIJ));
196:   PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_WRONG, "Matrix must be AIJ for greedy coloring");
197:   if (isMPIAIJ) {
198:     Mat_SeqAIJ *dseq;
199:     Mat_SeqAIJ *oseq;
200:     md    = aij->A;
201:     dseq  = (Mat_SeqAIJ *)md->data;
202:     mo    = aij->B;
203:     oseq  = (Mat_SeqAIJ *)mo->data;
204:     md_i  = dseq->i;
205:     md_j  = dseq->j;
206:     mo_i  = oseq->i;
207:     mo_j  = oseq->j;
208:     rmd_i = dseq->i;
209:     rmd_j = dseq->j;
210:     rmo_i = oseq->i;
211:     rmo_j = oseq->j;
212:   } else {
213:     Mat_SeqAIJ *dseq;
214:     /* no off-processor nodes */
215:     md    = m;
216:     dseq  = (Mat_SeqAIJ *)md->data;
217:     md_i  = dseq->i;
218:     md_j  = dseq->j;
219:     mo_i  = NULL;
220:     mo_j  = NULL;
221:     rmd_i = dseq->i;
222:     rmd_j = dseq->j;
223:     rmo_i = NULL;
224:     rmo_j = NULL;
225:   }
226:   if (!gr->symmetric) {
227:     Mat_SeqAIJ *dseq = NULL;

229:     PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)mc), PETSC_ERR_SUP, "Nonsymmetric greedy coloring only works in serial");
230:     PetscCall(MatTranspose(m, MAT_INITIAL_MATRIX, &mt));
231:     dseq  = (Mat_SeqAIJ *)mt->data;
232:     rmd_i = dseq->i;
233:     rmd_j = dseq->j;
234:     rmo_i = NULL;
235:     rmo_j = NULL;
236:   }
237:   /* create the vectors and communication structures if necessary */
238:   no = 0;
239:   if (mo) {
240:     PetscCall(VecGetLocalSize(aij->lvec, &no));
241:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)m), &sf));
242:     PetscCall(MatGetLayouts(m, &layout, NULL));
243:     PetscCall(PetscSFSetGraphLayout(sf, layout, no, NULL, PETSC_COPY_VALUES, aij->garray));
244:   }
245:   PetscCall(MatColoringGetMaxColors(mc, &maxcolors));
246:   masksize = n;
247:   nbad     = 0;
248:   badsize  = n;
249:   PetscCall(PetscMalloc1(masksize, &mask));
250:   PetscCall(PetscMalloc4(n, &d1cols, n, &dcolors, n, &conf, n, &bad));
251:   PetscCall(PetscMalloc2(badsize, &badidx, badsize, &badnext));
252:   for (i = 0; i < masksize; i++) mask[i] = -1;
253:   for (i = 0; i < n; i++) {
254:     dcolors[i] = maxcolors;
255:     bad[i]     = -1;
256:   }
257:   for (i = 0; i < badsize; i++) badnext[i] = -1;
258:   if (mo) {
259:     PetscCall(PetscMalloc3(no, &owts, no, &oconf, no, &ocolors));
260:     PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
261:     PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
262:     for (i = 0; i < no; i++) ocolors[i] = maxcolors;
263:   } else { /* Appease overzealous -Wmaybe-initialized */
264:     owts    = NULL;
265:     oconf   = NULL;
266:     ocolors = NULL;
267:   }
268:   mcol = 0;
269:   while (nd_global < n_global) {
270:     nd = n;
271:     /* assign lowest possible color to each local vertex */
272:     mcol_global = 0;
273:     PetscCall(PetscLogEventBegin(MATCOLORING_Local, mc, 0, 0, 0));
274:     for (i = 0; i < n; i++) {
275:       idx = lperm[i];
276:       if (dcolors[idx] == maxcolors) {
277:         /* entries in bad */
278:         cbad = bad[idx];
279:         while (cbad >= 0) {
280:           ccol = badidx[cbad];
281:           if (ccol >= masksize) {
282:             PetscInt *newmask;
283:             PetscCall(PetscMalloc1(masksize * 2, &newmask));
284:             for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
285:             for (k = 0; k < masksize; k++) newmask[k] = mask[k];
286:             PetscCall(PetscFree(mask));
287:             mask = newmask;
288:             masksize *= 2;
289:           }
290:           mask[ccol] = idx;
291:           cbad       = badnext[cbad];
292:         }
293:         /* diagonal distance-one rows */
294:         nd1cols = 0;
295:         ncols   = rmd_i[idx + 1] - rmd_i[idx];
296:         cidx    = &(rmd_j[rmd_i[idx]]);
297:         for (j = 0; j < ncols; j++) {
298:           d1cols[nd1cols] = cidx[j];
299:           nd1cols++;
300:           ccol = dcolors[cidx[j]];
301:           if (ccol != maxcolors) {
302:             if (ccol >= masksize) {
303:               PetscInt *newmask;
304:               PetscCall(PetscMalloc1(masksize * 2, &newmask));
305:               for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
306:               for (k = 0; k < masksize; k++) newmask[k] = mask[k];
307:               PetscCall(PetscFree(mask));
308:               mask = newmask;
309:               masksize *= 2;
310:             }
311:             mask[ccol] = idx;
312:           }
313:         }
314:         /* off-diagonal distance-one rows */
315:         if (mo) {
316:           ncols = rmo_i[idx + 1] - rmo_i[idx];
317:           cidx  = &(rmo_j[rmo_i[idx]]);
318:           for (j = 0; j < ncols; j++) {
319:             ccol = ocolors[cidx[j]];
320:             if (ccol != maxcolors) {
321:               if (ccol >= masksize) {
322:                 PetscInt *newmask;
323:                 PetscCall(PetscMalloc1(masksize * 2, &newmask));
324:                 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
325:                 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
326:                 PetscCall(PetscFree(mask));
327:                 mask = newmask;
328:                 masksize *= 2;
329:               }
330:               mask[ccol] = idx;
331:             }
332:           }
333:         }
334:         /* diagonal distance-two rows */
335:         for (j = 0; j < nd1cols; j++) {
336:           ncols = md_i[d1cols[j] + 1] - md_i[d1cols[j]];
337:           cidx  = &(md_j[md_i[d1cols[j]]]);
338:           for (l = 0; l < ncols; l++) {
339:             ccol = dcolors[cidx[l]];
340:             if (ccol != maxcolors) {
341:               if (ccol >= masksize) {
342:                 PetscInt *newmask;
343:                 PetscCall(PetscMalloc1(masksize * 2, &newmask));
344:                 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
345:                 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
346:                 PetscCall(PetscFree(mask));
347:                 mask = newmask;
348:                 masksize *= 2;
349:               }
350:               mask[ccol] = idx;
351:             }
352:           }
353:         }
354:         /* off-diagonal distance-two rows */
355:         if (mo) {
356:           for (j = 0; j < nd1cols; j++) {
357:             ncols = mo_i[d1cols[j] + 1] - mo_i[d1cols[j]];
358:             cidx  = &(mo_j[mo_i[d1cols[j]]]);
359:             for (l = 0; l < ncols; l++) {
360:               ccol = ocolors[cidx[l]];
361:               if (ccol != maxcolors) {
362:                 if (ccol >= masksize) {
363:                   PetscInt *newmask;
364:                   PetscCall(PetscMalloc1(masksize * 2, &newmask));
365:                   for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
366:                   for (k = 0; k < masksize; k++) newmask[k] = mask[k];
367:                   PetscCall(PetscFree(mask));
368:                   mask = newmask;
369:                   masksize *= 2;
370:                 }
371:                 mask[ccol] = idx;
372:               }
373:             }
374:           }
375:         }
376:         /* assign this one the lowest color possible by seeing if there's a gap in the sequence of sorted neighbor colors */
377:         for (j = 0; j < masksize; j++) {
378:           if (mask[j] != idx) break;
379:         }
380:         pcol = j;
381:         if (pcol > maxcolors) pcol = maxcolors;
382:         dcolors[idx] = pcol;
383:         if (pcol > mcol) mcol = pcol;
384:       }
385:     }
386:     PetscCall(PetscLogEventEnd(MATCOLORING_Local, mc, 0, 0, 0));
387:     if (mo) {
388:       /* transfer neighbor colors */
389:       PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dcolors, ocolors, MPI_REPLACE));
390:       PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dcolors, ocolors, MPI_REPLACE));
391:       /* find the maximum color assigned locally and allocate a mask */
392:       PetscCallMPI(MPIU_Allreduce(&mcol, &mcol_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
393:       PetscCall(PetscMalloc1(mcol_global + 1, &colorweights));
394:       /* check for conflicts */
395:       for (i = 0; i < n; i++) conf[i] = PETSC_FALSE;
396:       for (i = 0; i < no; i++) oconf[i] = PETSC_FALSE;
397:       for (i = 0; i < n; i++) {
398:         ncols = mo_i[i + 1] - mo_i[i];
399:         cidx  = &(mo_j[mo_i[i]]);
400:         if (ncols > 0) {
401:           /* fill in the mask */
402:           for (j = 0; j < mcol_global + 1; j++) colorweights[j] = 0;
403:           colorweights[dcolors[i]] = wts[i];
404:           /* fill in the off-diagonal part of the mask */
405:           for (j = 0; j < ncols; j++) {
406:             ccol = ocolors[cidx[j]];
407:             if (ccol < maxcolors) {
408:               if (colorweights[ccol] < owts[cidx[j]]) colorweights[ccol] = owts[cidx[j]];
409:             }
410:           }
411:           /* fill in the on-diagonal part of the mask */
412:           ncols = md_i[i + 1] - md_i[i];
413:           cidx  = &(md_j[md_i[i]]);
414:           for (j = 0; j < ncols; j++) {
415:             ccol = dcolors[cidx[j]];
416:             if (ccol < maxcolors) {
417:               if (colorweights[ccol] < wts[cidx[j]]) colorweights[ccol] = wts[cidx[j]];
418:             }
419:           }
420:           /* go back through and set up on and off-diagonal conflict vectors */
421:           ncols = md_i[i + 1] - md_i[i];
422:           cidx  = &(md_j[md_i[i]]);
423:           for (j = 0; j < ncols; j++) {
424:             ccol = dcolors[cidx[j]];
425:             if (ccol < maxcolors) {
426:               if (colorweights[ccol] > wts[cidx[j]]) conf[cidx[j]] = PETSC_TRUE;
427:             }
428:           }
429:           ncols = mo_i[i + 1] - mo_i[i];
430:           cidx  = &(mo_j[mo_i[i]]);
431:           for (j = 0; j < ncols; j++) {
432:             ccol = ocolors[cidx[j]];
433:             if (ccol < maxcolors) {
434:               if (colorweights[ccol] > owts[cidx[j]]) oconf[cidx[j]] = PETSC_TRUE;
435:             }
436:           }
437:         }
438:       }
439:       nd_global = 0;
440:       PetscCall(PetscFree(colorweights));
441:       PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
442:       PetscCall(PetscSFReduceBegin(sf, MPIU_INT, oconf, conf, MPI_SUM));
443:       PetscCall(PetscSFReduceEnd(sf, MPIU_INT, oconf, conf, MPI_SUM));
444:       PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
445:       /* go through and unset local colors that have conflicts */
446:       for (i = 0; i < n; i++) {
447:         if (conf[i] > 0) {
448:           /* push this color onto the bad stack */
449:           PetscCall(ISColoringValueCast(dcolors[i], &badidx[nbad]));
450:           badnext[nbad] = bad[i];
451:           bad[i]        = nbad;
452:           nbad++;
453:           if (nbad >= badsize) {
454:             PetscInt        *newbadnext;
455:             ISColoringValue *newbadidx;
456:             PetscCall(PetscMalloc2(badsize * 2, &newbadidx, badsize * 2, &newbadnext));
457:             for (k = 0; k < 2 * badsize; k++) newbadnext[k] = -1;
458:             for (k = 0; k < badsize; k++) {
459:               newbadidx[k]  = badidx[k];
460:               newbadnext[k] = badnext[k];
461:             }
462:             PetscCall(PetscFree2(badidx, badnext));
463:             badidx  = newbadidx;
464:             badnext = newbadnext;
465:             badsize *= 2;
466:           }
467:           dcolors[i] = maxcolors;
468:           nd--;
469:         }
470:       }
471:     }
472:     PetscCallMPI(MPIU_Allreduce(&nd, &nd_global, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
473:   }
474:   if (mo) {
475:     PetscCall(PetscSFDestroy(&sf));
476:     PetscCall(PetscFree3(owts, oconf, ocolors));
477:   }
478:   for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(dcolors[i], colors + i));
479:   PetscCall(PetscFree(mask));
480:   PetscCall(PetscFree4(d1cols, dcolors, conf, bad));
481:   PetscCall(PetscFree2(badidx, badnext));
482:   if (!gr->symmetric) PetscCall(MatDestroy(&mt));
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

486: static PetscErrorCode MatColoringApply_Greedy(MatColoring mc, ISColoring *iscoloring)
487: {
488:   PetscInt         finalcolor, finalcolor_global;
489:   ISColoringValue *colors;
490:   PetscInt         ncolstotal, ncols;
491:   PetscReal       *wts;
492:   PetscInt         i, *lperm;

494:   PetscFunctionBegin;
495:   PetscCall(MatGetSize(mc->mat, NULL, &ncolstotal));
496:   PetscCall(MatGetLocalSize(mc->mat, NULL, &ncols));
497:   if (!mc->user_weights) {
498:     PetscCall(MatColoringCreateWeights(mc, &wts, &lperm));
499:   } else {
500:     wts   = mc->user_weights;
501:     lperm = mc->user_lperm;
502:   }
503:   PetscCheck(mc->dist == 1 || mc->dist == 2, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_OUTOFRANGE, "Only distance 1 and distance 2 supported by MatColoringGreedy");
504:   PetscCall(PetscMalloc1(ncols, &colors));
505:   if (mc->dist == 1) {
506:     PetscCall(GreedyColoringLocalDistanceOne_Private(mc, wts, lperm, colors));
507:   } else {
508:     PetscCall(GreedyColoringLocalDistanceTwo_Private(mc, wts, lperm, colors));
509:   }
510:   finalcolor = 0;
511:   for (i = 0; i < ncols; i++) {
512:     if (colors[i] > finalcolor) finalcolor = colors[i];
513:   }
514:   finalcolor_global = 0;
515:   PetscCallMPI(MPIU_Allreduce(&finalcolor, &finalcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
516:   PetscCall(PetscLogEventBegin(MATCOLORING_ISCreate, mc, 0, 0, 0));
517:   PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)mc), finalcolor_global + 1, ncols, colors, PETSC_OWN_POINTER, iscoloring));
518:   PetscCall(PetscLogEventEnd(MATCOLORING_ISCreate, mc, 0, 0, 0));
519:   if (!mc->user_weights) {
520:     PetscCall(PetscFree(wts));
521:     PetscCall(PetscFree(lperm));
522:   }
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: static PetscErrorCode MatColoringSetFromOptions_Greedy(MatColoring mc, PetscOptionItems *PetscOptionsObject)
527: {
528:   MC_Greedy *gr = (MC_Greedy *)mc->data;

530:   PetscFunctionBegin;
531:   PetscOptionsHeadBegin(PetscOptionsObject, "Greedy options");
532:   PetscCall(PetscOptionsBool("-mat_coloring_greedy_symmetric", "Flag for assuming a symmetric matrix", "", gr->symmetric, &gr->symmetric, NULL));
533:   PetscOptionsHeadEnd();
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

537: /*MC
538:   MATCOLORINGGREEDY - Greedy-with-conflict correction based matrix coloring for distance 1 and 2 {cite}`bozdaug2005parallel`

540:    Level: beginner

542:    Notes:
543:    These algorithms proceed in two phases -- local coloring and conflict resolution.  The local coloring
544:    tentatively colors all vertices at the distance required given what's known of the global coloring.  Then,
545:    the updated colors are transferred to different processors at distance one.  In the distance one case, each
546:    vertex with nonlocal neighbors is then checked to see if it conforms, with the vertex being
547:    marked for recoloring if its lower weight than its same colored neighbor.  In the distance two case,
548:    each boundary vertex's immediate star is checked for validity of the coloring.  Lower-weight conflict
549:    vertices are marked, and then the conflicts are gathered back on owning processors.  In both cases
550:    this is done until each column has received a valid color.

552:    Supports both distance one and distance two colorings.

554: .seealso: `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MatColoringType`
555: M*/
556: PETSC_EXTERN PetscErrorCode MatColoringCreate_Greedy(MatColoring mc)
557: {
558:   MC_Greedy *gr;

560:   PetscFunctionBegin;
561:   PetscCall(PetscNew(&gr));
562:   mc->data                = gr;
563:   mc->ops->apply          = MatColoringApply_Greedy;
564:   mc->ops->view           = NULL;
565:   mc->ops->destroy        = MatColoringDestroy_Greedy;
566:   mc->ops->setfromoptions = MatColoringSetFromOptions_Greedy;

568:   gr->symmetric = PETSC_TRUE;
569:   PetscFunctionReturn(PETSC_SUCCESS);
570: }