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