Actual source code: agg.c

  1: /*
  2:  GAMG geometric-algebric multigrid PC - Mark Adams 2011
  3:  */

  5: #include <../src/ksp/pc/impls/gamg/gamg.h>
  6: #include <petscblaslapack.h>
  7: #include <petscdm.h>
  8: #include <petsc/private/kspimpl.h>

 10: typedef struct {
 11:   PetscInt   nsmooths;
 12:   PetscInt   aggressive_coarsening_levels; // number of aggressive coarsening levels (square or MISk)
 13:   PetscInt   aggressive_mis_k;             // the k in MIS-k
 14:   PetscBool  use_aggressive_square_graph;
 15:   PetscBool  use_minimum_degree_ordering;
 16:   MatCoarsen crs;
 17: } PC_GAMG_AGG;

 19: /*@
 20:   PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used for multigrid on all the levels

 22:   Logically Collective

 24:   Input Parameters:
 25: + pc - the preconditioner context
 26: - n  - the number of smooths

 28:   Options Database Key:
 29: . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation

 31:   Level: intermediate

 33: .seealso: [](ch_ksp), `PCMG`, `PCGAMG`
 34: @*/
 35: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
 36: {
 37:   PetscFunctionBegin;
 40:   PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
 41:   PetscFunctionReturn(PETSC_SUCCESS);
 42: }

 44: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
 45: {
 46:   PC_MG       *mg          = (PC_MG *)pc->data;
 47:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
 48:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

 50:   PetscFunctionBegin;
 51:   pc_gamg_agg->nsmooths = n;
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 53: }

 55: /*@
 56:   PCGAMGSetAggressiveLevels -  Use aggressive coarsening on first n levels

 58:   Logically Collective

 60:   Input Parameters:
 61: + pc - the preconditioner context
 62: - n  - 0, 1 or more

 64:   Options Database Key:
 65: . -pc_gamg_aggressive_coarsening <n,default = 1> - Number of levels to square the graph on before aggregating it

 67:   Level: intermediate

 69: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
 70: @*/
 71: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
 72: {
 73:   PetscFunctionBegin;
 76:   PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: /*@
 81:   PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (>2 is 'aggressive')

 83:   Logically Collective

 85:   Input Parameters:
 86: + pc - the preconditioner context
 87: - n  - 1 or more (default = 2)

 89:   Options Database Key:
 90: . -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive')

 92:   Level: intermediate

 94: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
 95: @*/
 96: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
 97: {
 98:   PetscFunctionBegin;
101:   PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: /*@
106:   PCGAMGSetAggressiveSquareGraph - Use graph square A'A for aggressive coarsening, old method

108:   Logically Collective

110:   Input Parameters:
111: + pc - the preconditioner context
112: - b  - default false - MIS-k is faster

114:   Options Database Key:
115: . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening

117:   Level: intermediate

119: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`
120: @*/
121: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
122: {
123:   PetscFunctionBegin;
126:   PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*@
131:   PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm

133:   Logically Collective

135:   Input Parameters:
136: + pc - the preconditioner context
137: - b  - default true

139:   Options Database Key:
140: . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm

142:   Level: intermediate

144: .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`
145: @*/
146: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
147: {
148:   PetscFunctionBegin;
151:   PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
156: {
157:   PC_MG       *mg          = (PC_MG *)pc->data;
158:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
159:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

161:   PetscFunctionBegin;
162:   pc_gamg_agg->aggressive_coarsening_levels = n;
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
167: {
168:   PC_MG       *mg          = (PC_MG *)pc->data;
169:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
170:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

172:   PetscFunctionBegin;
173:   pc_gamg_agg->aggressive_mis_k = n;
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
178: {
179:   PC_MG       *mg          = (PC_MG *)pc->data;
180:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
181:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

183:   PetscFunctionBegin;
184:   pc_gamg_agg->use_aggressive_square_graph = b;
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

188: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
189: {
190:   PC_MG       *mg          = (PC_MG *)pc->data;
191:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
192:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

194:   PetscFunctionBegin;
195:   pc_gamg_agg->use_minimum_degree_ordering = b;
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems *PetscOptionsObject)
200: {
201:   PC_MG       *mg          = (PC_MG *)pc->data;
202:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
203:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

205:   PetscFunctionBegin;
206:   PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
207:   {
208:     PetscBool flg;
209:     PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "smoothing steps for smoothed aggregation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL));
210:     PetscCall(PetscOptionsInt("-pc_gamg_aggressive_coarsening", "Number of aggressive coarsening (MIS-2) levels from finest", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &flg));
211:     if (!flg) {
212:       PetscCall(
213:         PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (deprecated alias for -pc_gamg_aggressive_coarsening)", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, NULL));
214:     } else {
215:       PetscCall(
216:         PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (alias for -pc_gamg_aggressive_coarsening, deprecated)", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &flg));
217:       if (flg) PetscCall(PetscInfo(pc, "Warning: both -pc_gamg_square_graph and -pc_gamg_aggressive_coarsening are used. -pc_gamg_square_graph is deprecated, Number of aggressive levels is %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels));
218:     }
219:     if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
220:       PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening", "PCGAMGSetAggressiveSquareGraph", pc_gamg_agg->use_aggressive_square_graph, &pc_gamg_agg->use_aggressive_square_graph, NULL));
221:     }
222:     PetscCall(PetscOptionsBool("-pc_gamg_mis_k_minimum_degree_ordering", "Use minimum degree ordering for greedy MIS", "PCGAMGMISkSetMinDegreeOrdering", pc_gamg_agg->use_minimum_degree_ordering, &pc_gamg_agg->use_minimum_degree_ordering, NULL));
223:     PetscCall(PetscOptionsInt("-pc_gamg_aggressive_mis_k", "Number of levels of multigrid to use.", "PCGAMGMISkSetAggressive", pc_gamg_agg->aggressive_mis_k, &pc_gamg_agg->aggressive_mis_k, NULL));
224:   }
225:   PetscOptionsHeadEnd();
226:   PetscFunctionReturn(PETSC_SUCCESS);
227: }

229: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
230: {
231:   PC_MG   *mg      = (PC_MG *)pc->data;
232:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

234:   PetscFunctionBegin;
235:   PetscCall(PetscFree(pc_gamg->subctx));
236:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
237:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
238:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
239:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
240:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
241:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
242:   PetscFunctionReturn(PETSC_SUCCESS);
243: }

245: /*
246:    PCSetCoordinates_AGG

248:    Collective

250:    Input Parameter:
251:    . pc - the preconditioner context
252:    . ndm - dimension of data (used for dof/vertex for Stokes)
253:    . a_nloc - number of vertices local
254:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
255: */

257: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
258: {
259:   PC_MG   *mg      = (PC_MG *)pc->data;
260:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
261:   PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
262:   Mat      mat = pc->pmat;

264:   PetscFunctionBegin;
267:   nloc = a_nloc;

269:   /* SA: null space vectors */
270:   PetscCall(MatGetBlockSize(mat, &ndf));               /* this does not work for Stokes */
271:   if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
272:   else if (coords) {
273:     PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
274:     pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
275:     if (ndm != ndf) PetscCheck(pc_gamg->data_cell_cols == ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Don't know how to create null space for ndm=%" PetscInt_FMT ", ndf=%" PetscInt_FMT ".  Use MatSetNearNullSpace().", ndm, ndf);
276:   } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
277:   pc_gamg->data_cell_rows = ndatarows = ndf;
278:   PetscCheck(pc_gamg->data_cell_cols > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "pc_gamg->data_cell_cols %" PetscInt_FMT " <= 0", pc_gamg->data_cell_cols);
279:   arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;

281:   if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
282:     PetscCall(PetscFree(pc_gamg->data));
283:     PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
284:   }
285:   /* copy data in - column oriented */
286:   for (kk = 0; kk < nloc; kk++) {
287:     const PetscInt M    = nloc * pc_gamg->data_cell_rows; /* stride into data */
288:     PetscReal     *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */
289:     if (pc_gamg->data_cell_cols == 1) *data = 1.0;
290:     else {
291:       /* translational modes */
292:       for (ii = 0; ii < ndatarows; ii++) {
293:         for (jj = 0; jj < ndatarows; jj++) {
294:           if (ii == jj) data[ii * M + jj] = 1.0;
295:           else data[ii * M + jj] = 0.0;
296:         }
297:       }

299:       /* rotational modes */
300:       if (coords) {
301:         if (ndm == 2) {
302:           data += 2 * M;
303:           data[0] = -coords[2 * kk + 1];
304:           data[1] = coords[2 * kk];
305:         } else {
306:           data += 3 * M;
307:           data[0]         = 0.0;
308:           data[M + 0]     = coords[3 * kk + 2];
309:           data[2 * M + 0] = -coords[3 * kk + 1];
310:           data[1]         = -coords[3 * kk + 2];
311:           data[M + 1]     = 0.0;
312:           data[2 * M + 1] = coords[3 * kk];
313:           data[2]         = coords[3 * kk + 1];
314:           data[M + 2]     = -coords[3 * kk];
315:           data[2 * M + 2] = 0.0;
316:         }
317:       }
318:     }
319:   }
320:   pc_gamg->data_sz = arrsz;
321:   PetscFunctionReturn(PETSC_SUCCESS);
322: }

324: /*
325:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
326:       Looks in Mat for near null space.
327:       Does not work for Stokes

329:   Input Parameter:
330:    . pc -
331:    . a_A - matrix to get (near) null space out of.
332: */
333: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
334: {
335:   PC_MG       *mg      = (PC_MG *)pc->data;
336:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
337:   MatNullSpace mnull;

339:   PetscFunctionBegin;
340:   PetscCall(MatGetNearNullSpace(a_A, &mnull));
341:   if (!mnull) {
342:     DM dm;
343:     PetscCall(PCGetDM(pc, &dm));
344:     if (!dm) PetscCall(MatGetDM(a_A, &dm));
345:     if (dm) {
346:       PetscObject deformation;
347:       PetscInt    Nf;

349:       PetscCall(DMGetNumFields(dm, &Nf));
350:       if (Nf) {
351:         PetscCall(DMGetField(dm, 0, NULL, &deformation));
352:         PetscCall(PetscObjectQuery((PetscObject)deformation, "nearnullspace", (PetscObject *)&mnull));
353:         if (!mnull) PetscCall(PetscObjectQuery((PetscObject)deformation, "nullspace", (PetscObject *)&mnull));
354:       }
355:     }
356:   }

358:   if (!mnull) {
359:     PetscInt bs, NN, MM;
360:     PetscCall(MatGetBlockSize(a_A, &bs));
361:     PetscCall(MatGetLocalSize(a_A, &MM, &NN));
362:     PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
363:     PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
364:   } else {
365:     PetscReal         *nullvec;
366:     PetscBool          has_const;
367:     PetscInt           i, j, mlocal, nvec, bs;
368:     const Vec         *vecs;
369:     const PetscScalar *v;

371:     PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
372:     PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
373:     for (i = 0; i < nvec; i++) {
374:       PetscCall(VecGetLocalSize(vecs[i], &j));
375:       PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
376:     }
377:     pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
378:     PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
379:     if (has_const)
380:       for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
381:     for (i = 0; i < nvec; i++) {
382:       PetscCall(VecGetArrayRead(vecs[i], &v));
383:       for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
384:       PetscCall(VecRestoreArrayRead(vecs[i], &v));
385:     }
386:     pc_gamg->data           = nullvec;
387:     pc_gamg->data_cell_cols = (nvec + !!has_const);
388:     PetscCall(MatGetBlockSize(a_A, &bs));
389:     pc_gamg->data_cell_rows = bs;
390:   }
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

394: /*
395:   formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0

397:   Input Parameter:
398:    . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
399:    . bs - row block size
400:    . nSAvec - column bs of new P
401:    . my0crs - global index of start of locals
402:    . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
403:    . data_in[data_stride*nSAvec] - local data on fine grid
404:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'

406:   Output Parameter:
407:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
408:    . a_Prol - prolongation operator
409: */
410: static PetscErrorCode formProl0(PetscCoarsenData *agg_llists, PetscInt bs, PetscInt nSAvec, PetscInt my0crs, PetscInt data_stride, PetscReal data_in[], const PetscInt flid_fgid[], PetscReal **a_data_out, Mat a_Prol)
411: {
412:   PetscInt        Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
413:   MPI_Comm        comm;
414:   PetscReal      *out_data;
415:   PetscCDIntNd   *pos;
416:   PCGAMGHashTable fgid_flid;

418:   PetscFunctionBegin;
419:   PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
420:   PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
421:   nloc = (Iend - Istart) / bs;
422:   my0  = Istart / bs;
423:   PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, Iend, Istart, bs);
424:   Iend /= bs;
425:   nghosts = data_stride / bs - nloc;

427:   PetscCall(PCGAMGHashTableCreate(2 * nghosts + 1, &fgid_flid));
428:   for (kk = 0; kk < nghosts; kk++) PetscCall(PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc + kk], nloc + kk));

430:   /* count selected -- same as number of cols of P */
431:   for (nSelected = mm = 0; mm < nloc; mm++) {
432:     PetscBool ise;
433:     PetscCall(PetscCDEmptyAt(agg_llists, mm, &ise));
434:     if (!ise) nSelected++;
435:   }
436:   PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
437:   PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT "  != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
438:   PetscCheck(nSelected == (jj - ii) / nSAvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nSelected %" PetscInt_FMT " != (jj %" PetscInt_FMT " - ii %" PetscInt_FMT ")/nSAvec %" PetscInt_FMT, nSelected, jj, ii, nSAvec);

440:   /* aloc space for coarse point data (output) */
441:   out_data_stride = nSelected * nSAvec;

443:   PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data));
444:   for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL;
445:   *a_data_out = out_data; /* output - stride nSelected*nSAvec */

447:   /* find points and set prolongation */
448:   minsz = 100;
449:   for (mm = clid = 0; mm < nloc; mm++) {
450:     PetscCall(PetscCDSizeAt(agg_llists, mm, &jj));
451:     if (jj > 0) {
452:       const PetscInt lid = mm, cgid = my0crs + clid;
453:       PetscInt       cids[100]; /* max bs */
454:       PetscBLASInt   asz = jj, M = asz * bs, N = nSAvec, INFO;
455:       PetscBLASInt   Mdata = M + ((N - M > 0) ? N - M : 0), LDA = Mdata, LWORK = N * bs;
456:       PetscScalar   *qqc, *qqr, *TAU, *WORK;
457:       PetscInt      *fids;
458:       PetscReal     *data;

460:       /* count agg */
461:       if (asz < minsz) minsz = asz;

463:       /* get block */
464:       PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids));

466:       aggID = 0;
467:       PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
468:       while (pos) {
469:         PetscInt gid1;
470:         PetscCall(PetscCDIntNdGetID(pos, &gid1));
471:         PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));

473:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
474:         else {
475:           PetscCall(PCGAMGHashTableFind(&fgid_flid, gid1, &flid));
476:           PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
477:         }
478:         /* copy in B_i matrix - column oriented */
479:         data = &data_in[flid * bs];
480:         for (ii = 0; ii < bs; ii++) {
481:           for (jj = 0; jj < N; jj++) {
482:             PetscReal d                       = data[jj * data_stride + ii];
483:             qqc[jj * Mdata + aggID * bs + ii] = d;
484:           }
485:         }
486:         /* set fine IDs */
487:         for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
488:         aggID++;
489:       }

491:       /* pad with zeros */
492:       for (ii = asz * bs; ii < Mdata; ii++) {
493:         for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
494:       }

496:       /* QR */
497:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
498:       PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
499:       PetscCall(PetscFPTrapPop());
500:       PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
501:       /* get R - column oriented - output B_{i+1} */
502:       {
503:         PetscReal *data = &out_data[clid * nSAvec];
504:         for (jj = 0; jj < nSAvec; jj++) {
505:           for (ii = 0; ii < nSAvec; ii++) {
506:             PetscCheck(data[jj * out_data_stride + ii] == PETSC_MAX_REAL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "data[jj*out_data_stride + ii] != %e", (double)PETSC_MAX_REAL);
507:             if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
508:             else data[jj * out_data_stride + ii] = 0.;
509:           }
510:         }
511:       }

513:       /* get Q - row oriented */
514:       PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
515:       PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xORGQR error arg %" PetscBLASInt_FMT, -INFO);

517:       for (ii = 0; ii < M; ii++) {
518:         for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
519:       }

521:       /* add diagonal block of P0 */
522:       for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ }
523:       PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
524:       PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
525:       clid++;
526:     } /* coarse agg */
527:   }   /* for all fine nodes */
528:   PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
529:   PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
530:   PetscCall(PCGAMGHashTableDestroy(&fgid_flid));
531:   PetscFunctionReturn(PETSC_SUCCESS);
532: }

534: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
535: {
536:   PC_MG       *mg          = (PC_MG *)pc->data;
537:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
538:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

540:   PetscFunctionBegin;
541:   PetscCall(PetscViewerASCIIPrintf(viewer, "      AGG specific options\n"));
542:   PetscCall(PetscViewerASCIIPrintf(viewer, "        Number of levels of aggressive coarsening %d\n", (int)pc_gamg_agg->aggressive_coarsening_levels));
543:   if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
544:     PetscCall(PetscViewerASCIIPrintf(viewer, "        %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
545:     if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(PetscViewerASCIIPrintf(viewer, "        MIS-%d coarsening on aggressive levels\n", (int)pc_gamg_agg->aggressive_mis_k));
546:   }
547:   PetscCall(PetscViewerASCIIPrintf(viewer, "        Number smoothing steps %d\n", (int)pc_gamg_agg->nsmooths));
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

551: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
552: {
553:   PC_MG          *mg          = (PC_MG *)pc->data;
554:   PC_GAMG        *pc_gamg     = (PC_GAMG *)mg->innerctx;
555:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
556:   const PetscReal vfilter     = pc_gamg->threshold[pc_gamg->current_level];
557:   PetscBool       ishem;
558:   const char     *prefix;
559:   MatInfo         info0, info1;
560:   PetscInt        bs;

562:   PetscFunctionBegin;
563:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
564:   /* Note: depending on the algorithm that will be used for computing the coarse grid points this should pass PETSC_TRUE or PETSC_FALSE as the first argument */

566:   /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
567:   PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
568:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
569:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
570:   PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
571:   PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
572:   if (ishem) pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense
573:   PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0));           /* global reduction */
574:   PetscCall(MatCreateGraph(Amat, PETSC_TRUE, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, a_Gmat));
575:   PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
576:   PetscCall(MatGetBlockSize(Amat, &bs));
577:   if (info0.nz_used > 0) PetscCall(PetscInfo(pc, "Filtering left %g %% edges in graph (%e %e)\n", 100.0 * info1.nz_used * (double)(bs * bs) / info0.nz_used, info0.nz_used, info1.nz_used));
578:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
579:   PetscFunctionReturn(PETSC_SUCCESS);
580: }

582: typedef PetscInt    NState;
583: static const NState NOT_DONE = -2;
584: static const NState DELETED  = -1;
585: static const NState REMOVED  = -3;
586: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)

588: /*
589:    fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD
590:      - AGG-MG specific: clears singletons out of 'selected_2'

592:    Input Parameter:
593:    . Gmat_2 - global matrix of squared graph (data not defined)
594:    . Gmat_1 - base graph to grab with base graph
595:    Input/Output Parameter:
596:    . aggs_2 - linked list of aggs with gids)
597: */
598: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
599: {
600:   PetscBool      isMPI;
601:   Mat_SeqAIJ    *matA_1, *matB_1 = NULL;
602:   MPI_Comm       comm;
603:   PetscInt       lid, *ii, *idx, ix, Iend, my0, kk, n, j;
604:   Mat_MPIAIJ    *mpimat_2 = NULL, *mpimat_1 = NULL;
605:   const PetscInt nloc = Gmat_2->rmap->n;
606:   PetscScalar   *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
607:   PetscInt      *lid_cprowID_1 = NULL;
608:   NState        *lid_state;
609:   Vec            ghost_par_orig2;
610:   PetscMPIInt    rank;

612:   PetscFunctionBegin;
613:   PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
614:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
615:   PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));

617:   /* get submatrices */
618:   PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
619:   PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
620:   PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
621:   for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
622:   if (isMPI) {
623:     /* grab matrix objects */
624:     mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
625:     mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
626:     matA_1   = (Mat_SeqAIJ *)mpimat_1->A->data;
627:     matB_1   = (Mat_SeqAIJ *)mpimat_1->B->data;

629:     /* force compressed row storage for B matrix in AuxMat */
630:     PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
631:     for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
632:       PetscInt lid = matB_1->compressedrow.rindex[ix];
633:       PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %d out of range. nloc = %d", (int)lid, (int)nloc);
634:       if (lid != -1) lid_cprowID_1[lid] = ix;
635:     }
636:   } else {
637:     PetscBool isAIJ;
638:     PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
639:     PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
640:     matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
641:   }
642:   if (nloc > 0) { PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???"); }
643:   /* get state of locals and selected gid for deleted */
644:   for (lid = 0; lid < nloc; lid++) {
645:     lid_parent_gid[lid] = -1.0;
646:     lid_state[lid]      = DELETED;
647:   }

649:   /* set lid_state */
650:   for (lid = 0; lid < nloc; lid++) {
651:     PetscCDIntNd *pos;
652:     PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
653:     if (pos) {
654:       PetscInt gid1;

656:       PetscCall(PetscCDIntNdGetID(pos, &gid1));
657:       PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %d != lid %d + my0 %d", (int)gid1, (int)lid, (int)my0);
658:       lid_state[lid] = gid1;
659:     }
660:   }

662:   /* map local to selected local, DELETED means a ghost owns it */
663:   for (lid = kk = 0; lid < nloc; lid++) {
664:     NState state = lid_state[lid];
665:     if (IS_SELECTED(state)) {
666:       PetscCDIntNd *pos;
667:       PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
668:       while (pos) {
669:         PetscInt gid1;
670:         PetscCall(PetscCDIntNdGetID(pos, &gid1));
671:         PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
672:         if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
673:       }
674:     }
675:   }
676:   /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
677:   if (isMPI) {
678:     Vec tempVec;
679:     /* get 'cpcol_1_state' */
680:     PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
681:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
682:       PetscScalar v = (PetscScalar)lid_state[kk];
683:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
684:     }
685:     PetscCall(VecAssemblyBegin(tempVec));
686:     PetscCall(VecAssemblyEnd(tempVec));
687:     PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
688:     PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
689:     PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
690:     /* get 'cpcol_2_state' */
691:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
692:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
693:     PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
694:     /* get 'cpcol_2_par_orig' */
695:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
696:       PetscScalar v = (PetscScalar)lid_parent_gid[kk];
697:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
698:     }
699:     PetscCall(VecAssemblyBegin(tempVec));
700:     PetscCall(VecAssemblyEnd(tempVec));
701:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
702:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
703:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
704:     PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));

706:     PetscCall(VecDestroy(&tempVec));
707:   } /* ismpi */
708:   for (lid = 0; lid < nloc; lid++) {
709:     NState state = lid_state[lid];
710:     if (IS_SELECTED(state)) {
711:       /* steal locals */
712:       ii  = matA_1->i;
713:       n   = ii[lid + 1] - ii[lid];
714:       idx = matA_1->j + ii[lid];
715:       for (j = 0; j < n; j++) {
716:         PetscInt lidj   = idx[j], sgid;
717:         NState   statej = lid_state[lidj];
718:         if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
719:           lid_parent_gid[lidj] = (PetscScalar)(lid + my0);                                              /* send this if sgid is not local */
720:           if (sgid >= my0 && sgid < Iend) {                                                             /* I'm stealing this local from a local sgid */
721:             PetscInt      hav = 0, slid = sgid - my0, gidj = lidj + my0;
722:             PetscCDIntNd *pos, *last = NULL;
723:             /* looking for local from local so id_llist_2 works */
724:             PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
725:             while (pos) {
726:               PetscInt gid;
727:               PetscCall(PetscCDIntNdGetID(pos, &gid));
728:               if (gid == gidj) {
729:                 PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
730:                 PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
731:                 PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
732:                 hav = 1;
733:                 break;
734:               } else last = pos;
735:               PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
736:             }
737:             if (hav != 1) {
738:               PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
739:               SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %d times???", (int)hav);
740:             }
741:           } else { /* I'm stealing this local, owned by a ghost */
742:             PetscCheck(sgid == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mat has an un-symmetric graph. Use '-%spc_gamg_sym_graph true' to symmetrize the graph or '-%spc_gamg_threshold -1' if the matrix is structurally symmetric.",
743:                        ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
744:             PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
745:           }
746:         }
747:       } /* local neighbors */
748:     } else if (state == DELETED /* && lid_cprowID_1 */) {
749:       PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
750:       /* see if I have a selected ghost neighbor that will steal me */
751:       if ((ix = lid_cprowID_1[lid]) != -1) {
752:         ii  = matB_1->compressedrow.i;
753:         n   = ii[ix + 1] - ii[ix];
754:         idx = matB_1->j + ii[ix];
755:         for (j = 0; j < n; j++) {
756:           PetscInt cpid   = idx[j];
757:           NState   statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
758:           if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
759:             lid_parent_gid[lid] = (PetscScalar)statej;              /* send who selected */
760:             if (sgidold >= my0 && sgidold < Iend) {                 /* this was mine */
761:               PetscInt      hav = 0, oldslidj = sgidold - my0;
762:               PetscCDIntNd *pos, *last        = NULL;
763:               /* remove from 'oldslidj' list */
764:               PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
765:               while (pos) {
766:                 PetscInt gid;
767:                 PetscCall(PetscCDIntNdGetID(pos, &gid));
768:                 if (lid + my0 == gid) {
769:                   /* id_llist_2[lastid] = id_llist_2[flid];   /\* remove lid from oldslidj list *\/ */
770:                   PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
771:                   PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
772:                   /* ghost (PetscScalar)statej will add this later */
773:                   hav = 1;
774:                   break;
775:                 } else last = pos;
776:                 PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
777:               }
778:               if (hav != 1) {
779:                 PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%d) adj in 'selected' lists - structurally unsymmetric matrix", (int)hav);
780:                 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %d times???", (int)hav);
781:               }
782:             } else {
783:               /* TODO: ghosts remove this later */
784:             }
785:           }
786:         }
787:       }
788:     } /* selected/deleted */
789:   }   /* node loop */

791:   if (isMPI) {
792:     PetscScalar    *cpcol_2_parent, *cpcol_2_gid;
793:     Vec             tempVec, ghostgids2, ghostparents2;
794:     PetscInt        cpid, nghost_2;
795:     PCGAMGHashTable gid_cpid;

797:     PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2));
798:     PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL));

800:     /* get 'cpcol_2_parent' */
801:     for (kk = 0, j = my0; kk < nloc; kk++, j++) { PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES)); }
802:     PetscCall(VecAssemblyBegin(tempVec));
803:     PetscCall(VecAssemblyEnd(tempVec));
804:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
805:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
806:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
807:     PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));

809:     /* get 'cpcol_2_gid' */
810:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
811:       PetscScalar v = (PetscScalar)j;
812:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
813:     }
814:     PetscCall(VecAssemblyBegin(tempVec));
815:     PetscCall(VecAssemblyEnd(tempVec));
816:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
817:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
818:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
819:     PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
820:     PetscCall(VecDestroy(&tempVec));

822:     /* look for deleted ghosts and add to table */
823:     PetscCall(PCGAMGHashTableCreate(2 * nghost_2 + 1, &gid_cpid));
824:     for (cpid = 0; cpid < nghost_2; cpid++) {
825:       NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
826:       if (state == DELETED) {
827:         PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
828:         PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
829:         if (sgid_old == -1 && sgid_new != -1) {
830:           PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
831:           PetscCall(PCGAMGHashTableAdd(&gid_cpid, gid, cpid));
832:         }
833:       }
834:     }

836:     /* look for deleted ghosts and see if they moved - remove it */
837:     for (lid = 0; lid < nloc; lid++) {
838:       NState state = lid_state[lid];
839:       if (IS_SELECTED(state)) {
840:         PetscCDIntNd *pos, *last = NULL;
841:         /* look for deleted ghosts and see if they moved */
842:         PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
843:         while (pos) {
844:           PetscInt gid;
845:           PetscCall(PetscCDIntNdGetID(pos, &gid));

847:           if (gid < my0 || gid >= Iend) {
848:             PetscCall(PCGAMGHashTableFind(&gid_cpid, gid, &cpid));
849:             if (cpid != -1) {
850:               /* a moved ghost - */
851:               /* id_llist_2[lastid] = id_llist_2[flid];    /\* remove 'flid' from list *\/ */
852:               PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
853:             } else last = pos;
854:           } else last = pos;

856:           PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
857:         } /* loop over list of deleted */
858:       }   /* selected */
859:     }
860:     PetscCall(PCGAMGHashTableDestroy(&gid_cpid));

862:     /* look at ghosts, see if they changed - and it */
863:     for (cpid = 0; cpid < nghost_2; cpid++) {
864:       PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
865:       if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
866:         PetscInt      gid      = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
867:         PetscInt      slid_new = sgid_new - my0, hav = 0;
868:         PetscCDIntNd *pos;

870:         /* search for this gid to see if I have it */
871:         PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
872:         while (pos) {
873:           PetscInt gidj;
874:           PetscCall(PetscCDIntNdGetID(pos, &gidj));
875:           PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));

877:           if (gidj == gid) {
878:             hav = 1;
879:             break;
880:           }
881:         }
882:         if (hav != 1) {
883:           /* insert 'flidj' into head of llist */
884:           PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
885:         }
886:       }
887:     }
888:     PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
889:     PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
890:     PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
891:     PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
892:     PetscCall(VecDestroy(&ghostgids2));
893:     PetscCall(VecDestroy(&ghostparents2));
894:     PetscCall(VecDestroy(&ghost_par_orig2));
895:   }
896:   PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
897:   PetscFunctionReturn(PETSC_SUCCESS);
898: }

900: /*
901:    PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for
902:      communication of QR data used with HEM and MISk coarsening

904:   Input Parameter:
905:    . a_pc - this

907:   Input/Output Parameter:
908:    . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out)

910:   Output Parameter:
911:    . agg_lists - list of aggregates

913: */
914: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
915: {
916:   PC_MG       *mg          = (PC_MG *)a_pc->data;
917:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
918:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
919:   Mat          mat, Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
920:   IS           perm;
921:   PetscInt     Istart, Iend, Ii, nloc, bs, nn;
922:   PetscInt    *permute, *degree;
923:   PetscBool   *bIndexSet;
924:   PetscReal    hashfact;
925:   PetscInt     iSwapIndex;
926:   PetscRandom  random;
927:   MPI_Comm     comm;

929:   PetscFunctionBegin;
930:   PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
931:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
932:   PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
933:   PetscCall(MatGetBlockSize(Gmat1, &bs));
934:   PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
935:   nloc = nn / bs;
936:   /* get MIS aggs - randomize */
937:   PetscCall(PetscMalloc2(nloc, &permute, nloc, &degree));
938:   PetscCall(PetscCalloc1(nloc, &bIndexSet));
939:   for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
940:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
941:   PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
942:   for (Ii = 0; Ii < nloc; Ii++) {
943:     PetscInt nc;
944:     PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
945:     degree[Ii] = nc;
946:     PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
947:   }
948:   for (Ii = 0; Ii < nloc; Ii++) {
949:     PetscCall(PetscRandomGetValueReal(random, &hashfact));
950:     iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
951:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
952:       PetscInt iTemp        = permute[iSwapIndex];
953:       permute[iSwapIndex]   = permute[Ii];
954:       permute[Ii]           = iTemp;
955:       iTemp                 = degree[iSwapIndex];
956:       degree[iSwapIndex]    = degree[Ii];
957:       degree[Ii]            = iTemp;
958:       bIndexSet[iSwapIndex] = PETSC_TRUE;
959:     }
960:   }
961:   // apply minimum degree ordering -- NEW
962:   if (pc_gamg_agg->use_minimum_degree_ordering) { PetscCall(PetscSortIntWithArray(nloc, degree, permute)); }
963:   PetscCall(PetscFree(bIndexSet));
964:   PetscCall(PetscRandomDestroy(&random));
965:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
966:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
967:   // square graph
968:   if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels && pc_gamg_agg->use_aggressive_square_graph) {
969:     PetscCall(PCGAMGSquareGraph_GAMG(a_pc, Gmat1, &Gmat2));
970:   } else Gmat2 = Gmat1;
971:   // switch to old MIS-1 for square graph
972:   if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
973:     if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
974:     else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS));                                                                   // old MIS -- side effect
975:   } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) {                                 // we reset the MIS
976:     const char *prefix;
977:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
978:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
979:     PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
980:   }
981:   PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
982:   PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
983:   PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
984:   PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
985:   PetscCall(MatCoarsenViewFromOptions(pc_gamg_agg->crs, NULL, "-mat_coarsen_view"));
986:   PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */
987:   PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));

989:   PetscCall(ISDestroy(&perm));
990:   PetscCall(PetscFree2(permute, degree));
991:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));

993:   if (Gmat2 != Gmat1) {
994:     PetscCoarsenData *llist = *agg_lists;
995:     PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
996:     PetscCall(MatDestroy(&Gmat1));
997:     *a_Gmat1 = Gmat2; /* output */
998:     PetscCall(PetscCDGetMat(llist, &mat));
999:     PetscCheck(!mat, comm, PETSC_ERR_ARG_WRONG, "Unexpected auxiliary matrix with squared graph");
1000:   } else {
1001:     PetscCoarsenData *llist = *agg_lists;
1002:     /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */
1003:     PetscCall(PetscCDGetMat(llist, &mat));
1004:     if (mat) {
1005:       PetscCall(MatDestroy(a_Gmat1));
1006:       *a_Gmat1 = mat; /* output */
1007:     }
1008:   }
1009:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1010:   PetscFunctionReturn(PETSC_SUCCESS);
1011: }

1013: /*
1014:  PCGAMGProlongator_AGG

1016:  Input Parameter:
1017:  . pc - this
1018:  . Amat - matrix on this fine level
1019:  . Graph - used to get ghost data for nodes in
1020:  . agg_lists - list of aggregates
1021:  Output Parameter:
1022:  . a_P_out - prolongation operator to the next level
1023:  */
1024: static PetscErrorCode PCGAMGProlongator_AGG(PC pc, Mat Amat, Mat Gmat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1025: {
1026:   PC_MG         *mg      = (PC_MG *)pc->data;
1027:   PC_GAMG       *pc_gamg = (PC_GAMG *)mg->innerctx;
1028:   const PetscInt col_bs  = pc_gamg->data_cell_cols;
1029:   PetscInt       Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1030:   Mat            Prol;
1031:   PetscMPIInt    size;
1032:   MPI_Comm       comm;
1033:   PetscReal     *data_w_ghost;
1034:   PetscInt       myCrs0, nbnodes = 0, *flid_fgid;
1035:   MatType        mtype;

1037:   PetscFunctionBegin;
1038:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1039:   PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1040:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1041:   PetscCallMPI(MPI_Comm_size(comm, &size));
1042:   PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1043:   PetscCall(MatGetBlockSize(Amat, &bs));
1044:   nloc = (Iend - Istart) / bs;
1045:   my0  = Istart / bs;
1046:   PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") not divisible by bs %" PetscInt_FMT, Iend, Istart, bs);

1048:   /* get 'nLocalSelected' */
1049:   for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1050:     PetscBool ise;
1051:     /* filter out singletons 0 or 1? */
1052:     PetscCall(PetscCDEmptyAt(agg_lists, ii, &ise));
1053:     if (!ise) nLocalSelected++;
1054:   }

1056:   /* create prolongator, create P matrix */
1057:   PetscCall(MatGetType(Amat, &mtype));
1058:   PetscCall(MatCreate(comm, &Prol));
1059:   PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1060:   PetscCall(MatSetBlockSizes(Prol, bs, col_bs));
1061:   PetscCall(MatSetType(Prol, mtype));
1062: #if PetscDefined(HAVE_DEVICE)
1063:   PetscBool flg;
1064:   PetscCall(MatBoundToCPU(Amat, &flg));
1065:   PetscCall(MatBindToCPU(Prol, flg));
1066:   if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1067: #endif
1068:   PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1069:   PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));

1071:   /* can get all points "removed" */
1072:   PetscCall(MatGetSize(Prol, &kk, &ii));
1073:   if (!ii) {
1074:     PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1075:     PetscCall(MatDestroy(&Prol));
1076:     *a_P_out = NULL; /* out */
1077:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1078:     PetscFunctionReturn(PETSC_SUCCESS);
1079:   }
1080:   PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1081:   PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));

1083:   PetscCheck((kk - myCrs0) % col_bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(kk %" PetscInt_FMT " -myCrs0 %" PetscInt_FMT ") not divisible by col_bs %" PetscInt_FMT, kk, myCrs0, col_bs);
1084:   myCrs0 = myCrs0 / col_bs;
1085:   PetscCheck((kk / col_bs - myCrs0) == nLocalSelected, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(kk %" PetscInt_FMT "/col_bs %" PetscInt_FMT " - myCrs0 %" PetscInt_FMT ") != nLocalSelected %" PetscInt_FMT ")", kk, col_bs, myCrs0, nLocalSelected);

1087:   /* create global vector of data in 'data_w_ghost' */
1088:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1089:   if (size > 1) { /* get ghost null space data */
1090:     PetscReal *tmp_gdata, *tmp_ldata, *tp2;
1091:     PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1092:     for (jj = 0; jj < col_bs; jj++) {
1093:       for (kk = 0; kk < bs; kk++) {
1094:         PetscInt         ii, stride;
1095:         const PetscReal *tp = pc_gamg->data + jj * bs * nloc + kk;
1096:         for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;

1098:         PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata));

1100:         if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1101:           PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1102:           nbnodes = bs * stride;
1103:         }
1104:         tp2 = data_w_ghost + jj * bs * stride + kk;
1105:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1106:         PetscCall(PetscFree(tmp_gdata));
1107:       }
1108:     }
1109:     PetscCall(PetscFree(tmp_ldata));
1110:   } else {
1111:     nbnodes      = bs * nloc;
1112:     data_w_ghost = (PetscReal *)pc_gamg->data;
1113:   }

1115:   /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */
1116:   if (size > 1) {
1117:     PetscReal *fid_glid_loc, *fiddata;
1118:     PetscInt   stride;

1120:     PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1121:     for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1122:     PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1123:     PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1124:     for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1125:     PetscCall(PetscFree(fiddata));

1127:     PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1128:     PetscCall(PetscFree(fid_glid_loc));
1129:   } else {
1130:     PetscCall(PetscMalloc1(nloc, &flid_fgid));
1131:     for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1132:   }
1133:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1134:   /* get P0 */
1135:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1136:   {
1137:     PetscReal *data_out = NULL;
1138:     PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1139:     PetscCall(PetscFree(pc_gamg->data));

1141:     pc_gamg->data           = data_out;
1142:     pc_gamg->data_cell_rows = col_bs;
1143:     pc_gamg->data_sz        = col_bs * col_bs * nLocalSelected;
1144:   }
1145:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1146:   if (size > 1) PetscCall(PetscFree(data_w_ghost));
1147:   PetscCall(PetscFree(flid_fgid));

1149:   *a_P_out = Prol; /* out */

1151:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1152:   PetscFunctionReturn(PETSC_SUCCESS);
1153: }

1155: /*
1156:    PCGAMGOptProlongator_AGG

1158:   Input Parameter:
1159:    . pc - this
1160:    . Amat - matrix on this fine level
1161:  In/Output Parameter:
1162:    . a_P - prolongation operator to the next level
1163: */
1164: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
1165: {
1166:   PC_MG       *mg          = (PC_MG *)pc->data;
1167:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
1168:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1169:   PetscInt     jj;
1170:   Mat          Prol = *a_P;
1171:   MPI_Comm     comm;
1172:   KSP          eksp;
1173:   Vec          bb, xx;
1174:   PC           epc;
1175:   PetscReal    alpha, emax, emin;

1177:   PetscFunctionBegin;
1178:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1179:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));

1181:   /* compute maximum singular value of operator to be used in smoother */
1182:   if (0 < pc_gamg_agg->nsmooths) {
1183:     /* get eigen estimates */
1184:     if (pc_gamg->emax > 0) {
1185:       emin = pc_gamg->emin;
1186:       emax = pc_gamg->emax;
1187:     } else {
1188:       const char *prefix;

1190:       PetscCall(MatCreateVecs(Amat, &bb, NULL));
1191:       PetscCall(MatCreateVecs(Amat, &xx, NULL));
1192:       PetscCall(KSPSetNoisy_Private(bb));

1194:       PetscCall(KSPCreate(comm, &eksp));
1195:       PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
1196:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
1197:       PetscCall(KSPSetOptionsPrefix(eksp, prefix));
1198:       PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
1199:       {
1200:         PetscBool isset, sflg;
1201:         PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
1202:         if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
1203:       }
1204:       PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
1205:       PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));

1207:       PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
1208:       PetscCall(KSPSetOperators(eksp, Amat, Amat));

1210:       PetscCall(KSPGetPC(eksp, &epc));
1211:       PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */

1213:       PetscCall(KSPSetTolerances(eksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 10)); // 10 is safer, but 5 is often fine, can override with -pc_gamg_esteig_ksp_max_it -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.2

1215:       PetscCall(KSPSetFromOptions(eksp));
1216:       PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
1217:       PetscCall(KSPSolve(eksp, bb, xx));
1218:       PetscCall(KSPCheckSolve(eksp, pc, xx));

1220:       PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
1221:       PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
1222:       PetscCall(VecDestroy(&xx));
1223:       PetscCall(VecDestroy(&bb));
1224:       PetscCall(KSPDestroy(&eksp));
1225:     }
1226:     if (pc_gamg->use_sa_esteig) {
1227:       mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
1228:       mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
1229:       PetscCall(PetscInfo(pc, "%s: Smooth P0: level %" PetscInt_FMT ", cache spectra %g %g\n", ((PetscObject)pc)->prefix, pc_gamg->current_level, (double)emin, (double)emax));
1230:     } else {
1231:       mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1232:       mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1233:     }
1234:   } else {
1235:     mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1236:     mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1237:   }

1239:   /* smooth P0 */
1240:   for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1241:     Mat tMat;
1242:     Vec diag;

1244:     PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));

1246:     /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1247:     PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1248:     PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat));
1249:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1250:     PetscCall(MatProductClear(tMat));
1251:     PetscCall(MatCreateVecs(Amat, &diag, NULL));
1252:     PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
1253:     PetscCall(VecReciprocal(diag));
1254:     PetscCall(MatDiagonalScale(tMat, diag, NULL));
1255:     PetscCall(VecDestroy(&diag));

1257:     /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
1258:     PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");
1259:     /* TODO: Document the 1.4 and don't hardwire it in this routine */
1260:     alpha = -1.4 / emax;

1262:     PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
1263:     PetscCall(MatDestroy(&Prol));
1264:     Prol = tMat;
1265:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1266:   }
1267:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1268:   *a_P = Prol;
1269:   PetscFunctionReturn(PETSC_SUCCESS);
1270: }

1272: /*
1273:    PCCreateGAMG_AGG

1275:   Input Parameter:
1276:    . pc -
1277: */
1278: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1279: {
1280:   PC_MG       *mg      = (PC_MG *)pc->data;
1281:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
1282:   PC_GAMG_AGG *pc_gamg_agg;

1284:   PetscFunctionBegin;
1285:   /* create sub context for SA */
1286:   PetscCall(PetscNew(&pc_gamg_agg));
1287:   pc_gamg->subctx = pc_gamg_agg;

1289:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1290:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
1291:   /* reset does not do anything; setup not virtual */

1293:   /* set internal function pointers */
1294:   pc_gamg->ops->creategraph       = PCGAMGCreateGraph_AGG;
1295:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_AGG;
1296:   pc_gamg->ops->prolongator       = PCGAMGProlongator_AGG;
1297:   pc_gamg->ops->optprolongator    = PCGAMGOptProlongator_AGG;
1298:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1299:   pc_gamg->ops->view              = PCView_GAMG_AGG;

1301:   pc_gamg_agg->nsmooths                     = 1;
1302:   pc_gamg_agg->aggressive_coarsening_levels = 1;
1303:   pc_gamg_agg->use_aggressive_square_graph  = PETSC_FALSE;
1304:   pc_gamg_agg->use_minimum_degree_ordering  = PETSC_FALSE;
1305:   pc_gamg_agg->aggressive_mis_k             = 2;

1307:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
1308:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
1309:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
1310:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
1311:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
1312:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
1313:   PetscFunctionReturn(PETSC_SUCCESS);
1314: }