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;                     // number of smoothing steps to construct prolongation
 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:   PetscBool  use_low_mem_filter;
 17:   PetscBool  graph_symmetrize;
 18:   MatCoarsen crs;
 19: } PC_GAMG_AGG;

 21: /*@
 22:   PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used to construct the prolongation operator

 24:   Logically Collective

 26:   Input Parameters:
 27: + pc - the preconditioner context
 28: - n  - the number of smooths

 30:   Options Database Key:
 31: . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use

 33:   Level: intermediate

 35:   Note:
 36:   This is a different concept from the number smoothing steps used during the linear solution process which
 37:   can be set with `-mg_levels_ksp_max_it`

 39:   Developer Note:
 40:   This should be named `PCGAMGAGGSetNSmooths()`.

 42: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCMG`, `PCGAMG`
 43: @*/
 44: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
 45: {
 46:   PetscFunctionBegin;
 49:   PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n));
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

 53: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
 54: {
 55:   PC_MG       *mg          = (PC_MG *)pc->data;
 56:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
 57:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

 59:   PetscFunctionBegin;
 60:   pc_gamg_agg->nsmooths = n;
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

 64: /*@
 65:   PCGAMGSetAggressiveLevels -  Use aggressive coarsening on first n levels

 67:   Logically Collective

 69:   Input Parameters:
 70: + pc - the preconditioner context
 71: - n  - 0, 1 or more

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

 76:   Level: intermediate

 78: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
 79: @*/
 80: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
 81: {
 82:   PetscFunctionBegin;
 85:   PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

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

 92:   Logically Collective

 94:   Input Parameters:
 95: + pc - the preconditioner context
 96: - n  - 1 or more (default = 2)

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

101:   Level: intermediate

103: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
104: @*/
105: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
106: {
107:   PetscFunctionBegin;
110:   PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

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

117:   Logically Collective

119:   Input Parameters:
120: + pc - the preconditioner context
121: - b  - default false - MIS-k is faster

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

126:   Level: intermediate

128: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()`
129: @*/
130: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
131: {
132:   PetscFunctionBegin;
135:   PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
136:   PetscFunctionReturn(PETSC_SUCCESS);
137: }

139: /*@
140:   PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm

142:   Logically Collective

144:   Input Parameters:
145: + pc - the preconditioner context
146: - b  - default true

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

151:   Level: intermediate

153: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGSetLowMemoryFilter()`
154: @*/
155: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
156: {
157:   PetscFunctionBegin;
160:   PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: /*@
165:   PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter

167:   Logically Collective

169:   Input Parameters:
170: + pc - the preconditioner context
171: - b  - default false

173:   Options Database Key:
174: . -pc_gamg_low_memory_threshold_filter <bool,default=false> - Use low memory graph/matrix filter

176:   Level: intermediate

178: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
179:   `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
180: @*/
181: PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b)
182: {
183:   PetscFunctionBegin;
186:   PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b));
187:   PetscFunctionReturn(PETSC_SUCCESS);
188: }

190: /*@
191:   PCGAMGSetGraphSymmetrize - Set the flag to symmetrize the graph used in coarsening

193:   Logically Collective

195:   Input Parameters:
196: + pc - the preconditioner context
197: - b  - default false

199:   Options Database Key:
200: . -pc_gamg_graph_symmetrize <bool,default=false> - Symmetrize the graph

202:   Level: intermediate

204: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
205:   `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
206: @*/
207: PetscErrorCode PCGAMGSetGraphSymmetrize(PC pc, PetscBool b)
208: {
209:   PetscFunctionBegin;
212:   PetscTryMethod(pc, "PCGAMGSetGraphSymmetrize_C", (PC, PetscBool), (pc, b));
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n)
217: {
218:   PC_MG       *mg          = (PC_MG *)pc->data;
219:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
220:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

222:   PetscFunctionBegin;
223:   pc_gamg_agg->aggressive_coarsening_levels = n;
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

227: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
228: {
229:   PC_MG       *mg          = (PC_MG *)pc->data;
230:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
231:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

233:   PetscFunctionBegin;
234:   pc_gamg_agg->aggressive_mis_k = n;
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b)
239: {
240:   PC_MG       *mg          = (PC_MG *)pc->data;
241:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
242:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

244:   PetscFunctionBegin;
245:   pc_gamg_agg->use_aggressive_square_graph = b;
246:   PetscFunctionReturn(PETSC_SUCCESS);
247: }

249: static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b)
250: {
251:   PC_MG       *mg          = (PC_MG *)pc->data;
252:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
253:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

255:   PetscFunctionBegin;
256:   pc_gamg_agg->use_low_mem_filter = b;
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: static PetscErrorCode PCGAMGSetGraphSymmetrize_AGG(PC pc, PetscBool b)
261: {
262:   PC_MG       *mg          = (PC_MG *)pc->data;
263:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
264:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

266:   PetscFunctionBegin;
267:   pc_gamg_agg->graph_symmetrize = b;
268:   PetscFunctionReturn(PETSC_SUCCESS);
269: }

271: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
272: {
273:   PC_MG       *mg          = (PC_MG *)pc->data;
274:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
275:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

277:   PetscFunctionBegin;
278:   pc_gamg_agg->use_minimum_degree_ordering = b;
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject)
283: {
284:   PC_MG       *mg          = (PC_MG *)pc->data;
285:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
286:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
287:   PetscBool    n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
288:   PetscInt     nsq_graph_old = 0;

290:   PetscFunctionBegin;
291:   PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
292:   PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "number of smoothing steps to construct prolongation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL));
293:   // aggressive coarsening logic with deprecated -pc_gamg_square_graph
294:   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, &n_aggressive_flg));
295:   if (!n_aggressive_flg)
296:     PetscCall(PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (deprecated alias for -pc_gamg_aggressive_coarsening)", "PCGAMGSetAggressiveLevels", nsq_graph_old, &nsq_graph_old, &old_sq_provided));
297:   PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening", "PCGAMGSetAggressiveSquareGraph", new_sqr_graph, &pc_gamg_agg->use_aggressive_square_graph, &new_sq_provided));
298:   if (!new_sq_provided && old_sq_provided) {
299:     pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
300:     pc_gamg_agg->use_aggressive_square_graph  = PETSC_TRUE;
301:   }
302:   if (new_sq_provided && old_sq_provided)
303:     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 %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
304:   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));
305:   PetscCall(PetscOptionsBool("-pc_gamg_low_memory_threshold_filter", "Use the (built-in) low memory graph/matrix filter", "PCGAMGSetLowMemoryFilter", pc_gamg_agg->use_low_mem_filter, &pc_gamg_agg->use_low_mem_filter, NULL));
306:   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));
307:   PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL));
308:   PetscOptionsHeadEnd();
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
313: {
314:   PC_MG       *mg          = (PC_MG *)pc->data;
315:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
316:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

318:   PetscFunctionBegin;
319:   PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
320:   PetscCall(PetscFree(pc_gamg->subctx));
321:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
322:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
323:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
324:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
325:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
326:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
327:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL));
328:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: /*
333:    PCSetCoordinates_AGG

335:    Collective

337:    Input Parameter:
338:    . pc - the preconditioner context
339:    . ndm - dimension of data (used for dof/vertex for Stokes)
340:    . a_nloc - number of vertices local
341:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
342: */

