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%D",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. %D nodes. Change new active set size %d --> %d (devCount=%d #nodes=%D)\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.")); 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 PetscCheckFalse(nactive%pc_gamg->level_reduction_factors[pc_gamg->current_level],PETSC_COMM_SELF,PETSC_ERR_PLIB,"odd number of active process %D wrt reduction factor %D",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/%D)\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 PetscCheckFalse(ncrs_eq % cr_bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"ncrs_eq %D not divisible by cr_bs %D",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)=%D\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, %D 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_%D.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 PetscCheckFalse(new_size==nactive,PETSC_COMM_SELF,PETSC_ERR_PLIB,"new_size==nactive. Should not happen"); 285 PetscCall(PetscInfo(pc,"%s: Number of equations (loc) %D 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 auxilary 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 %D 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 %D\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_%d_",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 %D\n",((PetscObject)pc)->prefix,level)); 558 } else { 559 PetscCall(PetscInfo(pc,"%s: RAP after first solve, new matrix level %D\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 %D) N=%D, n data rows=%D, n data cols=%D, nnz/row (ave)=%d, np=%D\n",((PetscObject)pc)->prefix,0,M,pc_gamg->data_cell_rows,pc_gamg->data_cell_cols,(int)(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 %D",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 %D\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 PetscCheckFalse(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: %D) N=%D, n data cols=%D, nnz/row (ave)=%d, %D active pes\n",((PetscObject)pc)->prefix,level1,M,pc_gamg->data_cell_cols,(int)(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 %D. Grid too small: %D 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: %D 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 %D 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 PetscErrorCode ierr; 777 ierr = PetscObjectOptionsBegin((PetscObject)pc);PetscCall(ierr); 778 PetscCall(PCSetFromOptions_MG(PetscOptionsObject,pc)); 779 ierr = PetscOptionsEnd();PetscCall(ierr); 780 PetscCall(PCMGSetGalerkin(pc,PC_MG_GALERKIN_EXTERNAL)); 781 782 /* setup cheby eigen estimates from SA */ 783 if (pc_gamg->use_sa_esteig) { 784 for (lidx = 1, level = pc_gamg->Nlevels-2; level >= 0 ; lidx++, level--) { 785 KSP smoother; 786 PetscBool ischeb; 787 788 PetscCall(PCMGGetSmoother(pc, lidx, &smoother)); 789 PetscCall(PetscObjectTypeCompare((PetscObject)smoother,KSPCHEBYSHEV,&ischeb)); 790 if (ischeb) { 791 KSP_Chebyshev *cheb = (KSP_Chebyshev*)smoother->data; 792 793 // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG 794 if (mg->max_eigen_DinvA[level] > 0) { 795 // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it. 796 // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.) 797 PetscReal emax,emin; 798 799 emin = mg->min_eigen_DinvA[level]; 800 emax = mg->max_eigen_DinvA[level]; 801 PetscCall(PetscInfo(pc,"%s: PCSetUp_GAMG: call KSPChebyshevSetEigenvalues on level %D (N=%D) with emax = %g emin = %g\n",((PetscObject)pc)->prefix,level,Aarr[level]->rmap->N,(double)emax,(double)emin)); 802 cheb->emin_provided = emin; 803 cheb->emax_provided = emax; 804 } 805 } 806 } 807 } 808 809 PetscCall(PCSetUp_MG(pc)); 810 811 /* clean up */ 812 for (level=1; level<pc_gamg->Nlevels; level++) { 813 PetscCall(MatDestroy(&Parr[level])); 814 PetscCall(MatDestroy(&Aarr[level])); 815 } 816 } else { 817 KSP smoother; 818 819 PetscCall(PetscInfo(pc,"%s: One level solver used (system is seen as DD). Using default solver.\n")); 820 PetscCall(PCMGGetSmoother(pc, 0, &smoother)); 821 PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0])); 822 PetscCall(KSPSetType(smoother, KSPPREONLY)); 823 PetscCall(PCSetUp_MG(pc)); 824 } 825 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP],0,0,0,0)); 826 PetscFunctionReturn(0); 827 } 828 829 /* ------------------------------------------------------------------------- */ 830 /* 831 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 832 that was created with PCCreate_GAMG(). 833 834 Input Parameter: 835 . pc - the preconditioner context 836 837 Application Interface Routine: PCDestroy() 838 */ 839 PetscErrorCode PCDestroy_GAMG(PC pc) 840 { 841 PC_MG *mg = (PC_MG*)pc->data; 842 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 843 844 PetscFunctionBegin; 845 PetscCall(PCReset_GAMG(pc)); 846 if (pc_gamg->ops->destroy) { 847 PetscCall((*pc_gamg->ops->destroy)(pc)); 848 } 849 PetscCall(PetscFree(pc_gamg->ops)); 850 PetscCall(PetscFree(pc_gamg->gamg_type_name)); 851 PetscCall(PetscFree(pc_gamg)); 852 PetscCall(PCDestroy_MG(pc)); 853 PetscFunctionReturn(0); 854 } 855 856 /*@ 857 PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction. 858 859 Logically Collective on PC 860 861 Input Parameters: 862 + pc - the preconditioner context 863 - n - the number of equations 864 865 Options Database Key: 866 . -pc_gamg_process_eq_limit <limit> - set the limit 867 868 Notes: 869 GAMG will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process 870 that has degrees of freedom 871 872 Level: intermediate 873 874 .seealso: PCGAMGSetCoarseEqLim(), PCGAMGSetRankReductionFactors() 875 @*/ 876 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 877 { 878 PetscFunctionBegin; 879 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 880 PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n)); 881 PetscFunctionReturn(0); 882 } 883 884 static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 885 { 886 PC_MG *mg = (PC_MG*)pc->data; 887 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 888 889 PetscFunctionBegin; 890 if (n>0) pc_gamg->min_eq_proc = n; 891 PetscFunctionReturn(0); 892 } 893 894 /*@ 895 PCGAMGSetCoarseEqLim - Set maximum number of equations on coarsest grid. 896 897 Collective on PC 898 899 Input Parameters: 900 + pc - the preconditioner context 901 - n - maximum number of equations to aim for 902 903 Options Database Key: 904 . -pc_gamg_coarse_eq_limit <limit> - set the limit 905 906 Notes: For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid 907 has less than 1000 unknowns. 908 909 Level: intermediate 910 911 .seealso: PCGAMGSetProcEqLim(), PCGAMGSetRankReductionFactors() 912 @*/ 913 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n) 914 { 915 PetscFunctionBegin; 916 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 917 PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n)); 918 PetscFunctionReturn(0); 919 } 920 921 static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n) 922 { 923 PC_MG *mg = (PC_MG*)pc->data; 924 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 925 926 PetscFunctionBegin; 927 if (n>0) pc_gamg->coarse_eq_limit = n; 928 PetscFunctionReturn(0); 929 } 930 931 /*@ 932 PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids 933 934 Collective on PC 935 936 Input Parameters: 937 + pc - the preconditioner context 938 - n - PETSC_TRUE or PETSC_FALSE 939 940 Options Database Key: 941 . -pc_gamg_repartition <true,false> - turn on the repartitioning 942 943 Notes: 944 this will generally improve the loading balancing of the work on each level 945 946 Level: intermediate 947 948 @*/ 949 PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n) 950 { 951 PetscFunctionBegin; 952 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 953 PetscTryMethod(pc,"PCGAMGSetRepartition_C",(PC,PetscBool),(pc,n)); 954 PetscFunctionReturn(0); 955 } 956 957 static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n) 958 { 959 PC_MG *mg = (PC_MG*)pc->data; 960 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 961 962 PetscFunctionBegin; 963 pc_gamg->repart = n; 964 PetscFunctionReturn(0); 965 } 966 967 /*@ 968 PCGAMGSetUseSAEstEig - Use eigen estimate from smoothed aggregation for Chebyshev smoother 969 970 Collective on PC 971 972 Input Parameters: 973 + pc - the preconditioner context 974 - n - number of its 975 976 Options Database Key: 977 . -pc_gamg_use_sa_esteig <true,false> - use the eigen estimate 978 979 Notes: 980 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$. 981 Eigenvalue estimates (based on a few CG or GMRES iterations) are computed to choose $\omega$ so that this is a stable smoothing operation. 982 If Chebyshev with Jacobi (diagonal) preconditioning is used for smoothing, then the eigenvalue estimates can be reused. 983 This option is only used when the smoother uses Jacobi, and should be turned off if a different PCJacobiType is used. 984 It became default in PETSc 3.17. 985 986 Level: advanced 987 988 .seealso: KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet() 989 @*/ 990 PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool n) 991 { 992 PetscFunctionBegin; 993 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 994 PetscTryMethod(pc,"PCGAMGSetUseSAEstEig_C",(PC,PetscBool),(pc,n)); 995 PetscFunctionReturn(0); 996 } 997 998 static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool n) 999 { 1000 PC_MG *mg = (PC_MG*)pc->data; 1001 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1002 1003 PetscFunctionBegin; 1004 pc_gamg->use_sa_esteig = n; 1005 PetscFunctionReturn(0); 1006 } 1007 1008 /*@ 1009 PCGAMGSetEigenvalues - Set eigenvalues 1010 1011 Collective on PC 1012 1013 Input Parameters: 1014 + pc - the preconditioner context 1015 - emax - max eigenvalue 1016 . emin - min eigenvalue 1017 1018 Options Database Key: 1019 . -pc_gamg_eigenvalues <emin,emax> - estimates of the eigenvalues 1020 1021 Level: intermediate 1022 1023 .seealso: PCGAMGSetUseSAEstEig() 1024 @*/ 1025 PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax,PetscReal emin) 1026 { 1027 PetscFunctionBegin; 1028 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1029 PetscTryMethod(pc,"PCGAMGSetEigenvalues_C",(PC,PetscReal,PetscReal),(pc,emax,emin)); 1030 PetscFunctionReturn(0); 1031 } 1032 1033 static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc,PetscReal emax,PetscReal emin) 1034 { 1035 PC_MG *mg = (PC_MG*)pc->data; 1036 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1037 1038 PetscFunctionBegin; 1039 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); 1040 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); 1041 pc_gamg->emax = emax; 1042 pc_gamg->emin = emin; 1043 PetscFunctionReturn(0); 1044 } 1045 1046 /*@ 1047 PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding algebraic multigrid preconditioner 1048 1049 Collective on PC 1050 1051 Input Parameters: 1052 + pc - the preconditioner context 1053 - n - PETSC_TRUE or PETSC_FALSE 1054 1055 Options Database Key: 1056 . -pc_gamg_reuse_interpolation <true,false> - reuse the previous interpolation 1057 1058 Level: intermediate 1059 1060 Notes: 1061 May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows 1062 rebuilding the preconditioner quicker. 1063 1064 @*/ 1065 PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n) 1066 { 1067 PetscFunctionBegin; 1068 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1069 PetscTryMethod(pc,"PCGAMGSetReuseInterpolation_C",(PC,PetscBool),(pc,n)); 1070 PetscFunctionReturn(0); 1071 } 1072 1073 static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n) 1074 { 1075 PC_MG *mg = (PC_MG*)pc->data; 1076 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1077 1078 PetscFunctionBegin; 1079 pc_gamg->reuse_prol = n; 1080 PetscFunctionReturn(0); 1081 } 1082 1083 /*@ 1084 PCGAMGASMSetUseAggs - Have the PCGAMG smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner. 1085 1086 Collective on PC 1087 1088 Input Parameters: 1089 + pc - the preconditioner context 1090 - flg - PETSC_TRUE to use aggregates, PETSC_FALSE to not 1091 1092 Options Database Key: 1093 . -pc_gamg_asm_use_agg <true,false> - use aggregates to define the additive Schwarz subdomains 1094 1095 Level: intermediate 1096 1097 @*/ 1098 PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg) 1099 { 1100 PetscFunctionBegin; 1101 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1102 PetscTryMethod(pc,"PCGAMGASMSetUseAggs_C",(PC,PetscBool),(pc,flg)); 1103 PetscFunctionReturn(0); 1104 } 1105 1106 static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg) 1107 { 1108 PC_MG *mg = (PC_MG*)pc->data; 1109 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1110 1111 PetscFunctionBegin; 1112 pc_gamg->use_aggs_in_asm = flg; 1113 PetscFunctionReturn(0); 1114 } 1115 1116 /*@ 1117 PCGAMGSetUseParallelCoarseGridSolve - allow a parallel coarse grid solver 1118 1119 Collective on PC 1120 1121 Input Parameters: 1122 + pc - the preconditioner context 1123 - flg - PETSC_TRUE to not force coarse grid onto one processor 1124 1125 Options Database Key: 1126 . -pc_gamg_use_parallel_coarse_grid_solver - use a parallel coarse grid direct solver 1127 1128 Level: intermediate 1129 1130 .seealso: PCGAMGSetCoarseGridLayoutType(), PCGAMGSetCpuPinCoarseGrids() 1131 @*/ 1132 PetscErrorCode PCGAMGSetUseParallelCoarseGridSolve(PC pc, PetscBool flg) 1133 { 1134 PetscFunctionBegin; 1135 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1136 PetscTryMethod(pc,"PCGAMGSetUseParallelCoarseGridSolve_C",(PC,PetscBool),(pc,flg)); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 static PetscErrorCode PCGAMGSetUseParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg) 1141 { 1142 PC_MG *mg = (PC_MG*)pc->data; 1143 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1144 1145 PetscFunctionBegin; 1146 pc_gamg->use_parallel_coarse_grid_solver = flg; 1147 PetscFunctionReturn(0); 1148 } 1149 1150 /*@ 1151 PCGAMGSetCpuPinCoarseGrids - pin reduced grids to CPU 1152 1153 Collective on PC 1154 1155 Input Parameters: 1156 + pc - the preconditioner context 1157 - flg - PETSC_TRUE to pin coarse grids to CPU 1158 1159 Options Database Key: 1160 . -pc_gamg_cpu_pin_coarse_grids - pin the coarse grids to the CPU 1161 1162 Level: intermediate 1163 1164 .seealso: PCGAMGSetCoarseGridLayoutType(), PCGAMGSetUseParallelCoarseGridSolve() 1165 @*/ 1166 PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg) 1167 { 1168 PetscFunctionBegin; 1169 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1170 PetscTryMethod(pc,"PCGAMGSetCpuPinCoarseGrids_C",(PC,PetscBool),(pc,flg)); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg) 1175 { 1176 PC_MG *mg = (PC_MG*)pc->data; 1177 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1178 1179 PetscFunctionBegin; 1180 pc_gamg->cpu_pin_coarse_grids = flg; 1181 PetscFunctionReturn(0); 1182 } 1183 1184 /*@ 1185 PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type) 1186 1187 Collective on PC 1188 1189 Input Parameters: 1190 + pc - the preconditioner context 1191 - flg - Layout type 1192 1193 Options Database Key: 1194 . -pc_gamg_coarse_grid_layout_type - place the coarse grids with natural ordering 1195 1196 Level: intermediate 1197 1198 .seealso: PCGAMGSetUseParallelCoarseGridSolve(), PCGAMGSetCpuPinCoarseGrids() 1199 @*/ 1200 PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg) 1201 { 1202 PetscFunctionBegin; 1203 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1204 PetscTryMethod(pc,"PCGAMGSetCoarseGridLayoutType_C",(PC,PCGAMGLayoutType),(pc,flg)); 1205 PetscFunctionReturn(0); 1206 } 1207 1208 static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg) 1209 { 1210 PC_MG *mg = (PC_MG*)pc->data; 1211 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1212 1213 PetscFunctionBegin; 1214 pc_gamg->layout_type = flg; 1215 PetscFunctionReturn(0); 1216 } 1217 1218 /*@ 1219 PCGAMGSetNlevels - Sets the maximum number of levels PCGAMG will use 1220 1221 Not collective on PC 1222 1223 Input Parameters: 1224 + pc - the preconditioner 1225 - n - the maximum number of levels to use 1226 1227 Options Database Key: 1228 . -pc_mg_levels <n> - set the maximum number of levels to allow 1229 1230 Level: intermediate 1231 1232 @*/ 1233 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n) 1234 { 1235 PetscFunctionBegin; 1236 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1237 PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n)); 1238 PetscFunctionReturn(0); 1239 } 1240 1241 static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n) 1242 { 1243 PC_MG *mg = (PC_MG*)pc->data; 1244 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1245 1246 PetscFunctionBegin; 1247 pc_gamg->Nlevels = n; 1248 PetscFunctionReturn(0); 1249 } 1250 1251 /*@ 1252 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 1253 1254 Not collective on PC 1255 1256 Input Parameters: 1257 + pc - the preconditioner context 1258 . 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 1259 - n - number of threshold values provided in array 1260 1261 Options Database Key: 1262 . -pc_gamg_threshold <threshold> - the threshold to drop edges 1263 1264 Notes: 1265 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. 1266 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. 1267 1268 If n is less than the total number of coarsenings (see PCGAMGSetNlevels()), then threshold scaling (see PCGAMGSetThresholdScale()) is used for each successive coarsening. 1269 In this case, PCGAMGSetThresholdScale() must be called before PCGAMGSetThreshold(). 1270 If n is greater than the total number of levels, the excess entries in threshold will not be used. 1271 1272 Level: intermediate 1273 1274 .seealso: PCGAMGFilterGraph(), PCGAMGSetSquareGraph() 1275 @*/ 1276 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n) 1277 { 1278 PetscFunctionBegin; 1279 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1280 if (n) PetscValidRealPointer(v,2); 1281 PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal[],PetscInt),(pc,v,n)); 1282 PetscFunctionReturn(0); 1283 } 1284 1285 static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n) 1286 { 1287 PC_MG *mg = (PC_MG*)pc->data; 1288 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1289 PetscInt i; 1290 PetscFunctionBegin; 1291 for (i=0; i<PetscMin(n,PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i]; 1292 for (; i<PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i-1]*pc_gamg->threshold_scale; 1293 PetscFunctionReturn(0); 1294 } 1295 1296 /*@ 1297 PCGAMGSetRankReductionFactors - Set manual schedule for process reduction on coarse grids 1298 1299 Collective on PC 1300 1301 Input Parameters: 1302 + pc - the preconditioner context 1303 . v - array of reduction factors. 0 for fist value forces a reduction to one process/device on first level in Cuda 1304 - n - number of values provided in array 1305 1306 Options Database Key: 1307 . -pc_gamg_rank_reduction_factors <factors> - provide the schedule 1308 1309 Level: intermediate 1310 1311 .seealso: PCGAMGSetProcEqLim(), PCGAMGSetCoarseEqLim() 1312 @*/ 1313 PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n) 1314 { 1315 PetscFunctionBegin; 1316 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1317 if (n) PetscValidIntPointer(v,2); 1318 PetscTryMethod(pc,"PCGAMGSetRankReductionFactors_C",(PC,PetscInt[],PetscInt),(pc,v,n)); 1319 PetscFunctionReturn(0); 1320 } 1321 1322 static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n) 1323 { 1324 PC_MG *mg = (PC_MG*)pc->data; 1325 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1326 PetscInt i; 1327 PetscFunctionBegin; 1328 for (i=0; i<PetscMin(n,PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i]; 1329 for (; i<PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */ 1330 PetscFunctionReturn(0); 1331 } 1332 1333 /*@ 1334 PCGAMGSetThresholdScale - Relative threshold reduction at each level 1335 1336 Not collective on PC 1337 1338 Input Parameters: 1339 + pc - the preconditioner context 1340 - scale - the threshold value reduction, usually < 1.0 1341 1342 Options Database Key: 1343 . -pc_gamg_threshold_scale <v> - set the relative threshold reduction on each level 1344 1345 Notes: 1346 The initial threshold (for an arbitrary number of levels starting from the finest) can be set with PCGAMGSetThreshold(). 1347 This scaling is used for each subsequent coarsening, but must be called before PCGAMGSetThreshold(). 1348 1349 Level: advanced 1350 1351 .seealso: PCGAMGSetThreshold() 1352 @*/ 1353 PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v) 1354 { 1355 PetscFunctionBegin; 1356 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1357 PetscTryMethod(pc,"PCGAMGSetThresholdScale_C",(PC,PetscReal),(pc,v)); 1358 PetscFunctionReturn(0); 1359 } 1360 1361 static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v) 1362 { 1363 PC_MG *mg = (PC_MG*)pc->data; 1364 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1365 PetscFunctionBegin; 1366 pc_gamg->threshold_scale = v; 1367 PetscFunctionReturn(0); 1368 } 1369 1370 /*@C 1371 PCGAMGSetType - Set solution method 1372 1373 Collective on PC 1374 1375 Input Parameters: 1376 + pc - the preconditioner context 1377 - type - PCGAMGAGG, PCGAMGGEO, or PCGAMGCLASSICAL 1378 1379 Options Database Key: 1380 . -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply 1381 1382 Level: intermediate 1383 1384 .seealso: PCGAMGGetType(), PCGAMG, PCGAMGType 1385 @*/ 1386 PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type) 1387 { 1388 PetscFunctionBegin; 1389 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1390 PetscTryMethod(pc,"PCGAMGSetType_C",(PC,PCGAMGType),(pc,type)); 1391 PetscFunctionReturn(0); 1392 } 1393 1394 /*@C 1395 PCGAMGGetType - Get solution method 1396 1397 Collective on PC 1398 1399 Input Parameter: 1400 . pc - the preconditioner context 1401 1402 Output Parameter: 1403 . type - the type of algorithm used 1404 1405 Level: intermediate 1406 1407 .seealso: PCGAMGSetType(), PCGAMGType 1408 @*/ 1409 PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type) 1410 { 1411 PetscFunctionBegin; 1412 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1413 PetscUseMethod(pc,"PCGAMGGetType_C",(PC,PCGAMGType*),(pc,type)); 1414 PetscFunctionReturn(0); 1415 } 1416 1417 static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type) 1418 { 1419 PC_MG *mg = (PC_MG*)pc->data; 1420 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1421 1422 PetscFunctionBegin; 1423 *type = pc_gamg->type; 1424 PetscFunctionReturn(0); 1425 } 1426 1427 static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type) 1428 { 1429 PC_MG *mg = (PC_MG*)pc->data; 1430 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1431 PetscErrorCode (*r)(PC); 1432 1433 PetscFunctionBegin; 1434 pc_gamg->type = type; 1435 PetscCall(PetscFunctionListFind(GAMGList,type,&r)); 1436 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown GAMG type %s given",type); 1437 if (pc_gamg->ops->destroy) { 1438 PetscCall((*pc_gamg->ops->destroy)(pc)); 1439 PetscCall(PetscMemzero(pc_gamg->ops,sizeof(struct _PCGAMGOps))); 1440 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1441 /* cleaning up common data in pc_gamg - this should disapear someday */ 1442 pc_gamg->data_cell_cols = 0; 1443 pc_gamg->data_cell_rows = 0; 1444 pc_gamg->orig_data_cell_cols = 0; 1445 pc_gamg->orig_data_cell_rows = 0; 1446 PetscCall(PetscFree(pc_gamg->data)); 1447 pc_gamg->data_sz = 0; 1448 } 1449 PetscCall(PetscFree(pc_gamg->gamg_type_name)); 1450 PetscCall(PetscStrallocpy(type,&pc_gamg->gamg_type_name)); 1451 PetscCall((*r)(pc)); 1452 PetscFunctionReturn(0); 1453 } 1454 1455 static PetscErrorCode PCView_GAMG(PC pc,PetscViewer viewer) 1456 { 1457 PC_MG *mg = (PC_MG*)pc->data; 1458 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1459 PetscReal gc=0, oc=0; 1460 1461 PetscFunctionBegin; 1462 PetscCall(PetscViewerASCIIPrintf(viewer," GAMG specific options\n")); 1463 PetscCall(PetscViewerASCIIPrintf(viewer," Threshold for dropping small values in graph on each level =")); 1464 for (PetscInt i=0;i<mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer," %g",(double)pc_gamg->threshold[i])); 1465 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 1466 PetscCall(PetscViewerASCIIPrintf(viewer," Threshold scaling factor for each level not specified = %g\n",(double)pc_gamg->threshold_scale)); 1467 if (pc_gamg->use_aggs_in_asm) { 1468 PetscCall(PetscViewerASCIIPrintf(viewer," Using aggregates from coarsening process to define subdomains for PCASM\n")); 1469 } 1470 if (pc_gamg->use_parallel_coarse_grid_solver) { 1471 PetscCall(PetscViewerASCIIPrintf(viewer," Using parallel coarse grid solver (all coarse grid equations not put on one process)\n")); 1472 } 1473 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 1474 if (pc_gamg->cpu_pin_coarse_grids) { 1475 /* PetscCall(PetscViewerASCIIPrintf(viewer," Pinning coarse grids to the CPU)\n")); */ 1476 } 1477 #endif 1478 /* if (pc_gamg->layout_type==PCGAMG_LAYOUT_COMPACT) { */ 1479 /* PetscCall(PetscViewerASCIIPrintf(viewer," Put reduced grids on processes in natural order (ie, 0,1,2...)\n")); */ 1480 /* } else { */ 1481 /* PetscCall(PetscViewerASCIIPrintf(viewer," Put reduced grids on whole machine (ie, 0,1*f,2*f...,np-f)\n")); */ 1482 /* } */ 1483 if (pc_gamg->ops->view) { 1484 PetscCall((*pc_gamg->ops->view)(pc,viewer)); 1485 } 1486 PetscCall(PCMGGetGridComplexity(pc,&gc,&oc)); 1487 PetscCall(PetscViewerASCIIPrintf(viewer," Complexity: grid = %g operator = %g\n",gc,oc)); 1488 PetscFunctionReturn(0); 1489 } 1490 1491 PetscErrorCode PCSetFromOptions_GAMG(PetscOptionItems *PetscOptionsObject,PC pc) 1492 { 1493 PC_MG *mg = (PC_MG*)pc->data; 1494 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1495 PetscBool flag; 1496 MPI_Comm comm; 1497 char prefix[256],tname[32]; 1498 PetscInt i,n; 1499 const char *pcpre; 1500 static const char *LayoutTypes[] = {"compact","spread","PCGAMGLayoutType","PC_GAMG_LAYOUT",NULL}; 1501 PetscFunctionBegin; 1502 PetscCall(PetscObjectGetComm((PetscObject)pc,&comm)); 1503 PetscCall(PetscOptionsHead(PetscOptionsObject,"GAMG options")); 1504 PetscCall(PetscOptionsFList("-pc_gamg_type","Type of AMG method","PCGAMGSetType",GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag)); 1505 if (flag) { 1506 PetscCall(PCGAMGSetType(pc,tname)); 1507 } 1508 PetscCall(PetscOptionsBool("-pc_gamg_repartition","Repartion coarse grids","PCGAMGSetRepartition",pc_gamg->repart,&pc_gamg->repart,NULL)); 1509 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)); 1510 PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation","Reuse prolongation operator","PCGAMGReuseInterpolation",pc_gamg->reuse_prol,&pc_gamg->reuse_prol,NULL)); 1511 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)); 1512 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)); 1513 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)); 1514 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)); 1515 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)); 1516 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)); 1517 PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale","Scaling of threshold for each level not specified","PCGAMGSetThresholdScale",pc_gamg->threshold_scale,&pc_gamg->threshold_scale,NULL)); 1518 n = PETSC_MG_MAXLEVELS; 1519 PetscCall(PetscOptionsRealArray("-pc_gamg_threshold","Relative threshold to use for dropping edges in aggregation graph","PCGAMGSetThreshold",pc_gamg->threshold,&n,&flag)); 1520 if (!flag || n < PETSC_MG_MAXLEVELS) { 1521 if (!flag) n = 1; 1522 i = n; 1523 do {pc_gamg->threshold[i] = pc_gamg->threshold[i-1]*pc_gamg->threshold_scale;} while (++i<PETSC_MG_MAXLEVELS); 1524 } 1525 n = PETSC_MG_MAXLEVELS; 1526 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)); 1527 if (!flag) i = 0; 1528 else i = n; 1529 do {pc_gamg->level_reduction_factors[i] = -1;} while (++i<PETSC_MG_MAXLEVELS); 1530 PetscCall(PetscOptionsInt("-pc_mg_levels","Set number of MG levels","PCGAMGSetNlevels",pc_gamg->Nlevels,&pc_gamg->Nlevels,NULL)); 1531 { 1532 PetscReal eminmax[2] = {0., 0.}; 1533 n = 2; 1534 PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues","extreme eigenvalues for smoothed aggregation","PCGAMGSetEigenvalues",eminmax,&n,&flag)); 1535 if (flag) { 1536 PetscCheckFalse(n != 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues"); 1537 PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0])); 1538 } 1539 } 1540 /* set options for subtype */ 1541 if (pc_gamg->ops->setfromoptions) PetscCall((*pc_gamg->ops->setfromoptions)(PetscOptionsObject,pc)); 1542 1543 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1544 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%spc_gamg_",pcpre ? pcpre : "")); 1545 PetscCall(PetscOptionsTail()); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 /* -------------------------------------------------------------------------- */ 1550 /*MC 1551 PCGAMG - Geometric algebraic multigrid (AMG) preconditioner 1552 1553 Options Database Keys: 1554 + -pc_gamg_type <type> - one of agg, geo, or classical 1555 . -pc_gamg_repartition <true,default=false> - repartition the degrees of freedom accross the coarse grids as they are determined 1556 . -pc_gamg_reuse_interpolation <true,default=false> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations 1557 . -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 1558 . -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> 1559 equations on each process that has degrees of freedom 1560 . -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for. 1561 . -pc_gamg_threshold[] <thresh,default=0> - Before aggregating the graph GAMG will remove small values from the graph on each level 1562 - -pc_gamg_threshold_scale <scale,default=1> - Scaling of threshold on each coarser grid if not specified 1563 1564 Options Database Keys for default Aggregation: 1565 + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation 1566 . -pc_gamg_sym_graph <true,default=false> - symmetrize the graph before computing the aggregation 1567 - -pc_gamg_square_graph <n,default=1> - number of levels to square the graph before aggregating it 1568 1569 Multigrid options: 1570 + -pc_mg_cycles <v> - v or w, see PCMGSetCycleType() 1571 . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see PCMGSetDistinctSmoothUp() 1572 . -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade 1573 - -pc_mg_levels <levels> - Number of levels of multigrid to use. 1574 1575 Notes: 1576 In order to obtain good performance for PCGAMG for vector valued problems you must 1577 Call MatSetBlockSize() to indicate the number of degrees of freedom per grid point 1578 Call MatSetNearNullSpace() (or PCSetCoordinates() if solving the equations of elasticity) to indicate the near null space of the operator 1579 See the Users Manual Chapter 4 for more details 1580 1581 Level: intermediate 1582 1583 .seealso: PCCreate(), PCSetType(), MatSetBlockSize(), PCMGType, PCSetCoordinates(), MatSetNearNullSpace(), PCGAMGSetType(), PCGAMGAGG, PCGAMGGEO, PCGAMGCLASSICAL, PCGAMGSetProcEqLim(), 1584 PCGAMGSetCoarseEqLim(), PCGAMGSetRepartition(), PCGAMGRegister(), PCGAMGSetReuseInterpolation(), PCGAMGASMSetUseAggs(), PCGAMGSetUseParallelCoarseGridSolve(), PCGAMGSetNlevels(), PCGAMGSetThreshold(), PCGAMGGetType(), PCGAMGSetReuseInterpolation(), PCGAMGSetUseSAEstEig() 1585 M*/ 1586 1587 PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc) 1588 { 1589 PC_GAMG *pc_gamg; 1590 PC_MG *mg; 1591 1592 PetscFunctionBegin; 1593 /* register AMG type */ 1594 PetscCall(PCGAMGInitializePackage()); 1595 1596 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1597 PetscCall(PCSetType(pc, PCMG)); 1598 PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG)); 1599 1600 /* create a supporting struct and attach it to pc */ 1601 PetscCall(PetscNewLog(pc,&pc_gamg)); 1602 PetscCall(PCMGSetGalerkin(pc,PC_MG_GALERKIN_EXTERNAL)); 1603 mg = (PC_MG*)pc->data; 1604 mg->innerctx = pc_gamg; 1605 1606 PetscCall(PetscNewLog(pc,&pc_gamg->ops)); 1607 1608 /* these should be in subctx but repartitioning needs simple arrays */ 1609 pc_gamg->data_sz = 0; 1610 pc_gamg->data = NULL; 1611 1612 /* overwrite the pointers of PCMG by the functions of base class PCGAMG */ 1613 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1614 pc->ops->setup = PCSetUp_GAMG; 1615 pc->ops->reset = PCReset_GAMG; 1616 pc->ops->destroy = PCDestroy_GAMG; 1617 mg->view = PCView_GAMG; 1618 1619 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG)); 1620 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG)); 1621 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetProcEqLim_C",PCGAMGSetProcEqLim_GAMG)); 1622 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseEqLim_C",PCGAMGSetCoarseEqLim_GAMG)); 1623 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRepartition_C",PCGAMGSetRepartition_GAMG)); 1624 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetEigenvalues_C",PCGAMGSetEigenvalues_GAMG)); 1625 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetUseSAEstEig_C",PCGAMGSetUseSAEstEig_GAMG)); 1626 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetReuseInterpolation_C",PCGAMGSetReuseInterpolation_GAMG)); 1627 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGASMSetUseAggs_C",PCGAMGASMSetUseAggs_GAMG)); 1628 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetUseParallelCoarseGridSolve_C",PCGAMGSetUseParallelCoarseGridSolve_GAMG)); 1629 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCpuPinCoarseGrids_C",PCGAMGSetCpuPinCoarseGrids_GAMG)); 1630 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseGridLayoutType_C",PCGAMGSetCoarseGridLayoutType_GAMG)); 1631 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThreshold_C",PCGAMGSetThreshold_GAMG)); 1632 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRankReductionFactors_C",PCGAMGSetRankReductionFactors_GAMG)); 1633 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThresholdScale_C",PCGAMGSetThresholdScale_GAMG)); 1634 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetType_C",PCGAMGSetType_GAMG)); 1635 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGGetType_C",PCGAMGGetType_GAMG)); 1636 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNlevels_C",PCGAMGSetNlevels_GAMG)); 1637 pc_gamg->repart = PETSC_FALSE; 1638 pc_gamg->reuse_prol = PETSC_FALSE; 1639 pc_gamg->use_aggs_in_asm = PETSC_FALSE; 1640 pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE; 1641 pc_gamg->cpu_pin_coarse_grids = PETSC_FALSE; 1642 pc_gamg->layout_type = PCGAMG_LAYOUT_SPREAD; 1643 pc_gamg->min_eq_proc = 50; 1644 pc_gamg->coarse_eq_limit = 50; 1645 PetscCall(PetscArrayzero(pc_gamg->threshold,PETSC_MG_MAXLEVELS)); 1646 pc_gamg->threshold_scale = 1.; 1647 pc_gamg->Nlevels = PETSC_MG_MAXLEVELS; 1648 pc_gamg->current_level = 0; /* don't need to init really */ 1649 pc_gamg->use_sa_esteig = PETSC_TRUE; 1650 pc_gamg->emin = 0; 1651 pc_gamg->emax = 0; 1652 1653 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1654 1655 /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */ 1656 PetscCall(PCGAMGSetType(pc,PCGAMGAGG)); 1657 PetscFunctionReturn(0); 1658 } 1659 1660 /*@C 1661 PCGAMGInitializePackage - This function initializes everything in the PCGAMG package. It is called 1662 from PCInitializePackage(). 1663 1664 Level: developer 1665 1666 .seealso: PetscInitialize() 1667 @*/ 1668 PetscErrorCode PCGAMGInitializePackage(void) 1669 { 1670 PetscInt l; 1671 1672 PetscFunctionBegin; 1673 if (PCGAMGPackageInitialized) PetscFunctionReturn(0); 1674 PCGAMGPackageInitialized = PETSC_TRUE; 1675 PetscCall(PetscFunctionListAdd(&GAMGList,PCGAMGGEO,PCCreateGAMG_GEO)); 1676 PetscCall(PetscFunctionListAdd(&GAMGList,PCGAMGAGG,PCCreateGAMG_AGG)); 1677 PetscCall(PetscFunctionListAdd(&GAMGList,PCGAMGCLASSICAL,PCCreateGAMG_Classical)); 1678 PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage)); 1679 1680 /* general events */ 1681 PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP])); 1682 PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH])); 1683 PetscCall(PetscLogEventRegister(" PCGAMGFilter", PC_CLASSID, &petsc_gamg_setup_events[GAMG_FILTER])); 1684 PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN])); 1685 PetscCall(PetscLogEventRegister(" GAMG Square G", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SQUARE])); 1686 PetscCall(PetscLogEventRegister(" GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS])); 1687 PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL])); 1688 PetscCall(PetscLogEventRegister(" GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA])); 1689 PetscCall(PetscLogEventRegister(" GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB])); 1690 PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT])); 1691 PetscCall(PetscLogEventRegister(" GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM])); 1692 PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL])); 1693 PetscCall(PetscLogEventRegister(" GAMG PtAP",PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP])); 1694 PetscCall(PetscLogEventRegister(" GAMG Reduce",PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE])); 1695 PetscCall(PetscLogEventRegister(" GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART])); 1696 /* PetscCall(PetscLogEventRegister(" GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */ 1697 /* PetscCall(PetscLogEventRegister(" GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */ 1698 /* PetscCall(PetscLogEventRegister(" GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */ 1699 for (l=0;l<PETSC_MG_MAXLEVELS;l++) { 1700 char ename[32]; 1701 1702 PetscCall(PetscSNPrintf(ename,sizeof(ename),"PCGAMG Squ l%02d",l)); 1703 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0])); 1704 PetscCall(PetscSNPrintf(ename,sizeof(ename),"PCGAMG Gal l%02d",l)); 1705 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1])); 1706 PetscCall(PetscSNPrintf(ename,sizeof(ename),"PCGAMG Opt l%02d",l)); 1707 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2])); 1708 } 1709 #if defined(GAMG_STAGES) 1710 { /* create timer stages */ 1711 char str[32]; 1712 sprintf(str,"GAMG Level %d",0); 1713 PetscCall(PetscLogStageRegister(str, &gamg_stages[0])); 1714 } 1715 #endif 1716 PetscFunctionReturn(0); 1717 } 1718 1719 /*@C 1720 PCGAMGFinalizePackage - This function frees everything from the PCGAMG package. It is 1721 called from PetscFinalize() automatically. 1722 1723 Level: developer 1724 1725 .seealso: PetscFinalize() 1726 @*/ 1727 PetscErrorCode PCGAMGFinalizePackage(void) 1728 { 1729 PetscFunctionBegin; 1730 PCGAMGPackageInitialized = PETSC_FALSE; 1731 PetscCall(PetscFunctionListDestroy(&GAMGList)); 1732 PetscFunctionReturn(0); 1733 } 1734 1735 /*@C 1736 PCGAMGRegister - Register a PCGAMG implementation. 1737 1738 Input Parameters: 1739 + type - string that will be used as the name of the GAMG type. 1740 - create - function for creating the gamg context. 1741 1742 Level: advanced 1743 1744 .seealso: PCGAMGType, PCGAMG, PCGAMGSetType() 1745 @*/ 1746 PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC)) 1747 { 1748 PetscFunctionBegin; 1749 PetscCall(PCGAMGInitializePackage()); 1750 PetscCall(PetscFunctionListAdd(&GAMGList,type,create)); 1751 PetscFunctionReturn(0); 1752 } 1753