1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/aij/seq/aij.h> 3 #include <../src/mat/impls/aij/mpi/mpiaij.h> 4 #include <petscsf.h> 5 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) 10 11 /* 12 MatCoarsenApply_MIS_private - parallel maximal independent set (MIS) with data locality info. MatAIJ specific!!! 13 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'; 18 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; 38 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 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isAIJ)); 53 PetscCheck(isAIJ, comm, PETSC_ERR_PLIB, "Require AIJ matrix."); 54 matA = (Mat_SeqAIJ *)Gmat->data; 55 } 56 PetscCall(MatGetOwnershipRange(Gmat, &my0, &Iend)); 57 PetscCall(PetscMalloc1(nloc, &lid_gid)); /* explicit array needed */ 58 if (mpimat) { 59 for (kk = 0, gid = my0; kk < nloc; kk++, gid++) lid_gid[kk] = gid; 60 PetscCall(VecGetLocalSize(mpimat->lvec, &num_fine_ghosts)); 61 PetscCall(PetscMalloc1(num_fine_ghosts, &cpcol_gid)); 62 PetscCall(PetscMalloc1(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; 70 71 PetscCall(PetscMalloc1(nloc, &lid_cprowID)); 72 PetscCall(PetscMalloc1(nloc, &lid_removed)); /* explicit array needed */ 73 if (strict_aggs) PetscCall(PetscMalloc1(nloc, &lid_parent_gid)); 74 PetscCall(PetscMalloc1(nloc, &lid_state)); 75 76 /* has ghost nodes for !strict and uses local indexing (yuck) */ 77 PetscCall(PetscCDCreate(strict_aggs ? nloc : num_fine_ghosts + nloc, &agg_lists)); 78 if (a_locals_llist) *a_locals_llist = agg_lists; 79 80 /* need an inverse map - locals */ 81 for (kk = 0; kk < nloc; kk++) { 82 lid_cprowID[kk] = -1; 83 lid_removed[kk] = PETSC_FALSE; 84 if (strict_aggs) lid_parent_gid[kk] = -1.0; 85 lid_state[kk] = MIS_NOT_DONE; 86 } 87 /* set index into cmpressed row 'lid_cprowID' */ 88 if (matB) { 89 for (ix = 0; ix < matB->compressedrow.nrows; ix++) { 90 lid = matB->compressedrow.rindex[ix]; 91 if (lid >= 0) lid_cprowID[lid] = ix; 92 } 93 } 94 /* MIS */ 95 nremoved = nDone = 0; 96 97 PetscCall(ISGetIndices(perm, &perm_ix)); 98 while (nDone < nloc || PETSC_TRUE) { /* asynchronous not implemented */ 99 /* check all vertices */ 100 for (kk = 0; kk < nloc; kk++) { 101 lid = perm_ix[kk]; 102 state = lid_state[lid]; 103 if (lid_removed[lid]) continue; 104 if (state == MIS_NOT_DONE) { 105 /* parallel test, delete if selected ghost */ 106 isOK = PETSC_TRUE; 107 if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */ 108 ii = matB->compressedrow.i; 109 n = ii[ix + 1] - ii[ix]; 110 idx = matB->j + ii[ix]; 111 for (j = 0; j < n; j++) { 112 cpid = idx[j]; /* compressed row ID in B mat */ 113 gid = cpcol_gid[cpid]; 114 statej = cpcol_state[cpid]; 115 PetscCheck(!MIS_IS_SELECTED(statej), PETSC_COMM_SELF, PETSC_ERR_SUP, "selected ghost: %d", (int)gid); 116 if (statej == MIS_NOT_DONE && gid >= Iend) { /* should be (pe>rank), use gid as pe proxy */ 117 isOK = PETSC_FALSE; /* can not delete */ 118 break; 119 } 120 } 121 } /* parallel test */ 122 if (isOK) { /* select or remove this vertex */ 123 nDone++; 124 /* check for singleton */ 125 ii = matA->i; 126 n = ii[lid + 1] - ii[lid]; 127 if (n < 2) { 128 /* if I have any ghost adj then not a sing */ 129 ix = lid_cprowID[lid]; 130 if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) { 131 nremoved++; 132 nrm_tot++; 133 lid_removed[lid] = PETSC_TRUE; 134 continue; 135 // lid_state[lidj] = MIS_REMOVED; /* add singleton to MIS (can cause low rank with elasticity on fine grid) */ 136 } 137 } 138 /* SELECTED state encoded with global index */ 139 lid_state[lid] = lid + my0; 140 nselected++; 141 if (strict_aggs) { 142 PetscCall(PetscCDAppendID(agg_lists, lid, lid + my0)); 143 } else { 144 PetscCall(PetscCDAppendID(agg_lists, lid, lid)); 145 } 146 /* delete local adj */ 147 idx = matA->j + ii[lid]; 148 for (j = 0; j < n; j++) { 149 lidj = idx[j]; 150 statej = lid_state[lidj]; 151 if (statej == MIS_NOT_DONE) { 152 nDone++; 153 if (strict_aggs) { 154 PetscCall(PetscCDAppendID(agg_lists, lid, lidj + my0)); 155 } else { 156 PetscCall(PetscCDAppendID(agg_lists, lid, lidj)); 157 } 158 lid_state[lidj] = MIS_DELETED; /* delete this */ 159 } 160 } 161 /* delete ghost adj of lid - deleted ghost done later for strict_aggs */ 162 if (!strict_aggs) { 163 if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */ 164 ii = matB->compressedrow.i; 165 n = ii[ix + 1] - ii[ix]; 166 idx = matB->j + ii[ix]; 167 for (j = 0; j < n; j++) { 168 cpid = idx[j]; /* compressed row ID in B mat */ 169 statej = cpcol_state[cpid]; 170 if (statej == MIS_NOT_DONE) PetscCall(PetscCDAppendID(agg_lists, lid, nloc + cpid)); 171 } 172 } 173 } 174 } /* selected */ 175 } /* not done vertex */ 176 } /* vertex loop */ 177 178 /* update ghost states and count todos */ 179 if (mpimat) { 180 /* scatter states, check for done */ 181 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE)); 182 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE)); 183 ii = matB->compressedrow.i; 184 for (ix = 0; ix < matB->compressedrow.nrows; ix++) { 185 lid = matB->compressedrow.rindex[ix]; /* local boundary node */ 186 state = lid_state[lid]; 187 if (state == MIS_NOT_DONE) { 188 /* look at ghosts */ 189 n = ii[ix + 1] - ii[ix]; 190 idx = matB->j + ii[ix]; 191 for (j = 0; j < n; j++) { 192 cpid = idx[j]; /* compressed row ID in B mat */ 193 statej = cpcol_state[cpid]; 194 if (MIS_IS_SELECTED(statej)) { /* lid is now deleted, do it */ 195 nDone++; 196 lid_state[lid] = MIS_DELETED; /* delete this */ 197 if (!strict_aggs) { 198 lidj = nloc + cpid; 199 PetscCall(PetscCDAppendID(agg_lists, lidj, lid)); 200 } else { 201 sgid = cpcol_gid[cpid]; 202 lid_parent_gid[lid] = sgid; /* keep track of proc that I belong to */ 203 } 204 break; 205 } 206 } 207 } 208 } 209 /* all done? */ 210 t1 = nloc - nDone; 211 PetscCall(MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm)); /* synchronous version */ 212 if (!t2) break; 213 } else break; /* all done */ 214 } /* outer parallel MIS loop */ 215 PetscCall(ISRestoreIndices(perm, &perm_ix)); 216 PetscCall(PetscInfo(info_is, "\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices. %" PetscInt_FMT " selected.\n", nremoved, nloc, nselected)); 217 218 /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */ 219 if (strict_aggs && matB) { 220 /* need to copy this to free buffer -- should do this globally */ 221 PetscCall(PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid)); 222 PetscCall(PetscMalloc1(num_fine_ghosts, &icpcol_gid)); 223 for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid]; 224 225 /* get proc of deleted ghost */ 226 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE)); 227 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE)); 228 for (cpid = 0; cpid < num_fine_ghosts; cpid++) { 229 sgid = cpcol_sel_gid[cpid]; 230 gid = icpcol_gid[cpid]; 231 if (sgid >= my0 && sgid < Iend) { /* I own this deleted */ 232 slid = sgid - my0; 233 PetscCall(PetscCDAppendID(agg_lists, slid, gid)); 234 } 235 } 236 PetscCall(PetscFree(icpcol_gid)); 237 PetscCall(PetscFree(cpcol_sel_gid)); 238 } 239 if (mpimat) { 240 PetscCall(PetscSFDestroy(&sf)); 241 PetscCall(PetscFree(cpcol_gid)); 242 PetscCall(PetscFree(cpcol_state)); 243 } 244 PetscCall(PetscFree(lid_cprowID)); 245 PetscCall(PetscFree(lid_gid)); 246 PetscCall(PetscFree(lid_removed)); 247 if (strict_aggs) PetscCall(PetscFree(lid_parent_gid)); 248 PetscCall(PetscFree(lid_state)); 249 if (strict_aggs) { 250 // check sizes -- all vertices must get in graph 251 PetscInt aa[2] = {0, nrm_tot}, bb[2], MM; 252 PetscCall(MatGetSize(Gmat, &MM, NULL)); 253 // check sizes -- all vertices must get in graph 254 PetscCall(PetscCDCount(agg_lists, &aa[0])); 255 PetscCall(MPIU_Allreduce(aa, bb, 2, MPIU_INT, MPI_SUM, comm)); 256 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])); 257 PetscCheck(MM >= bb[0], comm, PETSC_ERR_PLIB, "Sum of aggs too big"); 258 } 259 PetscCall(ISDestroy(&info_is)); 260 PetscFunctionReturn(PETSC_SUCCESS); 261 } 262 263 /* 264 MIS coarsen, simple greedy. 265 */ 266 static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse) 267 { 268 Mat mat = coarse->graph; 269 270 PetscFunctionBegin; 271 if (!coarse->perm) { 272 IS perm; 273 PetscInt n, m; 274 MPI_Comm comm; 275 276 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 277 PetscCall(MatGetLocalSize(mat, &m, &n)); 278 PetscCall(ISCreateStride(comm, m, 0, 1, &perm)); 279 PetscCall(MatCoarsenApply_MIS_private(perm, mat, coarse->strict_aggs, &coarse->agg_lists)); 280 PetscCall(ISDestroy(&perm)); 281 } else { 282 PetscCall(MatCoarsenApply_MIS_private(coarse->perm, mat, coarse->strict_aggs, &coarse->agg_lists)); 283 } 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 static PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse, PetscViewer viewer) 288 { 289 PetscMPIInt rank; 290 PetscBool iascii; 291 292 PetscFunctionBegin; 293 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank)); 294 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 295 if (iascii) { 296 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 297 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d] MIS aggregator\n", rank)); 298 if (!rank) { 299 PetscCDIntNd *pos, *pos2; 300 for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) { 301 PetscCall(PetscCDGetHeadPos(coarse->agg_lists, kk, &pos)); 302 if ((pos2 = pos)) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "selected %d: ", (int)kk)); 303 while (pos) { 304 PetscInt gid1; 305 PetscCall(PetscCDIntNdGetID(pos, &gid1)); 306 PetscCall(PetscCDGetNextPos(coarse->agg_lists, kk, &pos)); 307 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %d ", (int)gid1)); 308 } 309 if (pos2) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n")); 310 } 311 } 312 PetscCall(PetscViewerFlush(viewer)); 313 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 314 } 315 PetscFunctionReturn(PETSC_SUCCESS); 316 } 317 318 /*MC 319 MATCOARSENMIS - Creates a coarsening with a maximal independent set (MIS) algorithm 320 321 Collective 322 323 Input Parameter: 324 . coarse - the coarsen context 325 326 Level: beginner 327 328 .seealso: `MatCoarsen`, `MatCoarsenApply()`, `MatCoarsenGetData()`, `MatCoarsenSetType()`, `MatCoarsenType` 329 M*/ 330 PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MIS(MatCoarsen coarse) 331 { 332 PetscFunctionBegin; 333 coarse->ops->apply = MatCoarsenApply_MIS; 334 coarse->ops->view = MatCoarsenView_MIS; 335 PetscFunctionReturn(PETSC_SUCCESS); 336 } 337