344: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
345: {
346:   PC_MG   *mg      = (PC_MG *)pc->data;
347:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
348:   PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
349:   Mat      mat = pc->pmat;

351:   PetscFunctionBegin;
354:   nloc = a_nloc;

356:   /* SA: null space vectors */
357:   PetscCall(MatGetBlockSize(mat, &ndf));               /* this does not work for Stokes */
358:   if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
359:   else if (coords) {
360:     PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf);
361:     pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */
362:     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);
363:   } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
364:   pc_gamg->data_cell_rows = ndatarows = ndf;
365:   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);
366:   arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols;

368:   if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
369:     PetscCall(PetscFree(pc_gamg->data));
370:     PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
371:   }
372:   /* copy data in - column-oriented */
373:   for (kk = 0; kk < nloc; kk++) {
374:     const PetscInt M    = nloc * pc_gamg->data_cell_rows; /* stride into data */
375:     PetscReal     *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */

377:     if (pc_gamg->data_cell_cols == 1) *data = 1.0;
378:     else {
379:       /* translational modes */
380:       for (ii = 0; ii < ndatarows; ii++) {
381:         for (jj = 0; jj < ndatarows; jj++) {
382:           if (ii == jj) data[ii * M + jj] = 1.0;
383:           else data[ii * M + jj] = 0.0;
384:         }
385:       }

387:       /* rotational modes */
388:       if (coords) {
389:         if (ndm == 2) {
390:           data += 2 * M;
391:           data[0] = -coords[2 * kk + 1];
392:           data[1] = coords[2 * kk];
393:         } else {
394:           data += 3 * M;
395:           data[0]         = 0.0;
396:           data[M + 0]     = coords[3 * kk + 2];
397:           data[2 * M + 0] = -coords[3 * kk + 1];
398:           data[1]         = -coords[3 * kk + 2];
399:           data[M + 1]     = 0.0;
400:           data[2 * M + 1] = coords[3 * kk];
401:           data[2]         = coords[3 * kk + 1];
402:           data[M + 2]     = -coords[3 * kk];
403:           data[2 * M + 2] = 0.0;
404:         }
405:       }
406:     }
407:   }
408:   pc_gamg->data_sz = arrsz;
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: /*
413:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
414:       Looks in Mat for near null space.
415:       Does not work for Stokes

417:   Input Parameter:
418:    . pc -
419:    . a_A - matrix to get (near) null space out of.
420: */
421: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
422: {
423:   PC_MG       *mg      = (PC_MG *)pc->data;
424:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
425:   MatNullSpace mnull;

427:   PetscFunctionBegin;
428:   PetscCall(MatGetNearNullSpace(a_A, &mnull));
429:   if (!mnull) {
430:     DM dm;

432:     PetscCall(PCGetDM(pc, &dm));
433:     if (!dm) PetscCall(MatGetDM(a_A, &dm));
434:     if (dm) {
435:       PetscObject deformation;
436:       PetscInt    Nf;

438:       PetscCall(DMGetNumFields(dm, &Nf));
439:       if (Nf) {
440:         PetscCall(DMGetField(dm, 0, NULL, &deformation));
441:         if (deformation) {
442:           PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull));
443:           if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull));
444:         }
445:       }
446:     }
447:   }

449:   if (!mnull) {
450:     PetscInt bs, NN, MM;

452:     PetscCall(MatGetBlockSize(a_A, &bs));
453:     PetscCall(MatGetLocalSize(a_A, &MM, &NN));
454:     PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
455:     PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
456:   } else {
457:     PetscReal         *nullvec;
458:     PetscBool          has_const;
459:     PetscInt           i, j, mlocal, nvec, bs;
460:     const Vec         *vecs;
461:     const PetscScalar *v;

463:     PetscCall(MatGetLocalSize(a_A, &mlocal, NULL));
464:     PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
465:     for (i = 0; i < nvec; i++) {
466:       PetscCall(VecGetLocalSize(vecs[i], &j));
467:       PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal);
468:     }
469:     pc_gamg->data_sz = (nvec + !!has_const) * mlocal;
470:     PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec));
471:     if (has_const)
472:       for (i = 0; i < mlocal; i++) nullvec[i] = 1.0;
473:     for (i = 0; i < nvec; i++) {
474:       PetscCall(VecGetArrayRead(vecs[i], &v));
475:       for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]);
476:       PetscCall(VecRestoreArrayRead(vecs[i], &v));
477:     }
478:     pc_gamg->data           = nullvec;
479:     pc_gamg->data_cell_cols = (nvec + !!has_const);
480:     PetscCall(MatGetBlockSize(a_A, &bs));
481:     pc_gamg->data_cell_rows = bs;
482:   }
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

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

489:   Input Parameter:
490:    . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
491:    . bs - row block size
492:    . nSAvec - column bs of new P
493:    . my0crs - global index of start of locals
494:    . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
495:    . data_in[data_stride*nSAvec] - local data on fine grid
496:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'

