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> - the flag

     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> - the flag

     76:   Level: intermediate

     78:   Note:
     79:   By default, aggressive coarsening squares the matrix (computes $ A^T A$) before coarsening. Calling `PCGAMGSetAggressiveSquareGraph()` with a value of `PETSC_FALSE` changes the aggressive coarsening strategy to use MIS-k, see `PCGAMGMISkSetAggressive()`.

     81: .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()`
     82: @*/
     83: PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n)
     84: {
     85:   PetscFunctionBegin;
     88:   PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

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

     95:   Logically Collective

     97:   Input Parameters:
     98: + pc - the preconditioner context
     99: - n  - 1 or more (default = 2)

    101:   Options Database Key:
    102: . -pc_gamg_aggressive_mis_k <n,default=2> - the flag

    104:   Level: intermediate

    106: .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()`
    107: @*/
    108: PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n)
    109: {
    110:   PetscFunctionBegin;
    113:   PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n));
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: /*@
    118:   PCGAMGSetAggressiveSquareGraph - Use graph square, $A^T A$, for aggressive coarsening. Coarsening is slower than the alternative (MIS-2), which is faster and uses less memory

    120:   Logically Collective

    122:   Input Parameters:
    123: + pc - the preconditioner context
    124: - b  - default true

    126:   Options Database Key:
    127: . -pc_gamg_aggressive_square_graph <bool,default=true> - the flag

    129:   Level: intermediate

    131:   Notes:
    132:   If `b` is `PETSC_FALSE` then MIS-k is used for aggressive coarsening, see `PCGAMGMISkSetAggressive()`

    134:   Squaring the matrix to perform the aggressive coarsening is slower and requires more memory than using MIS-k, but may result in a better preconditioner
    135:   that converges faster.

    137: .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()`
    138: @*/
    139: PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b)
    140: {
    141:   PetscFunctionBegin;
    144:   PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b));
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: /*@
    149:   PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm

    151:   Logically Collective

    153:   Input Parameters:
    154: + pc - the preconditioner context
    155: - b  - default false

    157:   Options Database Key:
    158: . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=false> - the flag

    160:   Level: intermediate

    162: .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()`
    163: @*/
    164: PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b)
    165: {
    166:   PetscFunctionBegin;
    169:   PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: /*@
    174:   PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter

    176:   Logically Collective

    178:   Input Parameters:
    179: + pc - the preconditioner context
    180: - b  - default false

    182:   Options Database Key:
    183: . -pc_gamg_low_memory_threshold_filter <bool,default=false> - the flag

    185:   Level: intermediate

    187: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`,
    188:   `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
    189: @*/
    190: PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b)
    191: {
    192:   PetscFunctionBegin;
    195:   PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b));
    196:   PetscFunctionReturn(PETSC_SUCCESS);
    197: }

    199: /*@
    200:   PCGAMGSetGraphSymmetrize - Symmetrize graph used for coarsening. Defaults to true, but if matrix has symmetric attribute, then not needed since the graph is already known to be symmetric

    202:   Logically Collective

    204:   Input Parameters:
    205: + pc - the preconditioner context
    206: - b  - default true

    208:   Options Database Key:
    209: . -pc_gamg_graph_symmetrize <bool,default=true> - the flag

    211:   Level: intermediate

    213: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `MatCreateGraph()`,
    214:   `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`
    215: @*/
    216: PetscErrorCode PCGAMGSetGraphSymmetrize(PC pc, PetscBool b)
    217: {
    218:   PetscFunctionBegin;
    221:   PetscTryMethod(pc, "PCGAMGSetGraphSymmetrize_C", (PC, PetscBool), (pc, b));
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

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

    231:   PetscFunctionBegin;
    232:   pc_gamg_agg->aggressive_coarsening_levels = n;
    233:   PetscFunctionReturn(PETSC_SUCCESS);
    234: }

    236: static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n)
    237: {
    238:   PC_MG       *mg          = (PC_MG *)pc->data;
    239:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    240:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

    242:   PetscFunctionBegin;
    243:   pc_gamg_agg->aggressive_mis_k = n;
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }

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

    253:   PetscFunctionBegin;
    254:   pc_gamg_agg->use_aggressive_square_graph = b;
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

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

    264:   PetscFunctionBegin;
    265:   pc_gamg_agg->use_low_mem_filter = b;
    266:   PetscFunctionReturn(PETSC_SUCCESS);
    267: }

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

    275:   PetscFunctionBegin;
    276:   pc_gamg_agg->graph_symmetrize = b;
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b)
    281: {
    282:   PC_MG       *mg          = (PC_MG *)pc->data;
    283:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    284:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

    286:   PetscFunctionBegin;
    287:   pc_gamg_agg->use_minimum_degree_ordering = b;
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    289: }

    291: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject)
    292: {
    293:   PC_MG       *mg          = (PC_MG *)pc->data;
    294:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    295:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
    296:   PetscBool    n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph;
    297:   PetscInt     nsq_graph_old = 0;

    299:   PetscFunctionBegin;
    300:   PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options");
    301:   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));
    302:   // aggressive coarsening logic with deprecated -pc_gamg_square_graph
    303:   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));
    304:   if (!n_aggressive_flg)
    305:     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));
    306:   PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph $ (A^T A)$ for aggressive coarsening, if false, MIS-k (k=2) is used, see PCGAMGMISkSetAggressive()", "PCGAMGSetAggressiveSquareGraph", new_sqr_graph, &pc_gamg_agg->use_aggressive_square_graph, &new_sq_provided));
    307:   if (!new_sq_provided && old_sq_provided) {
    308:     pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero
    309:     pc_gamg_agg->use_aggressive_square_graph  = PETSC_TRUE;
    310:   }
    311:   if (new_sq_provided && old_sq_provided)
    312:     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));
    313:   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));
    314:   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));
    315:   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));
    316:   PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL));
    317:   PetscOptionsHeadEnd();
    318:   PetscFunctionReturn(PETSC_SUCCESS);
    319: }

    321: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
    322: {
    323:   PC_MG       *mg          = (PC_MG *)pc->data;
    324:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    325:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

    327:   PetscFunctionBegin;
    328:   PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
    329:   PetscCall(PetscFree(pc_gamg->subctx));
    330:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL));
    331:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL));
    332:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL));
    333:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL));
    334:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL));
    335:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL));
    336:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL));
    337:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
    338:   PetscFunctionReturn(PETSC_SUCCESS);
    339: }

    341: /*
    342:    PCSetCoordinates_AGG

    344:    Collective

    346:    Input Parameter:
    347:    . pc - the preconditioner context
    348:    . ndm - dimension of data (used for dof/vertex for Stokes)
    349:    . a_nloc - number of vertices local
    350:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
    351: */

    353: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
    354: {
    355:   PC_MG   *mg      = (PC_MG *)pc->data;
    356:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
    357:   PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf;
    358:   Mat      mat = pc->pmat;

    360:   PetscFunctionBegin;
    363:   nloc = a_nloc;

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

    377:   if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
    378:     PetscCall(PetscFree(pc_gamg->data));
    379:     PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data));
    380:   }
    381:   /* copy data in - column-oriented */
    382:   for (kk = 0; kk < nloc; kk++) {
    383:     const PetscInt M    = nloc * pc_gamg->data_cell_rows; /* stride into data */
    384:     PetscReal     *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */

    386:     if (pc_gamg->data_cell_cols == 1) *data = 1.0;
    387:     else {
    388:       /* translational modes */
    389:       for (ii = 0; ii < ndatarows; ii++) {
    390:         for (jj = 0; jj < ndatarows; jj++) {
    391:           if (ii == jj) data[ii * M + jj] = 1.0;
    392:           else data[ii * M + jj] = 0.0;
    393:         }
    394:       }

    396:       /* rotational modes */
    397:       if (coords) {
    398:         if (ndm == 2) {
    399:           data += 2 * M;
    400:           data[0] = -coords[2 * kk + 1];
    401:           data[1] = coords[2 * kk];
    402:         } else {
    403:           data += 3 * M;
    404:           data[0]         = 0.0;
    405:           data[M + 0]     = coords[3 * kk + 2];
    406:           data[2 * M + 0] = -coords[3 * kk + 1];
    407:           data[1]         = -coords[3 * kk + 2];
    408:           data[M + 1]     = 0.0;
    409:           data[2 * M + 1] = coords[3 * kk];
    410:           data[2]         = coords[3 * kk + 1];
    411:           data[M + 2]     = -coords[3 * kk];
    412:           data[2 * M + 2] = 0.0;
    413:         }
    414:       }
    415:     }
    416:   }
    417:   pc_gamg->data_sz = arrsz;
    418:   PetscFunctionReturn(PETSC_SUCCESS);
    419: }

    421: /*
    422:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
    423:       Looks in Mat for near null space.
    424:       Does not work for Stokes

    426:   Input Parameter:
    427:    . pc -
    428:    . a_A - matrix to get (near) null space out of.
    429: */
    430: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
    431: {
    432:   PC_MG       *mg      = (PC_MG *)pc->data;
    433:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
    434:   MatNullSpace mnull;

    436:   PetscFunctionBegin;
    437:   PetscCall(MatGetNearNullSpace(a_A, &mnull));
    438:   if (!mnull) {
    439:     DM dm;

    441:     PetscCall(PCGetDM(pc, &dm));
    442:     if (!dm) PetscCall(MatGetDM(a_A, &dm));
    443:     if (dm) {
    444:       PetscObject deformation;
    445:       PetscInt    Nf;

    447:       PetscCall(DMGetNumFields(dm, &Nf));
    448:       if (Nf) {
    449:         PetscCall(DMGetField(dm, 0, NULL, &deformation));
    450:         if (deformation) {
    451:           PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull));
    452:           if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull));
    453:         }
    454:       }
    455:     }
    456:   }

    458:   if (!mnull) {
    459:     PetscInt bs, NN, MM;

    461:     PetscCall(MatGetBlockSize(a_A, &bs));
    462:     PetscCall(MatGetLocalSize(a_A, &MM, &NN));
    463:     PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs);
    464:     PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL));
    465:   } else {
    466:     PetscReal         *nullvec;
    467:     PetscBool          has_const;
    468:     PetscInt           i, j, mlocal, nvec, bs;
    469:     const Vec         *vecs;
    470:     const PetscScalar *v;

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

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

    498:   Input Parameter:
    499:    . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices
    500:    . bs - row block size
    501:    . nSAvec - column bs of new P
    502:    . my0crs - global index of start of locals
    503:    . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
    504:    . data_in[data_stride*nSAvec] - local data on fine grid
    505:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'

    507:   Output Parameter:
    508:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
    509:    . a_Prol - prolongation operator
    510: */
    511: 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)
    512: {
    513:   PetscInt      Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride;
    514:   MPI_Comm      comm;
    515:   PetscReal    *out_data;
    516:   PetscCDIntNd *pos;
    517:   PetscHMapI    fgid_flid;

    519:   PetscFunctionBegin;
    520:   PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm));
    521:   PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend));
    522:   nloc = (Iend - Istart) / bs;
    523:   my0  = Istart / bs;
    524:   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);
    525:   Iend /= bs;
    526:   nghosts = data_stride / bs - nloc;

    528:   PetscCall(PetscHMapICreateWithSize(2 * nghosts + 1, &fgid_flid));

    530:   for (kk = 0; kk < nghosts; kk++) PetscCall(PetscHMapISet(fgid_flid, flid_fgid[nloc + kk], nloc + kk));

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

    536:     PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise));
    537:     if (!ise) nSelected++;
    538:   }
    539:   PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj));
    540:   PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT "  != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs);
    541:   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);

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

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

    550:   /* find points and set prolongation */
    551:   minsz = 100;
    552:   for (mm = clid = 0; mm < nloc; mm++) {
    553:     PetscCall(PetscCDCountAt(agg_llists, mm, &jj));
    554:     if (jj > 0) {
    555:       const PetscInt lid = mm, cgid = my0crs + clid;
    556:       PetscInt       cids[100]; /* max bs */
    557:       PetscBLASInt   asz, M, N, INFO;
    558:       PetscBLASInt   Mdata, LDA, LWORK;
    559:       PetscScalar   *qqc, *qqr, *TAU, *WORK;
    560:       PetscInt      *fids;
    561:       PetscReal     *data;

    563:       PetscCall(PetscBLASIntCast(jj, &asz));
    564:       PetscCall(PetscBLASIntCast(asz * bs, &M));
    565:       PetscCall(PetscBLASIntCast(nSAvec, &N));
    566:       PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata));
    567:       PetscCall(PetscBLASIntCast(Mdata, &LDA));
    568:       PetscCall(PetscBLASIntCast(N * bs, &LWORK));
    569:       /* count agg */
    570:       if (asz < minsz) minsz = asz;

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

    575:       aggID = 0;
    576:       PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
    577:       while (pos) {
    578:         PetscInt gid1;

    580:         PetscCall(PetscCDIntNdGetID(pos, &gid1));
    581:         PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));

    583:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
    584:         else {
    585:           PetscCall(PetscHMapIGet(fgid_flid, gid1, &flid));
    586:           PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table");
    587:         }
    588:         /* copy in B_i matrix - column-oriented */
    589:         data = &data_in[flid * bs];
    590:         for (ii = 0; ii < bs; ii++) {
    591:           for (jj = 0; jj < N; jj++) {
    592:             PetscReal d = data[jj * data_stride + ii];

    594:             qqc[jj * Mdata + aggID * bs + ii] = d;
    595:           }
    596:         }
    597:         /* set fine IDs */
    598:         for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk;
    599:         aggID++;
    600:       }

    602:       /* pad with zeros */
    603:       for (ii = asz * bs; ii < Mdata; ii++) {
    604:         for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0;
    605:       }

    607:       /* QR */
    608:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    609:       PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
    610:       PetscCall(PetscFPTrapPop());
    611:       PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error");
    612:       /* get R - column-oriented - output B_{i+1} */
    613:       {
    614:         PetscReal *data = &out_data[clid * nSAvec];

    616:         for (jj = 0; jj < nSAvec; jj++) {
    617:           for (ii = 0; ii < nSAvec; ii++) {
    618:             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);
    619:             if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]);
    620:             else data[jj * out_data_stride + ii] = 0.;
    621:           }
    622:         }
    623:       }

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

    629:       for (ii = 0; ii < M; ii++) {
    630:         for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii];
    631:       }

    633:       /* add diagonal block of P0 */
    634:       for (kk = 0; kk < N; kk++) cids[kk] = N * cgid + kk; /* global col IDs in P0 */
    635:       PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES));
    636:       PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids));
    637:       clid++;
    638:     } /* coarse agg */
    639:   } /* for all fine nodes */
    640:   PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY));
    641:   PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY));
    642:   PetscCall(PetscHMapIDestroy(&fgid_flid));
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }

    646: static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer)
    647: {
    648:   PC_MG       *mg          = (PC_MG *)pc->data;
    649:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    650:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;

    652:   PetscFunctionBegin;
    653:   PetscCall(PetscViewerASCIIPrintf(viewer, "      AGG specific options\n"));
    654:   PetscCall(PetscViewerASCIIPrintf(viewer, "        Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels));
    655:   if (pc_gamg_agg->aggressive_coarsening_levels > 0) {
    656:     PetscCall(PetscViewerASCIIPrintf(viewer, "        %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph"));
    657:     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));
    658:   }
    659:   PetscCall(PetscViewerASCIIPushTab(viewer));
    660:   PetscCall(PetscViewerASCIIPushTab(viewer));
    661:   PetscCall(PetscViewerASCIIPushTab(viewer));
    662:   PetscCall(PetscViewerASCIIPushTab(viewer));
    663:   if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer));
    664:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n"));
    665:   PetscCall(PetscViewerASCIIPopTab(viewer));
    666:   PetscCall(PetscViewerASCIIPopTab(viewer));
    667:   PetscCall(PetscViewerASCIIPopTab(viewer));
    668:   PetscCall(PetscViewerASCIIPopTab(viewer));
    669:   PetscCall(PetscViewerASCIIPrintf(viewer, "        Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths));
    670:   PetscFunctionReturn(PETSC_SUCCESS);
    671: }

    673: static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat)
    674: {
    675:   PC_MG          *mg          = (PC_MG *)pc->data;
    676:   PC_GAMG        *pc_gamg     = (PC_GAMG *)mg->innerctx;
    677:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
    678:   const PetscReal vfilter     = pc_gamg->threshold[pc_gamg->current_level];
    679:   PetscBool       ishem, ismis;
    680:   const char     *prefix;
    681:   MatInfo         info0, info1;
    682:   PetscInt        bs;

    684:   PetscFunctionBegin;
    685:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
    686:   /* 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 */
    687:   /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */
    688:   PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs));
    689:   PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs));
    690:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
    691:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
    692:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_"));
    693:   PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs));
    694:   PetscCall(MatGetBlockSize(Amat, &bs));
    695:   // check for valid indices wrt bs
    696:   for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) {
    697:     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",
    698:                pc_gamg_agg->crs->strength_index[ii], bs);
    699:   }
    700:   PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem));
    701:   if (ishem) {
    702:     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));
    703:     pc_gamg_agg->aggressive_coarsening_levels = 0;                                         // aggressive and HEM does not make sense
    704:     PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage
    705:     PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter));                          // for code coverage
    706:   } else {
    707:     PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis));
    708:     if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) {
    709:       PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n"));
    710:       pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE;
    711:     }
    712:   }
    713:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
    714:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
    715:   PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */

    717:   if (ishem || pc_gamg_agg->use_low_mem_filter) {
    718:     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));
    719:   } else {
    720:     // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive)
    721:     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));
    722:     if (vfilter >= 0) {
    723:       PetscInt           Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc;
    724:       Mat                tGmat, Gmat = *a_Gmat;
    725:       MPI_Comm           comm;
    726:       const PetscScalar *vals;
    727:       const PetscInt    *idx;
    728:       PetscInt          *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0;
    729:       MatScalar         *AA; // this is checked in graph
    730:       PetscBool          isseqaij;
    731:       Mat                a, b, c;
    732:       MatType            jtype;

    734:       PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm));
    735:       PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij));
    736:       PetscCall(MatGetType(Gmat, &jtype));
    737:       PetscCall(MatCreate(comm, &tGmat));
    738:       PetscCall(MatSetType(tGmat, jtype));

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

    744:       // global sizes
    745:       PetscCall(MatGetSize(Gmat, &MM, &NN));
    746:       PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend));
    747:       nloc = Iend - Istart;
    748:       PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz));
    749:       if (isseqaij) {
    750:         a = Gmat;
    751:         b = NULL;
    752:       } else {
    753:         Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data;

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

    790:               nnz1++;
    791:               if (c != a) cid = garray[idx[jj]];
    792:               AA[ncol_row] = vals[jj];
    793:               AJ[ncol_row] = cid;
    794:               ncol_row++;
    795:             }
    796:           }
    797:           PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals));
    798:           PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES));
    799:         }
    800:       }
    801:       PetscCall(PetscFree2(AA, AJ));
    802:       PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY));
    803:       PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY));
    804:       PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */
    805:       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));
    806:       PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view"));
    807:       PetscCall(MatDestroy(&Gmat));
    808:       *a_Gmat = tGmat;
    809:     }
    810:   }

    812:   PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */
    813:   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));
    814:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0));
    815:   PetscFunctionReturn(PETSC_SUCCESS);
    816: }

    818: typedef PetscInt    NState;
    819: static const NState NOT_DONE = -2;
    820: static const NState DELETED  = -1;
    821: static const NState REMOVED  = -3;
    822: #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED)

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

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

    848:   PetscFunctionBegin;
    849:   PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm));
    850:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    851:   PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend));

    853:   /* get submatrices */
    854:   PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI));
    855:   PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no"));
    856:   PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1));
    857:   for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
    858:   if (isMPI) {
    859:     /* grab matrix objects */
    860:     mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data;
    861:     mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data;
    862:     matA_1   = (Mat_SeqAIJ *)mpimat_1->A->data;
    863:     matB_1   = (Mat_SeqAIJ *)mpimat_1->B->data;

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

    870:       PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc);
    871:       if (lid != -1) lid_cprowID_1[lid] = ix;
    872:     }
    873:   } else {
    874:     PetscBool isAIJ;

    876:     PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ));
    877:     PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix.");
    878:     matA_1 = (Mat_SeqAIJ *)Gmat_1->data;
    879:   }
    880:   if (nloc > 0) PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???");
    881:   /* get state of locals and selected gid for deleted */
    882:   for (lid = 0; lid < nloc; lid++) {
    883:     lid_parent_gid[lid] = -1.0;
    884:     lid_state[lid]      = DELETED;
    885:   }

    887:   /* set lid_state */
    888:   for (lid = 0; lid < nloc; lid++) {
    889:     PetscCDIntNd *pos;

    891:     PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
    892:     if (pos) {
    893:       PetscInt gid1;

    895:       PetscCall(PetscCDIntNdGetID(pos, &gid1));
    896:       PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0);
    897:       lid_state[lid] = gid1;
    898:     }
    899:   }

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

    905:     if (IS_SELECTED(state)) {
    906:       PetscCDIntNd *pos;

    908:       PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
    909:       while (pos) {
    910:         PetscInt gid1;

    912:         PetscCall(PetscCDIntNdGetID(pos, &gid1));
    913:         PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
    914:         if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0);
    915:       }
    916:     }
    917:   }
    918:   /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
    919:   if (isMPI) {
    920:     Vec tempVec;

    922:     /* get 'cpcol_1_state' */
    923:     PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL));
    924:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
    925:       PetscScalar v = (PetscScalar)lid_state[kk];

    927:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
    928:     }
    929:     PetscCall(VecAssemblyBegin(tempVec));
    930:     PetscCall(VecAssemblyEnd(tempVec));
    931:     PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
    932:     PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD));
    933:     PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state));
    934:     /* get 'cpcol_2_state' */
    935:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
    936:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD));
    937:     PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state));
    938:     /* get 'cpcol_2_par_orig' */
    939:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
    940:       PetscScalar v = lid_parent_gid[kk];

    942:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
    943:     }
    944:     PetscCall(VecAssemblyBegin(tempVec));
    945:     PetscCall(VecAssemblyEnd(tempVec));
    946:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2));
    947:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
    948:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD));
    949:     PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig));

    951:     PetscCall(VecDestroy(&tempVec));
    952:   } /* ismpi */
    953:   for (lid = 0; lid < nloc; lid++) {
    954:     NState state = lid_state[lid];

    956:     if (IS_SELECTED(state)) {
    957:       /* steal locals */
    958:       ii  = matA_1->i;
    959:       n   = ii[lid + 1] - ii[lid];
    960:       idx = matA_1->j + ii[lid];
    961:       for (j = 0; j < n; j++) {
    962:         PetscInt lidj   = idx[j], sgid;
    963:         NState   statej = lid_state[lidj];

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

    971:             /* looking for local from local so id_llist_2 works */
    972:             PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos));
    973:             while (pos) {
    974:               PetscInt gid;

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

    1000:       /* see if I have a selected ghost neighbor that will steal me */
    1001:       if ((ix = lid_cprowID_1[lid]) != -1) {
    1002:         ii  = matB_1->compressedrow.i;
    1003:         n   = ii[ix + 1] - ii[ix];
    1004:         idx = matB_1->j + ii[ix];
    1005:         for (j = 0; j < n; j++) {
    1006:           PetscInt cpid   = idx[j];
    1007:           NState   statej = (NState)PetscRealPart(cpcol_1_state[cpid]);

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

    1015:               /* remove from 'oldslidj' list */
    1016:               PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos));
    1017:               while (pos) {
    1018:                 PetscInt gid;

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

    1044:   if (isMPI) {
    1045:     PetscScalar *cpcol_2_parent, *cpcol_2_gid;
    1046:     Vec          tempVec, ghostgids2, ghostparents2;
    1047:     PetscInt     cpid, nghost_2;
    1048:     PetscHMapI   gid_cpid;

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

    1053:     /* get 'cpcol_2_parent' */
    1054:     for (kk = 0, j = my0; kk < nloc; kk++, j++) PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES));
    1055:     PetscCall(VecAssemblyBegin(tempVec));
    1056:     PetscCall(VecAssemblyEnd(tempVec));
    1057:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2));
    1058:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
    1059:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD));
    1060:     PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent));

    1062:     /* get 'cpcol_2_gid' */
    1063:     for (kk = 0, j = my0; kk < nloc; kk++, j++) {
    1064:       PetscScalar v = (PetscScalar)j;

    1066:       PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES));
    1067:     }
    1068:     PetscCall(VecAssemblyBegin(tempVec));
    1069:     PetscCall(VecAssemblyEnd(tempVec));
    1070:     PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2));
    1071:     PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
    1072:     PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD));
    1073:     PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid));
    1074:     PetscCall(VecDestroy(&tempVec));

    1076:     /* look for deleted ghosts and add to table */
    1077:     PetscCall(PetscHMapICreateWithSize(2 * nghost_2 + 1, &gid_cpid));
    1078:     for (cpid = 0; cpid < nghost_2; cpid++) {
    1079:       NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);

    1081:       if (state == DELETED) {
    1082:         PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
    1083:         PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);

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

    1088:           PetscCall(PetscHMapISet(gid_cpid, gid, cpid));
    1089:         }
    1090:       }
    1091:     }

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

    1097:       if (IS_SELECTED(state)) {
    1098:         PetscCDIntNd *pos, *last = NULL;

    1100:         /* look for deleted ghosts and see if they moved */
    1101:         PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos));
    1102:         while (pos) {
    1103:           PetscInt gid;

    1105:           PetscCall(PetscCDIntNdGetID(pos, &gid));
    1106:           if (gid < my0 || gid >= Iend) {
    1107:             PetscCall(PetscHMapIGet(gid_cpid, gid, &cpid));
    1108:             if (cpid != -1) {
    1109:               /* a moved ghost - */
    1110:               /* id_llist_2[lastid] = id_llist_2[flid];    /\* remove 'flid' from list *\/ */
    1111:               PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last));
    1112:             } else last = pos;
    1113:           } else last = pos;

    1115:           PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos));
    1116:         } /* loop over list of deleted */
    1117:       } /* selected */
    1118:     }
    1119:     PetscCall(PetscHMapIDestroy(&gid_cpid));

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

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

    1130:         /* search for this gid to see if I have it */
    1131:         PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos));
    1132:         while (pos) {
    1133:           PetscInt gidj;

    1135:           PetscCall(PetscCDIntNdGetID(pos, &gidj));
    1136:           PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos));

    1138:           if (gidj == gid) {
    1139:             hav = 1;
    1140:             break;
    1141:           }
    1142:         }
    1143:         if (hav != 1) {
    1144:           /* insert 'flidj' into head of llist */
    1145:           PetscCall(PetscCDAppendID(aggs_2, slid_new, gid));
    1146:         }
    1147:       }
    1148:     }
    1149:     PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state));
    1150:     PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state));
    1151:     PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent));
    1152:     PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid));
    1153:     PetscCall(VecDestroy(&ghostgids2));
    1154:     PetscCall(VecDestroy(&ghostparents2));
    1155:     PetscCall(VecDestroy(&ghost_par_orig2));
    1156:   }
    1157:   PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1));
    1158:   PetscFunctionReturn(PETSC_SUCCESS);
    1159: }

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

    1165:   Input Parameter:
    1166:    . a_pc - this

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

    1171:   Output Parameter:
    1172:    . agg_lists - list of aggregates

    1174: */
    1175: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists)
    1176: {
    1177:   PC_MG       *mg          = (PC_MG *)a_pc->data;
    1178:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    1179:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
    1180:   Mat          Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */
    1181:   IS           perm;
    1182:   PetscInt     Istart, Iend, Ii, nloc, bs, nn;
    1183:   PetscInt    *permute, *degree;
    1184:   PetscBool   *bIndexSet;
    1185:   PetscReal    hashfact;
    1186:   PetscInt     iSwapIndex;
    1187:   PetscRandom  random;
    1188:   MPI_Comm     comm;

    1190:   PetscFunctionBegin;
    1191:   PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm));
    1192:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
    1193:   PetscCall(MatGetLocalSize(Gmat1, &nn, NULL));
    1194:   PetscCall(MatGetBlockSize(Gmat1, &bs));
    1195:   PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs);
    1196:   nloc = nn / bs;
    1197:   /* get MIS aggs - randomize */
    1198:   PetscCall(PetscMalloc2(nloc, &permute, nloc, &degree));
    1199:   PetscCall(PetscCalloc1(nloc, &bIndexSet));
    1200:   for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii;
    1201:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random));
    1202:   PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend));
    1203:   for (Ii = 0; Ii < nloc; Ii++) {
    1204:     PetscInt nc;

    1206:     PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
    1207:     degree[Ii] = nc;
    1208:     PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL));
    1209:   }
    1210:   for (Ii = 0; Ii < nloc; Ii++) {
    1211:     PetscCall(PetscRandomGetValueReal(random, &hashfact));
    1212:     iSwapIndex = (PetscInt)(hashfact * nloc) % nloc;
    1213:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
    1214:       PetscInt iTemp = permute[iSwapIndex];

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

    1240:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix));
    1241:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix));
    1242:     PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS
    1243:   }
    1244:   PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2));
    1245:   PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE));
    1246:   PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm));
    1247:   PetscCall(MatCoarsenApply(pc_gamg_agg->crs));
    1248:   PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */

    1250:   PetscCall(ISDestroy(&perm));
    1251:   PetscCall(PetscFree2(permute, degree));
    1252:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0));

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

    1257:     PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists));
    1258:     PetscCall(MatDestroy(&Gmat1));
    1259:     *a_Gmat1 = Gmat2;                          /* output */
    1260:     PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */
    1261:   }
    1262:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0));
    1263:   PetscFunctionReturn(PETSC_SUCCESS);
    1264: }

    1266: /*
    1267:  PCGAMGConstructProlongator_AGG

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

    1290:   PetscFunctionBegin;
    1291:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
    1292:   PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1");
    1293:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
    1294:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1295:   PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend));
    1296:   PetscCall(MatGetBlockSize(Amat, &bs));
    1297:   nloc = (Iend - Istart) / bs;
    1298:   my0  = Istart / bs;
    1299:   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);
    1300:   PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1

    1302:   /* get 'nLocalSelected' */
    1303:   for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) {
    1304:     PetscBool ise;

    1306:     /* filter out singletons 0 or 1? */
    1307:     PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise));
    1308:     if (!ise) nLocalSelected++;
    1309:   }

    1311:   /* create prolongator, create P matrix */
    1312:   PetscCall(MatGetType(Amat, &mtype));
    1313:   PetscCall(MatCreate(comm, &Prol));
    1314:   PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE));
    1315:   PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes?
    1316:   PetscCall(MatSetType(Prol, mtype));
    1317: #if PetscDefined(HAVE_DEVICE)
    1318:   PetscBool flg;
    1319:   PetscCall(MatBoundToCPU(Amat, &flg));
    1320:   PetscCall(MatBindToCPU(Prol, flg));
    1321:   if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
    1322: #endif
    1323:   PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL));
    1324:   PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL));

    1326:   /* can get all points "removed" */
    1327:   PetscCall(MatGetSize(Prol, &kk, &ii));
    1328:   if (!ii) {
    1329:     PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix));
    1330:     PetscCall(MatDestroy(&Prol));
    1331:     *a_P_out = NULL; /* out */
    1332:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
    1333:     PetscFunctionReturn(PETSC_SUCCESS);
    1334:   }
    1335:   PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs));
    1336:   PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk));

    1338:   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);
    1339:   myCrs0 = myCrs0 / col_bs;
    1340:   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);

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

    1347:     PetscCall(PetscMalloc1(nloc, &tmp_ldata));
    1348:     for (jj = 0; jj < col_bs; jj++) {
    1349:       for (kk = 0; kk < bs; kk++) {
    1350:         PetscInt         ii, stride;
    1351:         const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk);

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

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

    1357:         if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */
    1358:           PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost));
    1359:           nbnodes = bs * stride;
    1360:         }
    1361:         tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk);
    1362:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
    1363:         PetscCall(PetscFree(tmp_gdata));
    1364:       }
    1365:     }
    1366:     PetscCall(PetscFree(tmp_ldata));
    1367:   } else {
    1368:     nbnodes      = bs * nloc;
    1369:     data_w_ghost = pc_gamg->data;
    1370:   }

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

    1377:     PetscCall(PetscMalloc1(nloc, &fid_glid_loc));
    1378:     for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk);
    1379:     PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata));
    1380:     PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */
    1381:     for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
    1382:     PetscCall(PetscFree(fiddata));

    1384:     PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs);
    1385:     PetscCall(PetscFree(fid_glid_loc));
    1386:   } else {
    1387:     PetscCall(PetscMalloc1(nloc, &flid_fgid));
    1388:     for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk;
    1389:   }
    1390:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0));
    1391:   /* get P0 */
    1392:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
    1393:   {
    1394:     PetscReal *data_out = NULL;

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

    1399:     pc_gamg->data           = data_out;
    1400:     pc_gamg->data_cell_rows = col_bs;
    1401:     pc_gamg->data_sz        = col_bs * col_bs * nLocalSelected;
    1402:   }
    1403:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0));
    1404:   if (size > 1) PetscCall(PetscFree(data_w_ghost));
    1405:   PetscCall(PetscFree(flid_fgid));

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

    1410:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0));
    1411:   PetscFunctionReturn(PETSC_SUCCESS);
    1412: }

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

    1417:   Input Parameter:
    1418:    . pc - this
    1419:    . Amat - matrix on this fine level
    1420:  In/Output Parameter:
    1421:    . a_P - prolongation operator to the next level
    1422: */
    1423: static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P)
    1424: {
    1425:   PC_MG       *mg          = (PC_MG *)pc->data;
    1426:   PC_GAMG     *pc_gamg     = (PC_GAMG *)mg->innerctx;
    1427:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx;
    1428:   PetscInt     jj;
    1429:   Mat          Prol = *a_P;
    1430:   MPI_Comm     comm;
    1431:   KSP          eksp;
    1432:   Vec          bb, xx;
    1433:   PC           epc;
    1434:   PetscReal    alpha, emax, emin;

    1436:   PetscFunctionBegin;
    1437:   PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
    1438:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));

    1440:   /* compute maximum singular value of operator to be used in smoother */
    1441:   if (0 < pc_gamg_agg->nsmooths) {
    1442:     /* get eigen estimates */
    1443:     if (pc_gamg->emax > 0) {
    1444:       emin = pc_gamg->emin;
    1445:       emax = pc_gamg->emax;
    1446:     } else {
    1447:       const char *prefix;

    1449:       PetscCall(MatCreateVecs(Amat, &bb, NULL));
    1450:       PetscCall(MatCreateVecs(Amat, &xx, NULL));
    1451:       PetscCall(KSPSetNoisy_Private(Amat, bb));

    1453:       PetscCall(KSPCreate(comm, &eksp));
    1454:       PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel));
    1455:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
    1456:       PetscCall(KSPSetOptionsPrefix(eksp, prefix));
    1457:       PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_"));
    1458:       {
    1459:         PetscBool isset, sflg;

    1461:         PetscCall(MatIsSPDKnown(Amat, &isset, &sflg));
    1462:         if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG));
    1463:       }
    1464:       PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure));
    1465:       PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE));

    1467:       PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE));
    1468:       PetscCall(KSPSetOperators(eksp, Amat, Amat));

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

    1473:       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

    1475:       PetscCall(KSPSetFromOptions(eksp));
    1476:       PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE));
    1477:       PetscCall(KSPSolve(eksp, bb, xx));
    1478:       PetscCall(KSPCheckSolve(eksp, pc, xx));

    1480:       PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin));
    1481:       PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI));
    1482:       PetscCall(VecDestroy(&xx));
    1483:       PetscCall(VecDestroy(&bb));
    1484:       PetscCall(KSPDestroy(&eksp));
    1485:     }
    1486:     if (pc_gamg->use_sa_esteig) {
    1487:       mg->min_eigen_DinvA[pc_gamg->current_level] = emin;
    1488:       mg->max_eigen_DinvA[pc_gamg->current_level] = emax;
    1489:       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));
    1490:     } else {
    1491:       mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
    1492:       mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
    1493:     }
    1494:   } else {
    1495:     mg->min_eigen_DinvA[pc_gamg->current_level] = 0;
    1496:     mg->max_eigen_DinvA[pc_gamg->current_level] = 0;
    1497:   }

    1499:   /* smooth P0 */
    1500:   if (pc_gamg_agg->nsmooths > 0) {
    1501:     Vec diag;

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

    1506:     PetscCall(MatCreateVecs(Amat, &diag, NULL));
    1507:     PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */
    1508:     PetscCall(VecReciprocal(diag));

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

    1513:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
    1514:       /*
    1515:         Smooth aggregation on the prolongator

    1517:         P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1}
    1518:       */
    1519:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
    1520:       PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat));
    1521:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0));
    1522:       PetscCall(MatProductClear(tMat));
    1523:       PetscCall(MatDiagonalScale(tMat, diag, NULL));

    1525:       /* TODO: Document the 1.4 and don't hardwire it in this routine */
    1526:       alpha = -1.4 / emax;
    1527:       PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN));
    1528:       PetscCall(MatDestroy(&Prol));
    1529:       Prol = tMat;
    1530:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0));
    1531:     }
    1532:     PetscCall(VecDestroy(&diag));
    1533:   }
    1534:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0));
    1535:   PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation"));
    1536:   *a_P = Prol;
    1537:   PetscFunctionReturn(PETSC_SUCCESS);
    1538: }

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

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

    1551:   Level: intermediate

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

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

    1560: .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
    1561:           `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`,
    1562:           `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`,
    1563:           `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`,
    1564:           `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
    1565: M*/
    1566: PetscErrorCode PCCreateGAMG_AGG(PC pc)
    1567: {
    1568:   PC_MG       *mg      = (PC_MG *)pc->data;
    1569:   PC_GAMG     *pc_gamg = (PC_GAMG *)mg->innerctx;
    1570:   PC_GAMG_AGG *pc_gamg_agg;

    1572:   PetscFunctionBegin;
    1573:   /* create sub context for SA */
    1574:   PetscCall(PetscNew(&pc_gamg_agg));
    1575:   pc_gamg->subctx = pc_gamg_agg;

    1577:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
    1578:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
    1579:   /* reset does not do anything; setup not virtual */

    1581:   /* set internal function pointers */
    1582:   pc_gamg->ops->creategraph       = PCGAMGCreateGraph_AGG;
    1583:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_AGG;
    1584:   pc_gamg->ops->prolongator       = PCGAMGConstructProlongator_AGG;
    1585:   pc_gamg->ops->optprolongator    = PCGAMGOptimizeProlongator_AGG;
    1586:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;
    1587:   pc_gamg->ops->view              = PCView_GAMG_AGG;

    1589:   pc_gamg_agg->nsmooths                     = 1;
    1590:   pc_gamg_agg->aggressive_coarsening_levels = 1;
    1591:   pc_gamg_agg->use_aggressive_square_graph  = PETSC_TRUE;
    1592:   pc_gamg_agg->use_minimum_degree_ordering  = PETSC_FALSE;
    1593:   pc_gamg_agg->use_low_mem_filter           = PETSC_FALSE;
    1594:   pc_gamg_agg->aggressive_mis_k             = 2;
    1595:   pc_gamg_agg->graph_symmetrize             = PETSC_TRUE;

    1597:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG));
    1598:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG));
    1599:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG));
    1600:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG));
    1601:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG));
    1602:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG));
    1603:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG));
    1604:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG));
    1605:   PetscFunctionReturn(PETSC_SUCCESS);
    1606: }