xref: /petsc/src/mat/graphops/coarsen/impls/misk/misk.c (revision 8be712e46db5d855f641c6bd97b4543e0efe65bd)
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