498:   Output Parameter:
499:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
500:    . a_Prol - prolongation operator
501: */
502: 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)
503: {
504:   PetscInt        Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
505:   MPI_Comm        comm;
506:   PetscReal      *out_data;
507:   PetscCDIntNd   *pos;
508:   PCGAMGHashTable fgid_flid;

510:   PetscFunctionBegin;
511:   PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
512:   PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
513:   nloc = (Iend - Istart) / bs;
514:   my0  = Istart / bs;
515:   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);
516:   Iend /= bs;
517:   nghosts = data_stride / bs - nloc;

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

522:   /* count selected -- same as number of cols of P */
523:   for (nSelected = mm = 0; mm < nloc; mm++) {
524:     PetscBool ise;

526:     PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
527:     if (!ise) nSelected++;
528:   }
529:   PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
530:   PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT "  != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
531:   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);

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

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

540:   /* find points and set prolongation */
541:   minsz = 100;
542:   for (mm = clid = 0; mm < nloc; mm++) {
543:     PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
544:     if (jj > 0) {
545:       const PetscInt lid = mm, cgid = my0crs + clid;
546:       PetscInt       cids[100]; /* max bs */
547:       PetscBLASInt   asz, M, N, INFO;
548:       PetscBLASInt   Mdata, LDA, LWORK;
549:       PetscScalar   *qqc, *qqr, *TAU, *WORK;
550:       PetscInt      *fids;
551:       PetscReal     *data;

553:       PetscCall(PetscBLASIntCast(jj, &asz));
554:       PetscCall(PetscBLASIntCast(asz * bs, &M));
555:       PetscCall(PetscBLASIntCast(nSAvec, &N));
556:       PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata));
557:       PetscCall(PetscBLASIntCast(Mdata, &LDA));
558:       PetscCall(PetscBLASIntCast(N * bs, &LWORK));
559:       /* count agg */
560:       if (asz < minsz) minsz = asz;

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

565:       aggID = 0;
566:       PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
567:       while (pos) {
568:         PetscInt gid1;

570:         PetscCall(PetscCDIntNdGetID(pos, &gid1));
571:         PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));

573:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
574:         else {
575:           PetscCall(PCGAMGHashTableFind(&fgid_flid, gid1, &flid));
576:           PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
577:         }
578:         /* copy in B_i matrix - column-oriented */
579:         data = &data_in[flid * bs];
580:         for (ii = 0; ii < bs; ii++) {
581:           for (jj = 0; jj < N; jj++) {
582:             PetscReal d = data[jj * data_stride + ii];

584:             qqc[jj * Mdata + aggID * bs + ii] = d;
585:           }
586:         }
587:         /* set fine IDs */
588:         for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
589:         aggID++;
590:       }

592:       /* pad with zeros */
593:       for (ii = asz * bs; ii < Mdata; ii++) {
594:         for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
595:       }

597:       /* QR */
598:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
599:       PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
600:       PetscCall(PetscFPTrapPop());
601:       PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
602:       /* get R - column-oriented - output B_{i+1} */
603:       {
604:         PetscReal *data = &out_data[clid * nSAvec];

606:         for (jj = 0; jj < nSAvec; jj++) {
607:           for (ii = 0; ii < nSAvec; ii++) {
608:             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);
609:             if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
610:             else data[jj * out_data_stride + ii] = 0.;
611:           }
612:         }
613:       }

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

619:       for (ii = 0; ii < M; ii++) {
620:         for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
621:       }

623:       /* add diagonal block of P0 */
624:       for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ }
625:       PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
626:       PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
627:       clid++;
628:     } /* coarse agg */
629:   } /* for all fine nodes */
630:   PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
631:   PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
632:   PetscCall(PCGAMGHashTableDestroy(&fgid_flid));
633:   PetscFunctionReturn(PETSC_SUCCESS);
634: }

