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