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