636: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
637: {
638:   PC_MG       *mg          = (PC_MG *)pc->data;
639:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
640:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

642:   PetscFunctionBegin;
643:   PetscCall(PetscViewerASCIIPrintf(viewer, "      AGG specific options\n"));
644:   PetscCall(PetscViewerASCIIPrintf(viewer, "        Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
645:   if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
646:     PetscCall(PetscViewerASCIIPrintf(viewer, "        %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
647:     if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(PetscViewerASCIIPrintf(viewer, "        MIS-%" PetscInt_FMT " coarsening on aggressive levels\n", pc_gamg_agg->aggressive_mis_k));
648:   }
649:   PetscCall(PetscViewerASCIIPushTab(viewer));
650:   PetscCall(PetscViewerASCIIPushTab(viewer));
651:   PetscCall(PetscViewerASCIIPushTab(viewer));
652:   PetscCall(PetscViewerASCIIPushTab(viewer));
653:   if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer));
654:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n"));
655:   PetscCall(PetscViewerASCIIPopTab(viewer));
656:   PetscCall(PetscViewerASCIIPopTab(viewer));
657:   PetscCall(PetscViewerASCIIPopTab(viewer));
658:   PetscCall(PetscViewerASCIIPopTab(viewer));
659:   PetscCall(PetscViewerASCIIPrintf(viewer, "        Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
664: {
665:   PC_MG          *mg          = (PC_MG *)pc->data;
666:   PC_GAMG        *pc_gamg     = (PC_GAMG *)mg->innerctx;
667:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
668:   const PetscReal vfilter     = pc_gamg->threshold[pc_gamg->current_level];
669:   PetscBool       ishem, ismis;
670:   const char     *prefix;
671:   MatInfo         info0, info1;
672:   PetscInt        bs;

674:   PetscFunctionBegin;
675:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
676:   /* 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 */
677:   /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
678:   PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
679:   PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
680:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
681:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
682:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_"));
683:   PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
684:   PetscCall(MatGetBlockSize(Amat, &bs));
685:   // check for valid indices wrt bs
686:   for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
687:     PetscCheck(pc_gamg_agg->crs->strength_index[ii] >= 0 && pc_gamg_agg->crs->strength_index[ii] < bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Indices (%" PetscInt_FMT ") must be non-negative and < block size (%" PetscInt_FMT "), NB, can not use -mat_coarsen_strength_index with -mat_coarsen_strength_index",
688:                pc_gamg_agg->crs->strength_index[ii], bs);
689:   }
690:   PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
691:   if (ishem) {
692:     if (pc_gamg_agg->aggressive_coarsening_levels) PetscCall(PetscInfo(pc, "HEM and aggressive coarsening ignored: HEM using %" PetscInt_FMT " iterations\n", pc_gamg_agg->crs->max_it));
693:     pc_gamg_agg->aggressive_coarsening_levels = 0;                                         // aggressive and HEM does not make sense
694:     PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
695:     PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter));                          // for code coverage
696:   } else {
697:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
698:     if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
699:       PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
700:       pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
701:     }
702:   }
703:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
704:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
705:   PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */

707:   if (ishem || pc_gamg_agg->use_low_mem_filter) {
708:     PetscCall(MatCreateGraph(Amat, pc_gamg_agg->graph_symmetrize, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat));
709:   } else {
710:     // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive)
711:     PetscCall(MatCreateGraph(Amat, pc_gamg_agg->graph_symmetrize, PETSC_TRUE, -1, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat));
712:     if (vfilter >= 0) {
713:       PetscInt           Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
714:       Mat                tGmat, Gmat = *a_Gmat;
715:       MPI_Comm           comm;
716:       const PetscScalar *vals;
717:       const PetscInt    *idx;
718:       PetscInt          *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
719:       MatScalar         *AA; // this is checked in graph
720:       PetscBool          isseqaij;
721:       Mat                a, b, c;
722:       MatType            jtype;

724:       PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
725:       PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
726:       PetscCall(MatGetType(Gmat, &jtype));
727:       PetscCall(MatCreate(comm, &tGmat));
728:       PetscCall(MatSetType(tGmat, jtype));

730:       /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold?
731:         Also, if the matrix is symmetric, can we skip this
732:         operation? It can be very expensive on large matrices. */

734:       // global sizes
735:       PetscCall(MatGetSize(Gmat, &MM, &NN));
736:       PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
737:       nloc = Iend - Istart;
738:       PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
739:       if (isseqaij) {
740:         a = Gmat;
741:         b = NULL;
742:       } else {
743:         Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;

745:         a      = d->A;
746:         b      = d->B;
747:         garray = d->garray;
748:       }
749:       /* Determine upper bound on non-zeros needed in new filtered matrix */
750:       for (PetscInt row = 0; row < nloc; row++) {
751:         PetscCall(MatGetRow(a, row, &ncols, NULL, NULL));
752:         d_nnz[row] = ncols;
753:         if (ncols > maxcols) maxcols = ncols;
754:         PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL));
755:       }
756:       if (b) {
757:         for (PetscInt row = 0; row < nloc; row++) {
758:           PetscCall(MatGetRow(b, row, &ncols, NULL, NULL));
759:           o_nnz[row] = ncols;
760:           if (ncols > maxcols) maxcols = ncols;
761:           PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL));
762:         }
763:       }
764:       PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM));
765:       PetscCall(MatSetBlockSizes(tGmat, 1, 1));
766:       PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz));
767:       PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz));
768:       PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
769:       PetscCall(PetscFree2(d_nnz, o_nnz));
770:       PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ));
771:       nnz0 = nnz1 = 0;
772:       for (c = a, kk = 0; c && kk < 2; c = b, kk++) {
773:         for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) {
774:           PetscCall(MatGetRow(c, row, &ncols, &idx, &vals));
775:           for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) {
776:             PetscScalar sv = PetscAbs(PetscRealPart(vals[jj]));
777:             if (PetscRealPart(sv) > vfilter) {
778:               PetscInt cid = idx[jj] + Istart; //diag

780:               nnz1++;
781:               if (c != a) cid = garray[idx[jj]];
782:               AA[ncol_row] = vals[jj];
783:               AJ[ncol_row] = cid;
784:               ncol_row++;
785:             }
786:           }
787:           PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
788:           PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
789:         }
790:       }
791:       PetscCall(PetscFree2(AA, AJ));
792:       PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
793:       PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
794:       PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
795:       PetscCall(PetscInfo(pc, "\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%" PetscInt_FMT ", max row size %" PetscInt_FMT "\n", (!nnz0) ? 1. : 100. * (double)nnz1 / (double)nnz0, (double)vfilter, (!nloc) ? 1. : (double)nnz0 / (double)nloc, MM, maxcols));
796:       PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
797:       PetscCall(MatDestroy(&Gmat));
798:       *a_Gmat = tGmat;
799:     }
800:   }

802:   PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
803:   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));
804:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
805:   PetscFunctionReturn(PETSC_SUCCESS);
806: }

808: typedef PetscInt    NState;
809: static const NState NOT_DONE = -2;
810: static const NState DELETED  = -1;
811: static const NState REMOVED  = -3;
812: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)

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

818:    Input Parameter:
819:    . Gmat_2 - global matrix of squared graph (data not defined)
820:    . Gmat_1 - base graph to grab with base graph
821:    Input/Output Parameter:
822:    . aggs_2 - linked list of aggs with gids)
823: */
824: static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2)
825: {
826:   PetscBool      isMPI;
827:   Mat_SeqAIJ    *matA_1, *matB_1 = NULL;
828:   MPI_Comm       comm;
829:   PetscInt       lid, *ii, *idx, ix, Iend, my0, kk, n, j;
830:   Mat_MPIAIJ    *mpimat_2 = NULL, *mpimat_1 = NULL;
831:   const PetscInt nloc = Gmat_2->rmap->n;
832:   PetscScalar   *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid;
833:   PetscInt      *lid_cprowID_1 = NULL;
834:   NState        *lid_state;
835:   Vec            ghost_par_orig2;
836:   PetscMPIInt    rank;

838:   PetscFunctionBegin;
839:   PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
840:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
841:   PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));

