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: }