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