843:   /* get submatrices */
844:   PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
845:   PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
846:   PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
847:   for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
848:   if (isMPI) {
849:     /* grab matrix objects */
850:     mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
851:     mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
852:     matA_1   = (Mat_SeqAIJ *)mpimat_1->A->data;
853:     matB_1   = (Mat_SeqAIJ *)mpimat_1->B->data;

855:     /* force compressed row storage for B matrix in AuxMat */
856:     PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0));
857:     for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) {
858:       PetscInt lid = matB_1->compressedrow.rindex[ix];

860:       PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
861:       if (lid != -1) lid_cprowID_1[lid] = ix;
862:     }
863:   } else {
864:     PetscBool isAIJ;

866:     PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
867:     PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
868:     matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
869:   }
870:   if (nloc > 0) { PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???"); }
871:   /* get state of locals and selected gid for deleted */
872:   for (lid = 0; lid < nloc; lid++) {
873:     lid_parent_gid[lid] = -1.0;
874:     lid_state[lid]      = DELETED;
875:   }

877:   /* set lid_state */
878:   for (lid = 0; lid < nloc; lid++) {
879:     PetscCDIntNd *pos;

881:     PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
882:     if (pos) {
883:       PetscInt gid1;

885:       PetscCall(PetscCDIntNdGetID(pos, &gid1));
886:       PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
887:       lid_state[lid] = gid1;
888:     }
889:   }

891:   /* map local to selected local, DELETED means a ghost owns it */
892:   for (lid = 0; lid < nloc; lid++) {
893:     NState state = lid_state[lid];

895:     if (IS_SELECTED(state)) {
896:       PetscCDIntNd *pos;

898:       PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
899:       while (pos) {
900:         PetscInt gid1;

902:         PetscCall(PetscCDIntNdGetID(pos, &gid1));
903:         PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
904:         if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
905:       }
906:     }
907:   }
908:   /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
909:   if (isMPI) {
910:     Vec tempVec;

912:     /* get 'cpcol_1_state' */
913:     PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
914:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
915:       PetscScalar v = (PetscScalar)lid_state[kk];

917:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
918:     }
919:     PetscCall(VecAssemblyBegin(tempVec));
920:     PetscCall(VecAssemblyEnd(tempVec));
921:     PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
922:     PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
923:     PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
924:     /* get 'cpcol_2_state' */
925:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
926:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
927:     PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
928:     /* get 'cpcol_2_par_orig' */
929:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
930:       PetscScalar v = lid_parent_gid[kk];

932:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
933:     }
934:     PetscCall(VecAssemblyBegin(tempVec));
935:     PetscCall(VecAssemblyEnd(tempVec));
936:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
937:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
938:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
939:     PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));

941:     PetscCall(VecDestroy(&tempVec));
942:   } /* ismpi */
943:   for (lid = 0; lid < nloc; lid++) {
944:     NState state = lid_state[lid];

946:     if (IS_SELECTED(state)) {
947:       /* steal locals */
948:       ii  = matA_1->i;
949:       n   = ii[lid + 1] - ii[lid];
950:       idx = matA_1->j + ii[lid];
951:       for (j = 0; j < n; j++) {
952:         PetscInt lidj   = idx[j], sgid;
953:         NState   statej = lid_state[lidj];

955:         if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */
956:           lid_parent_gid[lidj] = (PetscScalar)(lid + my0);                                              /* send this if sgid is not local */
957:           if (sgid >= my0 && sgid < Iend) {                                                             /* I'm stealing this local from a local sgid */
958:             PetscInt      hav = 0, slid = sgid - my0, gidj = lidj + my0;
959:             PetscCDIntNd *pos, *last = NULL;

961:             /* looking for local from local so id_llist_2 works */
962:             PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
963:             while (pos) {
964:               PetscInt gid;

966:               PetscCall(PetscCDIntNdGetID(pos, &gid));
967:               if (gid == gidj) {
968:                 PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
969:                 PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last));
970:                 PetscCall(PetscCDAppendNode(aggs_2, lid, pos));
971:                 hav = 1;
972:                 break;
973:               } else last = pos;
974:               PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos));
975:             }
976:             if (hav != 1) {
977:               PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix");
978:               SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
979:             }
980:           } else { /* I'm stealing this local, owned by a ghost */
981:             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.",
982:                        ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "");
983:             PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0));
984:           }
985:         }
986:       } /* local neighbors */
987:     } else if (state == DELETED /* && lid_cprowID_1 */) {
988:       PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);

990:       /* see if I have a selected ghost neighbor that will steal me */
991:       if ((ix = lid_cprowID_1[lid]) != -1) {
992:         ii  = matB_1->compressedrow.i;
993:         n   = ii[ix + 1] - ii[ix];
994:         idx = matB_1->j + ii[ix];
995:         for (j = 0; j < n; j++) {
996:           PetscInt cpid   = idx[j];
997:           NState   statej = (NState)PetscRealPart(cpcol_1_state[cpid]);

999:           if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */
1000:             lid_parent_gid[lid] = (PetscScalar)statej;    /* send who selected */
1001:             if (sgidold >= my0 && sgidold < Iend) {       /* this was mine */
1002:               PetscInt      hav = 0, oldslidj = sgidold - my0;
1003:               PetscCDIntNd *pos, *last        = NULL;

1005:               /* remove from 'oldslidj' list */
1006:               PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
1007:               while (pos) {
1008:                 PetscInt gid;

1010:                 PetscCall(PetscCDIntNdGetID(pos, &gid));
1011:                 if (lid + my0 == gid) {
1012:                   /* id_llist_2[lastid] = id_llist_2[flid];   /\* remove lid from oldslidj list *\/ */
1013:                   PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null");
1014:                   PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last));
1015:                   /* ghost (PetscScalar)statej will add this later */
1016:                   hav = 1;
1017:                   break;
1018:                 } else last = pos;
1019:                 PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos));
1020:               }
1021:               if (hav != 1) {
1022:                 PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav);
1023:                 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav);
1024:               }
1025:             } else {
1026:               /* TODO: ghosts remove this later */
1027:             }
1028:           }
1029:         }
1030:       }
1031:     } /* selected/deleted */
1032:   } /* node loop */

