1 /* 2 GAMG geometric-algebric multigrid PC - Mark Adams 2011 3 */ 4 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 5 #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h> /*I "petscksp.h" I*/ 6 7 #if defined(PETSC_HAVE_CUDA) 8 #include <cuda_runtime.h> 9 #endif 10 11 #if defined(PETSC_HAVE_HIP) 12 #include <hip/hip_runtime.h> 13 #endif 14 15 PetscLogEvent petsc_gamg_setup_events[GAMG_NUM_SET]; 16 PetscLogEvent petsc_gamg_setup_matmat_events[PETSC_MG_MAXLEVELS][3]; 17 18 // #define GAMG_STAGES 19 #if defined(GAMG_STAGES) 20 static PetscLogStage gamg_stages[PETSC_MG_MAXLEVELS]; 21 #endif 22 23 static PetscFunctionList GAMGList = NULL; 24 static PetscBool PCGAMGPackageInitialized; 25 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; 30 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 PetscFunctionReturn(PETSC_SUCCESS); 43 } 44 45 /* 46 PCGAMGCreateLevel_GAMG: create coarse op with RAP. repartition and/or reduce number 47 of active processors. 48 49 Input Parameter: 50 . pc - parameters + side effect: coarse data in 'pc_gamg->data' and 51 'pc_gamg->data_sz' are changed via repartitioning/reduction. 52 . Amat_fine - matrix on this fine (k) level 53 . cr_bs - coarse block size 54 In/Output Parameter: 55 . a_P_inout - prolongation operator to the next level (k-->k-1) 56 . a_nactive_proc - number of active procs 57 Output Parameter: 58 . a_Amat_crs - coarse matrix that is created (k-1) 59 */ 60 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) 61 { 62 PC_MG *mg = (PC_MG *)pc->data; 63 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 64 Mat Cmat, Pold = *a_P_inout; 65 MPI_Comm comm; 66 PetscMPIInt rank, size, new_size, nactive = *a_nactive_proc; 67 PetscInt ncrs_eq, ncrs, f_bs; 68 69 PetscFunctionBegin; 70 PetscCall(PetscObjectGetComm((PetscObject)Amat_fine, &comm)); 71 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 72 PetscCallMPI(MPI_Comm_size(comm, &size)); 73 PetscCall(MatGetBlockSize(Amat_fine, &f_bs)); 74 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0)); 75 PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0)); 76 PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat)); 77 PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0)); 78 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0)); 79 80 if (Pcolumnperm) *Pcolumnperm = NULL; 81 82 /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/ 83 PetscCall(MatGetLocalSize(Cmat, &ncrs_eq, NULL)); 84 if (pc_gamg->data_cell_rows > 0) { 85 ncrs = pc_gamg->data_sz / pc_gamg->data_cell_cols / pc_gamg->data_cell_rows; 86 } else { 87 PetscInt bs; 88 PetscCall(MatGetBlockSize(Cmat, &bs)); 89 ncrs = ncrs_eq / bs; 90 } 91 /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */ 92 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. */ 93 #if defined(PETSC_HAVE_CUDA) 94 PetscShmComm pshmcomm; 95 PetscMPIInt locrank; 96 MPI_Comm loccomm; 97 PetscInt s_nnodes, r_nnodes, new_new_size; 98 cudaError_t cerr; 99 int devCount; 100 PetscCall(PetscShmCommGet(comm, &pshmcomm)); 101 PetscCall(PetscShmCommGetMpiShmComm(pshmcomm, &loccomm)); 102 PetscCallMPI(MPI_Comm_rank(loccomm, &locrank)); 103 s_nnodes = !locrank; 104 PetscCall(MPIU_Allreduce(&s_nnodes, &r_nnodes, 1, MPIU_INT, MPI_SUM, comm)); 105 PetscCheck((size % r_nnodes) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "odd number of nodes np=%d nnodes%" PetscInt_FMT, size, r_nnodes); 106 devCount = 0; 107 cerr = cudaGetDeviceCount(&devCount); 108 cudaGetLastError(); /* Reset the last error */ 109 if (cerr == cudaSuccess && devCount >= 1) { /* There are devices, else go to heuristic */ 110 new_new_size = r_nnodes * devCount; 111 new_size = new_new_size; 112 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)); 113 } else { 114 PetscCall(PetscInfo(pc, "%s: With Cuda but no device. Use heuristics.\n", ((PetscObject)pc)->prefix)); 115 goto HEURISTIC; 116 } 117 #else 118 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "should not be here"); 119 #endif 120 } else if (pc_gamg->level_reduction_factors[pc_gamg->current_level] > 0) { 121 if (nactive < pc_gamg->level_reduction_factors[pc_gamg->current_level]) { 122 new_size = 1; 123 PetscCall(PetscInfo(pc, "%s: reduction factor too small for %d active processes: reduce to one process\n", ((PetscObject)pc)->prefix, new_size)); 124 } else { 125 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]); 126 new_size = nactive / pc_gamg->level_reduction_factors[pc_gamg->current_level]; 127 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])); 128 } 129 } else if (is_last && !pc_gamg->use_parallel_coarse_grid_solver) { 130 new_size = 1; 131 PetscCall(PetscInfo(pc, "%s: Force coarsest grid reduction to %d active processes\n", ((PetscObject)pc)->prefix, new_size)); 132 } else { 133 PetscInt ncrs_eq_glob; 134 #if defined(PETSC_HAVE_CUDA) 135 HEURISTIC: 136 #endif 137 PetscCall(MatGetSize(Cmat, &ncrs_eq_glob, NULL)); 138 new_size = (PetscMPIInt)((float)ncrs_eq_glob / (float)pc_gamg->min_eq_proc + 0.5); /* hardwire min. number of eq/proc */ 139 if (!new_size) new_size = 1; /* not likely, possible? */ 140 else if (new_size >= nactive) new_size = nactive; /* no change, rare */ 141 PetscCall(PetscInfo(pc, "%s: Coarse grid reduction from %d to %d active processes\n", ((PetscObject)pc)->prefix, nactive, new_size)); 142 } 143 if (new_size == nactive) { 144 *a_Amat_crs = Cmat; /* output - no repartitioning or reduction - could bail here */ 145 if (new_size < size) { 146 /* odd case where multiple coarse grids are on one processor or no coarsening ... */ 147 PetscCall(PetscInfo(pc, "%s: reduced grid using same number of processors (%d) as last grid (use larger coarse grid)\n", ((PetscObject)pc)->prefix, nactive)); 148 if (pc_gamg->cpu_pin_coarse_grids) { 149 PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE)); 150 PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE)); 151 } 152 } 153 /* we know that the grid structure can be reused in MatPtAP */ 154 } else { /* reduce active processors - we know that the grid structure can NOT be reused in MatPtAP */ 155 PetscInt *counts, *newproc_idx, ii, jj, kk, strideNew, *tidx, ncrs_new, ncrs_eq_new, nloc_old, expand_factor = 1, rfactor = 1; 156 IS is_eq_newproc, is_eq_num, is_eq_num_prim, new_eq_indices; 157 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0)); 158 nloc_old = ncrs_eq / cr_bs; 159 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); 160 /* get new_size and rfactor */ 161 if (pc_gamg->layout_type == PCGAMG_LAYOUT_SPREAD || !pc_gamg->repart) { 162 /* find factor */ 163 if (new_size == 1) rfactor = size; /* don't modify */ 164 else { 165 PetscReal best_fact = 0.; 166 jj = -1; 167 for (kk = 1; kk <= size; kk++) { 168 if (!(size % kk)) { /* a candidate */ 169 PetscReal nactpe = (PetscReal)size / (PetscReal)kk, fact = nactpe / (PetscReal)new_size; 170 if (fact > 1.0) fact = 1. / fact; /* keep fact < 1 */ 171 if (fact > best_fact) { 172 best_fact = fact; 173 jj = kk; 174 } 175 } 176 } 177 if (jj != -1) rfactor = jj; 178 else rfactor = 1; /* a prime */ 179 if (pc_gamg->layout_type == PCGAMG_LAYOUT_COMPACT) expand_factor = 1; 180 else expand_factor = rfactor; 181 } 182 new_size = size / rfactor; /* make new size one that is factor */ 183 if (new_size == nactive) { /* no repartitioning or reduction, bail out because nested here (rare) */ 184 *a_Amat_crs = Cmat; 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 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 } 190 /* make 'is_eq_newproc' */ 191 PetscCall(PetscMalloc1(size, &counts)); 192 if (pc_gamg->repart) { /* Repartition Cmat_{k} and move columns of P^{k}_{k-1} and coordinates of primal part accordingly */ 193 Mat adj; 194 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0)); 195 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")); 196 /* get 'adj' */ 197 if (cr_bs == 1) { 198 PetscCall(MatConvert(Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj)); 199 } else { 200 /* make a scalar matrix to partition (no Stokes here) */ 201 Mat tMat; 202 PetscInt Istart_crs, Iend_crs, ncols, jj, Ii; 203 const PetscScalar *vals; 204 const PetscInt *idx; 205 PetscInt *d_nnz, *o_nnz, M, N; 206 static PetscInt llev = 0; /* ugly but just used for debugging */ 207 MatType mtype; 208 209 PetscCall(PetscMalloc2(ncrs, &d_nnz, ncrs, &o_nnz)); 210 PetscCall(MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs)); 211 PetscCall(MatGetSize(Cmat, &M, &N)); 212 for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) { 213 PetscCall(MatGetRow(Cmat, Ii, &ncols, NULL, NULL)); 214 d_nnz[jj] = ncols / cr_bs; 215 o_nnz[jj] = ncols / cr_bs; 216 PetscCall(MatRestoreRow(Cmat, Ii, &ncols, NULL, NULL)); 217 if (d_nnz[jj] > ncrs) d_nnz[jj] = ncrs; 218 if (o_nnz[jj] > (M / cr_bs - ncrs)) o_nnz[jj] = M / cr_bs - ncrs; 219 } 220 221 PetscCall(MatGetType(Amat_fine, &mtype)); 222 PetscCall(MatCreate(comm, &tMat)); 223 PetscCall(MatSetSizes(tMat, ncrs, ncrs, PETSC_DETERMINE, PETSC_DETERMINE)); 224 PetscCall(MatSetType(tMat, mtype)); 225 PetscCall(MatSeqAIJSetPreallocation(tMat, 0, d_nnz)); 226 PetscCall(MatMPIAIJSetPreallocation(tMat, 0, d_nnz, 0, o_nnz)); 227 PetscCall(PetscFree2(d_nnz, o_nnz)); 228 229 for (ii = Istart_crs; ii < Iend_crs; ii++) { 230 PetscInt dest_row = ii / cr_bs; 231 PetscCall(MatGetRow(Cmat, ii, &ncols, &idx, &vals)); 232 for (jj = 0; jj < ncols; jj++) { 233 PetscInt dest_col = idx[jj] / cr_bs; 234 PetscScalar v = 1.0; 235 PetscCall(MatSetValues(tMat, 1, &dest_row, 1, &dest_col, &v, ADD_VALUES)); 236 } 237 PetscCall(MatRestoreRow(Cmat, ii, &ncols, &idx, &vals)); 238 } 239 PetscCall(MatAssemblyBegin(tMat, MAT_FINAL_ASSEMBLY)); 240 PetscCall(MatAssemblyEnd(tMat, MAT_FINAL_ASSEMBLY)); 241 242 if (llev++ == -1) { 243 PetscViewer viewer; 244 char fname[32]; 245 PetscCall(PetscSNPrintf(fname, sizeof(fname), "part_mat_%" PetscInt_FMT ".mat", llev)); 246 PetscCall(PetscViewerBinaryOpen(comm, fname, FILE_MODE_WRITE, &viewer)); 247 PetscCall(MatView(tMat, viewer)); 248 PetscCall(PetscViewerDestroy(&viewer)); 249 } 250 PetscCall(MatConvert(tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj)); 251 PetscCall(MatDestroy(&tMat)); 252 } /* create 'adj' */ 253 254 { /* partition: get newproc_idx */ 255 char prefix[256]; 256 const char *pcpre; 257 const PetscInt *is_idx; 258 MatPartitioning mpart; 259 IS proc_is; 260 261 PetscCall(MatPartitioningCreate(comm, &mpart)); 262 PetscCall(MatPartitioningSetAdjacency(mpart, adj)); 263 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 264 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : "")); 265 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix)); 266 PetscCall(MatPartitioningSetFromOptions(mpart)); 267 PetscCall(MatPartitioningSetNParts(mpart, new_size)); 268 PetscCall(MatPartitioningApply(mpart, &proc_is)); 269 PetscCall(MatPartitioningDestroy(&mpart)); 270 271 /* collect IS info */ 272 PetscCall(PetscMalloc1(ncrs_eq, &newproc_idx)); 273 PetscCall(ISGetIndices(proc_is, &is_idx)); 274 for (kk = jj = 0; kk < nloc_old; kk++) { 275 for (ii = 0; ii < cr_bs; ii++, jj++) { newproc_idx[jj] = is_idx[kk] * expand_factor; /* distribution */ } 276 } 277 PetscCall(ISRestoreIndices(proc_is, &is_idx)); 278 PetscCall(ISDestroy(&proc_is)); 279 } 280 PetscCall(MatDestroy(&adj)); 281 282 PetscCall(ISCreateGeneral(comm, ncrs_eq, newproc_idx, PETSC_COPY_VALUES, &is_eq_newproc)); 283 PetscCall(PetscFree(newproc_idx)); 284 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0)); 285 } else { /* simple aggregation of parts -- 'is_eq_newproc' */ 286 PetscInt targetPE; 287 PetscCheck(new_size != nactive, PETSC_COMM_SELF, PETSC_ERR_PLIB, "new_size==nactive. Should not happen"); 288 PetscCall(PetscInfo(pc, "%s: Number of equations (loc) %" PetscInt_FMT " with simple aggregation\n", ((PetscObject)pc)->prefix, ncrs_eq)); 289 targetPE = (rank / rfactor) * expand_factor; 290 PetscCall(ISCreateStride(comm, ncrs_eq, targetPE, 0, &is_eq_newproc)); 291 } /* end simple 'is_eq_newproc' */ 292 293 /* 294 Create an index set from the is_eq_newproc index set to indicate the mapping TO 295 */ 296 PetscCall(ISPartitioningToNumbering(is_eq_newproc, &is_eq_num)); 297 is_eq_num_prim = is_eq_num; 298 /* 299 Determine how many equations/vertices are assigned to each processor 300 */ 301 PetscCall(ISPartitioningCount(is_eq_newproc, size, counts)); 302 ncrs_eq_new = counts[rank]; 303 PetscCall(ISDestroy(&is_eq_newproc)); 304 ncrs_new = ncrs_eq_new / cr_bs; 305 306 PetscCall(PetscFree(counts)); 307 /* 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 */ 308 { 309 Vec src_crd, dest_crd; 310 const PetscInt *idx, ndata_rows = pc_gamg->data_cell_rows, ndata_cols = pc_gamg->data_cell_cols, node_data_sz = ndata_rows * ndata_cols; 311 VecScatter vecscat; 312 PetscScalar *array; 313 IS isscat; 314 /* move data (for primal equations only) */ 315 /* Create a vector to contain the newly ordered element information */ 316 PetscCall(VecCreate(comm, &dest_crd)); 317 PetscCall(VecSetSizes(dest_crd, node_data_sz * ncrs_new, PETSC_DECIDE)); 318 PetscCall(VecSetType(dest_crd, VECSTANDARD)); /* this is needed! */ 319 /* 320 There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having 321 a block size of ...). Note, ISs are expanded into equation space by 'cr_bs'. 322 */ 323 PetscCall(PetscMalloc1(ncrs * node_data_sz, &tidx)); 324 PetscCall(ISGetIndices(is_eq_num_prim, &idx)); 325 for (ii = 0, jj = 0; ii < ncrs; ii++) { 326 PetscInt id = idx[ii * cr_bs] / cr_bs; /* get node back */ 327 for (kk = 0; kk < node_data_sz; kk++, jj++) tidx[jj] = id * node_data_sz + kk; 328 } 329 PetscCall(ISRestoreIndices(is_eq_num_prim, &idx)); 330 PetscCall(ISCreateGeneral(comm, node_data_sz * ncrs, tidx, PETSC_COPY_VALUES, &isscat)); 331 PetscCall(PetscFree(tidx)); 332 /* 333 Create a vector to contain the original vertex information for each element 334 */ 335 PetscCall(VecCreateSeq(PETSC_COMM_SELF, node_data_sz * ncrs, &src_crd)); 336 for (jj = 0; jj < ndata_cols; jj++) { 337 const PetscInt stride0 = ncrs * pc_gamg->data_cell_rows; 338 for (ii = 0; ii < ncrs; ii++) { 339 for (kk = 0; kk < ndata_rows; kk++) { 340 PetscInt ix = ii * ndata_rows + kk + jj * stride0, jx = ii * node_data_sz + kk * ndata_cols + jj; 341 PetscScalar tt = (PetscScalar)pc_gamg->data[ix]; 342 PetscCall(VecSetValues(src_crd, 1, &jx, &tt, INSERT_VALUES)); 343 } 344 } 345 } 346 PetscCall(VecAssemblyBegin(src_crd)); 347 PetscCall(VecAssemblyEnd(src_crd)); 348 /* 349 Scatter the element vertex information (still in the original vertex ordering) 350 to the correct processor 351 */ 352 PetscCall(VecScatterCreate(src_crd, NULL, dest_crd, isscat, &vecscat)); 353 PetscCall(ISDestroy(&isscat)); 354 PetscCall(VecScatterBegin(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD)); 355 PetscCall(VecScatterEnd(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD)); 356 PetscCall(VecScatterDestroy(&vecscat)); 357 PetscCall(VecDestroy(&src_crd)); 358 /* 359 Put the element vertex data into a new allocation of the gdata->ele 360 */ 361 PetscCall(PetscFree(pc_gamg->data)); 362 PetscCall(PetscMalloc1(node_data_sz * ncrs_new, &pc_gamg->data)); 363 364 pc_gamg->data_sz = node_data_sz * ncrs_new; 365 strideNew = ncrs_new * ndata_rows; 366 367 PetscCall(VecGetArray(dest_crd, &array)); 368 for (jj = 0; jj < ndata_cols; jj++) { 369 for (ii = 0; ii < ncrs_new; ii++) { 370 for (kk = 0; kk < ndata_rows; kk++) { 371 PetscInt ix = ii * ndata_rows + kk + jj * strideNew, jx = ii * node_data_sz + kk * ndata_cols + jj; 372 pc_gamg->data[ix] = PetscRealPart(array[jx]); 373 } 374 } 375 } 376 PetscCall(VecRestoreArray(dest_crd, &array)); 377 PetscCall(VecDestroy(&dest_crd)); 378 } 379 /* move A and P (columns) with new layout */ 380 /* 381 Invert for MatCreateSubMatrix 382 */ 383 PetscCall(ISInvertPermutation(is_eq_num, ncrs_eq_new, &new_eq_indices)); 384 PetscCall(ISSort(new_eq_indices)); /* is this needed? */ 385 PetscCall(ISSetBlockSize(new_eq_indices, cr_bs)); 386 if (is_eq_num != is_eq_num_prim) { PetscCall(ISDestroy(&is_eq_num_prim)); /* could be same as 'is_eq_num' */ } 387 if (Pcolumnperm) { 388 PetscCall(PetscObjectReference((PetscObject)new_eq_indices)); 389 *Pcolumnperm = new_eq_indices; 390 } 391 PetscCall(ISDestroy(&is_eq_num)); 392 393 /* 'a_Amat_crs' output */ 394 { 395 Mat mat; 396 PetscBool isset, isspd, isher; 397 #if !defined(PETSC_USE_COMPLEX) 398 PetscBool issym; 399 #endif 400 401 PetscCall(MatCreateSubMatrix(Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat)); 402 PetscCall(MatIsSPDKnown(Cmat, &isset, &isspd)); // like MatPropagateSymmetryOptions, but should set MAT_STRUCTURALLY_SYMMETRIC ? 403 if (isset) PetscCall(MatSetOption(mat, MAT_SPD, isspd)); 404 else { 405 PetscCall(MatIsHermitianKnown(Cmat, &isset, &isher)); 406 if (isset) PetscCall(MatSetOption(mat, MAT_HERMITIAN, isher)); 407 else { 408 #if !defined(PETSC_USE_COMPLEX) 409 PetscCall(MatIsSymmetricKnown(Cmat, &isset, &issym)); 410 if (isset) PetscCall(MatSetOption(mat, MAT_SYMMETRIC, issym)); 411 #endif 412 } 413 } 414 *a_Amat_crs = mat; 415 } 416 PetscCall(MatDestroy(&Cmat)); 417 418 /* prolongator */ 419 { 420 IS findices; 421 PetscInt Istart, Iend; 422 Mat Pnew; 423 424 PetscCall(MatGetOwnershipRange(Pold, &Istart, &Iend)); 425 PetscCall(ISCreateStride(comm, Iend - Istart, Istart, 1, &findices)); 426 PetscCall(ISSetBlockSize(findices, f_bs)); 427 PetscCall(MatCreateSubMatrix(Pold, findices, new_eq_indices, MAT_INITIAL_MATRIX, &Pnew)); 428 PetscCall(ISDestroy(&findices)); 429 PetscCall(MatSetOption(Pnew, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE)); 430 431 PetscCall(MatDestroy(a_P_inout)); 432 433 /* output - repartitioned */ 434 *a_P_inout = Pnew; 435 } 436 PetscCall(ISDestroy(&new_eq_indices)); 437 438 *a_nactive_proc = new_size; /* output */ 439 440 /* pinning on reduced grids, not a bad heuristic and optimization gets folded into process reduction optimization */ 441 if (pc_gamg->cpu_pin_coarse_grids) { 442 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 443 static PetscInt llev = 2; 444 PetscCall(PetscInfo(pc, "%s: Pinning level %" PetscInt_FMT " to the CPU\n", ((PetscObject)pc)->prefix, llev++)); 445 #endif 446 PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE)); 447 PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE)); 448 if (1) { /* HACK: move this to MatBindCPU_MPIAIJXXX; lvec is created, need to pin it, this is done in MatSetUpMultiply_MPIAIJ. Hack */ 449 Mat A = *a_Amat_crs, P = *a_P_inout; 450 PetscMPIInt size; 451 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 452 if (size > 1) { 453 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data; 454 PetscCall(VecBindToCPU(a->lvec, PETSC_TRUE)); 455 PetscCall(VecBindToCPU(p->lvec, PETSC_TRUE)); 456 } 457 } 458 } 459 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0)); 460 } // processor reduce 461 PetscFunctionReturn(PETSC_SUCCESS); 462 } 463 464 // used in GEO 465 PetscErrorCode PCGAMGSquareGraph_GAMG(PC a_pc, Mat Gmat1, Mat *Gmat2) 466 { 467 const char *prefix; 468 char addp[32]; 469 PC_MG *mg = (PC_MG *)a_pc->data; 470 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 471 472 PetscFunctionBegin; 473 PetscCall(PCGetOptionsPrefix(a_pc, &prefix)); 474 PetscCall(PetscInfo(a_pc, "%s: Square Graph on level %" PetscInt_FMT "\n", ((PetscObject)a_pc)->prefix, pc_gamg->current_level + 1)); 475 PetscCall(MatProductCreate(Gmat1, Gmat1, NULL, Gmat2)); 476 PetscCall(MatSetOptionsPrefix(*Gmat2, prefix)); 477 PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_square_%" PetscInt_FMT "_", pc_gamg->current_level)); 478 PetscCall(MatAppendOptionsPrefix(*Gmat2, addp)); 479 if ((*Gmat2)->structurally_symmetric == PETSC_BOOL3_TRUE) { 480 PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AB)); 481 } else { 482 PetscCall(MatSetOption(Gmat1, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE)); 483 PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AtB)); 484 } 485 PetscCall(MatProductSetFromOptions(*Gmat2)); 486 PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0)); 487 PetscCall(MatProductSymbolic(*Gmat2)); 488 PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0)); 489 PetscCall(MatProductClear(*Gmat2)); 490 /* we only need the sparsity, cheat and tell PETSc the matrix has been assembled */ 491 (*Gmat2)->assembled = PETSC_TRUE; 492 PetscFunctionReturn(PETSC_SUCCESS); 493 } 494 495 /* 496 PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner 497 by setting data structures and options. 498 499 Input Parameter: 500 . pc - the preconditioner context 501 502 */ 503 static PetscErrorCode PCSetUp_GAMG(PC pc) 504 { 505 PC_MG *mg = (PC_MG *)pc->data; 506 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 507 Mat Pmat = pc->pmat; 508 PetscInt fine_level, level, level1, bs, M, N, qq, lidx, nASMBlocksArr[PETSC_MG_MAXLEVELS], cr_bs; 509 MPI_Comm comm; 510 PetscMPIInt rank, size, nactivepe; 511 Mat Aarr[PETSC_MG_MAXLEVELS], Parr[PETSC_MG_MAXLEVELS]; 512 IS *ASMLocalIDsArr[PETSC_MG_MAXLEVELS]; 513 PetscBool is_last = PETSC_FALSE; 514 #if defined(PETSC_USE_INFO) 515 PetscLogDouble nnz0 = 0., nnztot = 0.; 516 MatInfo info; 517 #endif 518 519 PetscFunctionBegin; 520 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 521 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 522 PetscCallMPI(MPI_Comm_size(comm, &size)); 523 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0)); 524 if (pc->setupcalled) { 525 if (!pc_gamg->reuse_prol || pc->flag == DIFFERENT_NONZERO_PATTERN) { 526 /* reset everything */ 527 PetscCall(PCReset_MG(pc)); 528 pc->setupcalled = 0; 529 } else { 530 PC_MG_Levels **mglevels = mg->levels; 531 /* just do Galerkin grids */ 532 Mat B, dA, dB; 533 if (pc_gamg->Nlevels > 1) { 534 PetscInt gl; 535 /* currently only handle case where mat and pmat are the same on coarser levels */ 536 PetscCall(KSPGetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, &dA, &dB)); 537 /* (re)set to get dirty flag */ 538 PetscCall(KSPSetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, dA, dB)); 539 540 for (level = pc_gamg->Nlevels - 2, gl = 0; level >= 0; level--, gl++) { 541 MatReuse reuse = MAT_INITIAL_MATRIX; 542 #if defined(GAMG_STAGES) 543 PetscCall(PetscLogStagePush(gamg_stages[gl])); 544 #endif 545 /* matrix structure can change from repartitioning or process reduction but don't know if we have process reduction here. Should fix */ 546 PetscCall(KSPGetOperators(mglevels[level]->smoothd, NULL, &B)); 547 if (B->product) { 548 if (B->product->A == dB && B->product->B == mglevels[level + 1]->interpolate) reuse = MAT_REUSE_MATRIX; 549 } 550 if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDestroy(&mglevels[level]->A)); 551 if (reuse == MAT_REUSE_MATRIX) { 552 PetscCall(PetscInfo(pc, "%s: RAP after initial setup, reuse matrix level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level)); 553 } else { 554 PetscCall(PetscInfo(pc, "%s: RAP after initial setup, with repartitioning (new matrix) level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level)); 555 } 556 PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0)); 557 PetscCall(MatPtAP(dB, mglevels[level + 1]->interpolate, reuse, PETSC_DEFAULT, &B)); 558 PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0)); 559 if (reuse == MAT_INITIAL_MATRIX) mglevels[level]->A = B; 560 PetscCall(KSPSetOperators(mglevels[level]->smoothd, B, B)); 561 // check for redoing eigen estimates 562 if (pc_gamg->recompute_esteig) { 563 PetscBool ischeb; 564 KSP smoother; 565 PetscCall(PCMGGetSmoother(pc, level + 1, &smoother)); 566 PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb)); 567 if (ischeb) { 568 KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data; 569 cheb->emin_provided = 0; 570 cheb->emax_provided = 0; 571 } 572 /* we could call PetscCall(KSPChebyshevSetEigenvalues(smoother, 0, 0)); but the logic does not work properly */ 573 } 574 // inc 575 dB = B; 576 #if defined(GAMG_STAGES) 577 PetscCall(PetscLogStagePop()); 578 #endif 579 } 580 } 581 582 PetscCall(PCSetUp_MG(pc)); 583 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0)); 584 PetscFunctionReturn(PETSC_SUCCESS); 585 } 586 } 587 588 if (!pc_gamg->data) { 589 if (pc_gamg->orig_data) { 590 PetscCall(MatGetBlockSize(Pmat, &bs)); 591 PetscCall(MatGetLocalSize(Pmat, &qq, NULL)); 592 593 pc_gamg->data_sz = (qq / bs) * pc_gamg->orig_data_cell_rows * pc_gamg->orig_data_cell_cols; 594 pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows; 595 pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols; 596 597 PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data)); 598 for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq]; 599 } else { 600 PetscCheck(pc_gamg->ops->createdefaultdata, comm, PETSC_ERR_PLIB, "'createdefaultdata' not set(?) need to support NULL data"); 601 PetscCall(pc_gamg->ops->createdefaultdata(pc, Pmat)); 602 } 603 } 604 605 /* cache original data for reuse */ 606 if (!pc_gamg->orig_data && (PetscBool)(!pc_gamg->reuse_prol)) { 607 PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->orig_data)); 608 for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq]; 609 pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows; 610 pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols; 611 } 612 613 /* get basic dims */ 614 PetscCall(MatGetBlockSize(Pmat, &bs)); 615 PetscCall(MatGetSize(Pmat, &M, NULL)); 616 617 #if defined(PETSC_USE_INFO) 618 PetscCall(MatGetInfo(Pmat, MAT_GLOBAL_SUM, &info)); /* global reduction */ 619 nnz0 = info.nz_used; 620 nnztot = info.nz_used; 621 #endif 622 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)); 623 624 /* Get A_i and R_i */ 625 for (level = 0, Aarr[0] = Pmat, nactivepe = size; level < (pc_gamg->Nlevels - 1) && (level == 0 || M > pc_gamg->coarse_eq_limit); level++) { 626 pc_gamg->current_level = level; 627 PetscCheck(level < PETSC_MG_MAXLEVELS - 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many levels %" PetscInt_FMT, level + 1); 628 level1 = level + 1; 629 #if defined(GAMG_STAGES) 630 if (!gamg_stages[level]) { 631 char str[32]; 632 PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", (int)level)); 633 PetscCall(PetscLogStageRegister(str, &gamg_stages[level])); 634 } 635 PetscCall(PetscLogStagePush(gamg_stages[level])); 636 #endif 637 /* construct prolongator - Parr[level1] */ 638 if (level == 0 && pc_gamg->injection_index_size > 0) { 639 Mat Prol; 640 MatType mtype; 641 PetscInt prol_m, prol_n, Prol_N = (M / bs) * pc_gamg->injection_index_size, Istart, Iend, nn, row; 642 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" : "")); 643 PetscCall(MatGetOwnershipRange(Pmat, &Istart, &Iend)); 644 PetscCall(MatGetLocalSize(Pmat, &prol_m, NULL)); // rows m x n 645 prol_n = (prol_m / bs) * pc_gamg->injection_index_size; 646 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); 647 PetscCall(MatGetType(Pmat, &mtype)); 648 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &Prol)); 649 PetscCall(MatSetBlockSizes(Prol, bs, pc_gamg->injection_index_size)); 650 PetscCall(MatSetSizes(Prol, prol_m, prol_n, M, Prol_N)); 651 PetscCall(MatSetType(Prol, mtype)); 652 #if PetscDefined(HAVE_DEVICE) 653 PetscBool flg; 654 PetscCall(MatBoundToCPU(Pmat, &flg)); 655 PetscCall(MatBindToCPU(Prol, flg)); 656 if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE)); 657 #endif 658 PetscCall(MatSeqAIJSetPreallocation(Prol, 1, NULL)); 659 PetscCall(MatMPIAIJSetPreallocation(Prol, 1, NULL, 0, NULL)); 660 // set I \kron [1, 1, ... ]^T 661 for (PetscInt ii = Istart, col = (Istart / bs) * pc_gamg->injection_index_size; ii < Iend; ii += bs) { 662 const PetscScalar one = 1; 663 for (PetscInt jj = 0; jj < pc_gamg->injection_index_size; jj++, col++) { 664 PetscInt row = ii + pc_gamg->injection_index[jj]; 665 PetscCall(MatSetValues(Prol, 1, &row, 1, &col, &one, INSERT_VALUES)); 666 } 667 } 668 PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY)); 669 PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY)); 670 PetscCall(MatViewFromOptions(Prol, NULL, "-mat_view_injection")); 671 PetscCall(MatGetBlockSizes(Prol, NULL, &cr_bs)); // column size 672 Parr[level1] = Prol; 673 // can not deal with null space -- with array of 'injection cols' we could take 'injection rows and 'injection cols' to 'data' 674 if (pc_gamg->data) { 675 pc_gamg->data_cell_cols = pc_gamg->injection_index_size; 676 pc_gamg->data_cell_rows = pc_gamg->injection_index_size; 677 pc_gamg->orig_data_cell_cols = 0; 678 pc_gamg->orig_data_cell_rows = 0; 679 PetscCall(PetscFree(pc_gamg->data)); 680 pc_gamg->data_sz = pc_gamg->injection_index_size * prol_n; 681 PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data)); 682 for (row = nn = 0; row < prol_n; row += pc_gamg->injection_index_size) { 683 for (int jj = 0; jj < pc_gamg->injection_index_size; jj++) { 684 int idx = row * pc_gamg->injection_index_size + jj * pc_gamg->injection_index_size; 685 for (int kk = 0; kk < pc_gamg->injection_index_size; kk++, nn++) { pc_gamg->data[idx + kk] = (jj == kk) ? 1 : 0; } 686 } 687 } 688 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); 689 } 690 } else { 691 Mat Gmat, mat; 692 PetscCoarsenData *agg_lists; 693 Mat Prol11; 694 695 PetscCall(PCGAMGCreateGraph(pc, Aarr[level], &Gmat)); 696 PetscCall(pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists)); // Gmat may have ghosts for QR aggregates not in matrix 697 PetscCall(PetscCDGetMat(agg_lists, &mat)); 698 if (!mat) PetscCall(PetscCDSetMat(agg_lists, Gmat)); 699 PetscCall(pc_gamg->ops->prolongator(pc, Aarr[level], agg_lists, &Prol11)); 700 /* could have failed to create new level */ 701 if (Prol11) { 702 const char *prefix; 703 char addp[32]; 704 705 /* get new block size of coarse matrices */ 706 PetscCall(MatGetBlockSizes(Prol11, NULL, &cr_bs)); // column size 707 708 if (pc_gamg->ops->optprolongator) { 709 /* smooth */ 710 PetscCall(pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11)); 711 } 712 713 if (pc_gamg->use_aggs_in_asm) { 714 PetscInt bs; 715 PetscCall(MatGetBlockSizes(Prol11, &bs, NULL)); // row block size 716 PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level])); 717 PetscCall(PetscInfo(pc, "%d: %" PetscInt_FMT " ASM local domains, bs = %d\n", (int)level, nASMBlocksArr[level], (int)bs)); 718 } else if (pc_gamg->asm_hem_aggs) { 719 MatCoarsen crs; 720 const char *prefix; 721 PetscInt bs; 722 PetscCall(PetscCDGetMat(agg_lists, &mat)); 723 if (mat == Gmat) PetscCall(PetscCDClearMat(agg_lists)); // take the Mat away from the list (yuck) 724 PetscCall(PetscCDDestroy(agg_lists)); 725 PetscCall(PetscInfo(pc, "HEM ASM passes = %d\n", (int)pc_gamg->asm_hem_aggs)); 726 PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &crs)); 727 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix)); 728 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)crs, prefix)); 729 PetscCall(MatCoarsenSetFromOptions(crs)); // get strength args 730 PetscCall(MatCoarsenSetType(crs, MATCOARSENHEM)); 731 PetscCall(MatCoarsenSetMaximumIterations(crs, pc_gamg->asm_hem_aggs)); 732 PetscCall(MatCoarsenSetAdjacency(crs, Gmat)); 733 PetscCall(MatCoarsenSetStrictAggs(crs, PETSC_TRUE)); 734 PetscCall(MatCoarsenApply(crs)); 735 PetscCall(MatCoarsenViewFromOptions(crs, NULL, "-agg_hem_mat_coarsen_view")); 736 PetscCall(MatCoarsenGetData(crs, &agg_lists)); /* output */ 737 PetscCall(MatCoarsenDestroy(&crs)); 738 // create aggregates 739 PetscCall(MatGetBlockSizes(Aarr[level], &bs, NULL)); // row block size 740 PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level])); 741 } 742 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 743 PetscCall(MatSetOptionsPrefix(Prol11, prefix)); 744 PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_prolongator_%d_", (int)level)); 745 PetscCall(MatAppendOptionsPrefix(Prol11, addp)); 746 /* Always generate the transpose with CUDA 747 Such behaviour can be adapted with -pc_gamg_prolongator_ prefixed options */ 748 PetscCall(MatSetOption(Prol11, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE)); 749 PetscCall(MatSetFromOptions(Prol11)); 750 Parr[level1] = Prol11; 751 } else Parr[level1] = NULL; /* failed to coarsen */ 752 PetscCall(PetscCDGetMat(agg_lists, &mat)); 753 if (mat == Gmat) PetscCall(PetscCDClearMat(agg_lists)); // take the Mat away from the list (yuck) 754 PetscCall(MatDestroy(&Gmat)); 755 PetscCall(PetscCDDestroy(agg_lists)); 756 } /* construct prolongator scope */ 757 if (level == 0) Aarr[0] = Pmat; /* use Pmat for finest level setup */ 758 if (!Parr[level1]) { /* failed to coarsen */ 759 PetscCall(PetscInfo(pc, "%s: Stop gridding, level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level)); 760 #if defined(GAMG_STAGES) 761 PetscCall(PetscLogStagePop()); 762 #endif 763 break; 764 } 765 PetscCall(MatGetSize(Parr[level1], &M, &N)); /* N is next M, a loop test variables */ 766 PetscCheck(!is_last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Is last ?"); 767 if (N <= pc_gamg->coarse_eq_limit) is_last = PETSC_TRUE; 768 if (level1 == pc_gamg->Nlevels - 1) is_last = PETSC_TRUE; 769 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0)); 770 PetscCall(pc_gamg->ops->createlevel(pc, Aarr[level], cr_bs, &Parr[level1], &Aarr[level1], &nactivepe, NULL, is_last)); 771 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0)); 772 773 PetscCall(MatGetSize(Aarr[level1], &M, &N)); /* M is loop test variables */ 774 #if defined(PETSC_USE_INFO) 775 PetscCall(MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info)); 776 nnztot += info.nz_used; 777 #endif 778 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)); 779 780 #if defined(GAMG_STAGES) 781 PetscCall(PetscLogStagePop()); 782 #endif 783 /* stop if one node or one proc -- could pull back for singular problems */ 784 if ((pc_gamg->data_cell_cols && M / pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M / bs < 2)) { 785 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)); 786 level++; 787 break; 788 } 789 } /* levels */ 790 PetscCall(PetscFree(pc_gamg->data)); 791 792 PetscCall(PetscInfo(pc, "%s: %" PetscInt_FMT " levels, operator complexity = %g\n", ((PetscObject)pc)->prefix, level + 1, nnztot / nnz0)); 793 pc_gamg->Nlevels = level + 1; 794 fine_level = level; 795 PetscCall(PCMGSetLevels(pc, pc_gamg->Nlevels, NULL)); 796 797 if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */ 798 799 /* set default smoothers & set operators */ 800 for (lidx = 1, level = pc_gamg->Nlevels - 2; lidx <= fine_level; lidx++, level--) { 801 KSP smoother; 802 PC subpc; 803 804 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 805 PetscCall(KSPGetPC(smoother, &subpc)); 806 807 PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE)); 808 /* set ops */ 809 PetscCall(KSPSetOperators(smoother, Aarr[level], Aarr[level])); 810 PetscCall(PCMGSetInterpolation(pc, lidx, Parr[level + 1])); 811 812 /* set defaults */ 813 PetscCall(KSPSetType(smoother, KSPCHEBYSHEV)); 814 815 /* set blocks for ASM smoother that uses the 'aggregates' */ 816 if (pc_gamg->use_aggs_in_asm || pc_gamg->asm_hem_aggs) { 817 PetscInt sz; 818 IS *iss; 819 820 sz = nASMBlocksArr[level]; 821 iss = ASMLocalIDsArr[level]; 822 PetscCall(PCSetType(subpc, PCASM)); 823 PetscCall(PCASMSetOverlap(subpc, 0)); 824 PetscCall(PCASMSetType(subpc, PC_ASM_BASIC)); 825 if (!sz) { 826 IS is; 827 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is)); 828 PetscCall(PCASMSetLocalSubdomains(subpc, 1, NULL, &is)); 829 PetscCall(ISDestroy(&is)); 830 } else { 831 PetscInt kk; 832 PetscCall(PCASMSetLocalSubdomains(subpc, sz, iss, NULL)); 833 for (kk = 0; kk < sz; kk++) PetscCall(ISDestroy(&iss[kk])); 834 PetscCall(PetscFree(iss)); 835 } 836 ASMLocalIDsArr[level] = NULL; 837 nASMBlocksArr[level] = 0; 838 } else { 839 PetscCall(PCSetType(subpc, PCJACOBI)); 840 } 841 } 842 { 843 /* coarse grid */ 844 KSP smoother, *k2; 845 PC subpc, pc2; 846 PetscInt ii, first; 847 Mat Lmat = Aarr[(level = pc_gamg->Nlevels - 1)]; 848 lidx = 0; 849 850 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 851 PetscCall(KSPSetOperators(smoother, Lmat, Lmat)); 852 if (!pc_gamg->use_parallel_coarse_grid_solver) { 853 PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE)); 854 PetscCall(KSPGetPC(smoother, &subpc)); 855 PetscCall(PCSetType(subpc, PCBJACOBI)); 856 PetscCall(PCSetUp(subpc)); 857 PetscCall(PCBJacobiGetSubKSP(subpc, &ii, &first, &k2)); 858 PetscCheck(ii == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " is not one", ii); 859 PetscCall(KSPGetPC(k2[0], &pc2)); 860 PetscCall(PCSetType(pc2, PCLU)); 861 PetscCall(PCFactorSetShiftType(pc2, MAT_SHIFT_INBLOCKS)); 862 PetscCall(KSPSetTolerances(k2[0], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1)); 863 PetscCall(KSPSetType(k2[0], KSPPREONLY)); 864 } 865 } 866 867 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 868 PetscObjectOptionsBegin((PetscObject)pc); 869 PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject)); 870 PetscOptionsEnd(); 871 PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL)); 872 873 /* set cheby eigen estimates from SA to use in the solver */ 874 if (pc_gamg->use_sa_esteig) { 875 for (lidx = 1, level = pc_gamg->Nlevels - 2; level >= 0; lidx++, level--) { 876 KSP smoother; 877 PetscBool ischeb; 878 879 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 880 PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb)); 881 if (ischeb) { 882 KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data; 883 884 // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG 885 if (mg->max_eigen_DinvA[level] > 0) { 886 // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it. 887 // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.) 888 PetscReal emax, emin; 889 890 emin = mg->min_eigen_DinvA[level]; 891 emax = mg->max_eigen_DinvA[level]; 892 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)); 893 cheb->emin_provided = emin; 894 cheb->emax_provided = emax; 895 } 896 } 897 } 898 } 899 900 PetscCall(PCSetUp_MG(pc)); 901 902 /* clean up */ 903 for (level = 1; level < pc_gamg->Nlevels; level++) { 904 PetscCall(MatDestroy(&Parr[level])); 905 PetscCall(MatDestroy(&Aarr[level])); 906 } 907 } else { 908 KSP smoother; 909 910 PetscCall(PetscInfo(pc, "%s: One level solver used (system is seen as DD). Using default solver.\n", ((PetscObject)pc)->prefix)); 911 PetscCall(PCMGGetSmoother(pc, 0, &smoother)); 912 PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0])); 913 PetscCall(KSPSetType(smoother, KSPPREONLY)); 914 PetscCall(PCSetUp_MG(pc)); 915 } 916 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0)); 917 PetscFunctionReturn(PETSC_SUCCESS); 918 } 919 920 /* 921 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 922 that was created with PCCreate_GAMG(). 923 924 Input Parameter: 925 . pc - the preconditioner context 926 927 Application Interface Routine: PCDestroy() 928 */ 929 PetscErrorCode PCDestroy_GAMG(PC pc) 930 { 931 PC_MG *mg = (PC_MG *)pc->data; 932 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 933 934 PetscFunctionBegin; 935 PetscCall(PCReset_GAMG(pc)); 936 if (pc_gamg->ops->destroy) PetscCall((*pc_gamg->ops->destroy)(pc)); 937 PetscCall(PetscFree(pc_gamg->ops)); 938 PetscCall(PetscFree(pc_gamg->gamg_type_name)); 939 PetscCall(PetscFree(pc_gamg)); 940 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", NULL)); 941 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", NULL)); 942 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", NULL)); 943 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", NULL)); 944 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", NULL)); 945 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", NULL)); 946 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", NULL)); 947 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", NULL)); 948 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", NULL)); 949 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", NULL)); 950 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", NULL)); 951 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", NULL)); 952 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", NULL)); 953 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", NULL)); 954 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", NULL)); 955 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", NULL)); 956 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", NULL)); 957 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", NULL)); 958 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndices_C", NULL)); 959 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", NULL)); 960 PetscCall(PCDestroy_MG(pc)); 961 PetscFunctionReturn(PETSC_SUCCESS); 962 } 963 964 /*@ 965 PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction in `PCGAMG` 966 967 Logically Collective 968 969 Input Parameters: 970 + pc - the preconditioner context 971 - n - the number of equations 972 973 Options Database Key: 974 . -pc_gamg_process_eq_limit <limit> - set the limit 975 976 Level: intermediate 977 978 Note: 979 `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process 980 that has degrees of freedom 981 982 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()` 983 @*/ 984 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 985 { 986 PetscFunctionBegin; 987 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 988 PetscTryMethod(pc, "PCGAMGSetProcEqLim_C", (PC, PetscInt), (pc, n)); 989 PetscFunctionReturn(PETSC_SUCCESS); 990 } 991 992 static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 993 { 994 PC_MG *mg = (PC_MG *)pc->data; 995 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 996 997 PetscFunctionBegin; 998 if (n > 0) pc_gamg->min_eq_proc = n; 999 PetscFunctionReturn(PETSC_SUCCESS); 1000 } 1001 1002 /*@ 1003 PCGAMGSetCoarseEqLim - Set maximum number of equations on the coarsest grid of `PCGAMG` 1004 1005 Collective 1006 1007 Input Parameters: 1008 + pc - the preconditioner context 1009 - n - maximum number of equations to aim for 1010 1011 Options Database Key: 1012 . -pc_gamg_coarse_eq_limit <limit> - set the limit 1013 1014 Level: intermediate 1015 1016 Note: 1017 For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid 1018 has less than 1000 unknowns. 1019 1020 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()` 1021 @*/ 1022 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n) 1023 { 1024 PetscFunctionBegin; 1025 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1026 PetscTryMethod(pc, "PCGAMGSetCoarseEqLim_C", (PC, PetscInt), (pc, n)); 1027 PetscFunctionReturn(PETSC_SUCCESS); 1028 } 1029 1030 static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n) 1031 { 1032 PC_MG *mg = (PC_MG *)pc->data; 1033 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1034 1035 PetscFunctionBegin; 1036 if (n > 0) pc_gamg->coarse_eq_limit = n; 1037 PetscFunctionReturn(PETSC_SUCCESS); 1038 } 1039 1040 /*@ 1041 PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids when reducing the number of MPI ranks to use 1042 1043 Collective 1044 1045 Input Parameters: 1046 + pc - the preconditioner context 1047 - n - `PETSC_TRUE` or `PETSC_FALSE` 1048 1049 Options Database Key: 1050 . -pc_gamg_repartition <true,false> - turn on the repartitioning 1051 1052 Level: intermediate 1053 1054 Note: 1055 This will generally improve the loading balancing of the work on each level 1056 1057 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()` 1058 @*/ 1059 PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n) 1060 { 1061 PetscFunctionBegin; 1062 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1063 PetscTryMethod(pc, "PCGAMGSetRepartition_C", (PC, PetscBool), (pc, n)); 1064 PetscFunctionReturn(PETSC_SUCCESS); 1065 } 1066 1067 static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n) 1068 { 1069 PC_MG *mg = (PC_MG *)pc->data; 1070 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1071 1072 PetscFunctionBegin; 1073 pc_gamg->repart = n; 1074 PetscFunctionReturn(PETSC_SUCCESS); 1075 } 1076 1077 /*@ 1078 PCGAMGSetUseSAEstEig - Use the eigen estimate from smoothed aggregation for the Chebyshev smoother during the solution process 1079 1080 Collective 1081 1082 Input Parameters: 1083 + pc - the preconditioner context 1084 - b - flag 1085 1086 Options Database Key: 1087 . -pc_gamg_use_sa_esteig <true,false> - use the eigen estimate 1088 1089 Level: advanced 1090 1091 Notes: 1092 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$. 1093 Eigenvalue estimates (based on a few `PCCG` or `PCGMRES` iterations) are computed to choose $\omega$ so that this is a stable smoothing operation. 1094 If Chebyshev with Jacobi (diagonal) preconditioning is used for smoothing, then the eigenvalue estimates can be reused during the solution process 1095 This option is only used when the smoother uses Jacobi, and should be turned off if a different `PCJacobiType` is used. 1096 It became default in PETSc 3.17. 1097 1098 .seealso: [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `PCGAMGSetRecomputeEstEig()` 1099 @*/ 1100 PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool b) 1101 { 1102 PetscFunctionBegin; 1103 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1104 PetscTryMethod(pc, "PCGAMGSetUseSAEstEig_C", (PC, PetscBool), (pc, b)); 1105 PetscFunctionReturn(PETSC_SUCCESS); 1106 } 1107 1108 static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool b) 1109 { 1110 PC_MG *mg = (PC_MG *)pc->data; 1111 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1112 1113 PetscFunctionBegin; 1114 pc_gamg->use_sa_esteig = b; 1115 PetscFunctionReturn(PETSC_SUCCESS); 1116 } 1117 1118 /*@ 1119 PCGAMGSetRecomputeEstEig - Set flag for Chebyshev smoothers to recompute the eigen estimates 1120 1121 Collective 1122 1123 Input Parameters: 1124 + pc - the preconditioner context 1125 - b - flag 1126 1127 Options Database Key: 1128 . -pc_gamg_recompute_esteig <true> - use the eigen estimate 1129 1130 Level: advanced 1131 1132 .seealso: [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()` 1133 @*/ 1134 PetscErrorCode PCGAMGSetRecomputeEstEig(PC pc, PetscBool b) 1135 { 1136 PetscFunctionBegin; 1137 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1138 PetscTryMethod(pc, "PCGAMGSetRecomputeEstEig_C", (PC, PetscBool), (pc, b)); 1139 PetscFunctionReturn(PETSC_SUCCESS); 1140 } 1141 1142 static PetscErrorCode PCGAMGSetRecomputeEstEig_GAMG(PC pc, PetscBool b) 1143 { 1144 PC_MG *mg = (PC_MG *)pc->data; 1145 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1146 1147 PetscFunctionBegin; 1148 pc_gamg->recompute_esteig = b; 1149 PetscFunctionReturn(PETSC_SUCCESS); 1150 } 1151 1152 /*@ 1153 PCGAMGSetEigenvalues - Set WHAT eigenvalues WHY? 1154 1155 Collective 1156 1157 Input Parameters: 1158 + pc - the preconditioner context 1159 . emax - max eigenvalue 1160 - emin - min eigenvalue 1161 1162 Options Database Key: 1163 . -pc_gamg_eigenvalues <emin,emax> - estimates of the eigenvalues 1164 1165 Level: intermediate 1166 1167 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetUseSAEstEig()` 1168 @*/ 1169 PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax, PetscReal emin) 1170 { 1171 PetscFunctionBegin; 1172 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1173 PetscTryMethod(pc, "PCGAMGSetEigenvalues_C", (PC, PetscReal, PetscReal), (pc, emax, emin)); 1174 PetscFunctionReturn(PETSC_SUCCESS); 1175 } 1176 1177 static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc, PetscReal emax, PetscReal emin) 1178 { 1179 PC_MG *mg = (PC_MG *)pc->data; 1180 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1181 1182 PetscFunctionBegin; 1183 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); 1184 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); 1185 pc_gamg->emax = emax; 1186 pc_gamg->emin = emin; 1187 PetscFunctionReturn(PETSC_SUCCESS); 1188 } 1189 1190 /*@ 1191 PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding a `PCGAMG` algebraic multigrid preconditioner 1192 1193 Collective 1194 1195 Input Parameters: 1196 + pc - the preconditioner context 1197 - n - `PETSC_TRUE` or `PETSC_FALSE` 1198 1199 Options Database Key: 1200 . -pc_gamg_reuse_interpolation <true,false> - reuse the previous interpolation 1201 1202 Level: intermediate 1203 1204 Note: 1205 May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows 1206 rebuilding the preconditioner quicker. 1207 1208 .seealso: [](ch_ksp), `PCGAMG` 1209 @*/ 1210 PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n) 1211 { 1212 PetscFunctionBegin; 1213 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1214 PetscTryMethod(pc, "PCGAMGSetReuseInterpolation_C", (PC, PetscBool), (pc, n)); 1215 PetscFunctionReturn(PETSC_SUCCESS); 1216 } 1217 1218 static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n) 1219 { 1220 PC_MG *mg = (PC_MG *)pc->data; 1221 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1222 1223 PetscFunctionBegin; 1224 pc_gamg->reuse_prol = n; 1225 PetscFunctionReturn(PETSC_SUCCESS); 1226 } 1227 1228 /*@ 1229 PCGAMGASMSetUseAggs - Have the `PCGAMG` smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner 1230 used as the smoother 1231 1232 Collective 1233 1234 Input Parameters: 1235 + pc - the preconditioner context 1236 - flg - `PETSC_TRUE` to use aggregates, `PETSC_FALSE` to not 1237 1238 Options Database Key: 1239 . -pc_gamg_asm_use_agg <true,false> - use aggregates to define the additive Schwarz subdomains 1240 1241 Level: intermediate 1242 1243 .seealso: [](ch_ksp), `PCGAMG`, `PCASM`, `PCSetType` 1244 @*/ 1245 PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg) 1246 { 1247 PetscFunctionBegin; 1248 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1249 PetscTryMethod(pc, "PCGAMGASMSetUseAggs_C", (PC, PetscBool), (pc, flg)); 1250 PetscFunctionReturn(PETSC_SUCCESS); 1251 } 1252 1253 static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg) 1254 { 1255 PC_MG *mg = (PC_MG *)pc->data; 1256 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1257 1258 PetscFunctionBegin; 1259 pc_gamg->use_aggs_in_asm = flg; 1260 PetscFunctionReturn(PETSC_SUCCESS); 1261 } 1262 1263 /*@ 1264 PCGAMGSetParallelCoarseGridSolve - allow a parallel coarse grid solver 1265 1266 Collective 1267 1268 Input Parameters: 1269 + pc - the preconditioner context 1270 - flg - `PETSC_TRUE` to not force coarse grid onto one processor 1271 1272 Options Database Key: 1273 . -pc_gamg_parallel_coarse_grid_solver - use a parallel coarse grid direct solver 1274 1275 Level: intermediate 1276 1277 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetCpuPinCoarseGrids()` 1278 @*/ 1279 PetscErrorCode PCGAMGSetParallelCoarseGridSolve(PC pc, PetscBool flg) 1280 { 1281 PetscFunctionBegin; 1282 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1283 PetscTryMethod(pc, "PCGAMGSetParallelCoarseGridSolve_C", (PC, PetscBool), (pc, flg)); 1284 PetscFunctionReturn(PETSC_SUCCESS); 1285 } 1286 1287 static PetscErrorCode PCGAMGSetParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg) 1288 { 1289 PC_MG *mg = (PC_MG *)pc->data; 1290 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1291 1292 PetscFunctionBegin; 1293 pc_gamg->use_parallel_coarse_grid_solver = flg; 1294 PetscFunctionReturn(PETSC_SUCCESS); 1295 } 1296 1297 /*@ 1298 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 1299 1300 Collective 1301 1302 Input Parameters: 1303 + pc - the preconditioner context 1304 - flg - `PETSC_TRUE` to pin coarse grids to the CPU 1305 1306 Options Database Key: 1307 . -pc_gamg_cpu_pin_coarse_grids - pin the coarse grids to the CPU 1308 1309 Level: intermediate 1310 1311 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetParallelCoarseGridSolve()` 1312 @*/ 1313 PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg) 1314 { 1315 PetscFunctionBegin; 1316 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1317 PetscTryMethod(pc, "PCGAMGSetCpuPinCoarseGrids_C", (PC, PetscBool), (pc, flg)); 1318 PetscFunctionReturn(PETSC_SUCCESS); 1319 } 1320 1321 static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg) 1322 { 1323 PC_MG *mg = (PC_MG *)pc->data; 1324 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1325 1326 PetscFunctionBegin; 1327 pc_gamg->cpu_pin_coarse_grids = flg; 1328 PetscFunctionReturn(PETSC_SUCCESS); 1329 } 1330 1331 /*@ 1332 PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type) 1333 1334 Collective 1335 1336 Input Parameters: 1337 + pc - the preconditioner context 1338 - flg - `PCGAMGLayoutType` type, either `PCGAMG_LAYOUT_COMPACT` or `PCGAMG_LAYOUT_SPREAD` 1339 1340 Options Database Key: 1341 . -pc_gamg_coarse_grid_layout_type - place the coarse grids with natural ordering 1342 1343 Level: intermediate 1344 1345 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGLayoutType`, `PCGAMG_LAYOUT_COMPACT`, `PCGAMG_LAYOUT_SPREAD` 1346 @*/ 1347 PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg) 1348 { 1349 PetscFunctionBegin; 1350 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1351 PetscTryMethod(pc, "PCGAMGSetCoarseGridLayoutType_C", (PC, PCGAMGLayoutType), (pc, flg)); 1352 PetscFunctionReturn(PETSC_SUCCESS); 1353 } 1354 1355 static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg) 1356 { 1357 PC_MG *mg = (PC_MG *)pc->data; 1358 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1359 1360 PetscFunctionBegin; 1361 pc_gamg->layout_type = flg; 1362 PetscFunctionReturn(PETSC_SUCCESS); 1363 } 1364 1365 /*@ 1366 PCGAMGSetNlevels - Sets the maximum number of levels `PCGAMG` will use 1367 1368 Collective 1369 1370 Input Parameters: 1371 + pc - the preconditioner 1372 - n - the maximum number of levels to use 1373 1374 Options Database Key: 1375 . -pc_mg_levels <n> - set the maximum number of levels to allow 1376 1377 Level: intermediate 1378 1379 Developer Notes: 1380 Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG` 1381 1382 .seealso: [](ch_ksp), `PCGAMG` 1383 @*/ 1384 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n) 1385 { 1386 PetscFunctionBegin; 1387 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1388 PetscTryMethod(pc, "PCGAMGSetNlevels_C", (PC, PetscInt), (pc, n)); 1389 PetscFunctionReturn(PETSC_SUCCESS); 1390 } 1391 1392 static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n) 1393 { 1394 PC_MG *mg = (PC_MG *)pc->data; 1395 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1396 1397 PetscFunctionBegin; 1398 pc_gamg->Nlevels = n; 1399 PetscFunctionReturn(PETSC_SUCCESS); 1400 } 1401 1402 /*@ 1403 PCGAMGASMSetHEM - Sets the number of HEM matching passed 1404 1405 Collective 1406 1407 Input Parameters: 1408 + pc - the preconditioner 1409 - n - number of HEM matching passed to construct ASM subdomains 1410 1411 Options Database Key: 1412 . -pc_gamg_asm_hem <n> - set the number of HEM matching passed 1413 1414 Level: intermediate 1415 1416 Developer Notes: 1417 Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG` 1418 1419 .seealso: [](ch_ksp), `PCGAMG` 1420 @*/ 1421 PetscErrorCode PCGAMGASMSetHEM(PC pc, PetscInt n) 1422 { 1423 PetscFunctionBegin; 1424 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1425 PetscTryMethod(pc, "PCGAMGASMSetHEM_C", (PC, PetscInt), (pc, n)); 1426 PetscFunctionReturn(PETSC_SUCCESS); 1427 } 1428 1429 static PetscErrorCode PCGAMGASMSetHEM_GAMG(PC pc, PetscInt n) 1430 { 1431 PC_MG *mg = (PC_MG *)pc->data; 1432 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1433 1434 PetscFunctionBegin; 1435 pc_gamg->asm_hem_aggs = n; 1436 PetscFunctionReturn(PETSC_SUCCESS); 1437 } 1438 1439 /*@ 1440 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 1441 1442 Not Collective 1443 1444 Input Parameters: 1445 + pc - the preconditioner context 1446 . 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 1447 - n - number of threshold values provided in array 1448 1449 Options Database Key: 1450 . -pc_gamg_threshold <threshold> - the threshold to drop edges 1451 1452 Level: intermediate 1453 1454 Notes: 1455 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. 1456 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. 1457 1458 If `n` is less than the total number of coarsenings (see `PCGAMGSetNlevels()`), then threshold scaling (see `PCGAMGSetThresholdScale()`) is used for each successive coarsening. 1459 In this case, `PCGAMGSetThresholdScale()` must be called before `PCGAMGSetThreshold()`. 1460 If `n` is greater than the total number of levels, the excess entries in threshold will not be used. 1461 1462 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetMinDegreeOrderingMISk()`, `PCGAMGSetThresholdScale()` 1463 @*/ 1464 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n) 1465 { 1466 PetscFunctionBegin; 1467 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1468 if (n) PetscAssertPointer(v, 2); 1469 PetscTryMethod(pc, "PCGAMGSetThreshold_C", (PC, PetscReal[], PetscInt), (pc, v, n)); 1470 PetscFunctionReturn(PETSC_SUCCESS); 1471 } 1472 1473 static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n) 1474 { 1475 PC_MG *mg = (PC_MG *)pc->data; 1476 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1477 PetscInt i; 1478 1479 PetscFunctionBegin; 1480 for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i]; 1481 for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale; 1482 PetscFunctionReturn(PETSC_SUCCESS); 1483 } 1484 1485 /*@ 1486 PCGAMGSetRankReductionFactors - Set a manual schedule for MPI rank reduction on coarse grids 1487 1488 Collective 1489 1490 Input Parameters: 1491 + pc - the preconditioner context 1492 . v - array of reduction factors. 0 for first value forces a reduction to one process/device on first level in CUDA 1493 - n - number of values provided in array 1494 1495 Options Database Key: 1496 . -pc_gamg_rank_reduction_factors <factors> - provide the schedule 1497 1498 Level: intermediate 1499 1500 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()` 1501 @*/ 1502 PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n) 1503 { 1504 PetscFunctionBegin; 1505 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1506 if (n) PetscAssertPointer(v, 2); 1507 PetscTryMethod(pc, "PCGAMGSetRankReductionFactors_C", (PC, PetscInt[], PetscInt), (pc, v, n)); 1508 PetscFunctionReturn(PETSC_SUCCESS); 1509 } 1510 1511 static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n) 1512 { 1513 PC_MG *mg = (PC_MG *)pc->data; 1514 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1515 PetscInt i; 1516 1517 PetscFunctionBegin; 1518 for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i]; 1519 for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */ 1520 PetscFunctionReturn(PETSC_SUCCESS); 1521 } 1522 1523 /*@ 1524 PCGAMGSetThresholdScale - Relative threshold reduction at each level 1525 1526 Not Collective 1527 1528 Input Parameters: 1529 + pc - the preconditioner context 1530 - v - the threshold value reduction, usually < 1.0 1531 1532 Options Database Key: 1533 . -pc_gamg_threshold_scale <v> - set the relative threshold reduction on each level 1534 1535 Level: advanced 1536 1537 Note: 1538 The initial threshold (for an arbitrary number of levels starting from the finest) can be set with `PCGAMGSetThreshold()`. 1539 This scaling is used for each subsequent coarsening, but must be called before `PCGAMGSetThreshold()`. 1540 1541 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()` 1542 @*/ 1543 PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v) 1544 { 1545 PetscFunctionBegin; 1546 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1547 PetscTryMethod(pc, "PCGAMGSetThresholdScale_C", (PC, PetscReal), (pc, v)); 1548 PetscFunctionReturn(PETSC_SUCCESS); 1549 } 1550 1551 static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v) 1552 { 1553 PC_MG *mg = (PC_MG *)pc->data; 1554 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1555 1556 PetscFunctionBegin; 1557 pc_gamg->threshold_scale = v; 1558 PetscFunctionReturn(PETSC_SUCCESS); 1559 } 1560 1561 /*@C 1562 PCGAMGSetType - Set the type of algorithm `PCGAMG` should use 1563 1564 Collective 1565 1566 Input Parameters: 1567 + pc - the preconditioner context 1568 - type - `PCGAMGAGG`, `PCGAMGGEO`, or `PCGAMGCLASSICAL` 1569 1570 Options Database Key: 1571 . -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply - only agg supported 1572 1573 Level: intermediate 1574 1575 .seealso: [](ch_ksp), `PCGAMGGetType()`, `PCGAMG`, `PCGAMGType` 1576 @*/ 1577 PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type) 1578 { 1579 PetscFunctionBegin; 1580 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1581 PetscTryMethod(pc, "PCGAMGSetType_C", (PC, PCGAMGType), (pc, type)); 1582 PetscFunctionReturn(PETSC_SUCCESS); 1583 } 1584 1585 /*@C 1586 PCGAMGGetType - Get the type of algorithm `PCGAMG` will use 1587 1588 Collective 1589 1590 Input Parameter: 1591 . pc - the preconditioner context 1592 1593 Output Parameter: 1594 . type - the type of algorithm used 1595 1596 Level: intermediate 1597 1598 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetType()`, `PCGAMGType` 1599 @*/ 1600 PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type) 1601 { 1602 PetscFunctionBegin; 1603 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1604 PetscUseMethod(pc, "PCGAMGGetType_C", (PC, PCGAMGType *), (pc, type)); 1605 PetscFunctionReturn(PETSC_SUCCESS); 1606 } 1607 1608 static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type) 1609 { 1610 PC_MG *mg = (PC_MG *)pc->data; 1611 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1612 1613 PetscFunctionBegin; 1614 *type = pc_gamg->type; 1615 PetscFunctionReturn(PETSC_SUCCESS); 1616 } 1617 1618 static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type) 1619 { 1620 PC_MG *mg = (PC_MG *)pc->data; 1621 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1622 PetscErrorCode (*r)(PC); 1623 1624 PetscFunctionBegin; 1625 pc_gamg->type = type; 1626 PetscCall(PetscFunctionListFind(GAMGList, type, &r)); 1627 PetscCheck(r, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown GAMG type %s given", type); 1628 if (pc_gamg->ops->destroy) { 1629 PetscCall((*pc_gamg->ops->destroy)(pc)); 1630 PetscCall(PetscMemzero(pc_gamg->ops, sizeof(struct _PCGAMGOps))); 1631 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1632 /* cleaning up common data in pc_gamg - this should disappear someday */ 1633 pc_gamg->data_cell_cols = 0; 1634 pc_gamg->data_cell_rows = 0; 1635 pc_gamg->orig_data_cell_cols = 0; 1636 pc_gamg->orig_data_cell_rows = 0; 1637 PetscCall(PetscFree(pc_gamg->data)); 1638 pc_gamg->data_sz = 0; 1639 } 1640 PetscCall(PetscFree(pc_gamg->gamg_type_name)); 1641 PetscCall(PetscStrallocpy(type, &pc_gamg->gamg_type_name)); 1642 PetscCall((*r)(pc)); 1643 PetscFunctionReturn(PETSC_SUCCESS); 1644 } 1645 1646 static PetscErrorCode PCView_GAMG(PC pc, PetscViewer viewer) 1647 { 1648 PC_MG *mg = (PC_MG *)pc->data; 1649 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1650 PetscReal gc = 0, oc = 0; 1651 1652 PetscFunctionBegin; 1653 PetscCall(PetscViewerASCIIPrintf(viewer, " GAMG specific options\n")); 1654 PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for dropping small values in graph on each level =")); 1655 for (PetscInt i = 0; i < mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)pc_gamg->threshold[i])); 1656 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 1657 PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold scaling factor for each level not specified = %g\n", (double)pc_gamg->threshold_scale)); 1658 if (pc_gamg->use_aggs_in_asm) PetscCall(PetscViewerASCIIPrintf(viewer, " Using aggregates from coarsening process to define subdomains for PCASM\n")); // this take presedence 1659 else if (pc_gamg->asm_hem_aggs) 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)); 1660 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")); 1661 if (pc_gamg->injection_index_size) { 1662 PetscCall(PetscViewerASCIIPrintf(viewer, " Using injection restriction/prolongation on first level, dofs:")); 1663 for (int i = 0; i < pc_gamg->injection_index_size; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %d", (int)pc_gamg->injection_index[i])); 1664 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 1665 } 1666 if (pc_gamg->ops->view) PetscCall((*pc_gamg->ops->view)(pc, viewer)); 1667 PetscCall(PCMGGetGridComplexity(pc, &gc, &oc)); 1668 PetscCall(PetscViewerASCIIPrintf(viewer, " Complexity: grid = %g operator = %g\n", (double)gc, (double)oc)); 1669 PetscFunctionReturn(PETSC_SUCCESS); 1670 } 1671 1672 /*@ 1673 PCGAMGSetInjectionIndex - Array of subset of variables per vertex to inject into coarse grid space 1674 1675 Logically Collective 1676 1677 Input Parameters: 1678 + pc - the coarsen context 1679 . n - number of indices 1680 - idx - array of indices 1681 1682 Options Database Key: 1683 . -pc_gamg_injection_index - array of subset of variables per vertex to use for injection coarse grid space 1684 1685 Level: intermediate 1686 1687 .seealso: `PCGAMG` 1688 @*/ 1689 PetscErrorCode PCGAMGSetInjectionIndex(PC pc, PetscInt n, PetscInt idx[]) 1690 { 1691 PetscFunctionBegin; 1692 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1693 PetscValidLogicalCollectiveInt(pc, n, 2); 1694 PetscTryMethod(pc, "PCGAMGSetInjectionIndex_C", (PC, PetscInt, PetscInt[]), (pc, n, idx)); 1695 PetscFunctionReturn(PETSC_SUCCESS); 1696 } 1697 1698 static PetscErrorCode PCGAMGSetInjectionIndex_GAMG(PC pc, PetscInt n, PetscInt idx[]) 1699 { 1700 PC_MG *mg = (PC_MG *)pc->data; 1701 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1702 1703 PetscFunctionBegin; 1704 pc_gamg->injection_index_size = n; 1705 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); 1706 for (PetscInt i = 0; i < n; i++) pc_gamg->injection_index[i] = idx[i]; 1707 PetscFunctionReturn(PETSC_SUCCESS); 1708 } 1709 1710 static PetscErrorCode PCSetFromOptions_GAMG(PC pc, PetscOptionItems *PetscOptionsObject) 1711 { 1712 PC_MG *mg = (PC_MG *)pc->data; 1713 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1714 PetscBool flag; 1715 MPI_Comm comm; 1716 char prefix[256], tname[32]; 1717 PetscInt i, n; 1718 const char *pcpre; 1719 static const char *LayoutTypes[] = {"compact", "spread", "PCGAMGLayoutType", "PC_GAMG_LAYOUT", NULL}; 1720 1721 PetscFunctionBegin; 1722 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 1723 PetscOptionsHeadBegin(PetscOptionsObject, "GAMG options"); 1724 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)); 1725 if (flag) PetscCall(PCGAMGSetType(pc, tname)); 1726 PetscCall(PetscOptionsBool("-pc_gamg_repartition", "Repartion coarse grids", "PCGAMGSetRepartition", pc_gamg->repart, &pc_gamg->repart, NULL)); 1727 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)); 1728 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)); 1729 PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation", "Reuse prolongation operator", "PCGAMGReuseInterpolation", pc_gamg->reuse_prol, &pc_gamg->reuse_prol, NULL)); 1730 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)); 1731 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)); 1732 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)); 1733 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, 1734 (PetscEnum)pc_gamg->layout_type, (PetscEnum *)&pc_gamg->layout_type, NULL)); 1735 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)); 1736 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)); 1737 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)); 1738 PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale", "Scaling of threshold for each level not specified", "PCGAMGSetThresholdScale", pc_gamg->threshold_scale, &pc_gamg->threshold_scale, NULL)); 1739 n = PETSC_MG_MAXLEVELS; 1740 PetscCall(PetscOptionsRealArray("-pc_gamg_threshold", "Relative threshold to use for dropping edges in aggregation graph", "PCGAMGSetThreshold", pc_gamg->threshold, &n, &flag)); 1741 if (!flag || n < PETSC_MG_MAXLEVELS) { 1742 if (!flag) n = 1; 1743 i = n; 1744 do { 1745 pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale; 1746 } while (++i < PETSC_MG_MAXLEVELS); 1747 } 1748 PetscCall(PetscOptionsInt("-pc_mg_levels", "Set number of MG levels (should get from base class)", "PCGAMGSetNlevels", pc_gamg->Nlevels, &pc_gamg->Nlevels, NULL)); 1749 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); 1750 n = PETSC_MG_MAXLEVELS; 1751 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)); 1752 if (!flag) i = 0; 1753 else i = n; 1754 do { 1755 pc_gamg->level_reduction_factors[i] = -1; 1756 } while (++i < PETSC_MG_MAXLEVELS); 1757 { 1758 PetscReal eminmax[2] = {0., 0.}; 1759 n = 2; 1760 PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues", "extreme eigenvalues for smoothed aggregation", "PCGAMGSetEigenvalues", eminmax, &n, &flag)); 1761 if (flag) { 1762 PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues"); 1763 PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0])); 1764 } 1765 } 1766 pc_gamg->injection_index_size = MAT_COARSEN_STRENGTH_INDEX_SIZE; 1767 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)); 1768 /* set options for subtype */ 1769 PetscCall((*pc_gamg->ops->setfromoptions)(pc, PetscOptionsObject)); 1770 1771 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1772 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : "")); 1773 PetscOptionsHeadEnd(); 1774 PetscFunctionReturn(PETSC_SUCCESS); 1775 } 1776 1777 /*MC 1778 PCGAMG - Geometric algebraic multigrid (AMG) preconditioner 1779 1780 Options Database Keys: 1781 + -pc_gamg_type <type,default=agg> - one of agg, geo, or classical (only smoothed aggregation, agg, supported) 1782 . -pc_gamg_repartition <bool,default=false> - repartition the degrees of freedom across the coarse grids as they are determined 1783 . -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 1784 . -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> 1785 equations on each process that has degrees of freedom 1786 . -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for. 1787 . -pc_gamg_reuse_interpolation <bool,default=true> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations (should always be true) 1788 . -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) 1789 - -pc_gamg_threshold_scale <scale,default=1> - Scaling of threshold on each coarser grid if not specified 1790 1791 Options Database Keys for Aggregation: 1792 + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation 1793 . -pc_gamg_aggressive_coarsening <n,default=1> - number of aggressive coarsening (MIS-2) levels from finest. 1794 . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening 1795 . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm 1796 . -pc_gamg_pc_gamg_asm_hem_aggs <n,default=0> - Number of HEM aggregation steps for ASM smoother 1797 - -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive') 1798 1799 Options Database Keys for Multigrid: 1800 + -pc_mg_cycle_type <v> - v or w, see `PCMGSetCycleType()` 1801 . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see PCMGSetDistinctSmoothUp() 1802 . -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade 1803 - -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 1804 1805 Level: intermediate 1806 1807 Notes: 1808 To obtain good performance for `PCGAMG` for vector valued problems you must 1809 call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point 1810 call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator 1811 1812 The many options for `PCMG` also work directly for `PCGAMG` such as controlling the smoothers on each level etc. 1813 1814 See [the Users Manual section on PCGAMG](sec_amg) and [the Users Manual section on PCMG](sec_mg)for more details. 1815 1816 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, 1817 `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()` 1818 M*/ 1819 PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc) 1820 { 1821 PC_GAMG *pc_gamg; 1822 PC_MG *mg; 1823 1824 PetscFunctionBegin; 1825 /* register AMG type */ 1826 PetscCall(PCGAMGInitializePackage()); 1827 1828 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1829 PetscCall(PCSetType(pc, PCMG)); 1830 PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG)); 1831 1832 /* create a supporting struct and attach it to pc */ 1833 PetscCall(PetscNew(&pc_gamg)); 1834 PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL)); 1835 mg = (PC_MG *)pc->data; 1836 mg->innerctx = pc_gamg; 1837 1838 PetscCall(PetscNew(&pc_gamg->ops)); 1839 1840 /* these should be in subctx but repartitioning needs simple arrays */ 1841 pc_gamg->data_sz = 0; 1842 pc_gamg->data = NULL; 1843 1844 /* overwrite the pointers of PCMG by the functions of base class PCGAMG */ 1845 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1846 pc->ops->setup = PCSetUp_GAMG; 1847 pc->ops->reset = PCReset_GAMG; 1848 pc->ops->destroy = PCDestroy_GAMG; 1849 mg->view = PCView_GAMG; 1850 1851 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", PCGAMGSetProcEqLim_GAMG)); 1852 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", PCGAMGSetCoarseEqLim_GAMG)); 1853 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", PCGAMGSetRepartition_GAMG)); 1854 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", PCGAMGSetEigenvalues_GAMG)); 1855 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", PCGAMGSetUseSAEstEig_GAMG)); 1856 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", PCGAMGSetRecomputeEstEig_GAMG)); 1857 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", PCGAMGSetReuseInterpolation_GAMG)); 1858 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", PCGAMGASMSetUseAggs_GAMG)); 1859 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", PCGAMGSetParallelCoarseGridSolve_GAMG)); 1860 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", PCGAMGSetCpuPinCoarseGrids_GAMG)); 1861 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", PCGAMGSetCoarseGridLayoutType_GAMG)); 1862 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", PCGAMGSetThreshold_GAMG)); 1863 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", PCGAMGSetRankReductionFactors_GAMG)); 1864 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", PCGAMGSetThresholdScale_GAMG)); 1865 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", PCGAMGSetType_GAMG)); 1866 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", PCGAMGGetType_GAMG)); 1867 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", PCGAMGSetNlevels_GAMG)); 1868 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", PCGAMGASMSetHEM_GAMG)); 1869 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", PCGAMGSetInjectionIndex_GAMG)); 1870 pc_gamg->repart = PETSC_FALSE; 1871 pc_gamg->reuse_prol = PETSC_TRUE; 1872 pc_gamg->use_aggs_in_asm = PETSC_FALSE; 1873 pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE; 1874 pc_gamg->cpu_pin_coarse_grids = PETSC_FALSE; 1875 pc_gamg->layout_type = PCGAMG_LAYOUT_SPREAD; 1876 pc_gamg->min_eq_proc = 50; 1877 pc_gamg->asm_hem_aggs = 0; 1878 pc_gamg->coarse_eq_limit = 50; 1879 for (int i = 0; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = -1; 1880 pc_gamg->threshold_scale = 1.; 1881 pc_gamg->Nlevels = PETSC_MG_MAXLEVELS; 1882 pc_gamg->current_level = 0; /* don't need to init really */ 1883 pc_gamg->use_sa_esteig = PETSC_TRUE; 1884 pc_gamg->recompute_esteig = PETSC_TRUE; 1885 pc_gamg->emin = 0; 1886 pc_gamg->emax = 0; 1887 1888 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1889 1890 /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */ 1891 PetscCall(PCGAMGSetType(pc, PCGAMGAGG)); 1892 PetscFunctionReturn(PETSC_SUCCESS); 1893 } 1894 1895 /*@C 1896 PCGAMGInitializePackage - This function initializes everything in the `PCGAMG` package. It is called 1897 from `PCInitializePackage()`. 1898 1899 Level: developer 1900 1901 .seealso: [](ch_ksp), `PetscInitialize()` 1902 @*/ 1903 PetscErrorCode PCGAMGInitializePackage(void) 1904 { 1905 PetscInt l; 1906 1907 PetscFunctionBegin; 1908 if (PCGAMGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 1909 PCGAMGPackageInitialized = PETSC_TRUE; 1910 PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGGEO, PCCreateGAMG_GEO)); 1911 PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGAGG, PCCreateGAMG_AGG)); 1912 PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGCLASSICAL, PCCreateGAMG_Classical)); 1913 PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage)); 1914 1915 /* general events */ 1916 PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP])); 1917 PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH])); 1918 PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN])); 1919 PetscCall(PetscLogEventRegister(" GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS])); 1920 PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL])); 1921 PetscCall(PetscLogEventRegister(" GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA])); 1922 PetscCall(PetscLogEventRegister(" GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB])); 1923 PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT])); 1924 PetscCall(PetscLogEventRegister(" GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM])); 1925 PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL])); 1926 PetscCall(PetscLogEventRegister(" GAMG PtAP", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP])); 1927 PetscCall(PetscLogEventRegister(" GAMG Reduce", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE])); 1928 PetscCall(PetscLogEventRegister(" GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART])); 1929 /* PetscCall(PetscLogEventRegister(" GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */ 1930 /* PetscCall(PetscLogEventRegister(" GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */ 1931 /* PetscCall(PetscLogEventRegister(" GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */ 1932 for (l = 0; l < PETSC_MG_MAXLEVELS; l++) { 1933 char ename[32]; 1934 1935 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Squ l%02" PetscInt_FMT, l)); 1936 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0])); 1937 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Gal l%02" PetscInt_FMT, l)); 1938 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1])); 1939 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Opt l%02" PetscInt_FMT, l)); 1940 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2])); 1941 } 1942 #if defined(GAMG_STAGES) 1943 { /* create timer stages */ 1944 char str[32]; 1945 PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", 0)); 1946 PetscCall(PetscLogStageRegister(str, &gamg_stages[0])); 1947 } 1948 #endif 1949 PetscFunctionReturn(PETSC_SUCCESS); 1950 } 1951 1952 /*@C 1953 PCGAMGFinalizePackage - This function frees everything from the `PCGAMG` package. It is 1954 called from `PetscFinalize()` automatically. 1955 1956 Level: developer 1957 1958 .seealso: [](ch_ksp), `PetscFinalize()` 1959 @*/ 1960 PetscErrorCode PCGAMGFinalizePackage(void) 1961 { 1962 PetscFunctionBegin; 1963 PCGAMGPackageInitialized = PETSC_FALSE; 1964 PetscCall(PetscFunctionListDestroy(&GAMGList)); 1965 PetscFunctionReturn(PETSC_SUCCESS); 1966 } 1967 1968 /*@C 1969 PCGAMGRegister - Register a `PCGAMG` implementation. 1970 1971 Input Parameters: 1972 + type - string that will be used as the name of the `PCGAMG` type. 1973 - create - function for creating the gamg context. 1974 1975 Level: developer 1976 1977 .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()` 1978 @*/ 1979 PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC)) 1980 { 1981 PetscFunctionBegin; 1982 PetscCall(PCGAMGInitializePackage()); 1983 PetscCall(PetscFunctionListAdd(&GAMGList, type, create)); 1984 PetscFunctionReturn(PETSC_SUCCESS); 1985 } 1986 1987 /*@C 1988 PCGAMGCreateGraph - Creates a graph that is used by the ``PCGAMGType`` in the coarsening process 1989 1990 Input Parameters: 1991 + pc - the `PCGAMG` 1992 - A - the matrix, for any level 1993 1994 Output Parameter: 1995 . G - the graph 1996 1997 Level: advanced 1998 1999 .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()` 2000 @*/ 2001 PetscErrorCode PCGAMGCreateGraph(PC pc, Mat A, Mat *G) 2002 { 2003 PC_MG *mg = (PC_MG *)pc->data; 2004 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 2005 2006 PetscFunctionBegin; 2007 PetscCall(pc_gamg->ops->creategraph(pc, A, G)); 2008 PetscFunctionReturn(PETSC_SUCCESS); 2009 } 2010