1*8be712e4SBarry Smith #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2*8be712e4SBarry Smith #include <../src/mat/impls/aij/seq/aij.h> 3*8be712e4SBarry Smith #include <../src/mat/impls/aij/mpi/mpiaij.h> 4*8be712e4SBarry Smith #include <petscsf.h> 5*8be712e4SBarry Smith 6*8be712e4SBarry Smith #define MIS_NOT_DONE -2 7*8be712e4SBarry Smith #define MIS_DELETED -1 8*8be712e4SBarry Smith #define MIS_REMOVED -3 9*8be712e4SBarry Smith #define MIS_IS_SELECTED(s) (s >= 0) 10*8be712e4SBarry Smith 11*8be712e4SBarry Smith /* edge for priority queue */ 12*8be712e4SBarry Smith typedef struct edge_tag { 13*8be712e4SBarry Smith PetscReal weight; 14*8be712e4SBarry Smith PetscInt lid0, gid1, cpid1; 15*8be712e4SBarry Smith } Edge; 16*8be712e4SBarry Smith 17*8be712e4SBarry Smith static PetscErrorCode PetscCoarsenDataView_private(PetscCoarsenData *agg_lists, PetscViewer viewer) 18*8be712e4SBarry Smith { 19*8be712e4SBarry Smith PetscCDIntNd *pos, *pos2; 20*8be712e4SBarry Smith 21*8be712e4SBarry Smith PetscFunctionBegin; 22*8be712e4SBarry Smith for (PetscInt kk = 0; kk < agg_lists->size; kk++) { 23*8be712e4SBarry Smith PetscCall(PetscCDGetHeadPos(agg_lists, kk, &pos)); 24*8be712e4SBarry Smith if ((pos2 = pos)) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "selected local %d: ", (int)kk)); 25*8be712e4SBarry Smith while (pos) { 26*8be712e4SBarry Smith PetscInt gid1; 27*8be712e4SBarry Smith PetscCall(PetscCDIntNdGetID(pos, &gid1)); 28*8be712e4SBarry Smith PetscCall(PetscCDGetNextPos(agg_lists, kk, &pos)); 29*8be712e4SBarry Smith PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %d ", (int)gid1)); 30*8be712e4SBarry Smith } 31*8be712e4SBarry Smith if (pos2) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n")); 32*8be712e4SBarry Smith } 33*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 34*8be712e4SBarry Smith } 35*8be712e4SBarry Smith 36*8be712e4SBarry Smith /* 37*8be712e4SBarry Smith MatCoarsenApply_MISK_private - parallel heavy edge matching 38*8be712e4SBarry Smith 39*8be712e4SBarry Smith Input Parameter: 40*8be712e4SBarry Smith . perm - permutation 41*8be712e4SBarry Smith . Gmat - global matrix of graph (data not defined) 42*8be712e4SBarry Smith 43*8be712e4SBarry Smith Output Parameter: 44*8be712e4SBarry Smith . a_locals_llist - array of list of local nodes rooted at local node 45*8be712e4SBarry Smith */ 46*8be712e4SBarry Smith static PetscErrorCode MatCoarsenApply_MISK_private(IS perm, const PetscInt misk, Mat Gmat, PetscCoarsenData **a_locals_llist) 47*8be712e4SBarry Smith { 48*8be712e4SBarry Smith PetscBool isMPI; 49*8be712e4SBarry Smith MPI_Comm comm; 50*8be712e4SBarry Smith PetscMPIInt rank, size; 51*8be712e4SBarry Smith Mat cMat, Prols[5], Rtot; 52*8be712e4SBarry Smith PetscScalar one = 1; 53*8be712e4SBarry Smith 54*8be712e4SBarry Smith PetscFunctionBegin; 55*8be712e4SBarry Smith PetscValidHeaderSpecific(perm, IS_CLASSID, 1); 56*8be712e4SBarry Smith PetscValidHeaderSpecific(Gmat, MAT_CLASSID, 3); 57*8be712e4SBarry Smith PetscAssertPointer(a_locals_llist, 4); 58*8be712e4SBarry Smith PetscCheck(misk < 5 && misk > 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "too many/few levels: %d", (int)misk); 59*8be712e4SBarry Smith PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPI)); 60*8be712e4SBarry Smith PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm)); 61*8be712e4SBarry Smith PetscCallMPI(MPI_Comm_rank(comm, &rank)); 62*8be712e4SBarry Smith PetscCallMPI(MPI_Comm_size(comm, &size)); 63*8be712e4SBarry Smith PetscCall(PetscInfo(Gmat, "misk %d\n", (int)misk)); 64*8be712e4SBarry Smith /* make a copy of the graph, this gets destroyed in iterates */ 65*8be712e4SBarry Smith if (misk > 1) PetscCall(MatDuplicate(Gmat, MAT_COPY_VALUES, &cMat)); 66*8be712e4SBarry Smith else cMat = Gmat; 67*8be712e4SBarry Smith for (PetscInt iterIdx = 0; iterIdx < misk; iterIdx++) { 68*8be712e4SBarry Smith Mat_SeqAIJ *matA, *matB = NULL; 69*8be712e4SBarry Smith Mat_MPIAIJ *mpimat = NULL; 70*8be712e4SBarry Smith const PetscInt *perm_ix; 71*8be712e4SBarry Smith const PetscInt nloc_inner = cMat->rmap->n; 72*8be712e4SBarry Smith PetscCoarsenData *agg_lists; 73*8be712e4SBarry Smith PetscInt *cpcol_gid = NULL, *cpcol_state, *lid_cprowID, *lid_state, *lid_parent_gid = NULL; 74*8be712e4SBarry Smith PetscInt num_fine_ghosts, kk, n, ix, j, *idx, *ai, Iend, my0, nremoved, gid, cpid, lidj, sgid, t1, t2, slid, nDone, nselected = 0, state; 75*8be712e4SBarry Smith PetscBool *lid_removed, isOK; 76*8be712e4SBarry Smith PetscLayout layout; 77*8be712e4SBarry Smith PetscSF sf; 78*8be712e4SBarry Smith 79*8be712e4SBarry Smith if (isMPI) { 80*8be712e4SBarry Smith mpimat = (Mat_MPIAIJ *)cMat->data; 81*8be712e4SBarry Smith matA = (Mat_SeqAIJ *)mpimat->A->data; 82*8be712e4SBarry Smith matB = (Mat_SeqAIJ *)mpimat->B->data; 83*8be712e4SBarry Smith /* force compressed storage of B */ 84*8be712e4SBarry Smith PetscCall(MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, cMat->rmap->n, -1.0)); 85*8be712e4SBarry Smith } else { 86*8be712e4SBarry Smith PetscBool isAIJ; 87*8be712e4SBarry Smith PetscCall(PetscObjectBaseTypeCompare((PetscObject)cMat, MATSEQAIJ, &isAIJ)); 88*8be712e4SBarry Smith PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix."); 89*8be712e4SBarry Smith matA = (Mat_SeqAIJ *)cMat->data; 90*8be712e4SBarry Smith } 91*8be712e4SBarry Smith PetscCall(MatGetOwnershipRange(cMat, &my0, &Iend)); 92*8be712e4SBarry Smith if (mpimat) { 93*8be712e4SBarry Smith PetscInt *lid_gid; 94*8be712e4SBarry Smith PetscCall(PetscMalloc1(nloc_inner, &lid_gid)); /* explicit array needed */ 95*8be712e4SBarry Smith for (kk = 0, gid = my0; kk < nloc_inner; kk++, gid++) lid_gid[kk] = gid; 96*8be712e4SBarry Smith PetscCall(VecGetLocalSize(mpimat->lvec, &num_fine_ghosts)); 97*8be712e4SBarry Smith PetscCall(PetscMalloc1(num_fine_ghosts, &cpcol_gid)); 98*8be712e4SBarry Smith PetscCall(PetscMalloc1(num_fine_ghosts, &cpcol_state)); 99*8be712e4SBarry Smith PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)cMat), &sf)); 100*8be712e4SBarry Smith PetscCall(MatGetLayouts(cMat, &layout, NULL)); 101*8be712e4SBarry Smith PetscCall(PetscSFSetGraphLayout(sf, layout, num_fine_ghosts, NULL, PETSC_COPY_VALUES, mpimat->garray)); 102*8be712e4SBarry Smith PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE)); 103*8be712e4SBarry Smith PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_gid, cpcol_gid, MPI_REPLACE)); 104*8be712e4SBarry Smith for (kk = 0; kk < num_fine_ghosts; kk++) cpcol_state[kk] = MIS_NOT_DONE; 105*8be712e4SBarry Smith PetscCall(PetscFree(lid_gid)); 106*8be712e4SBarry Smith } else num_fine_ghosts = 0; 107*8be712e4SBarry Smith 108*8be712e4SBarry Smith PetscCall(PetscMalloc1(nloc_inner, &lid_cprowID)); 109*8be712e4SBarry Smith PetscCall(PetscMalloc1(nloc_inner, &lid_removed)); /* explicit array needed */ 110*8be712e4SBarry Smith PetscCall(PetscMalloc1(nloc_inner, &lid_parent_gid)); 111*8be712e4SBarry Smith PetscCall(PetscMalloc1(nloc_inner, &lid_state)); 112*8be712e4SBarry Smith 113*8be712e4SBarry Smith /* the data structure */ 114*8be712e4SBarry Smith PetscCall(PetscCDCreate(nloc_inner, &agg_lists)); 115*8be712e4SBarry Smith /* need an inverse map - locals */ 116*8be712e4SBarry Smith for (kk = 0; kk < nloc_inner; kk++) { 117*8be712e4SBarry Smith lid_cprowID[kk] = -1; 118*8be712e4SBarry Smith lid_removed[kk] = PETSC_FALSE; 119*8be712e4SBarry Smith lid_parent_gid[kk] = -1.0; 120*8be712e4SBarry Smith lid_state[kk] = MIS_NOT_DONE; 121*8be712e4SBarry Smith } 122*8be712e4SBarry Smith /* set index into cmpressed row 'lid_cprowID' */ 123*8be712e4SBarry Smith if (matB) { 124*8be712e4SBarry Smith for (ix = 0; ix < matB->compressedrow.nrows; ix++) { 125*8be712e4SBarry Smith const PetscInt lid = matB->compressedrow.rindex[ix]; 126*8be712e4SBarry Smith if (lid >= 0) lid_cprowID[lid] = ix; 127*8be712e4SBarry Smith } 128*8be712e4SBarry Smith } 129*8be712e4SBarry Smith /* MIS */ 130*8be712e4SBarry Smith nremoved = nDone = 0; 131*8be712e4SBarry Smith if (!iterIdx) PetscCall(ISGetIndices(perm, &perm_ix)); // use permutation on first MIS 132*8be712e4SBarry Smith else perm_ix = NULL; 133*8be712e4SBarry Smith while (nDone < nloc_inner || PETSC_TRUE) { /* asynchronous not implemented */ 134*8be712e4SBarry Smith /* check all vertices */ 135*8be712e4SBarry Smith for (kk = 0; kk < nloc_inner; kk++) { 136*8be712e4SBarry Smith const PetscInt lid = perm_ix ? perm_ix[kk] : kk; 137*8be712e4SBarry Smith state = lid_state[lid]; 138*8be712e4SBarry Smith if (iterIdx == 0 && lid_removed[lid]) continue; 139*8be712e4SBarry Smith if (state == MIS_NOT_DONE) { 140*8be712e4SBarry Smith /* parallel test, delete if selected ghost */ 141*8be712e4SBarry Smith isOK = PETSC_TRUE; 142*8be712e4SBarry Smith /* parallel test */ 143*8be712e4SBarry Smith if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */ 144*8be712e4SBarry Smith ai = matB->compressedrow.i; 145*8be712e4SBarry Smith n = ai[ix + 1] - ai[ix]; 146*8be712e4SBarry Smith idx = matB->j + ai[ix]; 147*8be712e4SBarry Smith for (j = 0; j < n; j++) { 148*8be712e4SBarry Smith cpid = idx[j]; /* compressed row ID in B mat */ 149*8be712e4SBarry Smith gid = cpcol_gid[cpid]; 150*8be712e4SBarry Smith if (cpcol_state[cpid] == MIS_NOT_DONE && gid >= Iend) { /* or pe>rank */ 151*8be712e4SBarry Smith isOK = PETSC_FALSE; /* can not delete */ 152*8be712e4SBarry Smith break; 153*8be712e4SBarry Smith } 154*8be712e4SBarry Smith } 155*8be712e4SBarry Smith } 156*8be712e4SBarry Smith if (isOK) { /* select or remove this vertex if it is a true singleton like a BC */ 157*8be712e4SBarry Smith nDone++; 158*8be712e4SBarry Smith /* check for singleton */ 159*8be712e4SBarry Smith ai = matA->i; 160*8be712e4SBarry Smith n = ai[lid + 1] - ai[lid]; 161*8be712e4SBarry Smith if (n < 2) { 162*8be712e4SBarry Smith /* if I have any ghost adj then not a singleton */ 163*8be712e4SBarry Smith ix = lid_cprowID[lid]; 164*8be712e4SBarry Smith if (ix == -1 || !(matB->compressedrow.i[ix + 1] - matB->compressedrow.i[ix])) { 165*8be712e4SBarry Smith if (iterIdx == 0) { 166*8be712e4SBarry Smith lid_removed[lid] = PETSC_TRUE; 167*8be712e4SBarry Smith nremoved++; // let it get selected 168*8be712e4SBarry Smith } 169*8be712e4SBarry Smith // PetscCall(PetscCDAppendID(agg_lists, lid, lid + my0)); 170*8be712e4SBarry Smith // lid_state[lid] = nselected; // >= 0 is selected, cache for ordering coarse grid 171*8be712e4SBarry Smith /* should select this because it is technically in the MIS but lets not */ 172*8be712e4SBarry Smith continue; /* one local adj (me) and no ghost - singleton */ 173*8be712e4SBarry Smith } 174*8be712e4SBarry Smith } 175*8be712e4SBarry Smith /* SELECTED state encoded with global index */ 176*8be712e4SBarry Smith lid_state[lid] = nselected; // >= 0 is selected, cache for ordering coarse grid 177*8be712e4SBarry Smith nselected++; 178*8be712e4SBarry Smith PetscCall(PetscCDAppendID(agg_lists, lid, lid + my0)); 179*8be712e4SBarry Smith /* delete local adj */ 180*8be712e4SBarry Smith idx = matA->j + ai[lid]; 181*8be712e4SBarry Smith for (j = 0; j < n; j++) { 182*8be712e4SBarry Smith lidj = idx[j]; 183*8be712e4SBarry Smith if (lid_state[lidj] == MIS_NOT_DONE) { 184*8be712e4SBarry Smith nDone++; 185*8be712e4SBarry Smith PetscCall(PetscCDAppendID(agg_lists, lid, lidj + my0)); 186*8be712e4SBarry Smith lid_state[lidj] = MIS_DELETED; /* delete this */ 187*8be712e4SBarry Smith } 188*8be712e4SBarry Smith } 189*8be712e4SBarry Smith } /* selected */ 190*8be712e4SBarry Smith } /* not done vertex */ 191*8be712e4SBarry Smith } /* vertex loop */ 192*8be712e4SBarry Smith 193*8be712e4SBarry Smith /* update ghost states and count todos */ 194*8be712e4SBarry Smith if (mpimat) { 195*8be712e4SBarry Smith /* scatter states, check for done */ 196*8be712e4SBarry Smith PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE)); 197*8be712e4SBarry Smith PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_state, cpcol_state, MPI_REPLACE)); 198*8be712e4SBarry Smith ai = matB->compressedrow.i; 199*8be712e4SBarry Smith for (ix = 0; ix < matB->compressedrow.nrows; ix++) { 200*8be712e4SBarry Smith const int lidj = matB->compressedrow.rindex[ix]; /* local boundary node */ 201*8be712e4SBarry Smith state = lid_state[lidj]; 202*8be712e4SBarry Smith if (state == MIS_NOT_DONE) { 203*8be712e4SBarry Smith /* look at ghosts */ 204*8be712e4SBarry Smith n = ai[ix + 1] - ai[ix]; 205*8be712e4SBarry Smith idx = matB->j + ai[ix]; 206*8be712e4SBarry Smith for (j = 0; j < n; j++) { 207*8be712e4SBarry Smith cpid = idx[j]; /* compressed row ID in B mat */ 208*8be712e4SBarry Smith if (MIS_IS_SELECTED(cpcol_state[cpid])) { /* lid is now deleted by ghost */ 209*8be712e4SBarry Smith nDone++; 210*8be712e4SBarry Smith lid_state[lidj] = MIS_DELETED; /* delete this */ 211*8be712e4SBarry Smith sgid = cpcol_gid[cpid]; 212*8be712e4SBarry Smith lid_parent_gid[lidj] = sgid; /* keep track of proc that I belong to */ 213*8be712e4SBarry Smith break; 214*8be712e4SBarry Smith } 215*8be712e4SBarry Smith } 216*8be712e4SBarry Smith } 217*8be712e4SBarry Smith } 218*8be712e4SBarry Smith /* all done? */ 219*8be712e4SBarry Smith t1 = nloc_inner - nDone; 220*8be712e4SBarry Smith PetscCall(MPIU_Allreduce(&t1, &t2, 1, MPIU_INT, MPI_SUM, comm)); /* synchronous version */ 221*8be712e4SBarry Smith if (!t2) break; 222*8be712e4SBarry Smith } else break; /* no mpi - all done */ 223*8be712e4SBarry Smith } /* outer parallel MIS loop */ 224*8be712e4SBarry Smith if (!iterIdx) PetscCall(ISRestoreIndices(perm, &perm_ix)); 225*8be712e4SBarry Smith PetscCall(PetscInfo(Gmat, "\t removed %" PetscInt_FMT " of %" PetscInt_FMT " vertices. %" PetscInt_FMT " selected.\n", nremoved, nloc_inner, nselected)); 226*8be712e4SBarry Smith 227*8be712e4SBarry Smith /* tell adj who my lid_parent_gid vertices belong to - fill in agg_lists selected ghost lists */ 228*8be712e4SBarry Smith if (matB) { 229*8be712e4SBarry Smith PetscInt *cpcol_sel_gid, *icpcol_gid; 230*8be712e4SBarry Smith /* need to copy this to free buffer -- should do this globally */ 231*8be712e4SBarry Smith PetscCall(PetscMalloc1(num_fine_ghosts, &cpcol_sel_gid)); 232*8be712e4SBarry Smith PetscCall(PetscMalloc1(num_fine_ghosts, &icpcol_gid)); 233*8be712e4SBarry Smith for (cpid = 0; cpid < num_fine_ghosts; cpid++) icpcol_gid[cpid] = cpcol_gid[cpid]; 234*8be712e4SBarry Smith /* get proc of deleted ghost */ 235*8be712e4SBarry Smith PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE)); 236*8be712e4SBarry Smith PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lid_parent_gid, cpcol_sel_gid, MPI_REPLACE)); 237*8be712e4SBarry Smith for (cpid = 0; cpid < num_fine_ghosts; cpid++) { 238*8be712e4SBarry Smith sgid = cpcol_sel_gid[cpid]; 239*8be712e4SBarry Smith gid = icpcol_gid[cpid]; 240*8be712e4SBarry Smith if (sgid >= my0 && sgid < Iend) { /* I own this deleted */ 241*8be712e4SBarry Smith slid = sgid - my0; 242*8be712e4SBarry Smith PetscCall(PetscCDAppendID(agg_lists, slid, gid)); 243*8be712e4SBarry Smith } 244*8be712e4SBarry Smith } 245*8be712e4SBarry Smith // done - cleanup 246*8be712e4SBarry Smith PetscCall(PetscFree(icpcol_gid)); 247*8be712e4SBarry Smith PetscCall(PetscFree(cpcol_sel_gid)); 248*8be712e4SBarry Smith PetscCall(PetscSFDestroy(&sf)); 249*8be712e4SBarry Smith PetscCall(PetscFree(cpcol_gid)); 250*8be712e4SBarry Smith PetscCall(PetscFree(cpcol_state)); 251*8be712e4SBarry Smith } 252*8be712e4SBarry Smith PetscCall(PetscFree(lid_cprowID)); 253*8be712e4SBarry Smith PetscCall(PetscFree(lid_removed)); 254*8be712e4SBarry Smith PetscCall(PetscFree(lid_parent_gid)); 255*8be712e4SBarry Smith PetscCall(PetscFree(lid_state)); 256*8be712e4SBarry Smith 257*8be712e4SBarry Smith /* MIS done - make projection matrix - P */ 258*8be712e4SBarry Smith MatType jtype; 259*8be712e4SBarry Smith PetscCall(MatGetType(Gmat, &jtype)); 260*8be712e4SBarry Smith PetscCall(MatCreate(comm, &Prols[iterIdx])); 261*8be712e4SBarry Smith PetscCall(MatSetType(Prols[iterIdx], jtype)); 262*8be712e4SBarry Smith PetscCall(MatSetSizes(Prols[iterIdx], nloc_inner, nselected, PETSC_DETERMINE, PETSC_DETERMINE)); 263*8be712e4SBarry Smith PetscCall(MatSeqAIJSetPreallocation(Prols[iterIdx], 1, NULL)); 264*8be712e4SBarry Smith PetscCall(MatMPIAIJSetPreallocation(Prols[iterIdx], 1, NULL, 1, NULL)); 265*8be712e4SBarry Smith { 266*8be712e4SBarry Smith PetscCDIntNd *pos, *pos2; 267*8be712e4SBarry Smith PetscInt colIndex, Iend, fgid; 268*8be712e4SBarry Smith PetscCall(MatGetOwnershipRangeColumn(Prols[iterIdx], &colIndex, &Iend)); 269*8be712e4SBarry Smith // TODO - order with permutation in lid_selected (reversed) 270*8be712e4SBarry Smith for (PetscInt lid = 0; lid < agg_lists->size; lid++) { 271*8be712e4SBarry Smith PetscCall(PetscCDGetHeadPos(agg_lists, lid, &pos)); 272*8be712e4SBarry Smith pos2 = pos; 273*8be712e4SBarry Smith while (pos) { 274*8be712e4SBarry Smith PetscCall(PetscCDIntNdGetID(pos, &fgid)); 275*8be712e4SBarry Smith PetscCall(PetscCDGetNextPos(agg_lists, lid, &pos)); 276*8be712e4SBarry Smith PetscCall(MatSetValues(Prols[iterIdx], 1, &fgid, 1, &colIndex, &one, INSERT_VALUES)); 277*8be712e4SBarry Smith } 278*8be712e4SBarry Smith if (pos2) colIndex++; 279*8be712e4SBarry Smith } 280*8be712e4SBarry Smith PetscCheck(Iend == colIndex, PETSC_COMM_SELF, PETSC_ERR_SUP, "Iend!=colIndex: %d %d", (int)Iend, (int)colIndex); 281*8be712e4SBarry Smith } 282*8be712e4SBarry Smith PetscCall(MatAssemblyBegin(Prols[iterIdx], MAT_FINAL_ASSEMBLY)); 283*8be712e4SBarry Smith PetscCall(MatAssemblyEnd(Prols[iterIdx], MAT_FINAL_ASSEMBLY)); 284*8be712e4SBarry Smith /* project to make new graph for next MIS, skip if last */ 285*8be712e4SBarry Smith if (iterIdx < misk - 1) { 286*8be712e4SBarry Smith Mat new_mat; 287*8be712e4SBarry Smith PetscCall(MatPtAP(cMat, Prols[iterIdx], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &new_mat)); 288*8be712e4SBarry Smith PetscCall(MatDestroy(&cMat)); 289*8be712e4SBarry Smith cMat = new_mat; // next iter 290*8be712e4SBarry Smith } else if (cMat != Gmat) PetscCall(MatDestroy(&cMat)); 291*8be712e4SBarry Smith // cleanup 292*8be712e4SBarry Smith PetscCall(PetscCDDestroy(agg_lists)); 293*8be712e4SBarry Smith } /* MIS-k iteration */ 294*8be712e4SBarry Smith /* make total prolongator Rtot = P_0 * P_1 * ... */ 295*8be712e4SBarry Smith Rtot = Prols[misk - 1]; // compose P then transpose to get R 296*8be712e4SBarry Smith for (PetscInt iterIdx = misk - 1; iterIdx > 0; iterIdx--) { 297*8be712e4SBarry Smith Mat P; 298*8be712e4SBarry Smith PetscCall(MatMatMult(Prols[iterIdx - 1], Rtot, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &P)); 299*8be712e4SBarry Smith PetscCall(MatDestroy(&Prols[iterIdx - 1])); 300*8be712e4SBarry Smith PetscCall(MatDestroy(&Rtot)); 301*8be712e4SBarry Smith Rtot = P; 302*8be712e4SBarry Smith } 303*8be712e4SBarry Smith PetscCall(MatTranspose(Rtot, MAT_INPLACE_MATRIX, &Rtot)); // R now 304*8be712e4SBarry Smith PetscCall(MatViewFromOptions(Rtot, NULL, "-misk_aggregation_view")); 305*8be712e4SBarry Smith /* make aggregates with Rtot - could use Rtot directly in theory but have to go through the aggregate list data structure */ 306*8be712e4SBarry Smith { 307*8be712e4SBarry Smith PetscInt Istart, Iend, ncols, NN, MM, jj = 0, max_osz = 0; 308*8be712e4SBarry Smith const PetscInt nloc = Gmat->rmap->n; 309*8be712e4SBarry Smith PetscCoarsenData *agg_lists; 310*8be712e4SBarry Smith Mat mat; 311*8be712e4SBarry Smith PetscCall(PetscCDCreate(nloc, &agg_lists)); 312*8be712e4SBarry Smith *a_locals_llist = agg_lists; // return 313*8be712e4SBarry Smith PetscCall(MatGetOwnershipRange(Rtot, &Istart, &Iend)); 314*8be712e4SBarry Smith for (int grow = Istart, lid = 0; grow < Iend; grow++, lid++) { 315*8be712e4SBarry Smith const PetscInt *idx; 316*8be712e4SBarry Smith PetscCall(MatGetRow(Rtot, grow, &ncols, &idx, NULL)); 317*8be712e4SBarry Smith for (int jj = 0; jj < ncols; jj++) { 318*8be712e4SBarry Smith PetscInt gcol = idx[jj]; 319*8be712e4SBarry Smith PetscCall(PetscCDAppendID(agg_lists, lid, gcol)); // local row, global column 320*8be712e4SBarry Smith } 321*8be712e4SBarry Smith PetscCall(MatRestoreRow(Rtot, grow, &ncols, &idx, NULL)); 322*8be712e4SBarry Smith } 323*8be712e4SBarry Smith PetscCall(MatDestroy(&Rtot)); 324*8be712e4SBarry Smith 325*8be712e4SBarry Smith /* make fake matrix, get largest nnz */ 326*8be712e4SBarry Smith for (int lid = 0; lid < nloc; lid++) { 327*8be712e4SBarry Smith PetscCall(PetscCDCountAt(agg_lists, lid, &jj)); 328*8be712e4SBarry Smith if (jj > max_osz) max_osz = jj; 329*8be712e4SBarry Smith } 330*8be712e4SBarry Smith PetscCall(MatGetSize(Gmat, &MM, &NN)); 331*8be712e4SBarry Smith if (max_osz > MM - nloc) max_osz = MM - nloc; 332*8be712e4SBarry Smith PetscCall(MatGetOwnershipRange(Gmat, &Istart, NULL)); 333*8be712e4SBarry Smith /* matrix of ghost adj for square graph */ 334*8be712e4SBarry Smith PetscCall(MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, max_osz, NULL, &mat)); 335*8be712e4SBarry Smith for (PetscInt lid = 0, gidi = Istart; lid < nloc; lid++, gidi++) { 336*8be712e4SBarry Smith PetscCDIntNd *pos; 337*8be712e4SBarry Smith PetscCall(PetscCDGetHeadPos(agg_lists, lid, &pos)); 338*8be712e4SBarry Smith while (pos) { 339*8be712e4SBarry Smith PetscInt gidj; 340*8be712e4SBarry Smith PetscCall(PetscCDIntNdGetID(pos, &gidj)); 341*8be712e4SBarry Smith PetscCall(PetscCDGetNextPos(agg_lists, lid, &pos)); 342*8be712e4SBarry Smith if (gidj < Istart || gidj >= Istart + nloc) PetscCall(MatSetValues(mat, 1, &gidi, 1, &gidj, &one, ADD_VALUES)); 343*8be712e4SBarry Smith } 344*8be712e4SBarry Smith } 345*8be712e4SBarry Smith PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 346*8be712e4SBarry Smith PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 347*8be712e4SBarry Smith PetscCall(PetscCDSetMat(agg_lists, mat)); 348*8be712e4SBarry Smith } 349*8be712e4SBarry Smith 350*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 351*8be712e4SBarry Smith } 352*8be712e4SBarry Smith 353*8be712e4SBarry Smith /* 354*8be712e4SBarry Smith Distance k MIS. k is in 'subctx' 355*8be712e4SBarry Smith */ 356*8be712e4SBarry Smith static PetscErrorCode MatCoarsenApply_MISK(MatCoarsen coarse) 357*8be712e4SBarry Smith { 358*8be712e4SBarry Smith Mat mat = coarse->graph; 359*8be712e4SBarry Smith PetscInt k; 360*8be712e4SBarry Smith 361*8be712e4SBarry Smith PetscFunctionBegin; 362*8be712e4SBarry Smith PetscCall(MatCoarsenMISKGetDistance(coarse, &k)); 363*8be712e4SBarry Smith PetscCheck(k > 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "too few levels: %d", (int)k); 364*8be712e4SBarry Smith if (!coarse->perm) { 365*8be712e4SBarry Smith IS perm; 366*8be712e4SBarry Smith PetscInt n, m; 367*8be712e4SBarry Smith 368*8be712e4SBarry Smith PetscCall(MatGetLocalSize(mat, &m, &n)); 369*8be712e4SBarry Smith PetscCall(ISCreateStride(PetscObjectComm((PetscObject)mat), m, 0, 1, &perm)); 370*8be712e4SBarry Smith PetscCall(MatCoarsenApply_MISK_private(perm, (PetscInt)k, mat, &coarse->agg_lists)); 371*8be712e4SBarry Smith PetscCall(ISDestroy(&perm)); 372*8be712e4SBarry Smith } else { 373*8be712e4SBarry Smith PetscCall(MatCoarsenApply_MISK_private(coarse->perm, (PetscInt)k, mat, &coarse->agg_lists)); 374*8be712e4SBarry Smith } 375*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 376*8be712e4SBarry Smith } 377*8be712e4SBarry Smith 378*8be712e4SBarry Smith static PetscErrorCode MatCoarsenView_MISK(MatCoarsen coarse, PetscViewer viewer) 379*8be712e4SBarry Smith { 380*8be712e4SBarry Smith PetscMPIInt rank; 381*8be712e4SBarry Smith PetscBool iascii; 382*8be712e4SBarry Smith 383*8be712e4SBarry Smith PetscFunctionBegin; 384*8be712e4SBarry Smith PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank)); 385*8be712e4SBarry Smith PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 386*8be712e4SBarry Smith if (iascii) { 387*8be712e4SBarry Smith PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 388*8be712e4SBarry Smith PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d] MISK aggregator\n", rank)); 389*8be712e4SBarry Smith if (!rank) PetscCall(PetscCoarsenDataView_private(coarse->agg_lists, viewer)); 390*8be712e4SBarry Smith PetscCall(PetscViewerFlush(viewer)); 391*8be712e4SBarry Smith PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 392*8be712e4SBarry Smith } 393*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 394*8be712e4SBarry Smith } 395*8be712e4SBarry Smith 396*8be712e4SBarry Smith static PetscErrorCode MatCoarsenSetFromOptions_MISK(MatCoarsen coarse, PetscOptionItems *PetscOptionsObject) 397*8be712e4SBarry Smith { 398*8be712e4SBarry Smith PetscInt k = 1; 399*8be712e4SBarry Smith PetscBool flg; 400*8be712e4SBarry Smith PetscFunctionBegin; 401*8be712e4SBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "MatCoarsen-MISk options"); 402*8be712e4SBarry Smith PetscCall(PetscOptionsInt("-mat_coarsen_misk_distance", "k distance for MIS", "", k, &k, &flg)); 403*8be712e4SBarry Smith if (flg) coarse->subctx = (void *)(size_t)k; 404*8be712e4SBarry Smith PetscOptionsHeadEnd(); 405*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 406*8be712e4SBarry Smith } 407*8be712e4SBarry Smith 408*8be712e4SBarry Smith /*MC 409*8be712e4SBarry Smith MATCOARSENMISK - A coarsener that uses MISK, a simple greedy coarsener 410*8be712e4SBarry Smith 411*8be712e4SBarry Smith Level: beginner 412*8be712e4SBarry Smith 413*8be712e4SBarry Smith Options Database Key: 414*8be712e4SBarry Smith . -mat_coarsen_misk_distance <k> - distance for MIS 415*8be712e4SBarry Smith 416*8be712e4SBarry Smith .seealso: `MatCoarsen`, `MatCoarsenMISKSetDistance()`, `MatCoarsenApply()`, `MatCoarsenSetType()`, `MatCoarsenType`, `MatCoarsenCreate()` 417*8be712e4SBarry Smith M*/ 418*8be712e4SBarry Smith 419*8be712e4SBarry Smith PETSC_EXTERN PetscErrorCode MatCoarsenCreate_MISK(MatCoarsen coarse) 420*8be712e4SBarry Smith { 421*8be712e4SBarry Smith PetscFunctionBegin; 422*8be712e4SBarry Smith coarse->ops->apply = MatCoarsenApply_MISK; 423*8be712e4SBarry Smith coarse->ops->view = MatCoarsenView_MISK; 424*8be712e4SBarry Smith coarse->subctx = (void *)(size_t)1; 425*8be712e4SBarry Smith coarse->ops->setfromoptions = MatCoarsenSetFromOptions_MISK; 426*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 427*8be712e4SBarry Smith } 428*8be712e4SBarry Smith 429*8be712e4SBarry Smith /*@ 430*8be712e4SBarry Smith MatCoarsenMISKSetDistance - the distance to be used by MISK 431*8be712e4SBarry Smith 432*8be712e4SBarry Smith Collective 433*8be712e4SBarry Smith 434*8be712e4SBarry Smith Input Parameters: 435*8be712e4SBarry Smith + crs - the coarsen 436*8be712e4SBarry Smith - k - the distance 437*8be712e4SBarry Smith 438*8be712e4SBarry Smith Options Database Key: 439*8be712e4SBarry Smith . -mat_coarsen_misk_distance <k> - distance for MIS 440*8be712e4SBarry Smith 441*8be712e4SBarry Smith Level: advanced 442*8be712e4SBarry Smith 443*8be712e4SBarry Smith .seealso: `MATCOARSENMISK`, `MatCoarsen`, `MatCoarseSetFromOptions()`, `MatCoarsenSetType()`, `MatCoarsenRegister()`, `MatCoarsenCreate()`, 444*8be712e4SBarry Smith `MatCoarsenDestroy()`, `MatCoarsenSetAdjacency()`, `MatCoarsenMISKGetDistance()` 445*8be712e4SBarry Smith `MatCoarsenGetData()` 446*8be712e4SBarry Smith @*/ 447*8be712e4SBarry Smith PetscErrorCode MatCoarsenMISKSetDistance(MatCoarsen crs, PetscInt k) 448*8be712e4SBarry Smith { 449*8be712e4SBarry Smith PetscFunctionBegin; 450*8be712e4SBarry Smith crs->subctx = (void *)(size_t)k; 451*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 452*8be712e4SBarry Smith } 453*8be712e4SBarry Smith 454*8be712e4SBarry Smith /*@ 455*8be712e4SBarry Smith MatCoarsenMISKGetDistance - gets the distance to be used by MISK 456*8be712e4SBarry Smith 457*8be712e4SBarry Smith Collective 458*8be712e4SBarry Smith 459*8be712e4SBarry Smith Input Parameter: 460*8be712e4SBarry Smith . crs - the coarsen 461*8be712e4SBarry Smith 462*8be712e4SBarry Smith Output Parameter: 463*8be712e4SBarry Smith . k - the distance 464*8be712e4SBarry Smith 465*8be712e4SBarry Smith Level: advanced 466*8be712e4SBarry Smith 467*8be712e4SBarry Smith .seealso: `MATCOARSENMISK`, `MatCoarsen`, `MatCoarseSetFromOptions()`, `MatCoarsenSetType()`, 468*8be712e4SBarry Smith `MatCoarsenRegister()`, `MatCoarsenCreate()`, `MatCoarsenDestroy()`, 469*8be712e4SBarry Smith `MatCoarsenSetAdjacency()`, `MatCoarsenGetData()` 470*8be712e4SBarry Smith @*/ 471*8be712e4SBarry Smith PetscErrorCode MatCoarsenMISKGetDistance(MatCoarsen crs, PetscInt *k) 472*8be712e4SBarry Smith { 473*8be712e4SBarry Smith PetscFunctionBegin; 474*8be712e4SBarry Smith *k = (PetscInt)(size_t)crs->subctx; 475*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS); 476*8be712e4SBarry Smith } 477