1034:   if (isMPI) {
1035:     PetscScalar    *cpcol_2_parent, *cpcol_2_gid;
1036:     Vec             tempVec, ghostgids2, ghostparents2;
1037:     PetscInt        cpid, nghost_2;
1038:     PCGAMGHashTable gid_cpid;

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

1043:     /* get 'cpcol_2_parent' */
1044:     for (kk = 0, j = my0; kk < nloc; kk++, j++) { PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES)); }
1045:     PetscCall(VecAssemblyBegin(tempVec));
1046:     PetscCall(VecAssemblyEnd(tempVec));
1047:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
1048:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1049:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
1050:     PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));

1052:     /* get 'cpcol_2_gid' */
1053:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
1054:       PetscScalar v = (PetscScalar)j;

1056:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
1057:     }
1058:     PetscCall(VecAssemblyBegin(tempVec));
1059:     PetscCall(VecAssemblyEnd(tempVec));
1060:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
1061:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1062:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
1063:     PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
1064:     PetscCall(VecDestroy(&tempVec));

1066:     /* look for deleted ghosts and add to table */
1067:     PetscCall(PCGAMGHashTableCreate(2 * nghost_2 + 1, &gid_cpid));
1068:     for (cpid = 0; cpid < nghost_2; cpid++) {
1069:       NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);

1071:       if (state == DELETED) {
1072:         PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
1073:         PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);

1075:         if (sgid_old == -1 && sgid_new != -1) {
1076:           PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);

1078:           PetscCall(PCGAMGHashTableAdd(&gid_cpid, gid, cpid));
1079:         }
1080:       }
1081:     }

1083:     /* look for deleted ghosts and see if they moved - remove it */
1084:     for (lid = 0; lid < nloc; lid++) {
1085:       NState state = lid_state[lid];

1087:       if (IS_SELECTED(state)) {
1088:         PetscCDIntNd *pos, *last = NULL;

1090:         /* look for deleted ghosts and see if they moved */
1091:         PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
1092:         while (pos) {
1093:           PetscInt gid;

1095:           PetscCall(PetscCDIntNdGetID(pos, &gid));
1096:           if (gid < my0 || gid >= Iend) {
1097:             PetscCall(PCGAMGHashTableFind(&gid_cpid, gid, &cpid));
1098:             if (cpid != -1) {
1099:               /* a moved ghost - */
1100:               /* id_llist_2[lastid] = id_llist_2[flid];    /\* remove 'flid' from list *\/ */
1101:               PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
1102:             } else last = pos;
1103:           } else last = pos;

1105:           PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
1106:         } /* loop over list of deleted */
1107:       } /* selected */
1108:     }
1109:     PetscCall(PCGAMGHashTableDestroy(&gid_cpid));

1111:     /* look at ghosts, see if they changed - and it */
1112:     for (cpid = 0; cpid < nghost_2; cpid++) {
1113:       PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);

1115:       if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
1116:         PetscInt      gid      = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
1117:         PetscInt      slid_new = sgid_new - my0, hav = 0;
1118:         PetscCDIntNd *pos;

1120:         /* search for this gid to see if I have it */
1121:         PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
1122:         while (pos) {
1123:           PetscInt gidj;

1125:           PetscCall(PetscCDIntNdGetID(pos, &gidj));
1126:           PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));

1128:           if (gidj == gid) {
1129:             hav = 1;
1130:             break;
1131:           }
1132:         }
1133:         if (hav != 1) {
1134:           /* insert 'flidj' into head of llist */
1135:           PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
1136:         }
1137:       }
1138:     }
1139:     PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
1140:     PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
1141:     PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
1142:     PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
1143:     PetscCall(VecDestroy(&ghostgids2));
1144:     PetscCall(VecDestroy(&ghostparents2));
1145:     PetscCall(VecDestroy(&ghost_par_orig2));
1146:   }
1147:   PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
1148:   PetscFunctionReturn(PETSC_SUCCESS);
1149: }

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

1155:   Input Parameter:
1156:    . a_pc - this

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

1161:   Output Parameter:
1162:    . agg_lists - list of aggregates

1164: */
1165: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
1166: {
1167:   PC_MG       *mg          = (PC_MG *)a_pc->data;
1168:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
1169:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1170:   Mat          Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
1171:   IS           perm;
1172:   PetscInt     Istart, Iend, Ii, nloc, bs, nn;
1173:   PetscInt    *permute, *degree;
1174:   PetscBool   *bIndexSet;
1175:   PetscReal    hashfact;
1176:   PetscInt     iSwapIndex;
1177:   PetscRandom  random;
1178:   MPI_Comm     comm;

1180:   PetscFunctionBegin;
1181:   PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
1182:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1183:   PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
1184:   PetscCall(MatGetBlockSize(Gmat1, &bs));
1185:   PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
1186:   nloc = nn / bs;
1187:   /* get MIS aggs - randomize */
1188:   PetscCall(PetscMalloc2(nloc, &permute, nloc, &degree));
1189:   PetscCall(PetscCalloc1(nloc, &bIndexSet));
1190:   for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
1191:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
1192:   PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
1193:   for (Ii = 0; Ii < nloc; Ii++) {
1194:     PetscInt nc;

1196:     PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1197:     degree[Ii] = nc;
1198:     PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
1199:   }
1200:   for (Ii = 0; Ii < nloc; Ii++) {
1201:     PetscCall(PetscRandomGetValueReal(random, &hashfact));
1202:     iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
1203:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1204:       PetscInt iTemp = permute[iSwapIndex];

1206:       permute[iSwapIndex]   = permute[Ii];
1207:       permute[Ii]           = iTemp;
1208:       iTemp                 = degree[iSwapIndex];
1209:       degree[iSwapIndex]    = degree[Ii];
1210:       degree[Ii]            = iTemp;
1211:       bIndexSet[iSwapIndex] = PETSC_TRUE;
1212:     }
1213:   }
1214:   // apply minimum degree ordering -- NEW
1215:   if (pc_gamg_agg->use_minimum_degree_ordering) { PetscCall(PetscSortIntWithArray(nloc, degree, permute)); }
1216:   PetscCall(PetscFree(bIndexSet));
1217:   PetscCall(PetscRandomDestroy(&random));
1218:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm));
1219:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));
1220:   // square graph
1221:   if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels && pc_gamg_agg->use_aggressive_square_graph) {
1222:     PetscCall(PCGAMGSquareGraph_GAMG(a_pc, Gmat1, &Gmat2));
1223:   } else Gmat2 = Gmat1;
1224:   // switch to old MIS-1 for square graph
1225:   if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) {
1226:     if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2
1227:     else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS));                                                                   // old MIS -- side effect
1228:   } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) {                                 // we reset the MIS
1229:     const char *prefix;

1231:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
1232:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
1233:     PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
1234:   }
1235:   PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
1236:   PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
1237:   PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
1238:   PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
1239:   PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */

