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;

 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 = 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--;
150:           }
151:         }
152:       }
153:       nd_global = 0;
154:     }
155:     PetscCallMPI(MPIU_Allreduce(&nd, &nd_global, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
156:   }
157:   for (i = 0; i < n; i++) colors[i] = (ISColoringValue)lcolors[i];
158:   PetscCall(PetscFree(mask));
159:   PetscCall(PetscFree(lcolors));
160:   if (mo) PetscCall(PetscFree2(ocolors, owts));
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: static PetscErrorCode GreedyColoringLocalDistanceTwo_Private(MatColoring mc, PetscReal *wts, PetscInt *lperm, ISColoringValue *colors)
165: {
166:   MC_Greedy       *gr = (MC_Greedy *)mc->data;
167:   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;
168:   Mat              m   = mc->mat, mt;
169:   Mat_MPIAIJ      *aij = (Mat_MPIAIJ *)m->data;
170:   Mat              md = NULL, mo = NULL;
171:   const PetscInt  *md_i, *mo_i, *md_j, *mo_j;
172:   const PetscInt  *rmd_i, *rmo_i, *rmd_j, *rmo_j;
173:   PetscBool        isMPIAIJ, isSEQAIJ;
174:   PetscInt         pcol, *dcolors, *ocolors;
175:   ISColoringValue *badidx;
176:   const PetscInt  *cidx;
177:   PetscReal       *owts, *colorweights;
178:   PetscInt        *oconf, *conf;
179:   PetscSF          sf;

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

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

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

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

514: static PetscErrorCode MatColoringSetFromOptions_Greedy(MatColoring mc, PetscOptionItems PetscOptionsObject)
515: {
516:   MC_Greedy *gr = (MC_Greedy *)mc->data;

518:   PetscFunctionBegin;
519:   PetscOptionsHeadBegin(PetscOptionsObject, "Greedy options");
520:   PetscCall(PetscOptionsBool("-mat_coloring_greedy_symmetric", "Flag for assuming a symmetric matrix", "", gr->symmetric, &gr->symmetric, NULL));
521:   PetscOptionsHeadEnd();
522:   PetscFunctionReturn(PETSC_SUCCESS);
523: }

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

528:    Level: beginner

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

540:    Supports both distance one and distance two colorings.

542: .seealso: [](sec_fdmatrix), [](sec_matfactor), `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`
543: M*/
544: PETSC_EXTERN PetscErrorCode MatColoringCreate_Greedy(MatColoring mc)
545: {
546:   MC_Greedy *gr;

548:   PetscFunctionBegin;
549:   PetscCall(PetscNew(&gr));
550:   mc->data                = gr;
551:   mc->ops->apply          = MatColoringApply_Greedy;
552:   mc->ops->view           = NULL;
553:   mc->ops->destroy        = MatColoringDestroy_Greedy;
554:   mc->ops->setfromoptions = MatColoringSetFromOptions_Greedy;

556:   gr->symmetric = PETSC_TRUE;
557:   PetscFunctionReturn(PETSC_SUCCESS);
558: }