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]; 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, &N)); 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 ", 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)); 623 624 /* Get A_i and R_i */ 625 for (level = 0, Aarr[0] = Pmat, nactivepe = size; level < (pc_gamg->Nlevels - 1) && (!level || 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 */ 638 Mat Gmat; 639 PetscCoarsenData *agg_lists; 640 Mat Prol11; 641 642 PetscCall(PCGAMGCreateGraph(pc, Aarr[level], &Gmat)); 643 PetscCall(pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists)); // Gmat may have ghosts for QR aggregates not in matrix 644 PetscCall(pc_gamg->ops->prolongator(pc, Aarr[level], Gmat, agg_lists, &Prol11)); 645 646 /* could have failed to create new level */ 647 if (Prol11) { 648 const char *prefix; 649 char addp[32]; 650 651 /* get new block size of coarse matrices */ 652 PetscCall(MatGetBlockSizes(Prol11, NULL, &bs)); 653 654 if (pc_gamg->ops->optprolongator) { 655 /* smooth */ 656 PetscCall(pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11)); 657 } 658 659 if (pc_gamg->use_aggs_in_asm) { 660 PetscInt bs, Istart; 661 PetscCall(MatGetBlockSizes(Prol11, &bs, NULL)); // not timed directly, ugly, could remove, but good ASM method 662 PetscCall(MatGetOwnershipRange(Prol11, &Istart, NULL)); 663 PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level])); 664 PetscCall(PetscInfo(pc, "%d: %" PetscInt_FMT " ASM local domains, bs = %d\n", (int)level, nASMBlocksArr[level], (int)bs)); 665 } 666 667 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 668 PetscCall(MatSetOptionsPrefix(Prol11, prefix)); 669 PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_prolongator_%d_", (int)level)); 670 PetscCall(MatAppendOptionsPrefix(Prol11, addp)); 671 /* Always generate the transpose with CUDA 672 Such behaviour can be adapted with -pc_gamg_prolongator_ prefixed options */ 673 PetscCall(MatSetOption(Prol11, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE)); 674 PetscCall(MatSetFromOptions(Prol11)); 675 Parr[level1] = Prol11; 676 } else Parr[level1] = NULL; /* failed to coarsen */ 677 678 PetscCall(MatDestroy(&Gmat)); 679 PetscCall(PetscCDDestroy(agg_lists)); 680 } /* construct prolongator scope */ 681 if (!level) Aarr[0] = Pmat; /* use Pmat for finest level setup */ 682 if (!Parr[level1]) { /* failed to coarsen */ 683 PetscCall(PetscInfo(pc, "%s: Stop gridding, level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level)); 684 #if defined(GAMG_STAGES) 685 PetscCall(PetscLogStagePop()); 686 #endif 687 break; 688 } 689 PetscCall(MatGetSize(Parr[level1], &M, &N)); /* N is next M, a loop test variables */ 690 PetscCheck(!is_last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Is last ?"); 691 if (N <= pc_gamg->coarse_eq_limit) is_last = PETSC_TRUE; 692 if (level1 == pc_gamg->Nlevels - 1) is_last = PETSC_TRUE; 693 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0)); 694 PetscCall(pc_gamg->ops->createlevel(pc, Aarr[level], bs, &Parr[level1], &Aarr[level1], &nactivepe, NULL, is_last)); 695 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0)); 696 697 PetscCall(MatGetSize(Aarr[level1], &M, &N)); /* M is loop test variables */ 698 #if defined(PETSC_USE_INFO) 699 PetscCall(MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info)); 700 nnztot += info.nz_used; 701 #endif 702 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)); 703 704 #if defined(GAMG_STAGES) 705 PetscCall(PetscLogStagePop()); 706 #endif 707 /* stop if one node or one proc -- could pull back for singular problems */ 708 if ((pc_gamg->data_cell_cols && M / pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M / bs < 2)) { 709 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)); 710 level++; 711 break; 712 } 713 } /* levels */ 714 PetscCall(PetscFree(pc_gamg->data)); 715 716 PetscCall(PetscInfo(pc, "%s: %" PetscInt_FMT " levels, operator complexity = %g\n", ((PetscObject)pc)->prefix, level + 1, nnztot / nnz0)); 717 pc_gamg->Nlevels = level + 1; 718 fine_level = level; 719 PetscCall(PCMGSetLevels(pc, pc_gamg->Nlevels, NULL)); 720 721 if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */ 722 723 /* set default smoothers & set operators */ 724 for (lidx = 1, level = pc_gamg->Nlevels - 2; lidx <= fine_level; lidx++, level--) { 725 KSP smoother; 726 PC subpc; 727 728 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 729 PetscCall(KSPGetPC(smoother, &subpc)); 730 731 PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE)); 732 /* set ops */ 733 PetscCall(KSPSetOperators(smoother, Aarr[level], Aarr[level])); 734 PetscCall(PCMGSetInterpolation(pc, lidx, Parr[level + 1])); 735 736 /* set defaults */ 737 PetscCall(KSPSetType(smoother, KSPCHEBYSHEV)); 738 739 /* set blocks for ASM smoother that uses the 'aggregates' */ 740 if (pc_gamg->use_aggs_in_asm) { 741 PetscInt sz; 742 IS *iss; 743 744 sz = nASMBlocksArr[level]; 745 iss = ASMLocalIDsArr[level]; 746 PetscCall(PCSetType(subpc, PCASM)); 747 PetscCall(PCASMSetOverlap(subpc, 0)); 748 PetscCall(PCASMSetType(subpc, PC_ASM_BASIC)); 749 if (!sz) { 750 IS is; 751 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is)); 752 PetscCall(PCASMSetLocalSubdomains(subpc, 1, NULL, &is)); 753 PetscCall(ISDestroy(&is)); 754 } else { 755 PetscInt kk; 756 PetscCall(PCASMSetLocalSubdomains(subpc, sz, iss, NULL)); 757 for (kk = 0; kk < sz; kk++) PetscCall(ISDestroy(&iss[kk])); 758 PetscCall(PetscFree(iss)); 759 } 760 ASMLocalIDsArr[level] = NULL; 761 nASMBlocksArr[level] = 0; 762 } else { 763 PetscCall(PCSetType(subpc, PCJACOBI)); 764 } 765 } 766 { 767 /* coarse grid */ 768 KSP smoother, *k2; 769 PC subpc, pc2; 770 PetscInt ii, first; 771 Mat Lmat = Aarr[(level = pc_gamg->Nlevels - 1)]; 772 lidx = 0; 773 774 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 775 PetscCall(KSPSetOperators(smoother, Lmat, Lmat)); 776 if (!pc_gamg->use_parallel_coarse_grid_solver) { 777 PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE)); 778 PetscCall(KSPGetPC(smoother, &subpc)); 779 PetscCall(PCSetType(subpc, PCBJACOBI)); 780 PetscCall(PCSetUp(subpc)); 781 PetscCall(PCBJacobiGetSubKSP(subpc, &ii, &first, &k2)); 782 PetscCheck(ii == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " is not one", ii); 783 PetscCall(KSPGetPC(k2[0], &pc2)); 784 PetscCall(PCSetType(pc2, PCLU)); 785 PetscCall(PCFactorSetShiftType(pc2, MAT_SHIFT_INBLOCKS)); 786 PetscCall(KSPSetTolerances(k2[0], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1)); 787 PetscCall(KSPSetType(k2[0], KSPPREONLY)); 788 } 789 } 790 791 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 792 PetscObjectOptionsBegin((PetscObject)pc); 793 PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject)); 794 PetscOptionsEnd(); 795 PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL)); 796 797 /* set cheby eigen estimates from SA to use in the solver */ 798 if (pc_gamg->use_sa_esteig) { 799 for (lidx = 1, level = pc_gamg->Nlevels - 2; level >= 0; lidx++, level--) { 800 KSP smoother; 801 PetscBool ischeb; 802 803 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 804 PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb)); 805 if (ischeb) { 806 KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data; 807 808 // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG 809 if (mg->max_eigen_DinvA[level] > 0) { 810 // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it. 811 // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.) 812 PetscReal emax, emin; 813 814 emin = mg->min_eigen_DinvA[level]; 815 emax = mg->max_eigen_DinvA[level]; 816 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)); 817 cheb->emin_provided = emin; 818 cheb->emax_provided = emax; 819 } 820 } 821 } 822 } 823 824 PetscCall(PCSetUp_MG(pc)); 825 826 /* clean up */ 827 for (level = 1; level < pc_gamg->Nlevels; level++) { 828 PetscCall(MatDestroy(&Parr[level])); 829 PetscCall(MatDestroy(&Aarr[level])); 830 } 831 } else { 832 KSP smoother; 833 834 PetscCall(PetscInfo(pc, "%s: One level solver used (system is seen as DD). Using default solver.\n", ((PetscObject)pc)->prefix)); 835 PetscCall(PCMGGetSmoother(pc, 0, &smoother)); 836 PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0])); 837 PetscCall(KSPSetType(smoother, KSPPREONLY)); 838 PetscCall(PCSetUp_MG(pc)); 839 } 840 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0)); 841 PetscFunctionReturn(PETSC_SUCCESS); 842 } 843 844 /* 845 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 846 that was created with PCCreate_GAMG(). 847 848 Input Parameter: 849 . pc - the preconditioner context 850 851 Application Interface Routine: PCDestroy() 852 */ 853 PetscErrorCode PCDestroy_GAMG(PC pc) 854 { 855 PC_MG *mg = (PC_MG *)pc->data; 856 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 857 858 PetscFunctionBegin; 859 PetscCall(PCReset_GAMG(pc)); 860 if (pc_gamg->ops->destroy) PetscCall((*pc_gamg->ops->destroy)(pc)); 861 PetscCall(PetscFree(pc_gamg->ops)); 862 PetscCall(PetscFree(pc_gamg->gamg_type_name)); 863 PetscCall(PetscFree(pc_gamg)); 864 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", NULL)); 865 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", NULL)); 866 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", NULL)); 867 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", NULL)); 868 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", NULL)); 869 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", NULL)); 870 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", NULL)); 871 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", NULL)); 872 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", NULL)); 873 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", NULL)); 874 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", NULL)); 875 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", NULL)); 876 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", NULL)); 877 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", NULL)); 878 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", NULL)); 879 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", NULL)); 880 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", NULL)); 881 PetscCall(PCDestroy_MG(pc)); 882 PetscFunctionReturn(PETSC_SUCCESS); 883 } 884 885 /*@ 886 PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction in `PCGAMG` 887 888 Logically Collective 889 890 Input Parameters: 891 + pc - the preconditioner context 892 - n - the number of equations 893 894 Options Database Key: 895 . -pc_gamg_process_eq_limit <limit> - set the limit 896 897 Level: intermediate 898 899 Note: 900 `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process 901 that has degrees of freedom 902 903 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()` 904 @*/ 905 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 906 { 907 PetscFunctionBegin; 908 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 909 PetscTryMethod(pc, "PCGAMGSetProcEqLim_C", (PC, PetscInt), (pc, n)); 910 PetscFunctionReturn(PETSC_SUCCESS); 911 } 912 913 static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 914 { 915 PC_MG *mg = (PC_MG *)pc->data; 916 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 917 918 PetscFunctionBegin; 919 if (n > 0) pc_gamg->min_eq_proc = n; 920 PetscFunctionReturn(PETSC_SUCCESS); 921 } 922 923 /*@ 924 PCGAMGSetCoarseEqLim - Set maximum number of equations on the coarsest grid of `PCGAMG` 925 926 Collective 927 928 Input Parameters: 929 + pc - the preconditioner context 930 - n - maximum number of equations to aim for 931 932 Options Database Key: 933 . -pc_gamg_coarse_eq_limit <limit> - set the limit 934 935 Level: intermediate 936 937 Note: 938 For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid 939 has less than 1000 unknowns. 940 941 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()` 942 @*/ 943 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n) 944 { 945 PetscFunctionBegin; 946 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 947 PetscTryMethod(pc, "PCGAMGSetCoarseEqLim_C", (PC, PetscInt), (pc, n)); 948 PetscFunctionReturn(PETSC_SUCCESS); 949 } 950 951 static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n) 952 { 953 PC_MG *mg = (PC_MG *)pc->data; 954 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 955 956 PetscFunctionBegin; 957 if (n > 0) pc_gamg->coarse_eq_limit = n; 958 PetscFunctionReturn(PETSC_SUCCESS); 959 } 960 961 /*@ 962 PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids when reducing the number of MPI ranks to use 963 964 Collective 965 966 Input Parameters: 967 + pc - the preconditioner context 968 - n - `PETSC_TRUE` or `PETSC_FALSE` 969 970 Options Database Key: 971 . -pc_gamg_repartition <true,false> - turn on the repartitioning 972 973 Level: intermediate 974 975 Note: 976 This will generally improve the loading balancing of the work on each level 977 978 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()` 979 @*/ 980 PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n) 981 { 982 PetscFunctionBegin; 983 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 984 PetscTryMethod(pc, "PCGAMGSetRepartition_C", (PC, PetscBool), (pc, n)); 985 PetscFunctionReturn(PETSC_SUCCESS); 986 } 987 988 static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n) 989 { 990 PC_MG *mg = (PC_MG *)pc->data; 991 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 992 993 PetscFunctionBegin; 994 pc_gamg->repart = n; 995 PetscFunctionReturn(PETSC_SUCCESS); 996 } 997 998 /*@ 999 PCGAMGSetUseSAEstEig - Use the eigen estimate from smoothed aggregation for the Chebyshev smoother during the solution process 1000 1001 Collective 1002 1003 Input Parameters: 1004 + pc - the preconditioner context 1005 - b - flag 1006 1007 Options Database Key: 1008 . -pc_gamg_use_sa_esteig <true,false> - use the eigen estimate 1009 1010 Level: advanced 1011 1012 Notes: 1013 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$. 1014 Eigenvalue estimates (based on a few `PCCG` or `PCGMRES` iterations) are computed to choose $\omega$ so that this is a stable smoothing operation. 1015 If Chebyshev with Jacobi (diagonal) preconditioning is used for smoothing, then the eigenvalue estimates can be reused during the solution process 1016 This option is only used when the smoother uses Jacobi, and should be turned off if a different `PCJacobiType` is used. 1017 It became default in PETSc 3.17. 1018 1019 .seealso: [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `PCGAMGSetRecomputeEstEig()` 1020 @*/ 1021 PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool b) 1022 { 1023 PetscFunctionBegin; 1024 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1025 PetscTryMethod(pc, "PCGAMGSetUseSAEstEig_C", (PC, PetscBool), (pc, b)); 1026 PetscFunctionReturn(PETSC_SUCCESS); 1027 } 1028 1029 static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool b) 1030 { 1031 PC_MG *mg = (PC_MG *)pc->data; 1032 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1033 1034 PetscFunctionBegin; 1035 pc_gamg->use_sa_esteig = b; 1036 PetscFunctionReturn(PETSC_SUCCESS); 1037 } 1038 1039 /*@ 1040 PCGAMGSetRecomputeEstEig - Set flag for Chebyshev smoothers to recompute the eigen estimates 1041 1042 Collective 1043 1044 Input Parameters: 1045 + pc - the preconditioner context 1046 - b - flag 1047 1048 Options Database Key: 1049 . -pc_gamg_recompute_esteig <true> - use the eigen estimate 1050 1051 Level: advanced 1052 1053 .seealso: [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()` 1054 @*/ 1055 PetscErrorCode PCGAMGSetRecomputeEstEig(PC pc, PetscBool b) 1056 { 1057 PetscFunctionBegin; 1058 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1059 PetscTryMethod(pc, "PCGAMGSetRecomputeEstEig_C", (PC, PetscBool), (pc, b)); 1060 PetscFunctionReturn(PETSC_SUCCESS); 1061 } 1062 1063 static PetscErrorCode PCGAMGSetRecomputeEstEig_GAMG(PC pc, PetscBool b) 1064 { 1065 PC_MG *mg = (PC_MG *)pc->data; 1066 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1067 1068 PetscFunctionBegin; 1069 pc_gamg->recompute_esteig = b; 1070 PetscFunctionReturn(PETSC_SUCCESS); 1071 } 1072 1073 /*@ 1074 PCGAMGSetEigenvalues - Set WHAT eigenvalues WHY? 1075 1076 Collective 1077 1078 Input Parameters: 1079 + pc - the preconditioner context 1080 . emax - max eigenvalue 1081 - emin - min eigenvalue 1082 1083 Options Database Key: 1084 . -pc_gamg_eigenvalues <emin,emax> - estimates of the eigenvalues 1085 1086 Level: intermediate 1087 1088 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetUseSAEstEig()` 1089 @*/ 1090 PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax, PetscReal emin) 1091 { 1092 PetscFunctionBegin; 1093 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1094 PetscTryMethod(pc, "PCGAMGSetEigenvalues_C", (PC, PetscReal, PetscReal), (pc, emax, emin)); 1095 PetscFunctionReturn(PETSC_SUCCESS); 1096 } 1097 1098 static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc, PetscReal emax, PetscReal emin) 1099 { 1100 PC_MG *mg = (PC_MG *)pc->data; 1101 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1102 1103 PetscFunctionBegin; 1104 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); 1105 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); 1106 pc_gamg->emax = emax; 1107 pc_gamg->emin = emin; 1108 PetscFunctionReturn(PETSC_SUCCESS); 1109 } 1110 1111 /*@ 1112 PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding a `PCGAMG` algebraic multigrid preconditioner 1113 1114 Collective 1115 1116 Input Parameters: 1117 + pc - the preconditioner context 1118 - n - `PETSC_TRUE` or `PETSC_FALSE` 1119 1120 Options Database Key: 1121 . -pc_gamg_reuse_interpolation <true,false> - reuse the previous interpolation 1122 1123 Level: intermediate 1124 1125 Note: 1126 May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows 1127 rebuilding the preconditioner quicker. 1128 1129 .seealso: [](ch_ksp), `PCGAMG` 1130 @*/ 1131 PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n) 1132 { 1133 PetscFunctionBegin; 1134 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1135 PetscTryMethod(pc, "PCGAMGSetReuseInterpolation_C", (PC, PetscBool), (pc, n)); 1136 PetscFunctionReturn(PETSC_SUCCESS); 1137 } 1138 1139 static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n) 1140 { 1141 PC_MG *mg = (PC_MG *)pc->data; 1142 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1143 1144 PetscFunctionBegin; 1145 pc_gamg->reuse_prol = n; 1146 PetscFunctionReturn(PETSC_SUCCESS); 1147 } 1148 1149 /*@ 1150 PCGAMGASMSetUseAggs - Have the `PCGAMG` smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner 1151 used as the smoother 1152 1153 Collective 1154 1155 Input Parameters: 1156 + pc - the preconditioner context 1157 - flg - `PETSC_TRUE` to use aggregates, `PETSC_FALSE` to not 1158 1159 Options Database Key: 1160 . -pc_gamg_asm_use_agg <true,false> - use aggregates to define the additive Schwarz subdomains 1161 1162 Level: intermediate 1163 1164 .seealso: [](ch_ksp), `PCGAMG`, `PCASM`, `PCSetType` 1165 @*/ 1166 PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg) 1167 { 1168 PetscFunctionBegin; 1169 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1170 PetscTryMethod(pc, "PCGAMGASMSetUseAggs_C", (PC, PetscBool), (pc, flg)); 1171 PetscFunctionReturn(PETSC_SUCCESS); 1172 } 1173 1174 static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg) 1175 { 1176 PC_MG *mg = (PC_MG *)pc->data; 1177 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1178 1179 PetscFunctionBegin; 1180 pc_gamg->use_aggs_in_asm = flg; 1181 PetscFunctionReturn(PETSC_SUCCESS); 1182 } 1183 1184 /*@ 1185 PCGAMGSetParallelCoarseGridSolve - allow a parallel coarse grid solver 1186 1187 Collective 1188 1189 Input Parameters: 1190 + pc - the preconditioner context 1191 - flg - `PETSC_TRUE` to not force coarse grid onto one processor 1192 1193 Options Database Key: 1194 . -pc_gamg_parallel_coarse_grid_solver - use a parallel coarse grid direct solver 1195 1196 Level: intermediate 1197 1198 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetCpuPinCoarseGrids()` 1199 @*/ 1200 PetscErrorCode PCGAMGSetParallelCoarseGridSolve(PC pc, PetscBool flg) 1201 { 1202 PetscFunctionBegin; 1203 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1204 PetscTryMethod(pc, "PCGAMGSetParallelCoarseGridSolve_C", (PC, PetscBool), (pc, flg)); 1205 PetscFunctionReturn(PETSC_SUCCESS); 1206 } 1207 1208 static PetscErrorCode PCGAMGSetParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg) 1209 { 1210 PC_MG *mg = (PC_MG *)pc->data; 1211 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1212 1213 PetscFunctionBegin; 1214 pc_gamg->use_parallel_coarse_grid_solver = flg; 1215 PetscFunctionReturn(PETSC_SUCCESS); 1216 } 1217 1218 /*@ 1219 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 1220 1221 Collective 1222 1223 Input Parameters: 1224 + pc - the preconditioner context 1225 - flg - `PETSC_TRUE` to pin coarse grids to the CPU 1226 1227 Options Database Key: 1228 . -pc_gamg_cpu_pin_coarse_grids - pin the coarse grids to the CPU 1229 1230 Level: intermediate 1231 1232 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetParallelCoarseGridSolve()` 1233 @*/ 1234 PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg) 1235 { 1236 PetscFunctionBegin; 1237 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1238 PetscTryMethod(pc, "PCGAMGSetCpuPinCoarseGrids_C", (PC, PetscBool), (pc, flg)); 1239 PetscFunctionReturn(PETSC_SUCCESS); 1240 } 1241 1242 static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg) 1243 { 1244 PC_MG *mg = (PC_MG *)pc->data; 1245 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1246 1247 PetscFunctionBegin; 1248 pc_gamg->cpu_pin_coarse_grids = flg; 1249 PetscFunctionReturn(PETSC_SUCCESS); 1250 } 1251 1252 /*@ 1253 PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type) 1254 1255 Collective 1256 1257 Input Parameters: 1258 + pc - the preconditioner context 1259 - flg - `PCGAMGLayoutType` type, either `PCGAMG_LAYOUT_COMPACT` or `PCGAMG_LAYOUT_SPREAD` 1260 1261 Options Database Key: 1262 . -pc_gamg_coarse_grid_layout_type - place the coarse grids with natural ordering 1263 1264 Level: intermediate 1265 1266 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGLayoutType`, `PCGAMG_LAYOUT_COMPACT`, `PCGAMG_LAYOUT_SPREAD` 1267 @*/ 1268 PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg) 1269 { 1270 PetscFunctionBegin; 1271 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1272 PetscTryMethod(pc, "PCGAMGSetCoarseGridLayoutType_C", (PC, PCGAMGLayoutType), (pc, flg)); 1273 PetscFunctionReturn(PETSC_SUCCESS); 1274 } 1275 1276 static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg) 1277 { 1278 PC_MG *mg = (PC_MG *)pc->data; 1279 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1280 1281 PetscFunctionBegin; 1282 pc_gamg->layout_type = flg; 1283 PetscFunctionReturn(PETSC_SUCCESS); 1284 } 1285 1286 /*@ 1287 PCGAMGSetNlevels - Sets the maximum number of levels `PCGAMG` will use 1288 1289 Not Collective 1290 1291 Input Parameters: 1292 + pc - the preconditioner 1293 - n - the maximum number of levels to use 1294 1295 Options Database Key: 1296 . -pc_mg_levels <n> - set the maximum number of levels to allow 1297 1298 Level: intermediate 1299 1300 Developer Notes: 1301 Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG` 1302 1303 .seealso: [](ch_ksp), `PCGAMG` 1304 @*/ 1305 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n) 1306 { 1307 PetscFunctionBegin; 1308 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1309 PetscTryMethod(pc, "PCGAMGSetNlevels_C", (PC, PetscInt), (pc, n)); 1310 PetscFunctionReturn(PETSC_SUCCESS); 1311 } 1312 1313 static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n) 1314 { 1315 PC_MG *mg = (PC_MG *)pc->data; 1316 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1317 1318 PetscFunctionBegin; 1319 pc_gamg->Nlevels = n; 1320 PetscFunctionReturn(PETSC_SUCCESS); 1321 } 1322 1323 /*@ 1324 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 1325 1326 Not Collective 1327 1328 Input Parameters: 1329 + pc - the preconditioner context 1330 . 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 1331 - n - number of threshold values provided in array 1332 1333 Options Database Key: 1334 . -pc_gamg_threshold <threshold> - the threshold to drop edges 1335 1336 Level: intermediate 1337 1338 Notes: 1339 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. 1340 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. 1341 1342 If `n` is less than the total number of coarsenings (see `PCGAMGSetNlevels()`), then threshold scaling (see `PCGAMGSetThresholdScale()`) is used for each successive coarsening. 1343 In this case, `PCGAMGSetThresholdScale()` must be called before `PCGAMGSetThreshold()`. 1344 If `n` is greater than the total number of levels, the excess entries in threshold will not be used. 1345 1346 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetMinDegreeOrderingMISk()`, `PCGAMGSetThresholdScale()` 1347 @*/ 1348 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n) 1349 { 1350 PetscFunctionBegin; 1351 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1352 if (n) PetscAssertPointer(v, 2); 1353 PetscTryMethod(pc, "PCGAMGSetThreshold_C", (PC, PetscReal[], PetscInt), (pc, v, n)); 1354 PetscFunctionReturn(PETSC_SUCCESS); 1355 } 1356 1357 static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n) 1358 { 1359 PC_MG *mg = (PC_MG *)pc->data; 1360 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1361 PetscInt i; 1362 PetscFunctionBegin; 1363 for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i]; 1364 for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale; 1365 PetscFunctionReturn(PETSC_SUCCESS); 1366 } 1367 1368 /*@ 1369 PCGAMGSetRankReductionFactors - Set a manual schedule for MPI rank reduction on coarse grids 1370 1371 Collective 1372 1373 Input Parameters: 1374 + pc - the preconditioner context 1375 . v - array of reduction factors. 0 for first value forces a reduction to one process/device on first level in CUDA 1376 - n - number of values provided in array 1377 1378 Options Database Key: 1379 . -pc_gamg_rank_reduction_factors <factors> - provide the schedule 1380 1381 Level: intermediate 1382 1383 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()` 1384 @*/ 1385 PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n) 1386 { 1387 PetscFunctionBegin; 1388 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1389 if (n) PetscAssertPointer(v, 2); 1390 PetscTryMethod(pc, "PCGAMGSetRankReductionFactors_C", (PC, PetscInt[], PetscInt), (pc, v, n)); 1391 PetscFunctionReturn(PETSC_SUCCESS); 1392 } 1393 1394 static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n) 1395 { 1396 PC_MG *mg = (PC_MG *)pc->data; 1397 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1398 PetscInt i; 1399 PetscFunctionBegin; 1400 for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i]; 1401 for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */ 1402 PetscFunctionReturn(PETSC_SUCCESS); 1403 } 1404 1405 /*@ 1406 PCGAMGSetThresholdScale - Relative threshold reduction at each level 1407 1408 Not Collective 1409 1410 Input Parameters: 1411 + pc - the preconditioner context 1412 - v - the threshold value reduction, usually < 1.0 1413 1414 Options Database Key: 1415 . -pc_gamg_threshold_scale <v> - set the relative threshold reduction on each level 1416 1417 Level: advanced 1418 1419 Note: 1420 The initial threshold (for an arbitrary number of levels starting from the finest) can be set with `PCGAMGSetThreshold()`. 1421 This scaling is used for each subsequent coarsening, but must be called before `PCGAMGSetThreshold()`. 1422 1423 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()` 1424 @*/ 1425 PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v) 1426 { 1427 PetscFunctionBegin; 1428 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1429 PetscTryMethod(pc, "PCGAMGSetThresholdScale_C", (PC, PetscReal), (pc, v)); 1430 PetscFunctionReturn(PETSC_SUCCESS); 1431 } 1432 1433 static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v) 1434 { 1435 PC_MG *mg = (PC_MG *)pc->data; 1436 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1437 PetscFunctionBegin; 1438 pc_gamg->threshold_scale = v; 1439 PetscFunctionReturn(PETSC_SUCCESS); 1440 } 1441 1442 /*@C 1443 PCGAMGSetType - Set the type of algorithm `PCGAMG` should use 1444 1445 Collective 1446 1447 Input Parameters: 1448 + pc - the preconditioner context 1449 - type - `PCGAMGAGG`, `PCGAMGGEO`, or `PCGAMGCLASSICAL` 1450 1451 Options Database Key: 1452 . -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply - only agg supported 1453 1454 Level: intermediate 1455 1456 .seealso: [](ch_ksp), `PCGAMGGetType()`, `PCGAMG`, `PCGAMGType` 1457 @*/ 1458 PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type) 1459 { 1460 PetscFunctionBegin; 1461 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1462 PetscTryMethod(pc, "PCGAMGSetType_C", (PC, PCGAMGType), (pc, type)); 1463 PetscFunctionReturn(PETSC_SUCCESS); 1464 } 1465 1466 /*@C 1467 PCGAMGGetType - Get the type of algorithm `PCGAMG` will use 1468 1469 Collective 1470 1471 Input Parameter: 1472 . pc - the preconditioner context 1473 1474 Output Parameter: 1475 . type - the type of algorithm used 1476 1477 Level: intermediate 1478 1479 .seealso: [](ch_ksp), `PCGAMG`, `PCGAMGSetType()`, `PCGAMGType` 1480 @*/ 1481 PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type) 1482 { 1483 PetscFunctionBegin; 1484 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1485 PetscUseMethod(pc, "PCGAMGGetType_C", (PC, PCGAMGType *), (pc, type)); 1486 PetscFunctionReturn(PETSC_SUCCESS); 1487 } 1488 1489 static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type) 1490 { 1491 PC_MG *mg = (PC_MG *)pc->data; 1492 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1493 1494 PetscFunctionBegin; 1495 *type = pc_gamg->type; 1496 PetscFunctionReturn(PETSC_SUCCESS); 1497 } 1498 1499 static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type) 1500 { 1501 PC_MG *mg = (PC_MG *)pc->data; 1502 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1503 PetscErrorCode (*r)(PC); 1504 1505 PetscFunctionBegin; 1506 pc_gamg->type = type; 1507 PetscCall(PetscFunctionListFind(GAMGList, type, &r)); 1508 PetscCheck(r, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown GAMG type %s given", type); 1509 if (pc_gamg->ops->destroy) { 1510 PetscCall((*pc_gamg->ops->destroy)(pc)); 1511 PetscCall(PetscMemzero(pc_gamg->ops, sizeof(struct _PCGAMGOps))); 1512 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1513 /* cleaning up common data in pc_gamg - this should disappear someday */ 1514 pc_gamg->data_cell_cols = 0; 1515 pc_gamg->data_cell_rows = 0; 1516 pc_gamg->orig_data_cell_cols = 0; 1517 pc_gamg->orig_data_cell_rows = 0; 1518 PetscCall(PetscFree(pc_gamg->data)); 1519 pc_gamg->data_sz = 0; 1520 } 1521 PetscCall(PetscFree(pc_gamg->gamg_type_name)); 1522 PetscCall(PetscStrallocpy(type, &pc_gamg->gamg_type_name)); 1523 PetscCall((*r)(pc)); 1524 PetscFunctionReturn(PETSC_SUCCESS); 1525 } 1526 1527 static PetscErrorCode PCView_GAMG(PC pc, PetscViewer viewer) 1528 { 1529 PC_MG *mg = (PC_MG *)pc->data; 1530 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1531 PetscReal gc = 0, oc = 0; 1532 1533 PetscFunctionBegin; 1534 PetscCall(PetscViewerASCIIPrintf(viewer, " GAMG specific options\n")); 1535 PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for dropping small values in graph on each level =")); 1536 for (PetscInt i = 0; i < mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)pc_gamg->threshold[i])); 1537 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 1538 PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold scaling factor for each level not specified = %g\n", (double)pc_gamg->threshold_scale)); 1539 if (pc_gamg->use_aggs_in_asm) PetscCall(PetscViewerASCIIPrintf(viewer, " Using aggregates from coarsening process to define subdomains for PCASM\n")); 1540 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")); 1541 if (pc_gamg->ops->view) PetscCall((*pc_gamg->ops->view)(pc, viewer)); 1542 PetscCall(PCMGGetGridComplexity(pc, &gc, &oc)); 1543 PetscCall(PetscViewerASCIIPrintf(viewer, " Complexity: grid = %g operator = %g\n", (double)gc, (double)oc)); 1544 PetscFunctionReturn(PETSC_SUCCESS); 1545 } 1546 1547 static PetscErrorCode PCSetFromOptions_GAMG(PC pc, PetscOptionItems *PetscOptionsObject) 1548 { 1549 PC_MG *mg = (PC_MG *)pc->data; 1550 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1551 PetscBool flag; 1552 MPI_Comm comm; 1553 char prefix[256], tname[32]; 1554 PetscInt i, n; 1555 const char *pcpre; 1556 static const char *LayoutTypes[] = {"compact", "spread", "PCGAMGLayoutType", "PC_GAMG_LAYOUT", NULL}; 1557 PetscFunctionBegin; 1558 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 1559 PetscOptionsHeadBegin(PetscOptionsObject, "GAMG options"); 1560 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)); 1561 if (flag) PetscCall(PCGAMGSetType(pc, tname)); 1562 PetscCall(PetscOptionsBool("-pc_gamg_repartition", "Repartion coarse grids", "PCGAMGSetRepartition", pc_gamg->repart, &pc_gamg->repart, NULL)); 1563 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)); 1564 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)); 1565 PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation", "Reuse prolongation operator", "PCGAMGReuseInterpolation", pc_gamg->reuse_prol, &pc_gamg->reuse_prol, NULL)); 1566 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)); 1567 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)); 1568 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)); 1569 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, 1570 (PetscEnum)pc_gamg->layout_type, (PetscEnum *)&pc_gamg->layout_type, NULL)); 1571 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)); 1572 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)); 1573 PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale", "Scaling of threshold for each level not specified", "PCGAMGSetThresholdScale", pc_gamg->threshold_scale, &pc_gamg->threshold_scale, NULL)); 1574 n = PETSC_MG_MAXLEVELS; 1575 PetscCall(PetscOptionsRealArray("-pc_gamg_threshold", "Relative threshold to use for dropping edges in aggregation graph", "PCGAMGSetThreshold", pc_gamg->threshold, &n, &flag)); 1576 if (!flag || n < PETSC_MG_MAXLEVELS) { 1577 if (!flag) n = 1; 1578 i = n; 1579 do { 1580 pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale; 1581 } while (++i < PETSC_MG_MAXLEVELS); 1582 } 1583 PetscCall(PetscOptionsInt("-pc_mg_levels", "Set number of MG levels (should get from base class)", "PCGAMGSetNlevels", pc_gamg->Nlevels, &pc_gamg->Nlevels, NULL)); 1584 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); 1585 n = PETSC_MG_MAXLEVELS; 1586 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)); 1587 if (!flag) i = 0; 1588 else i = n; 1589 do { 1590 pc_gamg->level_reduction_factors[i] = -1; 1591 } while (++i < PETSC_MG_MAXLEVELS); 1592 { 1593 PetscReal eminmax[2] = {0., 0.}; 1594 n = 2; 1595 PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues", "extreme eigenvalues for smoothed aggregation", "PCGAMGSetEigenvalues", eminmax, &n, &flag)); 1596 if (flag) { 1597 PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues"); 1598 PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0])); 1599 } 1600 } 1601 1602 /* set options for subtype */ 1603 PetscCall((*pc_gamg->ops->setfromoptions)(pc, PetscOptionsObject)); 1604 1605 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1606 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : "")); 1607 PetscOptionsHeadEnd(); 1608 PetscFunctionReturn(PETSC_SUCCESS); 1609 } 1610 1611 /*MC 1612 PCGAMG - Geometric algebraic multigrid (AMG) preconditioner 1613 1614 Options Database Keys: 1615 + -pc_gamg_type <type,default=agg> - one of agg, geo, or classical (only smoothed aggregation, agg, supported) 1616 . -pc_gamg_repartition <bool,default=false> - repartition the degrees of freedom across the coarse grids as they are determined 1617 . -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 1618 . -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> 1619 equations on each process that has degrees of freedom 1620 . -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for. 1621 . -pc_gamg_reuse_interpolation <bool,default=true> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations (should always be true) 1622 . -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) 1623 - -pc_gamg_threshold_scale <scale,default=1> - Scaling of threshold on each coarser grid if not specified 1624 1625 Options Database Keys for Aggregation: 1626 + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation 1627 . -pc_gamg_aggressive_coarsening <n,default=1> - number of aggressive coarsening (MIS-2) levels from finest. 1628 . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening 1629 . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm 1630 - -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive') 1631 1632 Options Database Keys for Multigrid: 1633 + -pc_mg_cycle_type <v> - v or w, see `PCMGSetCycleType()` 1634 . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see PCMGSetDistinctSmoothUp() 1635 . -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade 1636 - -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 1637 1638 Level: intermediate 1639 1640 Notes: 1641 To obtain good performance for `PCGAMG` for vector valued problems you must 1642 call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point 1643 call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator 1644 1645 The many options for `PCMG` also work directly for `PCGAMG` such as controlling the smoothers on each level etc. 1646 1647 See [the Users Manual section on PCGAMG](sec_amg) and [the Users Manual section on PCMG](sec_mg)for more details. 1648 1649 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, 1650 `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()` 1651 M*/ 1652 PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc) 1653 { 1654 PC_GAMG *pc_gamg; 1655 PC_MG *mg; 1656 1657 PetscFunctionBegin; 1658 /* register AMG type */ 1659 PetscCall(PCGAMGInitializePackage()); 1660 1661 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1662 PetscCall(PCSetType(pc, PCMG)); 1663 PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG)); 1664 1665 /* create a supporting struct and attach it to pc */ 1666 PetscCall(PetscNew(&pc_gamg)); 1667 PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL)); 1668 mg = (PC_MG *)pc->data; 1669 mg->innerctx = pc_gamg; 1670 1671 PetscCall(PetscNew(&pc_gamg->ops)); 1672 1673 /* these should be in subctx but repartitioning needs simple arrays */ 1674 pc_gamg->data_sz = 0; 1675 pc_gamg->data = NULL; 1676 1677 /* overwrite the pointers of PCMG by the functions of base class PCGAMG */ 1678 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1679 pc->ops->setup = PCSetUp_GAMG; 1680 pc->ops->reset = PCReset_GAMG; 1681 pc->ops->destroy = PCDestroy_GAMG; 1682 mg->view = PCView_GAMG; 1683 1684 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", PCGAMGSetProcEqLim_GAMG)); 1685 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", PCGAMGSetCoarseEqLim_GAMG)); 1686 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", PCGAMGSetRepartition_GAMG)); 1687 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", PCGAMGSetEigenvalues_GAMG)); 1688 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", PCGAMGSetUseSAEstEig_GAMG)); 1689 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", PCGAMGSetRecomputeEstEig_GAMG)); 1690 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", PCGAMGSetReuseInterpolation_GAMG)); 1691 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", PCGAMGASMSetUseAggs_GAMG)); 1692 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", PCGAMGSetParallelCoarseGridSolve_GAMG)); 1693 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", PCGAMGSetCpuPinCoarseGrids_GAMG)); 1694 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", PCGAMGSetCoarseGridLayoutType_GAMG)); 1695 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", PCGAMGSetThreshold_GAMG)); 1696 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", PCGAMGSetRankReductionFactors_GAMG)); 1697 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", PCGAMGSetThresholdScale_GAMG)); 1698 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", PCGAMGSetType_GAMG)); 1699 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", PCGAMGGetType_GAMG)); 1700 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", PCGAMGSetNlevels_GAMG)); 1701 pc_gamg->repart = PETSC_FALSE; 1702 pc_gamg->reuse_prol = PETSC_TRUE; 1703 pc_gamg->use_aggs_in_asm = PETSC_FALSE; 1704 pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE; 1705 pc_gamg->cpu_pin_coarse_grids = PETSC_FALSE; 1706 pc_gamg->layout_type = PCGAMG_LAYOUT_SPREAD; 1707 pc_gamg->min_eq_proc = 50; 1708 pc_gamg->coarse_eq_limit = 50; 1709 for (int i = 0; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = -1; 1710 pc_gamg->threshold_scale = 1.; 1711 pc_gamg->Nlevels = PETSC_MG_MAXLEVELS; 1712 pc_gamg->current_level = 0; /* don't need to init really */ 1713 pc_gamg->use_sa_esteig = PETSC_TRUE; 1714 pc_gamg->recompute_esteig = PETSC_TRUE; 1715 pc_gamg->emin = 0; 1716 pc_gamg->emax = 0; 1717 1718 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1719 1720 /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */ 1721 PetscCall(PCGAMGSetType(pc, PCGAMGAGG)); 1722 PetscFunctionReturn(PETSC_SUCCESS); 1723 } 1724 1725 /*@C 1726 PCGAMGInitializePackage - This function initializes everything in the `PCGAMG` package. It is called 1727 from `PCInitializePackage()`. 1728 1729 Level: developer 1730 1731 .seealso: [](ch_ksp), `PetscInitialize()` 1732 @*/ 1733 PetscErrorCode PCGAMGInitializePackage(void) 1734 { 1735 PetscInt l; 1736 1737 PetscFunctionBegin; 1738 if (PCGAMGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 1739 PCGAMGPackageInitialized = PETSC_TRUE; 1740 PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGGEO, PCCreateGAMG_GEO)); 1741 PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGAGG, PCCreateGAMG_AGG)); 1742 PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGCLASSICAL, PCCreateGAMG_Classical)); 1743 PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage)); 1744 1745 /* general events */ 1746 PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP])); 1747 PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH])); 1748 PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN])); 1749 PetscCall(PetscLogEventRegister(" GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS])); 1750 PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL])); 1751 PetscCall(PetscLogEventRegister(" GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA])); 1752 PetscCall(PetscLogEventRegister(" GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB])); 1753 PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT])); 1754 PetscCall(PetscLogEventRegister(" GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM])); 1755 PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL])); 1756 PetscCall(PetscLogEventRegister(" GAMG PtAP", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP])); 1757 PetscCall(PetscLogEventRegister(" GAMG Reduce", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE])); 1758 PetscCall(PetscLogEventRegister(" GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART])); 1759 /* PetscCall(PetscLogEventRegister(" GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */ 1760 /* PetscCall(PetscLogEventRegister(" GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */ 1761 /* PetscCall(PetscLogEventRegister(" GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */ 1762 for (l = 0; l < PETSC_MG_MAXLEVELS; l++) { 1763 char ename[32]; 1764 1765 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Squ l%02" PetscInt_FMT, l)); 1766 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0])); 1767 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Gal l%02" PetscInt_FMT, l)); 1768 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1])); 1769 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Opt l%02" PetscInt_FMT, l)); 1770 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2])); 1771 } 1772 #if defined(GAMG_STAGES) 1773 { /* create timer stages */ 1774 char str[32]; 1775 PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", 0)); 1776 PetscCall(PetscLogStageRegister(str, &gamg_stages[0])); 1777 } 1778 #endif 1779 PetscFunctionReturn(PETSC_SUCCESS); 1780 } 1781 1782 /*@C 1783 PCGAMGFinalizePackage - This function frees everything from the `PCGAMG` package. It is 1784 called from `PetscFinalize()` automatically. 1785 1786 Level: developer 1787 1788 .seealso: [](ch_ksp), `PetscFinalize()` 1789 @*/ 1790 PetscErrorCode PCGAMGFinalizePackage(void) 1791 { 1792 PetscFunctionBegin; 1793 PCGAMGPackageInitialized = PETSC_FALSE; 1794 PetscCall(PetscFunctionListDestroy(&GAMGList)); 1795 PetscFunctionReturn(PETSC_SUCCESS); 1796 } 1797 1798 /*@C 1799 PCGAMGRegister - Register a `PCGAMG` implementation. 1800 1801 Input Parameters: 1802 + type - string that will be used as the name of the `PCGAMG` type. 1803 - create - function for creating the gamg context. 1804 1805 Level: developer 1806 1807 .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()` 1808 @*/ 1809 PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC)) 1810 { 1811 PetscFunctionBegin; 1812 PetscCall(PCGAMGInitializePackage()); 1813 PetscCall(PetscFunctionListAdd(&GAMGList, type, create)); 1814 PetscFunctionReturn(PETSC_SUCCESS); 1815 } 1816 1817 /*@C 1818 PCGAMGCreateGraph - Creates a graph that is used by the ``PCGAMGType`` in the coarsening process 1819 1820 Input Parameters: 1821 + pc - the `PCGAMG` 1822 - A - the matrix, for any level 1823 1824 Output Parameter: 1825 . G - the graph 1826 1827 Level: advanced 1828 1829 .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()` 1830 @*/ 1831 PetscErrorCode PCGAMGCreateGraph(PC pc, Mat A, Mat *G) 1832 { 1833 PC_MG *mg = (PC_MG *)pc->data; 1834 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1835 1836 PetscFunctionBegin; 1837 PetscCall(pc_gamg->ops->creategraph(pc, A, G)); 1838 PetscFunctionReturn(PETSC_SUCCESS); 1839 } 1840