1241:   PetscCall(ISDestroy(&perm));
1242:   PetscCall(PetscFree2(permute, degree));
1243:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));

1245:   if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected
1246:     PetscCoarsenData *llist = *agg_lists;

1248:     PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
1249:     PetscCall(MatDestroy(&Gmat1));
1250:     *a_Gmat1 = Gmat2;                          /* output */
1251:     PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
1252:   }
1253:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
1254:   PetscFunctionReturn(PETSC_SUCCESS);
1255: }

1257: /*
1258:  PCGAMGConstructProlongator_AGG

1260:  Input Parameter:
1261:  . pc - this
1262:  . Amat - matrix on this fine level
1263:  . Graph - used to get ghost data for nodes in
1264:  . agg_lists - list of aggregates
1265:  Output Parameter:
1266:  . a_P_out - prolongation operator to the next level
1267:  */
1268: static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out)
1269: {
1270:   PC_MG         *mg      = (PC_MG *)pc->data;
1271:   PC_GAMG       *pc_gamg = (PC_GAMG *)mg->innerctx;
1272:   const PetscInt col_bs  = pc_gamg->data_cell_cols;
1273:   PetscInt       Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs;
1274:   Mat            Gmat, Prol;
1275:   PetscMPIInt    size;
1276:   MPI_Comm       comm;
1277:   PetscReal     *data_w_ghost;
1278:   PetscInt       myCrs0, nbnodes = 0, *flid_fgid;
1279:   MatType        mtype;

1281:   PetscFunctionBegin;
1282:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1283:   PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
1284:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1285:   PetscCallMPI(MPI_Comm_size(comm, &size));
1286:   PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
1287:   PetscCall(MatGetBlockSize(Amat, &bs));
1288:   nloc = (Iend - Istart) / bs;
1289:   my0  = Istart / bs;
1290:   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);
1291:   PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1

1293:   /* get 'nLocalSelected' */
1294:   for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
1295:     PetscBool ise;

1297:     /* filter out singletons 0 or 1? */
1298:     PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
1299:     if (!ise) nLocalSelected++;
1300:   }

1302:   /* create prolongator, create P matrix */
1303:   PetscCall(MatGetType(Amat, &mtype));
1304:   PetscCall(MatCreate(comm, &Prol));
1305:   PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
1306:   PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
1307:   PetscCall(MatSetType(Prol, mtype));
1308: #if PetscDefined(HAVE_DEVICE)
1309:   PetscBool flg;
1310:   PetscCall(MatBoundToCPU(Amat, &flg));
1311:   PetscCall(MatBindToCPU(Prol, flg));
1312:   if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
1313: #endif
1314:   PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
1315:   PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));

1317:   /* can get all points "removed" */
1318:   PetscCall(MatGetSize(Prol, &kk, &ii));
1319:   if (!ii) {
1320:     PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
1321:     PetscCall(MatDestroy(&Prol));
1322:     *a_P_out = NULL; /* out */
1323:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1324:     PetscFunctionReturn(PETSC_SUCCESS);
1325:   }
1326:   PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
1327:   PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));

1329:   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);
1330:   myCrs0 = myCrs0 / col_bs;
1331:   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);

1333:   /* create global vector of data in 'data_w_ghost' */
1334:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1335:   if (size > 1) { /* get ghost null space data */
1336:     PetscReal *tmp_gdata, *tmp_ldata, *tp2;

1338:     PetscCall(PetscMalloc1(nloc, &tmp_ldata));
1339:     for (jj = 0; jj < col_bs; jj++) {
1340:       for (kk = 0; kk < bs; kk++) {
1341:         PetscInt         ii, stride;
1342:         const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);

1344:         for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;

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

1348:         if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
1349:           PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
1350:           nbnodes = bs * stride;
1351:         }
1352:         tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
1353:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1354:         PetscCall(PetscFree(tmp_gdata));
1355:       }
1356:     }
1357:     PetscCall(PetscFree(tmp_ldata));
1358:   } else {
1359:     nbnodes      = bs * nloc;
1360:     data_w_ghost = pc_gamg->data;
1361:   }

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

1368:     PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
1369:     for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
1370:     PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
1371:     PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
1372:     for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1373:     PetscCall(PetscFree(fiddata));

1375:     PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
1376:     PetscCall(PetscFree(fid_glid_loc));
1377:   } else {
1378:     PetscCall(PetscMalloc1(nloc, &flid_fgid));
1379:     for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
1380:   }
1381:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
1382:   /* get P0 */
1383:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1384:   {
1385:     PetscReal *data_out = NULL;

1387:     PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol));
1388:     PetscCall(PetscFree(pc_gamg->data));

1390:     pc_gamg->data           = data_out;
1391:     pc_gamg->data_cell_rows = col_bs;
1392:     pc_gamg->data_sz        = col_bs * col_bs * nLocalSelected;
1393:   }
1394:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
1395:   if (size > 1) PetscCall(PetscFree(data_w_ghost));
1396:   PetscCall(PetscFree(flid_fgid));

1398:   *a_P_out = Prol; /* out */
1399:   PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation"));

1401:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
1402:   PetscFunctionReturn(PETSC_SUCCESS);
1403: }

1405: /*
1406:    PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times

1408:   Input Parameter:
1409:    . pc - this
1410:    . Amat - matrix on this fine level
1411:  In/Output Parameter:
1412:    . a_P - prolongation operator to the next level
1413: */
1414: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
1415: {
1416:   PC_MG       *mg          = (PC_MG *)pc->data;
1417:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
1418:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
1419:   PetscInt     jj;
1420:   Mat          Prol = *a_P;
1421:   MPI_Comm     comm;
1422:   KSP          eksp;
1423:   Vec          bb, xx;
1424:   PC           epc;
1425:   PetscReal    alpha, emax, emin;

1427:   PetscFunctionBegin;
1428:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
1429:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));

1431:   /* compute maximum singular value of operator to be used in smoother */
1432:   if (0 < pc_gamg_agg->nsmooths) {
1433:     /* get eigen estimates */
1434:     if (pc_gamg->emax > 0) {
1435:       emin = pc_gamg->emin;
1436:       emax = pc_gamg->emax;
1437:     } else {
1438:       const char *prefix;

1440:       PetscCall(MatCreateVecs(Amat, &bb, NULL));
1441:       PetscCall(MatCreateVecs(Amat, &xx, NULL));
1442:       PetscCall(KSPSetNoisy_Private(Amat, bb));

1444:       PetscCall(KSPCreate(comm, &eksp));
1445:       PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
1446:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
1447:       PetscCall(KSPSetOptionsPrefix(eksp, prefix));
1448:       PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
1449:       {
1450:         PetscBool isset, sflg;

1452:         PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
1453:         if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
1454:       }
1455:       PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
1456:       PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));

1458:       PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
1459:       PetscCall(KSPSetOperators(eksp, Amat, Amat));

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

1464:       PetscCall(KSPSetTolerances(eksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 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

1466:       PetscCall(KSPSetFromOptions(eksp));
1467:       PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
1468:       PetscCall(KSPSolve(eksp, bb, xx));
1469:       PetscCall(KSPCheckSolve(eksp, pc, xx));

1471:       PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
1472:       PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
1473:       PetscCall(VecDestroy(&xx));
1474:       PetscCall(VecDestroy(&bb));
1475:       PetscCall(KSPDestroy(&eksp));
1476:     }
1477:     if (pc_gamg->use_sa_esteig) {
1478:       mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
1479:       mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
1480:       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));
1481:     } else {
1482:       mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1483:       mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1484:     }
1485:   } else {
1486:     mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
1487:     mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
1488:   }

1490:   /* smooth P0 */
1491:   if (pc_gamg_agg->nsmooths > 0) {
1492:     Vec diag;

1494:     /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */
1495:     PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero");

1497:     PetscCall(MatCreateVecs(Amat, &diag, NULL));
1498:     PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
1499:     PetscCall(VecReciprocal(diag));

1501:     for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1502:       Mat tMat;

1504:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1505:       /*
1506:         Smooth aggregation on the prolongator

1508:         P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
1509:       */
1510:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1511:       PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
1512:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
1513:       PetscCall(MatProductClear(tMat));
1514:       PetscCall(MatDiagonalScale(tMat, diag, NULL));

1516:       /* TODO: Document the 1.4 and don't hardwire it in this routine */
1517:       alpha = -1.4 / emax;
1518:       PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
1519:       PetscCall(MatDestroy(&Prol));
1520:       Prol = tMat;
1521:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
1522:     }
1523:     PetscCall(VecDestroy(&diag));
1524:   }
1525:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
1526:   PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
1527:   *a_P = Prol;
1528:   PetscFunctionReturn(PETSC_SUCCESS);
1529: }

1531: /*MC
1532:   PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner

1534:   Options Database Keys:
1535: + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation to construct prolongation
1536: . -pc_gamg_aggressive_coarsening <n,default=1> - number of aggressive coarsening (MIS-2) levels from finest.
1537: . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening
1538: . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm
1539: . -pc_gamg_pc_gamg_asm_hem_aggs <n,default=0> - Number of HEM aggregation steps for ASM smoother
1540: - -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive')

1542:   Level: intermediate

1544:   Notes:
1545:   To obtain good performance for `PCGAMG` for vector valued problems you must
1546:   call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point.
1547:   Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator

1549:   The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG`

1551: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
1552:           `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
1553:           `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
1554:           `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
1555:           `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
1556: M*/
1557: PetscErrorCode PCCreateGAMG_AGG(PC pc)
1558: {
1559:   PC_MG       *mg      = (PC_MG *)pc->data;
1560:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
1561:   PC_GAMG_AGG *pc_gamg_agg;

1563:   PetscFunctionBegin;
1564:   /* create sub context for SA */
1565:   PetscCall(PetscNew(&pc_gamg_agg));
1566:   pc_gamg->subctx = pc_gamg_agg;

1568:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1569:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
1570:   /* reset does not do anything; setup not virtual */

1572:   /* set internal function pointers */
1573:   pc_gamg->ops->creategraph       = PCGAMGCreateGraph_AGG;
1574:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_AGG;
1575:   pc_gamg->ops->prolongator       = PCGAMGConstructProlongator_AGG;
1576:   pc_gamg->ops->optprolongator    = PCGAMGOptimizeProlongator_AGG;
1577:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1578:   pc_gamg->ops->view              = PCView_GAMG_AGG;

1580:   pc_gamg_agg->nsmooths                     = 1;
1581:   pc_gamg_agg->aggressive_coarsening_levels = 1;
1582:   pc_gamg_agg->use_aggressive_square_graph  = PETSC_TRUE;
1583:   pc_gamg_agg->use_minimum_degree_ordering  = PETSC_FALSE;
1584:   pc_gamg_agg->use_low_mem_filter           = PETSC_FALSE;
1585:   pc_gamg_agg->aggressive_mis_k             = 2;
1586:   pc_gamg_agg->graph_symmetrize             = PETSC_TRUE;

1588:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
1589:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
1590:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
1591:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
1592:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
1593:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
1594:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
1595:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
1596:   PetscFunctionReturn(PETSC_SUCCESS);
1597: }