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