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_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;

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

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

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

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

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

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

471: static PetscErrorCode MatColoringApply_Greedy(MatColoring mc, ISColoring *iscoloring)
472: {
473:   PetscInt         finalcolor_global;
474:   ISColoringValue *colors;
475:   PetscInt         ncolstotal, ncols;
476:   PetscReal       *wts;
477:   PetscInt         i, *lperm;

479:   PetscFunctionBegin;
480:   PetscCall(MatGetSize(mc->mat, NULL, &ncolstotal));
481:   PetscCall(MatGetLocalSize(mc->mat, NULL, &ncols));
482:   if (!mc->user_weights) {
483:     PetscCall(MatColoringCreateWeights(mc, &wts, &lperm));
484:   } else {
485:     wts   = mc->user_weights;
486:     lperm = mc->user_lperm;
487:   }
488:   PetscCheck(mc->dist == 1 || mc->dist == 2, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_OUTOFRANGE, "Only distance 1 and distance 2 supported by MatColoringGreedy");
489:   PetscCall(PetscMalloc1(ncols, &colors));
490:   if (mc->dist == 1) {
491:     PetscCall(GreedyColoringLocalDistanceOne_Private(mc, wts, lperm, colors));
492:   } else {
493:     PetscCall(GreedyColoringLocalDistanceTwo_Private(mc, wts, lperm, colors));
494:   }
495:   finalcolor_global = 0;
496:   for (i = 0; i < ncols; i++) {
497:     if (colors[i] > finalcolor_global) finalcolor_global = colors[i];
498:   }
499:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &finalcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
500:   PetscCall(PetscLogEventBegin(MATCOLORING_ISCreate, mc, 0, 0, 0));
501:   PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)mc), finalcolor_global + 1, ncols, colors, PETSC_OWN_POINTER, iscoloring));
502:   PetscCall(PetscLogEventEnd(MATCOLORING_ISCreate, mc, 0, 0, 0));
503:   if (!mc->user_weights) {
504:     PetscCall(PetscFree(wts));
505:     PetscCall(PetscFree(lperm));
506:   }
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }

510: static PetscErrorCode MatColoringSetFromOptions_Greedy(MatColoring mc, PetscOptionItems PetscOptionsObject)
511: {
512:   MC_Greedy *gr = (MC_Greedy *)mc->data;

514:   PetscFunctionBegin;
515:   PetscOptionsHeadBegin(PetscOptionsObject, "Greedy options");
516:   PetscCall(PetscOptionsBool("-mat_coloring_greedy_symmetric", "Flag for assuming a symmetric matrix", "", gr->symmetric, &gr->symmetric, NULL));
517:   PetscOptionsHeadEnd();
518:   PetscFunctionReturn(PETSC_SUCCESS);
519: }

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

524:    Level: beginner

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

536:    Supports both distance one and distance two colorings.

538: .seealso: [](sec_fdmatrix), [](sec_matfactor), `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`
539: M*/
540: PETSC_EXTERN PetscErrorCode MatColoringCreate_Greedy(MatColoring mc)
541: {
542:   MC_Greedy *gr;

544:   PetscFunctionBegin;
545:   PetscCall(PetscNew(&gr));
546:   mc->data                = gr;
547:   mc->ops->apply          = MatColoringApply_Greedy;
548:   mc->ops->view           = NULL;
549:   mc->ops->destroy        = MatColoringDestroy_Greedy;
550:   mc->ops->setfromoptions = MatColoringSetFromOptions_Greedy;

552:   gr->symmetric = PETSC_TRUE;
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }