Actual source code: mis.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: #define MIS_NOT_DONE       -2
      7: #define MIS_DELETED        -1
      8: #define MIS_REMOVED        -3
      9: #define MIS_IS_SELECTED(s) (s != MIS_DELETED && s != MIS_NOT_DONE && s != MIS_REMOVED)

     11: /*
     12:    MatCoarsenApply_MIS_private - parallel maximal independent set (MIS) with data locality info. MatAIJ specific!!!

     14:    Input Parameter:
     15:    . perm - serial permutation of rows of local to process in MIS
     16:    . Gmat - global matrix of graph (data not defined)
     17:    . strict_aggs - flag for whether to keep strict (non overlapping) aggregates in 'llist';

     19:    Output Parameter:
     20:    . a_selected - IS of selected vertices, includes 'ghost' nodes at end with natural local indices
     21:    . a_locals_llist - array of list of nodes rooted at selected nodes
     22: */
     23: static PetscErrorCode MatCoarsenApply_MIS_private(IS perm, Mat Gmat, PetscBool strict_aggs, PetscCoarsenData **a_locals_llist)
     24: {
     25:   Mat_SeqAIJ       *matA, *matB = NULL;
     26:   Mat_MPIAIJ       *mpimat = NULL;
     27:   MPI_Comm          comm;
     28:   PetscInt          num_fine_ghosts, kk, n, ix, j, *idx, *ii, Iend, my0, nremoved, gid, lid, cpid, lidj, sgid, t1, t2, slid, nDone, nselected = 0, state, statej;
     29:   PetscInt         *cpcol_gid, *cpcol_state, *lid_cprowID, *lid_gid, *cpcol_sel_gid, *icpcol_gid, *lid_state, *lid_parent_gid = NULL, nrm_tot = 0;
     30:   PetscBool        *lid_removed;
     31:   PetscBool         isMPI, isAIJ, isOK;
     32:   const PetscInt   *perm_ix;
     33:   const PetscInt    nloc = Gmat->rmap->n;
     34:   PetscCoarsenData *agg_lists;
     35:   PetscLayout       layout;
     36:   PetscSF           sf;
     37:   IS                info_is;

     39:   PetscFunctionBegin;
     40:   PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
     41:   PetscCall(ISCreate(comm, &info_is));
     42:   PetscCall(PetscInfo(info_is, "mis: nloc = %d\n", (int)nloc));
     43:   /* get submatrices */
     44:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPI));
     45:   if (isMPI) {
     46:     mpimat = (Mat_MPIAIJ *)Gmat->data;
     47:     matA   = (Mat_SeqAIJ *)mpimat->A->data;
     48:     matB   = (Mat_SeqAIJ *)mpimat->B->data;
     49:     /* force compressed storage of B */
     50:     PetscCall(MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, Gmat->rmap->n, -1.0));
     51:   } else {
     52:     matA = (Mat_SeqAIJ *)Gmat->data;
     53:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isAIJ));
     54:     PetscCheck(isAIJ, comm, PETSC_ERR_PLIB, "Require AIJ matrix.");
     55:   }
     56:   PetscCall(MatGetOwnershipRange(Gmat, &my0, &Iend));
     57:   PetscCall(PetscMalloc4(nloc, &lid_gid, nloc, &lid_cprowID, nloc, &lid_removed, nloc, &lid_state));
     58:   if (strict_aggs) PetscCall(PetscMalloc1(nloc, &lid_parent_gid));
     59:   if (isMPI) {
     60:     for (kk = 0, gid = my0; kk < nloc; kk++, gid++) lid_gid[kk] = gid;
     61:     PetscCall(VecGetLocalSize(mpimat->lvec, &num_fine_ghosts));
     62:     PetscCall(PetscMalloc2(num_fine_ghosts, &cpcol_gid, num_fine_ghosts, &cpcol_state));
     63:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)Gmat), &sf));
     64:     PetscCall(MatGetLayouts(Gmat, &layout, NULL));
     65:     PetscCall(PetscSFSetGraphLayout(sf, layout, num_fine_ghosts, NULL, PETSC_COPY_VALUES, mpimat->garray));
     66:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE));
     67:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE));
     68:     for (kk = 0; kk < num_fine_ghosts; kk++) cpcol_state[kk] = MIS_NOT_DONE;
     69:   } else num_fine_ghosts = 0;

     71:   /* has ghost nodes for !strict and uses local indexing (yuck) */
     72:   PetscCall(PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts + nloc, &agg_lists));
     73:   if (a_locals_llist) *a_locals_llist = agg_lists;

     75:   /* need an inverse map - locals */
     76:   for (kk = 0; kk < nloc; kk++) {
     77:     lid_cprowID[kk] = -1;
     78:     lid_removed[kk] = PETSC_FALSE;
     79:     if (strict_aggs) lid_parent_gid[kk] = -1.0;
     80:     lid_state[kk] = MIS_NOT_DONE;
     81:   }
     82:   /* set index into cmpressed row 'lid_cprowID' */
     83:   if (matB) {
     84:     for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
     85:       lid = matB->compressedrow.rindex[ix];
     86:       if (lid >= 0) lid_cprowID[lid] = ix;
     87:     }
     88:   }
     89:   /* MIS */
     90:   nremoved = nDone = 0;

     92:   PetscCall(ISGetIndices(perm, &perm_ix));
     93:   while (nDone < nloc || PETSC_TRUE) { /* asynchronous not implemented */
     94:     /* check all vertices */
     95:     for (kk = 0; kk < nloc; kk++) {
     96:       lid   = perm_ix[kk];
     97:       state = lid_state[lid];
     98:       if (lid_removed[lid]) continue;
     99:       if (state == MIS_NOT_DONE) {
    100:         /* parallel test, delete if selected ghost */
    101:         isOK = PETSC_TRUE;
    102:         if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
    103:           ii  = matB->compressedrow.i;
    104:           n   = ii[ix + 1] - ii[ix];
    105:           idx = matB->j + ii[ix];
    106:           for (j = 0; j < n; j++) {
    107:             cpid   = idx[j]; /* compressed row ID in B mat */
    108:             gid    = cpcol_gid[cpid];
    109:             statej = cpcol_state[cpid];
    110:             PetscCheck(!MIS_IS_SELECTED(statej), PETSC_COMM_SELF, PETSC_ERR_SUP, "selected ghost: %" PetscInt_FMT, gid);
    111:             if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
    112:               isOK = PETSC_FALSE;                        /* can not delete */
    113:               break;
    114:             }
    115:           }
    116:         } /* parallel test */
    117:         if (isOK) { /* select or remove this vertex */
    118:           nDone++;
    119:           /* check for singleton */
    120:           ii = matA->i;
    121:           n  = ii[lid + 1] - ii[lid];
    122:           if (n < 2) {
    123:             /* if I have any ghost adj then not a sing */
    124:             ix = lid_cprowID[lid];
    125:             if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) {
    126:               nremoved++;
    127:               nrm_tot++;
    128:               lid_removed[lid] = PETSC_TRUE;
    129:               continue;
    130:               // lid_state[lidj] = MIS_REMOVED; /* add singleton to MIS (can cause low rank with elasticity on fine grid) */
    131:             }
    132:           }
    133:           /* SELECTED state encoded with global index */
    134:           lid_state[lid] = lid + my0;
    135:           nselected++;
    136:           if (strict_aggs) {
    137:             PetscCall(PetscCDAppendID(agg_lists, lid, lid + my0));
    138:           } else {
    139:             PetscCall(PetscCDAppendID(agg_lists, lid, lid));
    140:           }
    141:           /* delete local adj */
    142:           idx = matA->j + ii[lid];
    143:           for (j = 0; j < n; j++) {
    144:             lidj   = idx[j];
    145:             statej = lid_state[lidj];
    146:             if (statej == MIS_NOT_DONE) {
    147:               nDone++;
    148:               if (strict_aggs) {
    149:                 PetscCall(PetscCDAppendID(agg_lists, lid, lidj + my0));
    150:               } else {
    151:                 PetscCall(PetscCDAppendID(agg_lists, lid, lidj));
    152:               }
    153:               lid_state[lidj] = MIS_DELETED; /* delete this */
    154:             }
    155:           }
    156:           /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
    157:           if (!strict_aggs) {
    158:             if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
    159:               ii  = matB->compressedrow.i;
    160:               n   = ii[ix + 1] - ii[ix];
    161:               idx = matB->j + ii[ix];
    162:               for (j = 0; j < n; j++) {
    163:                 cpid   = idx[j]; /* compressed row ID in B mat */
    164:                 statej = cpcol_state[cpid];
    165:                 if (statej == MIS_NOT_DONE) PetscCall(PetscCDAppendID(agg_lists, lid, nloc + cpid));
    166:               }
    167:             }
    168:           }
    169:         } /* selected */
    170:       } /* not done vertex */
    171:     } /* vertex loop */

    173:     /* update ghost states and count todos */
    174:     if (isMPI) {
    175:       /* scatter states, check for done */
    176:       PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE));
    177:       PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE));
    178:       ii = matB->compressedrow.i;
    179:       for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
    180:         lid   = matB->compressedrow.rindex[ix]; /* local boundary node */
    181:         state = lid_state[lid];
    182:         if (state == MIS_NOT_DONE) {
    183:           /* look at ghosts */
    184:           n   = ii[ix + 1] - ii[ix];
    185:           idx = matB->j + ii[ix];
    186:           for (j = 0; j < n; j++) {
    187:             cpid   = idx[j]; /* compressed row ID in B mat */
    188:             statej = cpcol_state[cpid];
    189:             if (MIS_IS_SELECTED(statej)) { /* lid is now deleted, do it */
    190:               nDone++;
    191:               lid_state[lid] = MIS_DELETED; /* delete this */
    192:               if (!strict_aggs) {
    193:                 lidj = nloc + cpid;
    194:                 PetscCall(PetscCDAppendID(agg_lists, lidj, lid));
    195:               } else {
    196:                 sgid                = cpcol_gid[cpid];
    197:                 lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */
    198:               }
    199:               break;
    200:             }
    201:           }
    202:         }
    203:       }
    204:       /* all done? */
    205:       t1 = nloc - nDone;
    206:       PetscCallMPI(MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm)); /* synchronous version */
    207:       if (!t2) break;
    208:     } else break; /* all done */
    209:   } /* outer parallel MIS loop */
    210:   PetscCall(ISRestoreIndices(perm, &perm_ix));
    211:   PetscCall(PetscInfo(info_is, "\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices.  %" PetscInt_FMT " selected.\n", nremoved, nloc, nselected));

    213:   /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
    214:   if (strict_aggs && matB) {
    215:     /* need to copy this to free buffer -- should do this globally */
    216:     PetscCall(PetscMalloc2(num_fine_ghosts, &cpcol_sel_gid, num_fine_ghosts, &icpcol_gid));
    217:     for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];

    219:     /* get proc of deleted ghost */
    220:     PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE));
    221:     PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE));
    222:     for (cpid = 0; cpid < num_fine_ghosts; cpid++) {
    223:       sgid = cpcol_sel_gid[cpid];
    224:       gid  = icpcol_gid[cpid];
    225:       if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
    226:         slid = sgid - my0;
    227:         PetscCall(PetscCDAppendID(agg_lists, slid, gid));
    228:       }
    229:     }
    230:     PetscCall(PetscFree2(cpcol_sel_gid, icpcol_gid));
    231:   }
    232:   if (isMPI) {
    233:     PetscCall(PetscSFDestroy(&sf));
    234:     PetscCall(PetscFree2(cpcol_gid, cpcol_state));
    235:   }
    236:   PetscCall(PetscFree4(lid_gid, lid_cprowID, lid_removed, lid_state));
    237:   if (strict_aggs) {
    238:     // check sizes -- all vertices must get in graph
    239:     PetscInt aa[2] = {0, nrm_tot}, bb[2], MM;

    241:     PetscCall(PetscFree(lid_parent_gid));
    242:     PetscCall(MatGetSize(Gmat, &MM, NULL));
    243:     // check sizes -- all vertices must get in graph
    244:     PetscCall(PetscCDCount(agg_lists, &aa[0]));
    245:     PetscCallMPI(MPIU_Allreduce(aa, bb, 2, MPIU_INT, MPI_SUM, comm));
    246:     if (MM != bb[0]) PetscCall(PetscInfo(info_is, "Warning: N = %" PetscInt_FMT ", sum of aggregates %" PetscInt_FMT ", %" PetscInt_FMT " removed total\n", MM, bb[0], bb[1]));
    247:     PetscCheck(MM >= bb[0], comm, PETSC_ERR_PLIB, "Sum of aggs is too large");
    248:   }
    249:   PetscCall(ISDestroy(&info_is));
    250:   PetscFunctionReturn(PETSC_SUCCESS);
    251: }

    253: /*
    254:    MIS coarsen, simple greedy.
    255: */
    256: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
    257: {
    258:   Mat mat = coarse->graph;

    260:   PetscFunctionBegin;
    261:   if (!coarse->perm) {
    262:     IS       perm;
    263:     PetscInt n, m;
    264:     MPI_Comm comm;

    266:     PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
    267:     PetscCall(MatGetLocalSize(mat, &m, &n));
    268:     PetscCall(ISCreateStride(comm, m, 0, 1, &perm));
    269:     PetscCall(MatCoarsenApply_MIS_private(perm, mat, coarse->strict_aggs, &coarse->agg_lists));
    270:     PetscCall(ISDestroy(&perm));
    271:   } else {
    272:     PetscCall(MatCoarsenApply_MIS_private(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists));
    273:   }
    274:   PetscFunctionReturn(PETSC_SUCCESS);
    275: }

    277: static PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse, PetscViewer viewer)
    278: {
    279:   PetscMPIInt       rank;
    280:   PetscBool         isascii;
    281:   PetscViewerFormat format;

    283:   PetscFunctionBegin;
    284:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
    285:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    286:   PetscCall(PetscViewerGetFormat(viewer, &format));
    287:   if (isascii && format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    288:     if (coarse->agg_lists) {
    289:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
    290:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] MIS aggregator\n", rank));
    291:       if (!rank) {
    292:         PetscCDIntNd *pos, *pos2;
    293:         for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) {
    294:           PetscCall(PetscCDGetHeadPos(coarse->agg_lists, kk, &pos));
    295:           if ((pos2 = pos)) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "selected %" PetscInt_FMT ": ", kk));
    296:           while (pos) {
    297:             PetscInt gid1;
    298:             PetscCall(PetscCDIntNdGetID(pos, &gid1));
    299:             PetscCall(PetscCDGetNextPos(coarse->agg_lists, kk, &pos));
    300:             PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " ", gid1));
    301:           }
    302:           if (pos2) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
    303:         }
    304:       }
    305:       PetscCall(PetscViewerFlush(viewer));
    306:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
    307:     } else {
    308:       PetscCall(PetscViewerASCIIPrintf(viewer, "  MIS aggregator lists are not available\n"));
    309:     }
    310:   }
    311:   PetscFunctionReturn(PETSC_SUCCESS);
    312: }

    314: /*MC
    315:    MATCOARSENMIS - Creates a coarsening object that uses a maximal independent set (MIS) algorithm

    317:    Collective

    319:    Input Parameter:
    320: .  coarse - the coarsen context

    322:    Level: beginner

    324: .seealso: `MatCoarsen`, `MatCoarsenApply()`, `MatCoarsenGetData()`, `MatCoarsenSetType()`, `MatCoarsenType`
    325: M*/
    326: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
    327: {
    328:   PetscFunctionBegin;
    329:   coarse->ops->apply = MatCoarsenApply_MIS;
    330:   coarse->ops->view  = MatCoarsenView_MIS;
    331:   PetscFunctionReturn(PETSC_SUCCESS);
    332: }