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: /*
     13:    maxIndSetAgg - parallel maximal independent set (MIS) with data locality info. MatAIJ specific!!!

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

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

     41:   PetscObjectGetComm((PetscObject)Gmat,&comm);

     43:   /* get submatrices */
     44:   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:     MatCheckCompressedRow(mpimat->B,matB->nonzerorowcnt,&matB->compressedrow,matB->i,Gmat->rmap->n,-1.0);
     51:   } else {
     52:     PetscObjectBaseTypeCompare((PetscObject)Gmat,MATSEQAIJ,&isAIJ);
     53:     if (!isAIJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Require AIJ matrix.");
     54:     matA = (Mat_SeqAIJ*)Gmat->data;
     55:   }
     56:   MatGetOwnershipRange(Gmat,&my0,&Iend);
     57:   PetscMalloc1(nloc,&lid_gid); /* explicit array needed */
     58:   if (mpimat) {
     59:     for (kk=0,gid=my0; kk<nloc; kk++,gid++) {
     60:       lid_gid[kk] = gid;
     61:     }
     62:     VecGetLocalSize(mpimat->lvec, &num_fine_ghosts);
     63:     PetscMalloc1(num_fine_ghosts,&cpcol_gid);
     64:     PetscMalloc1(num_fine_ghosts,&cpcol_state);
     65:     PetscSFCreate(PetscObjectComm((PetscObject)Gmat),&sf);
     66:     MatGetLayouts(Gmat,&layout,NULL);
     67:     PetscSFSetGraphLayout(sf,layout,num_fine_ghosts,NULL,PETSC_COPY_VALUES,mpimat->garray);
     68:     PetscSFBcastBegin(sf,MPIU_INT,lid_gid,cpcol_gid,MPI_REPLACE);
     69:     PetscSFBcastEnd(sf,MPIU_INT,lid_gid,cpcol_gid,MPI_REPLACE);
     70:     for (kk=0;kk<num_fine_ghosts;kk++) {
     71:       cpcol_state[kk]=MIS_NOT_DONE;
     72:     }
     73:   } else num_fine_ghosts = 0;

     75:   PetscMalloc1(nloc, &lid_cprowID);
     76:   PetscMalloc1(nloc, &lid_removed); /* explicit array needed */
     77:   if (strict_aggs) {
     78:     PetscMalloc1(nloc,&lid_parent_gid);
     79:   }
     80:   PetscMalloc1(nloc,&lid_state);

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

     86:   /* need an inverse map - locals */
     87:   for (kk=0; kk<nloc; kk++) {
     88:     lid_cprowID[kk] = -1; lid_removed[kk] = PETSC_FALSE;
     89:     if (strict_aggs) {
     90:       lid_parent_gid[kk] = -1.0;
     91:     }
     92:     lid_state[kk] = MIS_NOT_DONE;
     93:   }
     94:   /* set index into cmpressed row 'lid_cprowID' */
     95:   if (matB) {
     96:     for (ix=0; ix<matB->compressedrow.nrows; ix++) {
     97:       lid = matB->compressedrow.rindex[ix];
     98:       lid_cprowID[lid] = ix;
     99:     }
    100:   }
    101:   /* MIS */
    102:   iter = nremoved = nDone = 0;
    103:   ISGetIndices(perm, &perm_ix);
    104:   while (nDone < nloc || PETSC_TRUE) { /* asyncronous not implemented */
    105:     iter++;
    106:     /* check all vertices */
    107:     for (kk=0; kk<nloc; kk++) {
    108:       lid   = perm_ix[kk];
    109:       state = lid_state[lid];
    110:       if (lid_removed[lid]) continue;
    111:       if (state == MIS_NOT_DONE) {
    112:         /* parallel test, delete if selected ghost */
    113:         isOK = PETSC_TRUE;
    114:         if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
    115:           ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
    116:           idx = matB->j + ii[ix];
    117:           for (j=0; j<n; j++) {
    118:             cpid   = idx[j]; /* compressed row ID in B mat */
    119:             gid    = cpcol_gid[cpid];
    120:             statej = cpcol_state[cpid];
    121:             if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */
    122:               isOK = PETSC_FALSE; /* can not delete */
    123:               break;
    124:             }
    125:           }
    126:         } /* parallel test */
    127:         if (isOK) { /* select or remove this vertex */
    128:           nDone++;
    129:           /* check for singleton */
    130:           ii = matA->i; n = ii[lid+1] - ii[lid];
    131:           if (n < 2) {
    132:             /* if I have any ghost adj then not a sing */
    133:             ix = lid_cprowID[lid];
    134:             if (ix==-1 || !(matB->compressedrow.i[ix+1]-matB->compressedrow.i[ix])) {
    135:               nremoved++;
    136:               lid_removed[lid] = PETSC_TRUE;
    137:               /* should select this because it is technically in the MIS but lets not */
    138:               continue; /* one local adj (me) and no ghost - singleton */
    139:             }
    140:           }
    141:           /* SELECTED state encoded with global index */
    142:           lid_state[lid] = lid+my0; /* needed???? */
    143:           nselected++;
    144:           if (strict_aggs) {
    145:             PetscCDAppendID(agg_lists, lid, lid+my0);
    146:           } else {
    147:             PetscCDAppendID(agg_lists, lid, lid);
    148:           }
    149:           /* delete local adj */
    150:           idx = matA->j + ii[lid];
    151:           for (j=0; j<n; j++) {
    152:             lidj   = idx[j];
    153:             statej = lid_state[lidj];
    154:             if (statej == MIS_NOT_DONE) {
    155:               nDone++;
    156:               if (strict_aggs) {
    157:                 PetscCDAppendID(agg_lists, lid, lidj+my0);
    158:               } else {
    159:                 PetscCDAppendID(agg_lists, lid, lidj);
    160:               }
    161:               lid_state[lidj] = MIS_DELETED;  /* delete this */
    162:             }
    163:           }
    164:           /* delete ghost adj of lid - deleted ghost done later for strict_aggs */
    165:           if (!strict_aggs) {
    166:             if ((ix=lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
    167:               ii  = matB->compressedrow.i; n = ii[ix+1] - ii[ix];
    168:               idx = matB->j + ii[ix];
    169:               for (j=0; j<n; j++) {
    170:                 cpid   = idx[j]; /* compressed row ID in B mat */
    171:                 statej = cpcol_state[cpid];
    172:                 if (statej == MIS_NOT_DONE) {
    173:                   PetscCDAppendID(agg_lists, lid, nloc+cpid);
    174:                 }
    175:               }
    176:             }
    177:           }
    178:         } /* selected */
    179:       } /* not done vertex */
    180:     } /* vertex loop */

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

    222:   /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */
    223:   if (strict_aggs && matB) {
    224:     /* need to copy this to free buffer -- should do this globaly */
    225:     PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid);
    226:     PetscMalloc1(num_fine_ghosts, &icpcol_gid);
    227:     for (cpid=0; cpid<num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid];

    229:     /* get proc of deleted ghost */
    230:     PetscSFBcastBegin(sf,MPIU_INT,lid_parent_gid,cpcol_sel_gid,MPI_REPLACE);
    231:     PetscSFBcastEnd(sf,MPIU_INT,lid_parent_gid,cpcol_sel_gid,MPI_REPLACE);
    232:     for (cpid=0; cpid<num_fine_ghosts; cpid++) {
    233:       sgid = cpcol_sel_gid[cpid];
    234:       gid  = icpcol_gid[cpid];
    235:       if (sgid >= my0 && sgid < Iend) { /* I own this deleted */
    236:         slid = sgid - my0;
    237:         PetscCDAppendID(agg_lists, slid, gid);
    238:       }
    239:     }
    240:     PetscFree(icpcol_gid);
    241:     PetscFree(cpcol_sel_gid);
    242:   }
    243:   if (mpimat) {
    244:     PetscSFDestroy(&sf);
    245:     PetscFree(cpcol_gid);
    246:     PetscFree(cpcol_state);
    247:   }
    248:   PetscFree(lid_cprowID);
    249:   PetscFree(lid_gid);
    250:   PetscFree(lid_removed);
    251:   if (strict_aggs) {
    252:     PetscFree(lid_parent_gid);
    253:   }
    254:   PetscFree(lid_state);
    255:   return(0);
    256: }

    258: /*
    259:    MIS coarsen, simple greedy.
    260: */
    261: static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
    262: {
    264:   Mat            mat = coarse->graph;

    267:   if (!coarse->perm) {
    268:     IS       perm;
    269:     PetscInt n,m;
    270:     MPI_Comm comm;

    272:     PetscObjectGetComm((PetscObject)mat,&comm);
    273:     MatGetLocalSize(mat, &m, &n);
    274:     ISCreateStride(comm, m, 0, 1, &perm);
    275:     maxIndSetAgg(perm, mat, coarse->strict_aggs, &coarse->agg_lists);
    276:     ISDestroy(&perm);
    277:   } else {
    278:     maxIndSetAgg(coarse->perm, mat, coarse->strict_aggs,  &coarse->agg_lists);
    279:   }
    280:   return(0);
    281: }

    283: PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)
    284: {
    286:   PetscMPIInt    rank;
    287:   PetscBool      iascii;

    290:   MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
    291:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    292:   if (iascii) {
    293:     PetscViewerASCIIPushSynchronized(viewer);
    294:     PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] MIS aggregator\n",rank);
    295:     PetscViewerFlush(viewer);
    296:     PetscViewerASCIIPopSynchronized(viewer);
    297:   }
    298:   return(0);
    299: }

    301: /*MC
    302:    MATCOARSENMIS - Creates a coarsen context via the external package MIS.

    304:    Collective

    306:    Input Parameter:
    307: .  coarse - the coarsen context

    309:    Options Database Keys:
    310: .  -mat_coarsen_MIS_xxx -

    312:    Level: beginner

    314: .seealso: MatCoarsenSetType(), MatCoarsenType

    316: M*/

    318: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse)
    319: {
    321:   coarse->ops->apply = MatCoarsenApply_MIS;
    322:   coarse->ops->view  = MatCoarsenView_MIS;
    323:   /* coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MIS; */
    324:   return(0);
    325: }