Actual source code: gamg.c

      1: /*
      2:  GAMG geometric-algebric multigrid PC - Mark Adams 2011
      3:  */
      4: #include <../src/ksp/pc/impls/gamg/gamg.h>
      5: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>

      7: #if defined(PETSC_HAVE_CUDA)
      8:   #include <cuda_runtime.h>
      9: #endif

     11: #if defined(PETSC_HAVE_HIP)
     12:   #include <hip/hip_runtime.h>
     13: #endif

     15: PetscLogEvent petsc_gamg_setup_events[GAMG_NUM_SET];
     16: PetscLogEvent petsc_gamg_setup_matmat_events[PETSC_MG_MAXLEVELS][3];

     18: // #define GAMG_STAGES
     19: #if defined(GAMG_STAGES)
     20: static PetscLogStage gamg_stages[PETSC_MG_MAXLEVELS];
     21: #endif

     23: static PetscFunctionList GAMGList = NULL;
     24: static PetscBool         PCGAMGPackageInitialized;

     26: static PetscErrorCode PCReset_GAMG(PC pc)
     27: {
     28:   PC_MG   *mg      = (PC_MG *)pc->data;
     29:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

     31:   PetscFunctionBegin;
     32:   PetscCall(PetscFree(pc_gamg->data));
     33:   pc_gamg->data_sz = 0;
     34:   PetscCall(PetscFree(pc_gamg->orig_data));
     35:   for (PetscInt level = 0; level < PETSC_MG_MAXLEVELS; level++) {
     36:     mg->min_eigen_DinvA[level] = 0;
     37:     mg->max_eigen_DinvA[level] = 0;
     38:   }
     39:   pc_gamg->emin = 0;
     40:   pc_gamg->emax = 0;
     41:   PetscCall(PCReset_MG(pc));
     42:   PetscCall(MatCoarsenDestroy(&pc_gamg->asm_crs));
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }

     46: /*
     47:    PCGAMGCreateLevel_GAMG: create coarse op with RAP.  repartition and/or reduce number
     48:      of active processors.

     50:    Input Parameter:
     51:    . pc - parameters + side effect: coarse data in 'pc_gamg->data' and
     52:           'pc_gamg->data_sz' are changed via repartitioning/reduction.
     53:    . Amat_fine - matrix on this fine (k) level
     54:    . cr_bs - coarse block size
     55:    In/Output Parameter:
     56:    . a_P_inout - prolongation operator to the next level (k-->k-1)
     57:    . a_nactive_proc - number of active procs
     58:    Output Parameter:
     59:    . a_Amat_crs - coarse matrix that is created (k-1)
     60: */
     61: static PetscErrorCode PCGAMGCreateLevel_GAMG(PC pc, Mat Amat_fine, PetscInt cr_bs, Mat *a_P_inout, Mat *a_Amat_crs, PetscMPIInt *a_nactive_proc, IS *Pcolumnperm, PetscBool is_last)
     62: {
     63:   PC_MG      *mg      = (PC_MG *)pc->data;
     64:   PC_GAMG    *pc_gamg = (PC_GAMG *)mg->innerctx;
     65:   Mat         Cmat = NULL, Pold = *a_P_inout;
     66:   MPI_Comm    comm;
     67:   PetscMPIInt rank, size, new_size, nactive = *a_nactive_proc;
     68:   PetscInt    ncrs_eq, ncrs, f_bs;

     70:   PetscFunctionBegin;
     71:   PetscCall(PetscObjectGetComm((PetscObject)Amat_fine, &comm));
     72:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     73:   PetscCallMPI(MPI_Comm_size(comm, &size));
     74:   PetscCall(MatGetBlockSize(Amat_fine, &f_bs));

     76:   if (Pcolumnperm) *Pcolumnperm = NULL;

     78:   /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/
     79:   PetscCall(MatGetLocalSize(Pold, NULL, &ncrs_eq));
     80:   if (pc_gamg->data_cell_rows > 0) {
     81:     ncrs = pc_gamg->data_sz / pc_gamg->data_cell_cols / pc_gamg->data_cell_rows;
     82:   } else {
     83:     PetscInt bs;
     84:     PetscCall(MatGetBlockSizes(Pold, NULL, &bs));
     85:     ncrs = ncrs_eq / bs;
     86:   }
     87:   /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */
     88:   if (pc_gamg->level_reduction_factors[pc_gamg->current_level] == 0 && PetscDefined(HAVE_CUDA) && pc_gamg->current_level == 0) { /* 0 turns reducing to 1 process/device on; do for HIP, etc. */
     89: #if defined(PETSC_HAVE_CUDA)
     90:     PetscShmComm pshmcomm;
     91:     PetscMPIInt  locrank;
     92:     MPI_Comm     loccomm;
     93:     PetscInt     s_nnodes, r_nnodes, new_new_size;
     94:     cudaError_t  cerr;
     95:     int          devCount;
     96:     PetscCall(PetscShmCommGet(comm, &pshmcomm));
     97:     PetscCall(PetscShmCommGetMpiShmComm(pshmcomm, &loccomm));
     98:     PetscCallMPI(MPI_Comm_rank(loccomm, &locrank));
     99:     s_nnodes = !locrank;
    100:     PetscCallMPI(MPIU_Allreduce(&s_nnodes, &r_nnodes, 1, MPIU_INT, MPI_SUM, comm));
    101:     PetscCheck((size % r_nnodes) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "odd number of nodes np=%d nnodes%" PetscInt_FMT, size, r_nnodes);
    102:     devCount = 0;
    103:     cerr     = cudaGetDeviceCount(&devCount);
    104:     cudaGetLastError();                         /* Reset the last error */
    105:     if (cerr == cudaSuccess && devCount >= 1) { /* There are devices, else go to heuristic */
    106:       new_new_size = r_nnodes * devCount;
    107:       new_size     = new_new_size;
    108:       PetscCall(PetscInfo(pc, "%s: Fine grid with Cuda. %" PetscInt_FMT " nodes. Change new active set size %d --> %d (devCount=%d #nodes=%" PetscInt_FMT ")\n", ((PetscObject)pc)->prefix, r_nnodes, nactive, new_size, devCount, r_nnodes));
    109:     } else {
    110:       PetscCall(PetscInfo(pc, "%s: With Cuda but no device. Use heuristics.\n", ((PetscObject)pc)->prefix));
    111:       goto HEURISTIC;
    112:     }
    113: #else
    114:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "should not be here");
    115: #endif
    116:   } else if (pc_gamg->level_reduction_factors[pc_gamg->current_level] > 0) {
    117:     if (nactive < pc_gamg->level_reduction_factors[pc_gamg->current_level]) {
    118:       new_size = 1;
    119:       PetscCall(PetscInfo(pc, "%s: reduction factor too small for %d active processes: reduce to one process\n", ((PetscObject)pc)->prefix, new_size));
    120:     } else {
    121:       PetscCheck(nactive % pc_gamg->level_reduction_factors[pc_gamg->current_level] == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "odd number of active process %d wrt reduction factor %" PetscInt_FMT, nactive, pc_gamg->level_reduction_factors[pc_gamg->current_level]);
    122:       new_size = nactive / pc_gamg->level_reduction_factors[pc_gamg->current_level];
    123:       PetscCall(PetscInfo(pc, "%s: Manually setting reduction to %d active processes (%d/%" PetscInt_FMT ")\n", ((PetscObject)pc)->prefix, new_size, nactive, pc_gamg->level_reduction_factors[pc_gamg->current_level]));
    124:     }
    125:   } else if (is_last && !pc_gamg->use_parallel_coarse_grid_solver) {
    126:     new_size = 1;
    127:     PetscCall(PetscInfo(pc, "%s: Force coarsest grid reduction to %d active processes\n", ((PetscObject)pc)->prefix, new_size));
    128:   } else {
    129:     PetscInt ncrs_eq_glob;
    130: #if defined(PETSC_HAVE_CUDA)
    131:   HEURISTIC:
    132: #endif
    133:     PetscCall(MatGetSize(Pold, NULL, &ncrs_eq_glob));
    134:     new_size = (PetscMPIInt)((float)ncrs_eq_glob / (float)pc_gamg->min_eq_proc + 0.5); /* hardwire min. number of eq/proc */
    135:     if (!new_size) new_size = 1;                                                       /* not likely, possible? */
    136:     else if (new_size >= nactive) new_size = nactive;                                  /* no change, rare */
    137:     PetscCall(PetscInfo(pc, "%s: Coarse grid reduction from %d to %d active processes\n", ((PetscObject)pc)->prefix, nactive, new_size));
    138:   }
    139:   if (new_size == nactive) {
    140:     /* output - no repartitioning or reduction - could bail here
    141:        we know that the grid structure can be reused in MatPtAP */
    142:     PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    143:     PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    144:     PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, a_Amat_crs));
    145:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    146:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    147:     if (new_size < size) {
    148:       /* odd case where multiple coarse grids are on one processor or no coarsening ... */
    149:       PetscCall(PetscInfo(pc, "%s: reduced grid using same number of processors (%d) as last grid (use larger coarse grid)\n", ((PetscObject)pc)->prefix, nactive));
    150:       if (pc_gamg->cpu_pin_coarse_grids) {
    151:         PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE));
    152:         PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE));
    153:       }
    154:     }
    155:   } else { /* reduce active processors - we know that the grid structure can NOT be reused in MatPtAP */
    156:     PetscInt *counts, *newproc_idx, ii, jj, kk, strideNew, *tidx, ncrs_new, ncrs_eq_new, nloc_old, expand_factor = 1, rfactor = 1;
    157:     IS        is_eq_newproc, is_eq_num, new_eq_indices;
    158:     PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
    159:     nloc_old = ncrs_eq / cr_bs;
    160:     PetscCheck(ncrs_eq % cr_bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ncrs_eq %" PetscInt_FMT " not divisible by cr_bs %" PetscInt_FMT, ncrs_eq, cr_bs);
    161:     /* get new_size and rfactor */
    162:     if (pc_gamg->layout_type == PCGAMG_LAYOUT_SPREAD || !pc_gamg->repart) {
    163:       /* find factor */
    164:       if (new_size == 1) rfactor = size; /* don't modify */
    165:       else {
    166:         PetscReal best_fact = 0.;
    167:         jj                  = -1;
    168:         for (kk = 1; kk <= size; kk++) {
    169:           if (!(size % kk)) { /* a candidate */
    170:             PetscReal nactpe = (PetscReal)size / (PetscReal)kk, fact = nactpe / (PetscReal)new_size;
    171:             if (fact > 1.0) fact = 1. / fact; /* keep fact < 1 */
    172:             if (fact > best_fact) {
    173:               best_fact = fact;
    174:               jj        = kk;
    175:             }
    176:           }
    177:         }
    178:         if (jj != -1) rfactor = jj;
    179:         else rfactor = 1; /* a prime */
    180:         if (pc_gamg->layout_type == PCGAMG_LAYOUT_COMPACT) expand_factor = 1;
    181:         else expand_factor = rfactor;
    182:       }
    183:       new_size = size / rfactor; /* make new size one that is factor */
    184:       if (new_size == nactive) { /* no repartitioning or reduction, bail out because nested here (rare) */
    185:         PetscCall(PetscInfo(pc, "%s: Finding factorable processor set stopped reduction: new_size=%d, neq(loc)=%" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, new_size, ncrs_eq));
    186:         PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
    187:         PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    188:         PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    189:         PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, a_Amat_crs));
    190:         PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    191:         PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    192:         PetscFunctionReturn(PETSC_SUCCESS);
    193:       }
    194:     }
    195:     /* make 'is_eq_newproc' */
    196:     if (pc_gamg->repart) { /* Repartition Cmat_{k} and move columns of P^{k}_{k-1} and coordinates of primal part accordingly */
    197:       Mat adj;

    199:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    200:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    201:       PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat));
    202:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    203:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    204:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0));
    205:       PetscCall(PetscInfo(pc, "%s: Repartition: size (active): %d --> %d, %" PetscInt_FMT " local equations, using %s process layout\n", ((PetscObject)pc)->prefix, *a_nactive_proc, new_size, ncrs_eq, (pc_gamg->layout_type == PCGAMG_LAYOUT_COMPACT) ? "compact" : "spread"));
    206:       /* get 'adj' */
    207:       if (cr_bs == 1) {
    208:         PetscCall(MatConvert(Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
    209:       } else {
    210:         /* make a scalar matrix to partition (no Stokes here) */
    211:         Mat                tMat;
    212:         PetscInt           Istart_crs, Iend_crs, ncols, jj, Ii;
    213:         const PetscScalar *vals;
    214:         const PetscInt    *idx;
    215:         PetscInt          *d_nnz, *o_nnz, M, N, maxnnz = 0, *j_buf = NULL;
    216:         PetscScalar       *v_buff = NULL;
    217:         static PetscInt    llev   = 0; /* ugly but just used for debugging */
    218:         MatType            mtype;

    220:         PetscCall(PetscMalloc2(ncrs, &d_nnz, ncrs, &o_nnz));
    221:         PetscCall(MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs));
    222:         PetscCall(MatGetSize(Cmat, &M, &N));
    223:         for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) {
    224:           PetscCall(MatGetRow(Cmat, Ii, &ncols, NULL, NULL));
    225:           d_nnz[jj] = ncols / cr_bs;
    226:           o_nnz[jj] = ncols / cr_bs;
    227:           if (ncols > maxnnz) maxnnz = ncols;
    228:           PetscCall(MatRestoreRow(Cmat, Ii, &ncols, NULL, NULL));
    229:           if (d_nnz[jj] > ncrs) d_nnz[jj] = ncrs;
    230:           if (o_nnz[jj] > (M / cr_bs - ncrs)) o_nnz[jj] = M / cr_bs - ncrs;
    231:         }

    233:         PetscCall(MatGetType(Amat_fine, &mtype));
    234:         PetscCall(MatCreate(comm, &tMat));
    235:         PetscCall(MatSetSizes(tMat, ncrs, ncrs, PETSC_DETERMINE, PETSC_DETERMINE));
    236:         PetscCall(MatSetType(tMat, mtype));
    237:         PetscCall(MatSeqAIJSetPreallocation(tMat, 0, d_nnz));
    238:         PetscCall(MatMPIAIJSetPreallocation(tMat, 0, d_nnz, 0, o_nnz));
    239:         PetscCall(PetscFree2(d_nnz, o_nnz));
    240:         PetscCall(PetscMalloc2(maxnnz, &j_buf, maxnnz, &v_buff));
    241:         for (ii = 0; ii < maxnnz; ii++) v_buff[ii] = 1.;

    243:         for (ii = Istart_crs; ii < Iend_crs; ii++) {
    244:           PetscInt dest_row = ii / cr_bs;
    245:           PetscCall(MatGetRow(Cmat, ii, &ncols, &idx, &vals));
    246:           for (jj = 0; jj < ncols; jj++) j_buf[jj] = idx[jj] / cr_bs;
    247:           PetscCall(MatSetValues(tMat, 1, &dest_row, ncols, j_buf, v_buff, ADD_VALUES));
    248:           PetscCall(MatRestoreRow(Cmat, ii, &ncols, &idx, &vals));
    249:         }
    250:         PetscCall(MatAssemblyBegin(tMat, MAT_FINAL_ASSEMBLY));
    251:         PetscCall(MatAssemblyEnd(tMat, MAT_FINAL_ASSEMBLY));
    252:         PetscCall(PetscFree2(j_buf, v_buff));

    254:         if (llev++ == -1) {
    255:           PetscViewer viewer;
    256:           char        fname[32];
    257:           PetscCall(PetscSNPrintf(fname, sizeof(fname), "part_mat_%" PetscInt_FMT ".mat", llev));
    258:           PetscCall(PetscViewerBinaryOpen(comm, fname, FILE_MODE_WRITE, &viewer));
    259:           PetscCall(MatView(tMat, viewer));
    260:           PetscCall(PetscViewerDestroy(&viewer));
    261:         }
    262:         PetscCall(MatConvert(tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
    263:         PetscCall(MatDestroy(&tMat));
    264:       } /* create 'adj' */

    266:       { /* partition: get newproc_idx */
    267:         char            prefix[256];
    268:         const char     *pcpre;
    269:         const PetscInt *is_idx;
    270:         MatPartitioning mpart;
    271:         IS              proc_is;

    273:         PetscCall(MatPartitioningCreate(comm, &mpart));
    274:         PetscCall(MatPartitioningSetAdjacency(mpart, adj));
    275:         PetscCall(PCGetOptionsPrefix(pc, &pcpre));
    276:         PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
    277:         PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
    278:         PetscCall(MatPartitioningSetFromOptions(mpart));
    279:         PetscCall(MatPartitioningSetNParts(mpart, new_size));
    280:         PetscCall(MatPartitioningApply(mpart, &proc_is));
    281:         PetscCall(MatPartitioningDestroy(&mpart));

    283:         /* collect IS info */
    284:         PetscCall(PetscMalloc1(ncrs_eq, &newproc_idx));
    285:         PetscCall(ISGetIndices(proc_is, &is_idx));
    286:         for (kk = jj = 0; kk < nloc_old; kk++) {
    287:           for (ii = 0; ii < cr_bs; ii++, jj++) { newproc_idx[jj] = is_idx[kk] * expand_factor; /* distribution */ }
    288:         }
    289:         PetscCall(ISRestoreIndices(proc_is, &is_idx));
    290:         PetscCall(ISDestroy(&proc_is));
    291:       }
    292:       PetscCall(MatDestroy(&adj));

    294:       PetscCall(ISCreateGeneral(comm, ncrs_eq, newproc_idx, PETSC_OWN_POINTER, &is_eq_newproc));
    295:       /*
    296:         Create an index set from the is_eq_newproc index set to indicate the mapping TO
    297:       */
    298:       PetscCall(ISPartitioningToNumbering(is_eq_newproc, &is_eq_num));
    299:       /*
    300:         Determine how many equations/vertices are assigned to each processor
    301:       */
    302:       PetscCall(PetscMalloc1(size, &counts));
    303:       PetscCall(ISPartitioningCount(is_eq_newproc, size, counts));
    304:       ncrs_eq_new = counts[rank];
    305:       PetscCall(ISDestroy(&is_eq_newproc));
    306:       PetscCall(PetscFree(counts));
    307:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0));
    308:     } else { /* simple aggregation of parts -- 'is_eq_newproc' */
    309:       const PetscInt *ranges;
    310:       PetscInt        newstart = 0;
    311:       PetscLayout     ilay;

    313:       PetscCheck(new_size != nactive, PETSC_COMM_SELF, PETSC_ERR_PLIB, "new_size==nactive. Should not happen");
    314:       PetscCall(PetscInfo(pc, "%s: Number of equations (loc) %" PetscInt_FMT " with simple aggregation\n", ((PetscObject)pc)->prefix, ncrs_eq));
    315:       PetscCallMPI(MPI_Exscan(&ncrs_eq, &newstart, 1, MPIU_INT, MPI_SUM, comm));
    316:       PetscCall(ISCreateStride(comm, ncrs_eq, newstart, 1, &is_eq_num));
    317:       PetscCall(ISSetPermutation(is_eq_num));
    318:       PetscCall(ISGetLayout(is_eq_num, &ilay));
    319:       PetscCall(PetscLayoutGetRanges(ilay, &ranges));
    320:       ncrs_eq_new = 0;
    321:       for (PetscInt r = 0; r < size; r++)
    322:         if (rank == (r / rfactor) * expand_factor) ncrs_eq_new += ranges[r + 1] - ranges[r];
    323:       //targetPE = (rank / rfactor) * expand_factor;
    324:       //PetscCall(ISCreateStride(comm, ncrs_eq, targetPE, 0, &is_eq_newproc));
    325:       //PetscCall(ISPartitioningToNumbering(is_eq_newproc, &is_eq_num));
    326:       //PetscCall(PetscMalloc1(size, &counts));
    327:       //PetscCall(ISPartitioningCount(is_eq_newproc, size, counts));
    328:       //ncrs_eq_new = counts[rank];
    329:       //PetscCall(ISDestroy(&is_eq_newproc));
    330:       //PetscCall(PetscFree(counts));
    331:     } /* end simple 'is_eq_newproc' */

    333:     ncrs_new = ncrs_eq_new / cr_bs;

    335:     /* data movement scope -- this could be moved to subclasses so that we don't try to cram all auxiliary data into some complex abstracted thing */
    336:     {
    337:       Vec             src_crd, dest_crd;
    338:       const PetscInt *idx, ndata_rows = pc_gamg->data_cell_rows, ndata_cols = pc_gamg->data_cell_cols, node_data_sz = ndata_rows * ndata_cols;
    339:       VecScatter      vecscat;
    340:       PetscScalar    *array;
    341:       IS              isscat;
    342:       /* move data (for primal equations only) */
    343:       /* Create a vector to contain the newly ordered element information */
    344:       PetscCall(VecCreate(comm, &dest_crd));
    345:       PetscCall(VecSetSizes(dest_crd, node_data_sz * ncrs_new, PETSC_DECIDE));
    346:       PetscCall(VecSetType(dest_crd, VECSTANDARD)); /* this is needed! */
    347:       /*
    348:         There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having
    349:         a block size of ...).  Note, ISs are expanded into equation space by 'cr_bs'.
    350:       */
    351:       PetscCall(PetscMalloc1(ncrs * node_data_sz, &tidx));
    352:       PetscCall(ISGetIndices(is_eq_num, &idx));
    353:       for (ii = 0, jj = 0; ii < ncrs; ii++) {
    354:         PetscInt id = idx[ii * cr_bs] / cr_bs; /* get node back */
    355:         for (kk = 0; kk < node_data_sz; kk++, jj++) tidx[jj] = id * node_data_sz + kk;
    356:       }
    357:       PetscCall(ISRestoreIndices(is_eq_num, &idx));
    358:       PetscCall(ISCreateGeneral(comm, node_data_sz * ncrs, tidx, PETSC_COPY_VALUES, &isscat));
    359:       PetscCall(PetscFree(tidx));
    360:       /*
    361:         Create a vector to contain the original vertex information for each element
    362:       */
    363:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, node_data_sz * ncrs, &src_crd));
    364:       for (jj = 0; jj < ndata_cols; jj++) {
    365:         const PetscInt stride0 = ncrs * pc_gamg->data_cell_rows;
    366:         for (ii = 0; ii < ncrs; ii++) {
    367:           for (kk = 0; kk < ndata_rows; kk++) {
    368:             PetscInt    ix = ii * ndata_rows + kk + jj * stride0, jx = ii * node_data_sz + kk * ndata_cols + jj;
    369:             PetscScalar tt = (PetscScalar)pc_gamg->data[ix];
    370:             PetscCall(VecSetValues(src_crd, 1, &jx, &tt, INSERT_VALUES));
    371:           }
    372:         }
    373:       }
    374:       PetscCall(VecAssemblyBegin(src_crd));
    375:       PetscCall(VecAssemblyEnd(src_crd));
    376:       /*
    377:         Scatter the element vertex information (still in the original vertex ordering)
    378:         to the correct processor
    379:       */
    380:       PetscCall(VecScatterCreate(src_crd, NULL, dest_crd, isscat, &vecscat));
    381:       PetscCall(ISDestroy(&isscat));
    382:       PetscCall(VecScatterBegin(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD));
    383:       PetscCall(VecScatterEnd(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD));
    384:       PetscCall(VecScatterDestroy(&vecscat));
    385:       PetscCall(VecDestroy(&src_crd));
    386:       /*
    387:         Put the element vertex data into a new allocation of the gdata->ele
    388:       */
    389:       PetscCall(PetscFree(pc_gamg->data));
    390:       PetscCall(PetscMalloc1(node_data_sz * ncrs_new, &pc_gamg->data));

    392:       pc_gamg->data_sz = node_data_sz * ncrs_new;
    393:       strideNew        = ncrs_new * ndata_rows;

    395:       PetscCall(VecGetArray(dest_crd, &array));
    396:       for (jj = 0; jj < ndata_cols; jj++) {
    397:         for (ii = 0; ii < ncrs_new; ii++) {
    398:           for (kk = 0; kk < ndata_rows; kk++) {
    399:             PetscInt ix = ii * ndata_rows + kk + jj * strideNew, jx = ii * node_data_sz + kk * ndata_cols + jj;
    400:             pc_gamg->data[ix] = PetscRealPart(array[jx]);
    401:           }
    402:         }
    403:       }
    404:       PetscCall(VecRestoreArray(dest_crd, &array));
    405:       PetscCall(VecDestroy(&dest_crd));
    406:     }
    407:     /* move A and P (columns) with new layout */
    408:     /*
    409:       Invert for MatCreateSubMatrix
    410:     */
    411:     PetscCall(ISInvertPermutation(is_eq_num, ncrs_eq_new, &new_eq_indices));
    412:     PetscCall(ISSort(new_eq_indices));
    413:     PetscCall(ISSetBlockSize(new_eq_indices, cr_bs));
    414:     if (Pcolumnperm) {
    415:       PetscCall(PetscObjectReference((PetscObject)new_eq_indices));
    416:       *Pcolumnperm = new_eq_indices;
    417:     }
    418:     PetscCall(ISDestroy(&is_eq_num));

    420:     /* 'a_Amat_crs' output */
    421:     if (Cmat) { /* repartitioning from Cmat adjacency case */
    422:       Mat       mat;
    423:       PetscBool isset, isspd, isher;
    424: #if !defined(PETSC_USE_COMPLEX)
    425:       PetscBool issym;
    426: #endif

    428:       PetscCall(MatCreateSubMatrix(Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat));
    429:       PetscCall(MatIsSPDKnown(Cmat, &isset, &isspd)); // like MatPropagateSymmetryOptions, but should set MAT_STRUCTURALLY_SYMMETRIC ?
    430:       if (isset) PetscCall(MatSetOption(mat, MAT_SPD, isspd));
    431:       else {
    432:         PetscCall(MatIsHermitianKnown(Cmat, &isset, &isher));
    433:         if (isset) PetscCall(MatSetOption(mat, MAT_HERMITIAN, isher));
    434:         else {
    435: #if !defined(PETSC_USE_COMPLEX)
    436:           PetscCall(MatIsSymmetricKnown(Cmat, &isset, &issym));
    437:           if (isset) PetscCall(MatSetOption(mat, MAT_SYMMETRIC, issym));
    438: #endif
    439:         }
    440:       }
    441:       *a_Amat_crs = mat;
    442:     }

    444:     /* prolongator */
    445:     {
    446:       IS       findices;
    447:       PetscInt Istart, Iend;
    448:       Mat      Pnew;

    450:       PetscCall(MatGetOwnershipRange(Pold, &Istart, &Iend));
    451:       PetscCall(ISCreateStride(comm, Iend - Istart, Istart, 1, &findices));
    452:       PetscCall(ISSetBlockSize(findices, f_bs));
    453:       PetscCall(MatCreateSubMatrix(Pold, findices, new_eq_indices, MAT_INITIAL_MATRIX, &Pnew));
    454:       PetscCall(ISDestroy(&findices));
    455:       PetscCall(MatSetOption(Pnew, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));

    457:       PetscCall(MatDestroy(a_P_inout));

    459:       /* output - repartitioned */
    460:       *a_P_inout = Pnew;
    461:     }

    463:     if (!Cmat) { /* simple repartitioning case */
    464:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    465:       PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    466:       PetscCall(MatPtAP(Amat_fine, *a_P_inout, MAT_INITIAL_MATRIX, 2.0, a_Amat_crs));
    467:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
    468:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
    469:     }
    470:     PetscCall(MatDestroy(&Cmat));
    471:     PetscCall(ISDestroy(&new_eq_indices));

    473:     *a_nactive_proc = new_size; /* output */

    475:     /* pinning on reduced grids, not a bad heuristic and optimization gets folded into process reduction optimization */
    476:     if (pc_gamg->cpu_pin_coarse_grids) {
    477: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
    478:       static PetscInt llev = 2;
    479:       PetscCall(PetscInfo(pc, "%s: Pinning level %" PetscInt_FMT " to the CPU\n", ((PetscObject)pc)->prefix, llev++));
    480: #endif
    481:       PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE));
    482:       PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE));
    483:       if (1) { /* HACK: move this to MatBindCPU_MPIAIJXXX; lvec is created, need to pin it, this is done in MatSetUpMultiply_MPIAIJ. Hack */
    484:         Mat         A = *a_Amat_crs, P = *a_P_inout;
    485:         PetscMPIInt size;
    486:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    487:         if (size > 1) {
    488:           Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
    489:           PetscCall(VecBindToCPU(a->lvec, PETSC_TRUE));
    490:           PetscCall(VecBindToCPU(p->lvec, PETSC_TRUE));
    491:         }
    492:       }
    493:     }
    494:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
    495:   } // processor reduce
    496:   PetscFunctionReturn(PETSC_SUCCESS);
    497: }

    499: // used in GEO
    500: PetscErrorCode PCGAMGSquareGraph_GAMG(PC a_pc, Mat Gmat1, Mat *Gmat2)
    501: {
    502:   const char *prefix;
    503:   char        addp[32];
    504:   PC_MG      *mg      = (PC_MG *)a_pc->data;
    505:   PC_GAMG    *pc_gamg = (PC_GAMG *)mg->innerctx;

    507:   PetscFunctionBegin;
    508:   PetscCall(PCGetOptionsPrefix(a_pc, &prefix));
    509:   PetscCall(PetscInfo(a_pc, "%s: Square Graph on level %" PetscInt_FMT "\n", ((PetscObject)a_pc)->prefix, pc_gamg->current_level + 1));
    510:   PetscCall(MatProductCreate(Gmat1, Gmat1, NULL, Gmat2));
    511:   PetscCall(MatSetOptionsPrefix(*Gmat2, prefix));
    512:   PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_square_%" PetscInt_FMT "_", pc_gamg->current_level));
    513:   PetscCall(MatAppendOptionsPrefix(*Gmat2, addp));
    514:   if ((*Gmat2)->structurally_symmetric == PETSC_BOOL3_TRUE) {
    515:     PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AB));
    516:   } else {
    517:     PetscCall(MatSetOption(Gmat1, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
    518:     PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AtB));
    519:   }
    520:   PetscCall(MatProductSetFromOptions(*Gmat2));
    521:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0));
    522:   PetscCall(MatProductSymbolic(*Gmat2));
    523:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0));
    524:   PetscCall(MatProductClear(*Gmat2));
    525:   /* we only need the sparsity, cheat and tell PETSc the matrix has been assembled */
    526:   (*Gmat2)->assembled = PETSC_TRUE;
    527:   PetscFunctionReturn(PETSC_SUCCESS);
    528: }

    530: /*
    531:    PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner
    532:                     by setting data structures and options.

    534:    Input Parameter:
    535: .  pc - the preconditioner context

    537: */
    538: static PetscErrorCode PCSetUp_GAMG(PC pc)
    539: {
    540:   PC_MG      *mg      = (PC_MG *)pc->data;
    541:   PC_GAMG    *pc_gamg = (PC_GAMG *)mg->innerctx;
    542:   Mat         Pmat    = pc->pmat;
    543:   PetscInt    fine_level, level, level1, bs, M, N, qq, lidx, nASMBlocksArr[PETSC_MG_MAXLEVELS], cr_bs;
    544:   MPI_Comm    comm;
    545:   PetscMPIInt rank, size, nactivepe;
    546:   Mat         Aarr[PETSC_MG_MAXLEVELS], Parr[PETSC_MG_MAXLEVELS];
    547:   IS         *ASMLocalIDsArr[PETSC_MG_MAXLEVELS];
    548:   PetscBool   is_last = PETSC_FALSE;
    549: #if defined(PETSC_USE_INFO)
    550:   PetscLogDouble nnz0 = 0., nnztot = 0.;
    551:   MatInfo        info;
    552: #endif

    554:   PetscFunctionBegin;
    555:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
    556:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    557:   PetscCallMPI(MPI_Comm_size(comm, &size));
    558:   PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
    559:   if (pc->setupcalled) {
    560:     if (!pc_gamg->reuse_prol || pc->flag == DIFFERENT_NONZERO_PATTERN) {
    561:       /* reset everything */
    562:       PetscCall(PCReset_MG(pc));
    563:       pc->setupcalled = 0;
    564:     } else {
    565:       PC_MG_Levels **mglevels = mg->levels;
    566:       /* just do Galerkin grids */
    567:       Mat B, dA, dB;
    568:       if (pc_gamg->Nlevels > 1) {
    569:         PetscInt gl;
    570:         /* currently only handle case where mat and pmat are the same on coarser levels */
    571:         PetscCall(KSPGetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, &dA, &dB));
    572:         /* (re)set to get dirty flag */
    573:         PetscCall(KSPSetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, dA, dB));

    575:         for (level = pc_gamg->Nlevels - 2, gl = 0; level >= 0; level--, gl++) {
    576:           MatReuse reuse = MAT_INITIAL_MATRIX;
    577: #if defined(GAMG_STAGES)
    578:           PetscCall(PetscLogStagePush(gamg_stages[gl]));
    579: #endif
    580:           /* matrix nonzero structure can change from repartitioning or process reduction but don't know if we have process reduction here. Should fix */
    581:           PetscCall(KSPGetOperators(mglevels[level]->smoothd, NULL, &B));
    582:           if (B->product) {
    583:             if (B->product->A == dB && B->product->B == mglevels[level + 1]->interpolate) reuse = MAT_REUSE_MATRIX;
    584:           }
    585:           if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDestroy(&mglevels[level]->A));
    586:           if (reuse == MAT_REUSE_MATRIX) {
    587:             PetscCall(PetscInfo(pc, "%s: RAP after initial setup, reuse matrix level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
    588:           } else {
    589:             PetscCall(PetscInfo(pc, "%s: RAP after initial setup, with repartitioning (new matrix) level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
    590:           }
    591:           PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0));
    592:           PetscCall(MatPtAP(dB, mglevels[level + 1]->interpolate, reuse, PETSC_DETERMINE, &B));
    593:           PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0));
    594:           if (reuse == MAT_INITIAL_MATRIX) mglevels[level]->A = B;
    595:           PetscCall(KSPSetOperators(mglevels[level]->smoothd, B, B));
    596:           // check for redoing eigen estimates
    597:           if (pc_gamg->recompute_esteig) {
    598:             PetscBool ischeb;
    599:             KSP       smoother;
    600:             PetscCall(PCMGGetSmoother(pc, level + 1, &smoother));
    601:             PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb));
    602:             if (ischeb) {
    603:               KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data;
    604:               cheb->emin_provided = 0;
    605:               cheb->emax_provided = 0;
    606:             }
    607:             /* we could call PetscCall(KSPChebyshevSetEigenvalues(smoother, 0, 0)); but the logic does not work properly */
    608:           }
    609:           // inc
    610:           dB = B;
    611: #if defined(GAMG_STAGES)
    612:           PetscCall(PetscLogStagePop());
    613: #endif
    614:         }
    615:       }

    617:       PetscCall(PCSetUp_MG(pc));
    618:       PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
    619:       PetscFunctionReturn(PETSC_SUCCESS);
    620:     }
    621:   }

    623:   if (!pc_gamg->data) {
    624:     if (pc_gamg->orig_data) {
    625:       PetscCall(MatGetBlockSize(Pmat, &bs));
    626:       PetscCall(MatGetLocalSize(Pmat, &qq, NULL));

    628:       pc_gamg->data_sz        = (qq / bs) * pc_gamg->orig_data_cell_rows * pc_gamg->orig_data_cell_cols;
    629:       pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows;
    630:       pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols;

    632:       PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data));
    633:       for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq];
    634:     } else {
    635:       PetscCheck(pc_gamg->ops->createdefaultdata, comm, PETSC_ERR_PLIB, "'createdefaultdata' not set(?) need to support NULL data");
    636:       PetscCall(pc_gamg->ops->createdefaultdata(pc, Pmat));
    637:     }
    638:   }

    640:   /* cache original data for reuse */
    641:   if (!pc_gamg->orig_data && (PetscBool)(!pc_gamg->reuse_prol)) {
    642:     PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->orig_data));
    643:     for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq];
    644:     pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows;
    645:     pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols;
    646:   }

    648:   /* get basic dims */
    649:   PetscCall(MatGetBlockSize(Pmat, &bs));
    650:   PetscCall(MatGetSize(Pmat, &M, NULL));

    652: #if defined(PETSC_USE_INFO)
    653:   PetscCall(MatGetInfo(Pmat, MAT_GLOBAL_SUM, &info)); /* global reduction */
    654:   nnz0   = info.nz_used;
    655:   nnztot = info.nz_used;
    656: #endif
    657:   PetscCall(PetscInfo(pc, "%s: level %d) N=%" PetscInt_FMT ", n data rows=%" PetscInt_FMT ", n data cols=%" PetscInt_FMT ", nnz/row (ave)=%" PetscInt_FMT ", block size %d, np=%d\n", ((PetscObject)pc)->prefix, 0, M, pc_gamg->data_cell_rows, pc_gamg->data_cell_cols, (PetscInt)(nnz0 / (PetscReal)M + 0.5), (int)bs, size));

    659:   /* Get A_i and R_i */
    660:   for (level = 0, Aarr[0] = Pmat, nactivepe = size; level < (pc_gamg->Nlevels - 1) && (level == 0 || M > pc_gamg->coarse_eq_limit); level++) {
    661:     pc_gamg->current_level = level;
    662:     level1                 = level + 1;
    663: #if defined(GAMG_STAGES)
    664:     if (!gamg_stages[level]) {
    665:       char str[32];
    666:       PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", (int)level));
    667:       PetscCall(PetscLogStageRegister(str, &gamg_stages[level]));
    668:     }
    669:     PetscCall(PetscLogStagePush(gamg_stages[level]));
    670: #endif
    671:     /* construct prolongator - Parr[level1] */
    672:     if (level == 0 && pc_gamg->injection_index_size > 0) {
    673:       Mat      Prol;
    674:       MatType  mtype;
    675:       PetscInt prol_m, prol_n, Prol_N = (M / bs) * pc_gamg->injection_index_size, Istart, Iend, nn, row;
    676:       PetscCall(PetscInfo(pc, "Create fine grid injection space prolongation %" PetscInt_FMT " x %" PetscInt_FMT ". %s\n", M, Prol_N, pc_gamg->data ? "delete null space data" : ""));
    677:       PetscCall(MatGetOwnershipRange(Pmat, &Istart, &Iend));
    678:       PetscCall(MatGetLocalSize(Pmat, &prol_m, NULL)); // rows m x n
    679:       prol_n = (prol_m / bs) * pc_gamg->injection_index_size;
    680:       PetscCheck(pc_gamg->injection_index_size < bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Injection size %d must be less that block size %d", (int)pc_gamg->injection_index_size, (int)bs);
    681:       PetscCall(MatGetType(Pmat, &mtype));
    682:       PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &Prol));
    683:       PetscCall(MatSetBlockSizes(Prol, bs, pc_gamg->injection_index_size));
    684:       PetscCall(MatSetSizes(Prol, prol_m, prol_n, M, Prol_N));
    685:       PetscCall(MatSetType(Prol, mtype));
    686: #if PetscDefined(HAVE_DEVICE)
    687:       PetscBool flg;
    688:       PetscCall(MatBoundToCPU(Pmat, &flg));
    689:       PetscCall(MatBindToCPU(Prol, flg));
    690:       if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
    691: #endif
    692:       PetscCall(MatSeqAIJSetPreallocation(Prol, 1, NULL));
    693:       PetscCall(MatMPIAIJSetPreallocation(Prol, 1, NULL, 0, NULL));
    694:       // set I \kron [1, 1, ... ]^T
    695:       for (PetscInt ii = Istart, col = (Istart / bs) * pc_gamg->injection_index_size; ii < Iend; ii += bs) {
    696:         const PetscScalar one = 1;
    697:         for (PetscInt jj = 0; jj < pc_gamg->injection_index_size; jj++, col++) {
    698:           PetscInt row = ii + pc_gamg->injection_index[jj];
    699:           PetscCall(MatSetValues(Prol, 1, &row, 1, &col, &one, INSERT_VALUES));
    700:         }
    701:       }
    702:       PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
    703:       PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
    704:       PetscCall(MatViewFromOptions(Prol, NULL, "-mat_view_injection"));
    705:       PetscCall(MatGetBlockSizes(Prol, NULL, &cr_bs)); // column size
    706:       Parr[level1] = Prol;
    707:       // can not deal with null space -- with array of 'injection cols' we could take 'injection rows and 'injection cols' to 'data'
    708:       if (pc_gamg->data) {
    709:         pc_gamg->data_cell_cols      = pc_gamg->injection_index_size;
    710:         pc_gamg->data_cell_rows      = pc_gamg->injection_index_size;
    711:         pc_gamg->orig_data_cell_cols = 0;
    712:         pc_gamg->orig_data_cell_rows = 0;
    713:         PetscCall(PetscFree(pc_gamg->data));
    714:         pc_gamg->data_sz = pc_gamg->injection_index_size * prol_n;
    715:         PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data));
    716:         for (row = nn = 0; row < prol_n; row += pc_gamg->injection_index_size) {
    717:           for (PetscInt jj = 0; jj < pc_gamg->injection_index_size; jj++) {
    718:             PetscInt idx = row * pc_gamg->injection_index_size + jj * pc_gamg->injection_index_size;
    719:             for (PetscInt kk = 0; kk < pc_gamg->injection_index_size; kk++, nn++) { pc_gamg->data[idx + kk] = (jj == kk) ? 1 : 0; }
    720:           }
    721:         }
    722:         PetscCheck(nn == pc_gamg->data_sz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nn != pc_gamg->data_sz %" PetscInt_FMT " %" PetscInt_FMT, pc_gamg->data_sz, nn);
    723:       }
    724:     } else {
    725:       Mat               Gmat, mat;
    726:       PetscCoarsenData *agg_lists;
    727:       Mat               Prol11;

    729:       PetscCall(PCGAMGCreateGraph(pc, Aarr[level], &Gmat));
    730:       PetscCall(pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists)); // Gmat may have ghosts for QR aggregates not in matrix
    731:       PetscCall(PetscCDGetMat(agg_lists, &mat));
    732:       if (!mat) PetscCall(PetscCDSetMat(agg_lists, Gmat));
    733:       PetscCall(pc_gamg->ops->prolongator(pc, Aarr[level], agg_lists, &Prol11));
    734:       /* could have failed to create new level */
    735:       if (Prol11) {
    736:         const char *prefix;
    737:         char        addp[32];

    739:         /* get new block size of coarse matrices */
    740:         PetscCall(MatGetBlockSizes(Prol11, NULL, &cr_bs)); // column size

    742:         if (pc_gamg->ops->optprolongator) {
    743:           /* smooth */
    744:           PetscCall(pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11));
    745:         }

    747:         if (pc_gamg->use_aggs_in_asm) {
    748:           PetscInt bs;
    749:           PetscCall(MatGetBlockSizes(Prol11, &bs, NULL)); // row block size
    750:           PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level]));
    751:           PetscCall(PetscInfo(pc, "%d: %" PetscInt_FMT " ASM local domains,  bs = %d\n", (int)level, nASMBlocksArr[level], (int)bs));
    752:         } else if (pc_gamg->asm_hem_aggs) {
    753:           const char *prefix;
    754:           PetscInt    bs;

    756:           /*
    757:              Do not use aggs created for defining coarser problems, instead create aggs specifically to use
    758:              to define PCASM blocks.
    759:           */
    760:           PetscCall(PetscCDGetMat(agg_lists, &mat));
    761:           if (mat == Gmat) PetscCall(PetscCDClearMat(agg_lists)); // take the Mat away from the list (yuck)
    762:           PetscCall(PetscCDDestroy(agg_lists));
    763:           PetscCall(PetscInfo(pc, "HEM ASM passes = %d\n", (int)pc_gamg->asm_hem_aggs));
    764:           PetscCall(MatCoarsenDestroy(&pc_gamg->asm_crs));
    765:           PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg->asm_crs));
    766:           PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
    767:           PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg->asm_crs, prefix));
    768:           PetscCall(MatCoarsenSetFromOptions(pc_gamg->asm_crs)); // get strength args
    769:           PetscCall(MatCoarsenSetType(pc_gamg->asm_crs, MATCOARSENHEM));
    770:           PetscCall(MatCoarsenSetMaximumIterations(pc_gamg->asm_crs, pc_gamg->asm_hem_aggs));
    771:           PetscCall(MatCoarsenSetAdjacency(pc_gamg->asm_crs, Gmat));
    772:           PetscCall(MatCoarsenSetStrictAggs(pc_gamg->asm_crs, PETSC_TRUE));
    773:           PetscCall(MatCoarsenApply(pc_gamg->asm_crs));
    774:           PetscCall(MatCoarsenGetData(pc_gamg->asm_crs, &agg_lists)); /* output */
    775:           // create aggregates
    776:           PetscCall(MatGetBlockSizes(Aarr[level], &bs, NULL)); // row block size
    777:           PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level]));
    778:         }
    779:         PetscCall(PCGetOptionsPrefix(pc, &prefix));
    780:         PetscCall(MatSetOptionsPrefix(Prol11, prefix));
    781:         PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_prolongator_%d_", (int)level));
    782:         PetscCall(MatAppendOptionsPrefix(Prol11, addp));
    783:         /* Always generate the transpose with CUDA
    784:            Such behaviour can be adapted with -pc_gamg_prolongator_ prefixed options */
    785:         PetscCall(MatSetOption(Prol11, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
    786:         PetscCall(MatSetFromOptions(Prol11));
    787:         Parr[level1] = Prol11;
    788:       } else Parr[level1] = NULL; /* failed to coarsen */
    789:       PetscCall(PetscCDGetMat(agg_lists, &mat));
    790:       if (mat == Gmat) PetscCall(PetscCDClearMat(agg_lists)); // take the Mat away from the list (yuck)
    791:       PetscCall(MatDestroy(&Gmat));
    792:       PetscCall(PetscCDDestroy(agg_lists));
    793:     } /* construct prolongator scope */
    794:     if (level == 0) Aarr[0] = Pmat; /* use Pmat for finest level setup */
    795:     if (!Parr[level1]) {            /* failed to coarsen */
    796:       PetscCall(PetscInfo(pc, "%s: Stop gridding, level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
    797: #if defined(GAMG_STAGES)
    798:       PetscCall(PetscLogStagePop());
    799: #endif
    800:       break;
    801:     }
    802:     PetscCall(MatGetSize(Parr[level1], &M, &N)); /* N is next M, a loop test variables */
    803:     PetscCheck(!is_last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Is last ?");
    804:     if (N <= pc_gamg->coarse_eq_limit) is_last = PETSC_TRUE;
    805:     if (level1 == pc_gamg->Nlevels - 1) is_last = PETSC_TRUE;
    806:     if (level == PETSC_MG_MAXLEVELS - 2) is_last = PETSC_TRUE;
    807:     PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0));
    808:     PetscCall(pc_gamg->ops->createlevel(pc, Aarr[level], cr_bs, &Parr[level1], &Aarr[level1], &nactivepe, NULL, is_last));
    809:     PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0));

    811:     PetscCall(MatGetSize(Aarr[level1], &M, &N)); /* M is loop test variables */
    812: #if defined(PETSC_USE_INFO)
    813:     PetscCall(MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info));
    814:     nnztot += info.nz_used;
    815: #endif
    816:     PetscCall(PetscInfo(pc, "%s: %d) N=%" PetscInt_FMT ", n data cols=%" PetscInt_FMT ", nnz/row (ave)=%" PetscInt_FMT ", %d active pes\n", ((PetscObject)pc)->prefix, (int)level1, M, pc_gamg->data_cell_cols, (PetscInt)(info.nz_used / (PetscReal)M), nactivepe));

    818: #if defined(GAMG_STAGES)
    819:     PetscCall(PetscLogStagePop());
    820: #endif
    821:     /* stop if one node or one proc -- could pull back for singular problems */
    822:     if ((pc_gamg->data_cell_cols && M / pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M / bs < 2)) {
    823:       PetscCall(PetscInfo(pc, "%s: HARD stop of coarsening on level %" PetscInt_FMT ".  Grid too small: %" PetscInt_FMT " block nodes\n", ((PetscObject)pc)->prefix, level, M / bs));
    824:       level++;
    825:       break;
    826:     } else if (level == PETSC_MG_MAXLEVELS - 2) { /* stop if we are limited by PC_MG_MAXLEVELS */
    827:       PetscCall(PetscInfo(pc, "%s: HARD stop of coarsening on level %" PetscInt_FMT ".  PC_MG_MAXLEVELS reached\n", ((PetscObject)pc)->prefix, level));
    828:       level++;
    829:       break;
    830:     }
    831:   } /* levels */
    832:   PetscCall(PetscFree(pc_gamg->data));

    834:   PetscCall(PetscInfo(pc, "%s: %" PetscInt_FMT " levels, operator complexity = %g\n", ((PetscObject)pc)->prefix, level + 1, nnztot / nnz0));
    835:   pc_gamg->Nlevels = level + 1;
    836:   fine_level       = level;
    837:   PetscCall(PCMGSetLevels(pc, pc_gamg->Nlevels, NULL));

    839:   if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */

    841:     /* set default smoothers & set operators */
    842:     for (lidx = 1, level = pc_gamg->Nlevels - 2; lidx <= fine_level; lidx++, level--) {
    843:       KSP smoother;
    844:       PC  subpc;

    846:       PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
    847:       PetscCall(KSPGetPC(smoother, &subpc));

    849:       PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
    850:       /* set ops */
    851:       PetscCall(KSPSetOperators(smoother, Aarr[level], Aarr[level]));
    852:       PetscCall(PCMGSetInterpolation(pc, lidx, Parr[level + 1]));

    854:       /* set defaults */
    855:       PetscCall(KSPSetType(smoother, KSPCHEBYSHEV));

    857:       /* set blocks for ASM smoother that uses the 'aggregates' */
    858:       if (pc_gamg->use_aggs_in_asm || pc_gamg->asm_hem_aggs) {
    859:         PetscInt sz;
    860:         IS      *iss;

    862:         sz  = nASMBlocksArr[level];
    863:         iss = ASMLocalIDsArr[level];
    864:         PetscCall(PCSetType(subpc, PCASM));
    865:         PetscCall(PCASMSetOverlap(subpc, 0));
    866:         PetscCall(PCASMSetType(subpc, PC_ASM_BASIC));
    867:         if (!sz) {
    868:           IS is;
    869:           PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is));
    870:           PetscCall(PCASMSetLocalSubdomains(subpc, 1, NULL, &is));
    871:           PetscCall(ISDestroy(&is));
    872:         } else {
    873:           PetscInt kk;
    874:           PetscCall(PCASMSetLocalSubdomains(subpc, sz, iss, NULL));
    875:           for (kk = 0; kk < sz; kk++) PetscCall(ISDestroy(&iss[kk]));
    876:           PetscCall(PetscFree(iss));
    877:         }
    878:         ASMLocalIDsArr[level] = NULL;
    879:         nASMBlocksArr[level]  = 0;
    880:       } else {
    881:         PetscCall(PCSetType(subpc, PCJACOBI));
    882:       }
    883:     }
    884:     {
    885:       /* coarse grid */
    886:       KSP      smoother, *k2;
    887:       PC       subpc, pc2;
    888:       PetscInt ii, first;
    889:       Mat      Lmat = Aarr[pc_gamg->Nlevels - 1];
    890:       lidx          = 0;

    892:       PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
    893:       PetscCall(KSPSetOperators(smoother, Lmat, Lmat));
    894:       if (!pc_gamg->use_parallel_coarse_grid_solver) {
    895:         PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
    896:         PetscCall(KSPGetPC(smoother, &subpc));
    897:         PetscCall(PCSetType(subpc, PCBJACOBI));
    898:         PetscCall(PCSetUp(subpc));
    899:         PetscCall(PCBJacobiGetSubKSP(subpc, &ii, &first, &k2));
    900:         PetscCheck(ii == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " is not one", ii);
    901:         PetscCall(KSPGetPC(k2[0], &pc2));
    902:         PetscCall(PCSetType(pc2, PCLU));
    903:         PetscCall(PCFactorSetShiftType(pc2, MAT_SHIFT_INBLOCKS));
    904:         PetscCall(KSPSetTolerances(k2[0], PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1));
    905:         PetscCall(KSPSetType(k2[0], KSPPREONLY));
    906:       }
    907:     }

    909:     /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
    910:     PetscObjectOptionsBegin((PetscObject)pc);
    911:     PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject));
    912:     PetscOptionsEnd();
    913:     PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));

    915:     /* set cheby eigen estimates from SA to use in the solver */
    916:     if (pc_gamg->use_sa_esteig) {
    917:       for (lidx = 1, level = pc_gamg->Nlevels - 2; level >= 0; lidx++, level--) {
    918:         KSP       smoother;
    919:         PetscBool ischeb;

    921:         PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
    922:         PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb));
    923:         if (ischeb) {
    924:           KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data;

    926:           // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG
    927:           if (mg->max_eigen_DinvA[level] > 0) {
    928:             // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it.
    929:             // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.)
    930:             PetscReal emax, emin;

    932:             emin = mg->min_eigen_DinvA[level];
    933:             emax = mg->max_eigen_DinvA[level];
    934:             PetscCall(PetscInfo(pc, "%s: PCSetUp_GAMG: call KSPChebyshevSetEigenvalues on level %" PetscInt_FMT " (N=%" PetscInt_FMT ") with emax = %g emin = %g\n", ((PetscObject)pc)->prefix, level, Aarr[level]->rmap->N, (double)emax, (double)emin));
    935:             cheb->emin_provided = emin;
    936:             cheb->emax_provided = emax;
    937:           }
    938:         }
    939:       }
    940:     }

    942:     PetscCall(PCSetUp_MG(pc));

    944:     /* clean up */
    945:     for (level = 1; level < pc_gamg->Nlevels; level++) {
    946:       PetscCall(MatDestroy(&Parr[level]));
    947:       PetscCall(MatDestroy(&Aarr[level]));
    948:     }
    949:   } else {
    950:     KSP smoother;

    952:     PetscCall(PetscInfo(pc, "%s: One level solver used (system is seen as DD). Using default solver.\n", ((PetscObject)pc)->prefix));
    953:     PetscCall(PCMGGetSmoother(pc, 0, &smoother));
    954:     PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0]));
    955:     PetscCall(KSPSetType(smoother, KSPPREONLY));
    956:     PetscCall(PCSetUp_MG(pc));
    957:   }
    958:   PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
    959:   PetscFunctionReturn(PETSC_SUCCESS);
    960: }

    962: /*
    963:  PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
    964:    that was created with PCCreate_GAMG().

    966:    Input Parameter:
    967: .  pc - the preconditioner context

    969:    Application Interface Routine: PCDestroy()
    970: */
    971: PetscErrorCode PCDestroy_GAMG(PC pc)
    972: {
    973:   PC_MG   *mg      = (PC_MG *)pc->data;
    974:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    976:   PetscFunctionBegin;
    977:   PetscCall(PCReset_GAMG(pc));
    978:   if (pc_gamg->ops->destroy) PetscCall((*pc_gamg->ops->destroy)(pc));
    979:   PetscCall(PetscFree(pc_gamg->ops));
    980:   PetscCall(PetscFree(pc_gamg->gamg_type_name));
    981:   PetscCall(PetscFree(pc_gamg));
    982:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", NULL));
    983:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", NULL));
    984:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", NULL));
    985:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", NULL));
    986:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", NULL));
    987:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", NULL));
    988:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", NULL));
    989:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", NULL));
    990:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", NULL));
    991:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", NULL));
    992:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", NULL));
    993:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", NULL));
    994:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", NULL));
    995:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", NULL));
    996:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", NULL));
    997:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", NULL));
    998:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", NULL));
    999:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", NULL));
    1000:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndices_C", NULL));
    1001:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", NULL));
    1002:   PetscCall(PCDestroy_MG(pc));
    1003:   PetscFunctionReturn(PETSC_SUCCESS);
    1004: }

    1006: /*@
    1007:   PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction in `PCGAMG`

    1009:   Logically Collective

    1011:   Input Parameters:
    1012: + pc - the preconditioner context
    1013: - n  - the number of equations

    1015:   Options Database Key:
    1016: . -pc_gamg_process_eq_limit <limit> - set the limit

    1018:   Level: intermediate

    1020:   Note:
    1021:   `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process
    1022:   that has degrees of freedom

    1024:   Developer Note:
    1025:   Should be named `PCGAMGSetProcessEquationLimit()`.

    1027: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`
    1028: @*/
    1029: PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n)
    1030: {
    1031:   PetscFunctionBegin;
    1033:   PetscTryMethod(pc, "PCGAMGSetProcEqLim_C", (PC, PetscInt), (pc, n));
    1034:   PetscFunctionReturn(PETSC_SUCCESS);
    1035: }

    1037: static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
    1038: {
    1039:   PC_MG   *mg      = (PC_MG *)pc->data;
    1040:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1042:   PetscFunctionBegin;
    1043:   if (n > 0) pc_gamg->min_eq_proc = n;
    1044:   PetscFunctionReturn(PETSC_SUCCESS);
    1045: }

    1047: /*@
    1048:   PCGAMGSetCoarseEqLim - Set maximum number of equations on the coarsest grid of `PCGAMG`

    1050:   Collective

    1052:   Input Parameters:
    1053: + pc - the preconditioner context
    1054: - n  - maximum number of equations to aim for

    1056:   Options Database Key:
    1057: . -pc_gamg_coarse_eq_limit <limit> - set the limit

    1059:   Level: intermediate

    1061:   Note:
    1062:   For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid
    1063:   has less than 1000 unknowns.

    1065: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`,
    1066:           `PCGAMGSetParallelCoarseGridSolve()`
    1067: @*/
    1068: PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
    1069: {
    1070:   PetscFunctionBegin;
    1072:   PetscTryMethod(pc, "PCGAMGSetCoarseEqLim_C", (PC, PetscInt), (pc, n));
    1073:   PetscFunctionReturn(PETSC_SUCCESS);
    1074: }

    1076: static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
    1077: {
    1078:   PC_MG   *mg      = (PC_MG *)pc->data;
    1079:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1081:   PetscFunctionBegin;
    1082:   if (n > 0) pc_gamg->coarse_eq_limit = n;
    1083:   PetscFunctionReturn(PETSC_SUCCESS);
    1084: }

    1086: /*@
    1087:   PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids when reducing the number of MPI ranks to use

    1089:   Collective

    1091:   Input Parameters:
    1092: + pc - the preconditioner context
    1093: - n  - `PETSC_TRUE` or `PETSC_FALSE`

    1095:   Options Database Key:
    1096: . -pc_gamg_repartition <true,false> - turn on the repartitioning

    1098:   Level: intermediate

    1100:   Note:
    1101:   This will generally improve the loading balancing of the work on each level so the solves will be faster but it adds to the
    1102:   preconditioner setup time.

    1104: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`
    1105: @*/
    1106: PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n)
    1107: {
    1108:   PetscFunctionBegin;
    1110:   PetscTryMethod(pc, "PCGAMGSetRepartition_C", (PC, PetscBool), (pc, n));
    1111:   PetscFunctionReturn(PETSC_SUCCESS);
    1112: }

    1114: static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n)
    1115: {
    1116:   PC_MG   *mg      = (PC_MG *)pc->data;
    1117:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1119:   PetscFunctionBegin;
    1120:   pc_gamg->repart = n;
    1121:   PetscFunctionReturn(PETSC_SUCCESS);
    1122: }

    1124: /*@
    1125:   PCGAMGSetUseSAEstEig - Use the eigen estimate from smoothed aggregation for the Chebyshev smoother during the solution process

    1127:   Collective

    1129:   Input Parameters:
    1130: + pc - the preconditioner context
    1131: - b  - flag

    1133:   Options Database Key:
    1134: . -pc_gamg_use_sa_esteig <true,false> - use the eigen estimate

    1136:   Level: advanced

    1138:   Notes:
    1139:   Smoothed aggregation constructs the smoothed prolongator $P = (I - \omega D^{-1} A) T$ where $T$ is the tentative prolongator and $D$ is the diagonal of $A$.
    1140:   Eigenvalue estimates (based on a few `PCCG` or `PCGMRES` iterations) are computed to choose $\omega$ so that this is a stable smoothing operation.
    1141:   If `KSPCHEBYSHEV` with `PCJACOBI` (diagonal) preconditioning is used for smoothing on the finest level, then the eigenvalue estimates
    1142:   can be reused during the solution process.
    1143:   This option is only used when the smoother uses `PCJACOBI`, and should be turned off when a different `PCJacobiType` is used.
    1144:   It became default in PETSc 3.17.

    1146: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `PCGAMGSetRecomputeEstEig()`
    1147: @*/
    1148: PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool b)
    1149: {
    1150:   PetscFunctionBegin;
    1152:   PetscTryMethod(pc, "PCGAMGSetUseSAEstEig_C", (PC, PetscBool), (pc, b));
    1153:   PetscFunctionReturn(PETSC_SUCCESS);
    1154: }

    1156: static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool b)
    1157: {
    1158:   PC_MG   *mg      = (PC_MG *)pc->data;
    1159:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1161:   PetscFunctionBegin;
    1162:   pc_gamg->use_sa_esteig = b;
    1163:   PetscFunctionReturn(PETSC_SUCCESS);
    1164: }

    1166: /*@
    1167:   PCGAMGSetRecomputeEstEig - Set flag for Chebyshev smoothers to recompute the eigen estimates when a new matrix is used

    1169:   Collective

    1171:   Input Parameters:
    1172: + pc - the preconditioner context
    1173: - b  - flag, default is `PETSC_TRUE`

    1175:   Options Database Key:
    1176: . -pc_gamg_recompute_esteig <true> - use the eigen estimate

    1178:   Level: advanced

    1180:   Note:
    1181:   If the matrix changes only slightly in a new solve using ``PETSC_FALSE`` will save time in the setting up of the preconditioner
    1182:   and may not affect the solution time much.

    1184: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`
    1185: @*/
    1186: PetscErrorCode PCGAMGSetRecomputeEstEig(PC pc, PetscBool b)
    1187: {
    1188:   PetscFunctionBegin;
    1190:   PetscTryMethod(pc, "PCGAMGSetRecomputeEstEig_C", (PC, PetscBool), (pc, b));
    1191:   PetscFunctionReturn(PETSC_SUCCESS);
    1192: }

    1194: static PetscErrorCode PCGAMGSetRecomputeEstEig_GAMG(PC pc, PetscBool b)
    1195: {
    1196:   PC_MG   *mg      = (PC_MG *)pc->data;
    1197:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1199:   PetscFunctionBegin;
    1200:   pc_gamg->recompute_esteig = b;
    1201:   PetscFunctionReturn(PETSC_SUCCESS);
    1202: }

    1204: /*@
    1205:   PCGAMGSetEigenvalues - Set WHAT eigenvalues WHY?

    1207:   Collective

    1209:   Input Parameters:
    1210: + pc   - the preconditioner context
    1211: . emax - max eigenvalue
    1212: - emin - min eigenvalue

    1214:   Options Database Key:
    1215: . -pc_gamg_eigenvalues <emin,emax> - estimates of the eigenvalues

    1217:   Level: intermediate

    1219: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetUseSAEstEig()`
    1220: @*/
    1221: PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax, PetscReal emin)
    1222: {
    1223:   PetscFunctionBegin;
    1225:   PetscTryMethod(pc, "PCGAMGSetEigenvalues_C", (PC, PetscReal, PetscReal), (pc, emax, emin));
    1226:   PetscFunctionReturn(PETSC_SUCCESS);
    1227: }

    1229: static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc, PetscReal emax, PetscReal emin)
    1230: {
    1231:   PC_MG   *mg      = (PC_MG *)pc->data;
    1232:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1234:   PetscFunctionBegin;
    1235:   PetscCheck(emax > emin, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Maximum eigenvalue must be larger than minimum: max %g min %g", (double)emax, (double)emin);
    1236:   PetscCheck(emax * emin > 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Both eigenvalues must be of the same sign: max %g min %g", (double)emax, (double)emin);
    1237:   pc_gamg->emax = emax;
    1238:   pc_gamg->emin = emin;
    1239:   PetscFunctionReturn(PETSC_SUCCESS);
    1240: }

    1242: /*@
    1243:   PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding a `PCGAMG` algebraic multigrid preconditioner

    1245:   Collective

    1247:   Input Parameters:
    1248: + pc - the preconditioner context
    1249: - n  - `PETSC_TRUE` or `PETSC_FALSE`

    1251:   Options Database Key:
    1252: . -pc_gamg_reuse_interpolation <true,false> - reuse the previous interpolation

    1254:   Level: intermediate

    1256:   Note:
    1257:   May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows
    1258:   rebuilding the preconditioner quicker.

    1260: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
    1261: @*/
    1262: PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n)
    1263: {
    1264:   PetscFunctionBegin;
    1266:   PetscTryMethod(pc, "PCGAMGSetReuseInterpolation_C", (PC, PetscBool), (pc, n));
    1267:   PetscFunctionReturn(PETSC_SUCCESS);
    1268: }

    1270: static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n)
    1271: {
    1272:   PC_MG   *mg      = (PC_MG *)pc->data;
    1273:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1275:   PetscFunctionBegin;
    1276:   pc_gamg->reuse_prol = n;
    1277:   PetscFunctionReturn(PETSC_SUCCESS);
    1278: }

    1280: /*@
    1281:   PCGAMGASMSetUseAggs - Have the `PCGAMG` smoother on each level use `PCASM` where the aggregates defined by the coarsening process are
    1282:   the subdomains for the additive Schwarz preconditioner used as the smoother

    1284:   Collective

    1286:   Input Parameters:
    1287: + pc  - the preconditioner context
    1288: - flg - `PETSC_TRUE` to use aggregates, `PETSC_FALSE` to not

    1290:   Options Database Key:
    1291: . -pc_gamg_asm_use_agg <true,false> - use aggregates to define the additive Schwarz subdomains

    1293:   Level: intermediate

    1295:   Note:
    1296:   This option automatically sets the preconditioner on the levels to be `PCASM`.

    1298: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCASM`, `PCSetType`
    1299: @*/
    1300: PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg)
    1301: {
    1302:   PetscFunctionBegin;
    1304:   PetscTryMethod(pc, "PCGAMGASMSetUseAggs_C", (PC, PetscBool), (pc, flg));
    1305:   PetscFunctionReturn(PETSC_SUCCESS);
    1306: }

    1308: static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg)
    1309: {
    1310:   PC_MG   *mg      = (PC_MG *)pc->data;
    1311:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1313:   PetscFunctionBegin;
    1314:   pc_gamg->use_aggs_in_asm = flg;
    1315:   PetscFunctionReturn(PETSC_SUCCESS);
    1316: }

    1318: /*@
    1319:   PCGAMGSetParallelCoarseGridSolve - allow a parallel coarse grid solver

    1321:   Collective

    1323:   Input Parameters:
    1324: + pc  - the preconditioner context
    1325: - flg - `PETSC_TRUE` to not force coarse grid onto one processor

    1327:   Options Database Key:
    1328: . -pc_gamg_parallel_coarse_grid_solver - use a parallel coarse grid solver

    1330:   Level: intermediate

    1332: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGSetRankReductionFactors()`
    1333: @*/
    1334: PetscErrorCode PCGAMGSetParallelCoarseGridSolve(PC pc, PetscBool flg)
    1335: {
    1336:   PetscFunctionBegin;
    1338:   PetscTryMethod(pc, "PCGAMGSetParallelCoarseGridSolve_C", (PC, PetscBool), (pc, flg));
    1339:   PetscFunctionReturn(PETSC_SUCCESS);
    1340: }

    1342: static PetscErrorCode PCGAMGSetParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg)
    1343: {
    1344:   PC_MG   *mg      = (PC_MG *)pc->data;
    1345:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1347:   PetscFunctionBegin;
    1348:   pc_gamg->use_parallel_coarse_grid_solver = flg;
    1349:   PetscFunctionReturn(PETSC_SUCCESS);
    1350: }

    1352: /*@
    1353:   PCGAMGSetCpuPinCoarseGrids - pin the coarse grids created in `PCGAMG` to run only on the CPU since the problems may be too small to run efficiently on the GPUs

    1355:   Collective

    1357:   Input Parameters:
    1358: + pc  - the preconditioner context
    1359: - flg - `PETSC_TRUE` to pin coarse grids to the CPU

    1361:   Options Database Key:
    1362: . -pc_gamg_cpu_pin_coarse_grids - pin the coarse grids to the CPU

    1364:   Level: intermediate

    1366: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetParallelCoarseGridSolve()`
    1367: @*/
    1368: PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg)
    1369: {
    1370:   PetscFunctionBegin;
    1372:   PetscTryMethod(pc, "PCGAMGSetCpuPinCoarseGrids_C", (PC, PetscBool), (pc, flg));
    1373:   PetscFunctionReturn(PETSC_SUCCESS);
    1374: }

    1376: static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg)
    1377: {
    1378:   PC_MG   *mg      = (PC_MG *)pc->data;
    1379:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1381:   PetscFunctionBegin;
    1382:   pc_gamg->cpu_pin_coarse_grids = flg;
    1383:   PetscFunctionReturn(PETSC_SUCCESS);
    1384: }

    1386: /*@
    1387:   PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type)

    1389:   Collective

    1391:   Input Parameters:
    1392: + pc  - the preconditioner context
    1393: - flg - `PCGAMGLayoutType` type, either `PCGAMG_LAYOUT_COMPACT` or `PCGAMG_LAYOUT_SPREAD`

    1395:   Options Database Key:
    1396: . -pc_gamg_coarse_grid_layout_type - place the coarse grids with natural ordering

    1398:   Level: intermediate

    1400: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGLayoutType`, `PCGAMG_LAYOUT_COMPACT`, `PCGAMG_LAYOUT_SPREAD`
    1401: @*/
    1402: PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg)
    1403: {
    1404:   PetscFunctionBegin;
    1406:   PetscTryMethod(pc, "PCGAMGSetCoarseGridLayoutType_C", (PC, PCGAMGLayoutType), (pc, flg));
    1407:   PetscFunctionReturn(PETSC_SUCCESS);
    1408: }

    1410: static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg)
    1411: {
    1412:   PC_MG   *mg      = (PC_MG *)pc->data;
    1413:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1415:   PetscFunctionBegin;
    1416:   pc_gamg->layout_type = flg;
    1417:   PetscFunctionReturn(PETSC_SUCCESS);
    1418: }

    1420: /*@
    1421:   PCGAMGSetNlevels -  Sets the maximum number of levels `PCGAMG` will use

    1423:   Collective

    1425:   Input Parameters:
    1426: + pc - the preconditioner
    1427: - n  - the maximum number of levels to use

    1429:   Options Database Key:
    1430: . -pc_mg_levels <n> - set the maximum number of levels to allow

    1432:   Level: intermediate

    1434:   Developer Notes:
    1435:   Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG`

    1437: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
    1438: @*/
    1439: PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
    1440: {
    1441:   PetscFunctionBegin;
    1443:   PetscTryMethod(pc, "PCGAMGSetNlevels_C", (PC, PetscInt), (pc, n));
    1444:   PetscFunctionReturn(PETSC_SUCCESS);
    1445: }

    1447: static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
    1448: {
    1449:   PC_MG   *mg      = (PC_MG *)pc->data;
    1450:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1452:   PetscFunctionBegin;
    1453:   pc_gamg->Nlevels = n;
    1454:   PetscFunctionReturn(PETSC_SUCCESS);
    1455: }

    1457: /*@
    1458:   PCGAMGASMSetHEM -  Sets the number of HEM matching passed

    1460:   Collective

    1462:   Input Parameters:
    1463: + pc - the preconditioner
    1464: - n  - number of HEM matching passed to construct ASM subdomains

    1466:   Options Database Key:
    1467: . -pc_gamg_asm_hem <n> - set the number of HEM matching passed

    1469:   Level: intermediate

    1471:   Developer Notes:
    1472:   Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG`

    1474: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
    1475: @*/
    1476: PetscErrorCode PCGAMGASMSetHEM(PC pc, PetscInt n)
    1477: {
    1478:   PetscFunctionBegin;
    1480:   PetscTryMethod(pc, "PCGAMGASMSetHEM_C", (PC, PetscInt), (pc, n));
    1481:   PetscFunctionReturn(PETSC_SUCCESS);
    1482: }

    1484: static PetscErrorCode PCGAMGASMSetHEM_GAMG(PC pc, PetscInt n)
    1485: {
    1486:   PC_MG   *mg      = (PC_MG *)pc->data;
    1487:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1489:   PetscFunctionBegin;
    1490:   pc_gamg->asm_hem_aggs = n;
    1491:   PetscFunctionReturn(PETSC_SUCCESS);
    1492: }

    1494: /*@
    1495:   PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph

    1497:   Not Collective

    1499:   Input Parameters:
    1500: + pc - the preconditioner context
    1501: . v  - array of threshold values for finest n levels; 0.0 means keep all nonzero entries in the graph; negative means keep even zero entries in the graph
    1502: - n  - number of threshold values provided in array

    1504:   Options Database Key:
    1505: . -pc_gamg_threshold <threshold> - the threshold to drop edges

    1507:   Level: intermediate

    1509:   Notes:
    1510:   Increasing the threshold decreases the rate of coarsening. Conversely reducing the threshold increases the rate of coarsening (aggressive coarsening) and thereby reduces the complexity of the coarse grids, and generally results in slower solver converge rates. Reducing coarse grid complexity reduced the complexity of Galerkin coarse grid construction considerably.
    1511:   Before coarsening or aggregating the graph, `PCGAMG` removes small values from the graph with this threshold, and thus reducing the coupling in the graph and a different (perhaps better) coarser set of points.

    1513:   If `n` is less than the total number of coarsenings (see `PCGAMGSetNlevels()`), then threshold scaling (see `PCGAMGSetThresholdScale()`) is used for each successive coarsening.
    1514:   In this case, `PCGAMGSetThresholdScale()` must be called before `PCGAMGSetThreshold()`.
    1515:   If `n` is greater than the total number of levels, the excess entries in threshold will not be used.

    1517: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetMinDegreeOrderingMISk()`, `PCGAMGSetThresholdScale()`
    1518: @*/
    1519: PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n)
    1520: {
    1521:   PetscFunctionBegin;
    1523:   if (n) PetscAssertPointer(v, 2);
    1524:   PetscTryMethod(pc, "PCGAMGSetThreshold_C", (PC, PetscReal[], PetscInt), (pc, v, n));
    1525:   PetscFunctionReturn(PETSC_SUCCESS);
    1526: }

    1528: static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n)
    1529: {
    1530:   PC_MG   *mg      = (PC_MG *)pc->data;
    1531:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
    1532:   PetscInt i;

    1534:   PetscFunctionBegin;
    1535:   for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i];
    1536:   for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
    1537:   PetscFunctionReturn(PETSC_SUCCESS);
    1538: }

    1540: /*@
    1541:   PCGAMGSetRankReductionFactors - Set a manual schedule for MPI rank reduction on coarse grids

    1543:   Collective

    1545:   Input Parameters:
    1546: + pc - the preconditioner context
    1547: . v  - array of reduction factors. 0 for first value forces a reduction to one process/device on first level in CUDA
    1548: - n  - number of values provided in array

    1550:   Options Database Key:
    1551: . -pc_gamg_rank_reduction_factors <factors> - provide the schedule

    1553:   Level: intermediate

    1555: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetParallelCoarseGridSolve()`
    1556: @*/
    1557: PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n)
    1558: {
    1559:   PetscFunctionBegin;
    1561:   if (n) PetscAssertPointer(v, 2);
    1562:   PetscTryMethod(pc, "PCGAMGSetRankReductionFactors_C", (PC, PetscInt[], PetscInt), (pc, v, n));
    1563:   PetscFunctionReturn(PETSC_SUCCESS);
    1564: }

    1566: static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n)
    1567: {
    1568:   PC_MG   *mg      = (PC_MG *)pc->data;
    1569:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
    1570:   PetscInt i;

    1572:   PetscFunctionBegin;
    1573:   for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i];
    1574:   for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */
    1575:   PetscFunctionReturn(PETSC_SUCCESS);
    1576: }

    1578: /*@
    1579:   PCGAMGSetThresholdScale - Relative threshold reduction at each level

    1581:   Not Collective

    1583:   Input Parameters:
    1584: + pc - the preconditioner context
    1585: - v  - the threshold value reduction, usually < 1.0

    1587:   Options Database Key:
    1588: . -pc_gamg_threshold_scale <v> - set the relative threshold reduction on each level

    1590:   Level: advanced

    1592:   Note:
    1593:   The initial threshold (for an arbitrary number of levels starting from the finest) can be set with `PCGAMGSetThreshold()`.
    1594:   This scaling is used for each subsequent coarsening, but must be called before `PCGAMGSetThreshold()`.

    1596: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`
    1597: @*/
    1598: PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v)
    1599: {
    1600:   PetscFunctionBegin;
    1602:   PetscTryMethod(pc, "PCGAMGSetThresholdScale_C", (PC, PetscReal), (pc, v));
    1603:   PetscFunctionReturn(PETSC_SUCCESS);
    1604: }

    1606: static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v)
    1607: {
    1608:   PC_MG   *mg      = (PC_MG *)pc->data;
    1609:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1611:   PetscFunctionBegin;
    1612:   pc_gamg->threshold_scale = v;
    1613:   PetscFunctionReturn(PETSC_SUCCESS);
    1614: }

    1616: /*@
    1617:   PCGAMGSetType - Set the type of algorithm `PCGAMG` should use

    1619:   Collective

    1621:   Input Parameters:
    1622: + pc   - the preconditioner context
    1623: - type - `PCGAMGAGG`, `PCGAMGGEO`, or `PCGAMGCLASSICAL`

    1625:   Options Database Key:
    1626: . -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply - only agg is supported

    1628:   Level: intermediate

    1630: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMGGetType()`, `PCGAMG`, `PCGAMGType`
    1631: @*/
    1632: PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
    1633: {
    1634:   PetscFunctionBegin;
    1636:   PetscTryMethod(pc, "PCGAMGSetType_C", (PC, PCGAMGType), (pc, type));
    1637:   PetscFunctionReturn(PETSC_SUCCESS);
    1638: }

    1640: /*@
    1641:   PCGAMGGetType - Get the type of algorithm `PCGAMG` will use

    1643:   Collective

    1645:   Input Parameter:
    1646: . pc - the preconditioner context

    1648:   Output Parameter:
    1649: . type - the type of algorithm used

    1651:   Level: intermediate

    1653: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetType()`, `PCGAMGType`
    1654: @*/
    1655: PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type)
    1656: {
    1657:   PetscFunctionBegin;
    1659:   PetscUseMethod(pc, "PCGAMGGetType_C", (PC, PCGAMGType *), (pc, type));
    1660:   PetscFunctionReturn(PETSC_SUCCESS);
    1661: }

    1663: static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type)
    1664: {
    1665:   PC_MG   *mg      = (PC_MG *)pc->data;
    1666:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1668:   PetscFunctionBegin;
    1669:   *type = pc_gamg->type;
    1670:   PetscFunctionReturn(PETSC_SUCCESS);
    1671: }

    1673: static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
    1674: {
    1675:   PC_MG   *mg      = (PC_MG *)pc->data;
    1676:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
    1677:   PetscErrorCode (*r)(PC);

    1679:   PetscFunctionBegin;
    1680:   pc_gamg->type = type;
    1681:   PetscCall(PetscFunctionListFind(GAMGList, type, &r));
    1682:   PetscCheck(r, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown GAMG type %s given", type);
    1683:   if (pc_gamg->ops->destroy) {
    1684:     PetscCall((*pc_gamg->ops->destroy)(pc));
    1685:     PetscCall(PetscMemzero(pc_gamg->ops, sizeof(struct _PCGAMGOps)));
    1686:     pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
    1687:     /* cleaning up common data in pc_gamg - this should disappear someday */
    1688:     pc_gamg->data_cell_cols      = 0;
    1689:     pc_gamg->data_cell_rows      = 0;
    1690:     pc_gamg->orig_data_cell_cols = 0;
    1691:     pc_gamg->orig_data_cell_rows = 0;
    1692:     PetscCall(PetscFree(pc_gamg->data));
    1693:     pc_gamg->data_sz = 0;
    1694:   }
    1695:   PetscCall(PetscFree(pc_gamg->gamg_type_name));
    1696:   PetscCall(PetscStrallocpy(type, &pc_gamg->gamg_type_name));
    1697:   PetscCall((*r)(pc));
    1698:   PetscFunctionReturn(PETSC_SUCCESS);
    1699: }

    1701: static PetscErrorCode PCView_GAMG(PC pc, PetscViewer viewer)
    1702: {
    1703:   PC_MG         *mg       = (PC_MG *)pc->data;
    1704:   PC_MG_Levels **mglevels = mg->levels;
    1705:   PC_GAMG       *pc_gamg  = (PC_GAMG *)mg->innerctx;
    1706:   PetscReal      gc, oc;

    1708:   PetscFunctionBegin;
    1709:   PetscCall(PetscViewerASCIIPrintf(viewer, "    GAMG specific options\n"));
    1710:   PetscCall(PetscViewerASCIIPrintf(viewer, "      Threshold for dropping small values in graph on each level ="));
    1711:   for (PetscInt i = 0; i < mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)pc_gamg->threshold[i]));
    1712:   PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    1713:   PetscCall(PetscViewerASCIIPrintf(viewer, "      Threshold scaling factor for each level not specified = %g\n", (double)pc_gamg->threshold_scale));
    1714:   if (pc_gamg->use_aggs_in_asm) PetscCall(PetscViewerASCIIPrintf(viewer, "      Using aggregates from GAMG coarsening process to define subdomains for PCASM\n")); // this take presedence
    1715:   else if (pc_gamg->asm_hem_aggs) {
    1716:     PetscCall(PetscViewerASCIIPrintf(viewer, "      Using aggregates made with %d applications of heavy edge matching (HEM) to define subdomains for PCASM\n", (int)pc_gamg->asm_hem_aggs));
    1717:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1718:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1719:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1720:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1721:     PetscCall(MatCoarsenView(pc_gamg->asm_crs, viewer));
    1722:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1723:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1724:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1725:   }
    1726:   if (pc_gamg->use_parallel_coarse_grid_solver) PetscCall(PetscViewerASCIIPrintf(viewer, "      Using parallel coarse grid solver (all coarse grid equations not put on one process)\n"));
    1727:   if (pc_gamg->injection_index_size) {
    1728:     PetscCall(PetscViewerASCIIPrintf(viewer, "      Using injection restriction/prolongation on first level, dofs:"));
    1729:     for (int i = 0; i < pc_gamg->injection_index_size; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %d", (int)pc_gamg->injection_index[i]));
    1730:     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    1731:   }
    1732:   if (pc_gamg->ops->view) PetscCall((*pc_gamg->ops->view)(pc, viewer));
    1733:   gc = oc = 0;
    1734:   PetscCall(PCMGGetGridComplexity(pc, &gc, &oc));
    1735:   PetscCall(PetscViewerASCIIPrintf(viewer, "      Complexity:    grid = %g    operator = %g\n", (double)gc, (double)oc));
    1736:   PetscCall(PetscViewerASCIIPrintf(viewer, "      Per-level complexity: op = operator, int = interpolation\n"));
    1737:   PetscCall(PetscViewerASCIIPrintf(viewer, "          #equations  | #active PEs | avg nnz/row op | avg nnz/row int\n"));
    1738:   for (PetscInt i = 0; i < mg->nlevels; i++) {
    1739:     MatInfo   info;
    1740:     Mat       A;
    1741:     PetscReal rd[3];
    1742:     PetscInt  rst, ren, N;

    1744:     PetscCall(KSPGetOperators(mglevels[i]->smoothd, NULL, &A));
    1745:     PetscCall(MatGetOwnershipRange(A, &rst, &ren));
    1746:     PetscCall(MatGetSize(A, &N, NULL));
    1747:     PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
    1748:     rd[0] = (ren - rst > 0) ? 1 : 0;
    1749:     rd[1] = info.nz_used;
    1750:     rd[2] = 0;
    1751:     if (i) {
    1752:       Mat P;
    1753:       PetscCall(PCMGGetInterpolation(pc, i, &P));
    1754:       PetscCall(MatGetInfo(P, MAT_LOCAL, &info));
    1755:       rd[2] = info.nz_used;
    1756:     }
    1757:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, rd, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)pc)));
    1758:     PetscCall(PetscViewerASCIIPrintf(viewer, "     %12" PetscInt_FMT " %12" PetscInt_FMT "   %12" PetscInt_FMT "     %12" PetscInt_FMT "\n", N, (PetscInt)rd[0], (PetscInt)PetscCeilReal(rd[1] / N), (PetscInt)PetscCeilReal(rd[2] / N)));
    1759:   }
    1760:   PetscFunctionReturn(PETSC_SUCCESS);
    1761: }

    1763: /*@
    1764:   PCGAMGSetInjectionIndex - Array of subset of variables per vertex to inject into coarse grid space

    1766:   Logically Collective

    1768:   Input Parameters:
    1769: + pc  - the coarsen context
    1770: . n   - number of indices
    1771: - idx - array of indices

    1773:   Options Database Key:
    1774: . -pc_gamg_injection_index - array of subset of variables per vertex to use for injection coarse grid space

    1776:   Level: intermediate

    1778: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`
    1779: @*/
    1780: PetscErrorCode PCGAMGSetInjectionIndex(PC pc, PetscInt n, PetscInt idx[])
    1781: {
    1782:   PetscFunctionBegin;
    1785:   PetscTryMethod(pc, "PCGAMGSetInjectionIndex_C", (PC, PetscInt, PetscInt[]), (pc, n, idx));
    1786:   PetscFunctionReturn(PETSC_SUCCESS);
    1787: }

    1789: static PetscErrorCode PCGAMGSetInjectionIndex_GAMG(PC pc, PetscInt n, PetscInt idx[])
    1790: {
    1791:   PC_MG   *mg      = (PC_MG *)pc->data;
    1792:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    1794:   PetscFunctionBegin;
    1795:   pc_gamg->injection_index_size = n;
    1796:   PetscCheck(n < MAT_COARSEN_STRENGTH_INDEX_SIZE, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "array size %d larger than max %d", (int)n, MAT_COARSEN_STRENGTH_INDEX_SIZE);
    1797:   for (PetscInt i = 0; i < n; i++) pc_gamg->injection_index[i] = idx[i];
    1798:   PetscFunctionReturn(PETSC_SUCCESS);
    1799: }

    1801: static PetscErrorCode PCSetFromOptions_GAMG(PC pc, PetscOptionItems *PetscOptionsObject)
    1802: {
    1803:   PC_MG             *mg      = (PC_MG *)pc->data;
    1804:   PC_GAMG           *pc_gamg = (PC_GAMG *)mg->innerctx;
    1805:   PetscBool          flag;
    1806:   MPI_Comm           comm;
    1807:   char               prefix[256], tname[32];
    1808:   PetscInt           i, n;
    1809:   const char        *pcpre;
    1810:   static const char *LayoutTypes[] = {"compact", "spread", "PCGAMGLayoutType", "PC_GAMG_LAYOUT", NULL};

    1812:   PetscFunctionBegin;
    1813:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
    1814:   PetscOptionsHeadBegin(PetscOptionsObject, "GAMG options");
    1815:   PetscCall(PetscOptionsInt("-pc_mg_levels", "Number of Levels", "PCMGSetLevels", -1, &n, &flag));
    1816:   PetscCheck(!flag, comm, PETSC_ERR_ARG_WRONG, "Invalid flag -pc_mg_levels. GAMG does not allow the number of levels to be set.");
    1817:   PetscCall(PetscOptionsFList("-pc_gamg_type", "Type of AMG method (only 'agg' supported and useful)", "PCGAMGSetType", GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag));
    1818:   if (flag) PetscCall(PCGAMGSetType(pc, tname));
    1819:   PetscCall(PetscOptionsBool("-pc_gamg_repartition", "Repartion coarse grids", "PCGAMGSetRepartition", pc_gamg->repart, &pc_gamg->repart, NULL));
    1820:   PetscCall(PetscOptionsBool("-pc_gamg_use_sa_esteig", "Use eigen estimate from smoothed aggregation for smoother", "PCGAMGSetUseSAEstEig", pc_gamg->use_sa_esteig, &pc_gamg->use_sa_esteig, NULL));
    1821:   PetscCall(PetscOptionsBool("-pc_gamg_recompute_esteig", "Set flag to recompute eigen estimates for Chebyshev when matrix changes", "PCGAMGSetRecomputeEstEig", pc_gamg->recompute_esteig, &pc_gamg->recompute_esteig, NULL));
    1822:   PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation", "Reuse prolongation operator", "PCGAMGReuseInterpolation", pc_gamg->reuse_prol, &pc_gamg->reuse_prol, NULL));
    1823:   PetscCall(PetscOptionsBool("-pc_gamg_asm_use_agg", "Use aggregation aggregates for ASM smoother", "PCGAMGASMSetUseAggs", pc_gamg->use_aggs_in_asm, &pc_gamg->use_aggs_in_asm, NULL));
    1824:   PetscCall(PetscOptionsBool("-pc_gamg_parallel_coarse_grid_solver", "Use parallel coarse grid solver (otherwise put last grid on one process)", "PCGAMGSetParallelCoarseGridSolve", pc_gamg->use_parallel_coarse_grid_solver, &pc_gamg->use_parallel_coarse_grid_solver, NULL));
    1825:   PetscCall(PetscOptionsBool("-pc_gamg_cpu_pin_coarse_grids", "Pin coarse grids to the CPU", "PCGAMGSetCpuPinCoarseGrids", pc_gamg->cpu_pin_coarse_grids, &pc_gamg->cpu_pin_coarse_grids, NULL));
    1826:   PetscCall(PetscOptionsEnum("-pc_gamg_coarse_grid_layout_type", "compact: place reduced grids on processes in natural order; spread: distribute to whole machine for more memory bandwidth", "PCGAMGSetCoarseGridLayoutType", LayoutTypes,
    1827:                              (PetscEnum)pc_gamg->layout_type, (PetscEnum *)&pc_gamg->layout_type, NULL));
    1828:   PetscCall(PetscOptionsInt("-pc_gamg_process_eq_limit", "Limit (goal) on number of equations per process on coarse grids", "PCGAMGSetProcEqLim", pc_gamg->min_eq_proc, &pc_gamg->min_eq_proc, NULL));
    1829:   PetscCall(PetscOptionsInt("-pc_gamg_coarse_eq_limit", "Limit on number of equations for the coarse grid", "PCGAMGSetCoarseEqLim", pc_gamg->coarse_eq_limit, &pc_gamg->coarse_eq_limit, NULL));
    1830:   PetscCall(PetscOptionsInt("-pc_gamg_asm_hem_aggs", "Number of HEM matching passed in aggregates for ASM smoother", "PCGAMGASMSetHEM", pc_gamg->asm_hem_aggs, &pc_gamg->asm_hem_aggs, NULL));
    1831:   PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale", "Scaling of threshold for each level not specified", "PCGAMGSetThresholdScale", pc_gamg->threshold_scale, &pc_gamg->threshold_scale, NULL));
    1832:   n = PETSC_MG_MAXLEVELS;
    1833:   PetscCall(PetscOptionsRealArray("-pc_gamg_threshold", "Relative threshold to use for dropping edges in aggregation graph", "PCGAMGSetThreshold", pc_gamg->threshold, &n, &flag));
    1834:   if (!flag || n < PETSC_MG_MAXLEVELS) {
    1835:     if (!flag) n = 1;
    1836:     i = n;
    1837:     do {
    1838:       pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
    1839:     } while (++i < PETSC_MG_MAXLEVELS);
    1840:   }
    1841:   PetscCall(PetscOptionsInt("-pc_mg_levels", "Set number of MG levels (should get from base class)", "PCGAMGSetNlevels", pc_gamg->Nlevels, &pc_gamg->Nlevels, NULL));
    1842:   PetscCheck(pc_gamg->Nlevels <= PETSC_MG_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_mg_levels (%d) >= PETSC_MG_MAXLEVELS (%d)", (int)pc_gamg->Nlevels, (int)PETSC_MG_MAXLEVELS);
    1843:   n = PETSC_MG_MAXLEVELS;
    1844:   PetscCall(PetscOptionsIntArray("-pc_gamg_rank_reduction_factors", "Manual schedule of coarse grid reduction factors that overrides internal heuristics (0 for first reduction puts one process/device)", "PCGAMGSetRankReductionFactors", pc_gamg->level_reduction_factors, &n, &flag));
    1845:   if (!flag) i = 0;
    1846:   else i = n;
    1847:   do {
    1848:     pc_gamg->level_reduction_factors[i] = -1;
    1849:   } while (++i < PETSC_MG_MAXLEVELS);
    1850:   {
    1851:     PetscReal eminmax[2] = {0., 0.};
    1852:     n                    = 2;
    1853:     PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues", "extreme eigenvalues for smoothed aggregation", "PCGAMGSetEigenvalues", eminmax, &n, &flag));
    1854:     if (flag) {
    1855:       PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues");
    1856:       PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0]));
    1857:     }
    1858:   }
    1859:   pc_gamg->injection_index_size = MAT_COARSEN_STRENGTH_INDEX_SIZE;
    1860:   PetscCall(PetscOptionsIntArray("-pc_gamg_injection_index", "Array of indices to use to use injection coarse grid space", "PCGAMGSetInjectionIndex", pc_gamg->injection_index, &pc_gamg->injection_index_size, NULL));
    1861:   /* set options for subtype */
    1862:   PetscCall((*pc_gamg->ops->setfromoptions)(pc, PetscOptionsObject));

    1864:   PetscCall(PCGetOptionsPrefix(pc, &pcpre));
    1865:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
    1866:   PetscOptionsHeadEnd();
    1867:   PetscFunctionReturn(PETSC_SUCCESS);
    1868: }

    1870: /*MC
    1871:   PCGAMG - Geometric algebraic multigrid (AMG) preconditioner

    1873:   Options Database Keys:
    1874: + -pc_gamg_type <type,default=agg>                 - one of agg, geo, or classical (only smoothed aggregation, agg, supported)
    1875: . -pc_gamg_repartition  <bool,default=false>       - repartition the degrees of freedom across the coarse grids as they are determined
    1876: . -pc_gamg_asm_use_agg <bool,default=false>        - use the aggregates from the coasening process to defined the subdomains on each level for the `PCASM` smoother.
    1877:                                                      That is using `-mg_levels_pc_type asm`
    1878: . -pc_gamg_process_eq_limit <limit, default=50>    - `PCGAMG` will reduce the number of MPI ranks used directly on the coarse grids so that there are around <limit>
    1879:                                                      equations on each process that has degrees of freedom
    1880: . -pc_gamg_coarse_eq_limit <limit, default=50>     - Set maximum number of equations on coarsest grid to aim for.
    1881: . -pc_gamg_reuse_interpolation <bool,default=true> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations (should always be true)
    1882: . -pc_gamg_threshold[] <thresh,default=[-1,...]>   - Before aggregating the graph `PCGAMG` will remove small values from the graph on each level (< 0 does no filtering)
    1883: - -pc_gamg_threshold_scale <scale,default=1>       - Scaling of threshold on each coarser grid if not specified

    1885:   Options Database Keys for Aggregation:
    1886: + -pc_gamg_agg_nsmooths <nsmooth, default=1>                 - number of smoothing steps to use with smooth aggregation to construct prolongation
    1887: . -pc_gamg_aggressive_coarsening <n,default=1>               - number of aggressive coarsening (MIS-2) levels from finest.
    1888: . -pc_gamg_aggressive_square_graph <bool,default=false>      - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening
    1889: . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm
    1890: . -pc_gamg_pc_gamg_asm_hem_aggs <n,default=0>                - Number of HEM aggregation steps for `PCASM` smoother
    1891: - -pc_gamg_aggressive_mis_k <n,default=2>                    - Number (k) distance in MIS coarsening (>2 is 'aggressive')

    1893:   Options Database Keys for Multigrid:
    1894: + -pc_mg_cycle_type <v>        - v or w, see `PCMGSetCycleType()`
    1895: . -pc_mg_distinct_smoothup     - configure the up and down (pre and post) smoothers separately, see PCMGSetDistinctSmoothUp()
    1896: . -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade
    1897: - -pc_mg_levels <levels>       - Number of levels of multigrid to use. GAMG has a heuristic so pc_mg_levels is not usually used with GAMG

    1899:   Level: intermediate

    1901:   Notes:
    1902:   To obtain good performance for `PCGAMG` for vector valued problems you must
    1903:   call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point
    1904:   call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator

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

    1908: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
    1909:           `MatSetBlockSize()`,
    1910:           `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`,
    1911:           `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`,
    1912:           `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
    1913: M*/
    1914: PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
    1915: {
    1916:   PC_GAMG *pc_gamg;
    1917:   PC_MG   *mg;

    1919:   PetscFunctionBegin;
    1920:   /* register AMG type */
    1921:   PetscCall(PCGAMGInitializePackage());

    1923:   /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
    1924:   PetscCall(PCSetType(pc, PCMG));
    1925:   PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG));

    1927:   /* create a supporting struct and attach it to pc */
    1928:   PetscCall(PetscNew(&pc_gamg));
    1929:   PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
    1930:   mg           = (PC_MG *)pc->data;
    1931:   mg->innerctx = pc_gamg;

    1933:   PetscCall(PetscNew(&pc_gamg->ops));

    1935:   /* these should be in subctx but repartitioning needs simple arrays */
    1936:   pc_gamg->data_sz = 0;
    1937:   pc_gamg->data    = NULL;

    1939:   /* overwrite the pointers of PCMG by the functions of base class PCGAMG */
    1940:   pc->ops->setfromoptions = PCSetFromOptions_GAMG;
    1941:   pc->ops->setup          = PCSetUp_GAMG;
    1942:   pc->ops->reset          = PCReset_GAMG;
    1943:   pc->ops->destroy        = PCDestroy_GAMG;
    1944:   mg->view                = PCView_GAMG;

    1946:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", PCGAMGSetProcEqLim_GAMG));
    1947:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", PCGAMGSetCoarseEqLim_GAMG));
    1948:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", PCGAMGSetRepartition_GAMG));
    1949:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", PCGAMGSetEigenvalues_GAMG));
    1950:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", PCGAMGSetUseSAEstEig_GAMG));
    1951:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", PCGAMGSetRecomputeEstEig_GAMG));
    1952:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", PCGAMGSetReuseInterpolation_GAMG));
    1953:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", PCGAMGASMSetUseAggs_GAMG));
    1954:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", PCGAMGSetParallelCoarseGridSolve_GAMG));
    1955:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", PCGAMGSetCpuPinCoarseGrids_GAMG));
    1956:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", PCGAMGSetCoarseGridLayoutType_GAMG));
    1957:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", PCGAMGSetThreshold_GAMG));
    1958:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", PCGAMGSetRankReductionFactors_GAMG));
    1959:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", PCGAMGSetThresholdScale_GAMG));
    1960:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", PCGAMGSetType_GAMG));
    1961:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", PCGAMGGetType_GAMG));
    1962:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", PCGAMGSetNlevels_GAMG));
    1963:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", PCGAMGASMSetHEM_GAMG));
    1964:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", PCGAMGSetInjectionIndex_GAMG));
    1965:   pc_gamg->repart                          = PETSC_FALSE;
    1966:   pc_gamg->reuse_prol                      = PETSC_TRUE;
    1967:   pc_gamg->use_aggs_in_asm                 = PETSC_FALSE;
    1968:   pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE;
    1969:   pc_gamg->cpu_pin_coarse_grids            = PETSC_FALSE;
    1970:   pc_gamg->layout_type                     = PCGAMG_LAYOUT_SPREAD;
    1971:   pc_gamg->min_eq_proc                     = 50;
    1972:   pc_gamg->asm_hem_aggs                    = 0;
    1973:   pc_gamg->coarse_eq_limit                 = 50;
    1974:   for (int i = 0; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = -1;
    1975:   pc_gamg->threshold_scale  = 1.;
    1976:   pc_gamg->Nlevels          = PETSC_MG_MAXLEVELS;
    1977:   pc_gamg->current_level    = 0; /* don't need to init really */
    1978:   pc_gamg->use_sa_esteig    = PETSC_TRUE;
    1979:   pc_gamg->recompute_esteig = PETSC_TRUE;
    1980:   pc_gamg->emin             = 0;
    1981:   pc_gamg->emax             = 0;

    1983:   pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;

    1985:   /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */
    1986:   PetscCall(PCGAMGSetType(pc, PCGAMGAGG));
    1987:   PetscFunctionReturn(PETSC_SUCCESS);
    1988: }

    1990: /*@C
    1991:   PCGAMGInitializePackage - This function initializes everything in the `PCGAMG` package. It is called
    1992:   from `PCInitializePackage()`.

    1994:   Level: developer

    1996: .seealso: [](ch_ksp), `PetscInitialize()`
    1997: @*/
    1998: PetscErrorCode PCGAMGInitializePackage(void)
    1999: {
    2000:   PetscInt l;

    2002:   PetscFunctionBegin;
    2003:   if (PCGAMGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
    2004:   PCGAMGPackageInitialized = PETSC_TRUE;
    2005:   PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGGEO, PCCreateGAMG_GEO));
    2006:   PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGAGG, PCCreateGAMG_AGG));
    2007:   PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGCLASSICAL, PCCreateGAMG_Classical));
    2008:   PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage));

    2010:   /* general events */
    2011:   PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP]));
    2012:   PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH]));
    2013:   PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN]));
    2014:   PetscCall(PetscLogEventRegister("  GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS]));
    2015:   PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL]));
    2016:   PetscCall(PetscLogEventRegister("  GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA]));
    2017:   PetscCall(PetscLogEventRegister("  GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB]));
    2018:   PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT]));
    2019:   PetscCall(PetscLogEventRegister("  GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM]));
    2020:   PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL]));
    2021:   PetscCall(PetscLogEventRegister("  GAMG PtAP", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP]));
    2022:   PetscCall(PetscLogEventRegister("  GAMG Reduce", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE]));
    2023:   PetscCall(PetscLogEventRegister("   GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART]));
    2024:   /* PetscCall(PetscLogEventRegister("   GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */
    2025:   /* PetscCall(PetscLogEventRegister("   GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */
    2026:   /* PetscCall(PetscLogEventRegister("   GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */
    2027:   for (l = 0; l < PETSC_MG_MAXLEVELS; l++) {
    2028:     char ename[32];

    2030:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Squ l%02" PetscInt_FMT, l));
    2031:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0]));
    2032:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Gal l%02" PetscInt_FMT, l));
    2033:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1]));
    2034:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Opt l%02" PetscInt_FMT, l));
    2035:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2]));
    2036:   }
    2037: #if defined(GAMG_STAGES)
    2038:   { /* create timer stages */
    2039:     char str[32];
    2040:     PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", 0));
    2041:     PetscCall(PetscLogStageRegister(str, &gamg_stages[0]));
    2042:   }
    2043: #endif
    2044:   PetscFunctionReturn(PETSC_SUCCESS);
    2045: }

    2047: /*@C
    2048:   PCGAMGFinalizePackage - This function frees everything from the `PCGAMG` package. It is
    2049:   called from `PetscFinalize()` automatically.

    2051:   Level: developer

    2053: .seealso: [](ch_ksp), `PetscFinalize()`
    2054: @*/
    2055: PetscErrorCode PCGAMGFinalizePackage(void)
    2056: {
    2057:   PetscFunctionBegin;
    2058:   PCGAMGPackageInitialized = PETSC_FALSE;
    2059:   PetscCall(PetscFunctionListDestroy(&GAMGList));
    2060:   PetscFunctionReturn(PETSC_SUCCESS);
    2061: }

    2063: /*@C
    2064:   PCGAMGRegister - Register a `PCGAMG` implementation.

    2066:   Input Parameters:
    2067: + type   - string that will be used as the name of the `PCGAMG` type.
    2068: - create - function for creating the gamg context.

    2070:   Level: developer

    2072: .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
    2073: @*/
    2074: PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC))
    2075: {
    2076:   PetscFunctionBegin;
    2077:   PetscCall(PCGAMGInitializePackage());
    2078:   PetscCall(PetscFunctionListAdd(&GAMGList, type, create));
    2079:   PetscFunctionReturn(PETSC_SUCCESS);
    2080: }

    2082: /*@
    2083:   PCGAMGCreateGraph - Creates a graph that is used by the ``PCGAMGType`` in the coarsening process

    2085:   Input Parameters:
    2086: + pc - the `PCGAMG`
    2087: - A  - the matrix, for any level

    2089:   Output Parameter:
    2090: . G - the graph

    2092:   Level: advanced

    2094: .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
    2095: @*/
    2096: PetscErrorCode PCGAMGCreateGraph(PC pc, Mat A, Mat *G)
    2097: {
    2098:   PC_MG   *mg      = (PC_MG *)pc->data;
    2099:   PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;

    2101:   PetscFunctionBegin;
    2102:   PetscCall(pc_gamg->ops->creategraph(pc, A, G));
    2103:   PetscFunctionReturn(PETSC_SUCCESS);
    2104: }