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);CHKERRMPI(ierr); 114 s_nnodes = !locrank; 115 ierr = MPI_Allreduce(&s_nnodes,&r_nnodes,1,MPIU_INT,MPI_SUM,comm);CHKERRMPI(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 ierr = MatSetOption(Pnew,MAT_FORM_EXPLICIT_TRANSPOSE,PETSC_TRUE);CHKERRQ(ierr); 423 424 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET15],0,0,0,0);CHKERRQ(ierr); 425 ierr = MatDestroy(a_P_inout);CHKERRQ(ierr); 426 427 /* output - repartitioned */ 428 *a_P_inout = Pnew; 429 } 430 ierr = ISDestroy(&new_eq_indices);CHKERRQ(ierr); 431 432 *a_nactive_proc = new_size; /* output */ 433 434 /* pinning on reduced grids, not a bad heuristic and optimization gets folded into process reduction optimization */ 435 if (pc_gamg->cpu_pin_coarse_grids) { 436 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 437 static PetscInt llev = 2; 438 ierr = PetscInfo1(pc,"Pinning level %D to the CPU\n",llev++);CHKERRQ(ierr); 439 #endif 440 ierr = MatBindToCPU(*a_Amat_crs,PETSC_TRUE);CHKERRQ(ierr); 441 ierr = MatBindToCPU(*a_P_inout,PETSC_TRUE);CHKERRQ(ierr); 442 if (1) { /* HACK: move this to MatBindCPU_MPIAIJXXX; lvec is created, need to pin it, this is done in MatSetUpMultiply_MPIAIJ. Hack */ 443 Mat A = *a_Amat_crs, P = *a_P_inout; 444 PetscMPIInt size; 445 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRMPI(ierr); 446 if (size > 1) { 447 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data, *p = (Mat_MPIAIJ*)P->data; 448 ierr = VecBindToCPU(a->lvec,PETSC_TRUE);CHKERRQ(ierr); 449 ierr = VecBindToCPU(p->lvec,PETSC_TRUE);CHKERRQ(ierr); 450 } 451 } 452 } 453 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET12],0,0,0,0);CHKERRQ(ierr); 454 } 455 PetscFunctionReturn(0); 456 } 457 458 PetscErrorCode PCGAMGSquareGraph_GAMG(PC a_pc, Mat Gmat1, Mat* Gmat2) 459 { 460 PetscErrorCode ierr; 461 const char *prefix; 462 char addp[32]; 463 PC_MG *mg = (PC_MG*)a_pc->data; 464 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 465 466 PetscFunctionBegin; 467 ierr = PCGetOptionsPrefix(a_pc,&prefix);CHKERRQ(ierr); 468 ierr = PetscInfo1(a_pc,"Square Graph on level %D\n",pc_gamg->current_level+1);CHKERRQ(ierr); 469 ierr = MatProductCreate(Gmat1,Gmat1,NULL,Gmat2);CHKERRQ(ierr); 470 ierr = MatSetOptionsPrefix(*Gmat2,prefix);CHKERRQ(ierr); 471 ierr = PetscSNPrintf(addp,sizeof(addp),"pc_gamg_square_%d_",pc_gamg->current_level);CHKERRQ(ierr); 472 ierr = MatAppendOptionsPrefix(*Gmat2,addp);CHKERRQ(ierr); 473 if ((*Gmat2)->structurally_symmetric) { 474 ierr = MatProductSetType(*Gmat2,MATPRODUCT_AB);CHKERRQ(ierr); 475 } else { 476 ierr = MatSetOption(Gmat1,MAT_FORM_EXPLICIT_TRANSPOSE,PETSC_TRUE);CHKERRQ(ierr); 477 ierr = MatProductSetType(*Gmat2,MATPRODUCT_AtB);CHKERRQ(ierr); 478 } 479 ierr = MatProductSetFromOptions(*Gmat2);CHKERRQ(ierr); 480 ierr = PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0],0,0,0,0);CHKERRQ(ierr); 481 ierr = MatProductSymbolic(*Gmat2);CHKERRQ(ierr); 482 ierr = PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0],0,0,0,0);CHKERRQ(ierr); 483 ierr = MatProductClear(*Gmat2);CHKERRQ(ierr); 484 /* we only need the sparsity, cheat and tell PETSc the matrix has been assembled */ 485 (*Gmat2)->assembled = PETSC_TRUE; 486 PetscFunctionReturn(0); 487 } 488 489 /* -------------------------------------------------------------------------- */ 490 /* 491 PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner 492 by setting data structures and options. 493 494 Input Parameter: 495 . pc - the preconditioner context 496 497 */ 498 PetscErrorCode PCSetUp_GAMG(PC pc) 499 { 500 PetscErrorCode ierr; 501 PC_MG *mg = (PC_MG*)pc->data; 502 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 503 Mat Pmat = pc->pmat; 504 PetscInt fine_level,level,level1,bs,M,N,qq,lidx,nASMBlocksArr[PETSC_MG_MAXLEVELS]; 505 MPI_Comm comm; 506 PetscMPIInt rank,size,nactivepe; 507 Mat Aarr[PETSC_MG_MAXLEVELS],Parr[PETSC_MG_MAXLEVELS]; 508 IS *ASMLocalIDsArr[PETSC_MG_MAXLEVELS]; 509 PetscLogDouble nnz0=0.,nnztot=0.; 510 MatInfo info; 511 PetscBool is_last = PETSC_FALSE; 512 513 PetscFunctionBegin; 514 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 515 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 516 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 517 518 if (pc->setupcalled) { 519 if (!pc_gamg->reuse_prol || pc->flag == DIFFERENT_NONZERO_PATTERN) { 520 /* reset everything */ 521 ierr = PCReset_MG(pc);CHKERRQ(ierr); 522 pc->setupcalled = 0; 523 } else { 524 PC_MG_Levels **mglevels = mg->levels; 525 /* just do Galerkin grids */ 526 Mat B,dA,dB; 527 528 if (pc_gamg->Nlevels > 1) { 529 PetscInt gl; 530 /* currently only handle case where mat and pmat are the same on coarser levels */ 531 ierr = KSPGetOperators(mglevels[pc_gamg->Nlevels-1]->smoothd,&dA,&dB);CHKERRQ(ierr); 532 /* (re)set to get dirty flag */ 533 ierr = KSPSetOperators(mglevels[pc_gamg->Nlevels-1]->smoothd,dA,dB);CHKERRQ(ierr); 534 535 for (level=pc_gamg->Nlevels-2,gl=0; level>=0; level--,gl++) { 536 MatReuse reuse = MAT_INITIAL_MATRIX ; 537 538 /* matrix structure can change from repartitioning or process reduction but don't know if we have process reduction here. Should fix */ 539 ierr = KSPGetOperators(mglevels[level]->smoothd,NULL,&B);CHKERRQ(ierr); 540 if (B->product) { 541 if (B->product->A == dB && B->product->B == mglevels[level+1]->interpolate) { 542 reuse = MAT_REUSE_MATRIX; 543 } 544 } 545 if (reuse == MAT_INITIAL_MATRIX) { ierr = MatDestroy(&mglevels[level]->A);CHKERRQ(ierr); } 546 if (reuse == MAT_REUSE_MATRIX) { 547 ierr = PetscInfo1(pc,"RAP after first solve, reuse matrix level %D\n",level);CHKERRQ(ierr); 548 } else { 549 ierr = PetscInfo1(pc,"RAP after first solve, new matrix level %D\n",level);CHKERRQ(ierr); 550 } 551 ierr = PetscLogEventBegin(petsc_gamg_setup_matmat_events[gl][1],0,0,0,0);CHKERRQ(ierr); 552 ierr = MatPtAP(dB,mglevels[level+1]->interpolate,reuse,PETSC_DEFAULT,&B);CHKERRQ(ierr); 553 ierr = PetscLogEventEnd(petsc_gamg_setup_matmat_events[gl][1],0,0,0,0);CHKERRQ(ierr); 554 mglevels[level]->A = B; 555 ierr = KSPSetOperators(mglevels[level]->smoothd,B,B);CHKERRQ(ierr); 556 dB = B; 557 } 558 } 559 560 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563 } 564 565 if (!pc_gamg->data) { 566 if (pc_gamg->orig_data) { 567 ierr = MatGetBlockSize(Pmat, &bs);CHKERRQ(ierr); 568 ierr = MatGetLocalSize(Pmat, &qq, NULL);CHKERRQ(ierr); 569 570 pc_gamg->data_sz = (qq/bs)*pc_gamg->orig_data_cell_rows*pc_gamg->orig_data_cell_cols; 571 pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows; 572 pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols; 573 574 ierr = PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data);CHKERRQ(ierr); 575 for (qq=0; qq<pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq]; 576 } else { 577 if (!pc_gamg->ops->createdefaultdata) SETERRQ(comm,PETSC_ERR_PLIB,"'createdefaultdata' not set(?) need to support NULL data"); 578 ierr = pc_gamg->ops->createdefaultdata(pc,Pmat);CHKERRQ(ierr); 579 } 580 } 581 582 /* cache original data for reuse */ 583 if (!pc_gamg->orig_data && (PetscBool)(!pc_gamg->reuse_prol)) { 584 ierr = PetscMalloc1(pc_gamg->data_sz, &pc_gamg->orig_data);CHKERRQ(ierr); 585 for (qq=0; qq<pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq]; 586 pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows; 587 pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols; 588 } 589 590 /* get basic dims */ 591 ierr = MatGetBlockSize(Pmat, &bs);CHKERRQ(ierr); 592 ierr = MatGetSize(Pmat, &M, &N);CHKERRQ(ierr); 593 594 ierr = MatGetInfo(Pmat,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); /* global reduction */ 595 nnz0 = info.nz_used; 596 nnztot = info.nz_used; 597 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); 598 599 /* Get A_i and R_i */ 600 for (level=0, Aarr[0]=Pmat, nactivepe = size; level < (pc_gamg->Nlevels-1) && (!level || M>pc_gamg->coarse_eq_limit); level++) { 601 pc_gamg->current_level = level; 602 if (level >= PETSC_MG_MAXLEVELS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Too many levels %D",level); 603 level1 = level + 1; 604 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr); 605 #if defined(GAMG_STAGES) 606 ierr = PetscLogStagePush(gamg_stages[level]);CHKERRQ(ierr); 607 #endif 608 { /* construct prolongator */ 609 Mat Gmat; 610 PetscCoarsenData *agg_lists; 611 Mat Prol11; 612 613 ierr = pc_gamg->ops->graph(pc,Aarr[level], &Gmat);CHKERRQ(ierr); 614 ierr = pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists);CHKERRQ(ierr); 615 ierr = pc_gamg->ops->prolongator(pc,Aarr[level],Gmat,agg_lists,&Prol11);CHKERRQ(ierr); 616 617 /* could have failed to create new level */ 618 if (Prol11) { 619 const char *prefix; 620 char addp[32]; 621 622 /* get new block size of coarse matrices */ 623 ierr = MatGetBlockSizes(Prol11, NULL, &bs);CHKERRQ(ierr); 624 625 if (pc_gamg->ops->optprolongator) { 626 /* smooth */ 627 ierr = pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11);CHKERRQ(ierr); 628 } 629 630 if (pc_gamg->use_aggs_in_asm) { 631 PetscInt bs; 632 ierr = MatGetBlockSizes(Prol11, &bs, NULL);CHKERRQ(ierr); 633 ierr = PetscCDGetASMBlocks(agg_lists, bs, Gmat, &nASMBlocksArr[level], &ASMLocalIDsArr[level]);CHKERRQ(ierr); 634 } 635 636 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 637 ierr = MatSetOptionsPrefix(Prol11,prefix);CHKERRQ(ierr); 638 ierr = PetscSNPrintf(addp,sizeof(addp),"pc_gamg_prolongator_%d_",(int)level);CHKERRQ(ierr); 639 ierr = MatAppendOptionsPrefix(Prol11,addp);CHKERRQ(ierr); 640 /* Always generate the transpose with CUDA 641 Such behaviour can be adapted with -pc_gamg_prolongator_ prefixed options */ 642 ierr = MatSetOption(Prol11,MAT_FORM_EXPLICIT_TRANSPOSE,PETSC_TRUE);CHKERRQ(ierr); 643 ierr = MatSetFromOptions(Prol11);CHKERRQ(ierr); 644 Parr[level1] = Prol11; 645 } else Parr[level1] = NULL; /* failed to coarsen */ 646 647 ierr = MatDestroy(&Gmat);CHKERRQ(ierr); 648 ierr = PetscCDDestroy(agg_lists);CHKERRQ(ierr); 649 } /* construct prolongator scope */ 650 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr); 651 if (!level) Aarr[0] = Pmat; /* use Pmat for finest level setup */ 652 if (!Parr[level1]) { /* failed to coarsen */ 653 ierr = PetscInfo1(pc,"Stop gridding, level %D\n",level);CHKERRQ(ierr); 654 #if defined(GAMG_STAGES) 655 ierr = PetscLogStagePop();CHKERRQ(ierr); 656 #endif 657 break; 658 } 659 ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 660 ierr = MatGetSize(Parr[level1], &M, &N);CHKERRQ(ierr); /* N is next M, a loop test variables */ 661 if (is_last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Is last ?"); 662 if (N <= pc_gamg->coarse_eq_limit) is_last = PETSC_TRUE; 663 if (level1 == pc_gamg->Nlevels-1) is_last = PETSC_TRUE; 664 ierr = pc_gamg->ops->createlevel(pc, Aarr[level], bs, &Parr[level1], &Aarr[level1], &nactivepe, NULL, is_last);CHKERRQ(ierr); 665 666 ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr); 667 ierr = MatGetSize(Aarr[level1], &M, &N);CHKERRQ(ierr); /* M is loop test variables */ 668 ierr = MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info);CHKERRQ(ierr); 669 nnztot += info.nz_used; 670 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); 671 672 #if defined(GAMG_STAGES) 673 ierr = PetscLogStagePop();CHKERRQ(ierr); 674 #endif 675 /* stop if one node or one proc -- could pull back for singular problems */ 676 if ((pc_gamg->data_cell_cols && M/pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M/bs < 2)) { 677 ierr = PetscInfo2(pc,"HARD stop of coarsening on level %D. Grid too small: %D block nodes\n",level,M/bs);CHKERRQ(ierr); 678 level++; 679 break; 680 } 681 } /* levels */ 682 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 683 684 ierr = PetscInfo2(pc,"%D levels, grid complexity = %g\n",level+1,nnztot/nnz0);CHKERRQ(ierr); 685 pc_gamg->Nlevels = level + 1; 686 fine_level = level; 687 ierr = PCMGSetLevels(pc,pc_gamg->Nlevels,NULL);CHKERRQ(ierr); 688 689 if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */ 690 PetscErrorCode (*savesetfromoptions[PETSC_MG_MAXLEVELS])(PetscOptionItems*,KSP); 691 692 /* set default smoothers & set operators */ 693 for (lidx = 1, level = pc_gamg->Nlevels-2; lidx <= fine_level; lidx++, level--) { 694 KSP smoother; 695 PC subpc; 696 697 ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr); 698 ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr); 699 700 ierr = KSPSetNormType(smoother, KSP_NORM_NONE);CHKERRQ(ierr); 701 /* set ops */ 702 ierr = KSPSetOperators(smoother, Aarr[level], Aarr[level]);CHKERRQ(ierr); 703 ierr = PCMGSetInterpolation(pc, lidx, Parr[level+1]);CHKERRQ(ierr); 704 705 /* set defaults */ 706 ierr = KSPSetType(smoother, KSPCHEBYSHEV);CHKERRQ(ierr); 707 708 /* set blocks for ASM smoother that uses the 'aggregates' */ 709 if (pc_gamg->use_aggs_in_asm) { 710 PetscInt sz; 711 IS *iss; 712 713 sz = nASMBlocksArr[level]; 714 iss = ASMLocalIDsArr[level]; 715 ierr = PCSetType(subpc, PCASM);CHKERRQ(ierr); 716 ierr = PCASMSetOverlap(subpc, 0);CHKERRQ(ierr); 717 ierr = PCASMSetType(subpc,PC_ASM_BASIC);CHKERRQ(ierr); 718 if (!sz) { 719 IS is; 720 ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is);CHKERRQ(ierr); 721 ierr = PCASMSetLocalSubdomains(subpc, 1, NULL, &is);CHKERRQ(ierr); 722 ierr = ISDestroy(&is);CHKERRQ(ierr); 723 } else { 724 PetscInt kk; 725 ierr = PCASMSetLocalSubdomains(subpc, sz, NULL, iss);CHKERRQ(ierr); 726 for (kk=0; kk<sz; kk++) { 727 ierr = ISDestroy(&iss[kk]);CHKERRQ(ierr); 728 } 729 ierr = PetscFree(iss);CHKERRQ(ierr); 730 } 731 ASMLocalIDsArr[level] = NULL; 732 nASMBlocksArr[level] = 0; 733 } else { 734 ierr = PCSetType(subpc, PCSOR);CHKERRQ(ierr); 735 } 736 } 737 { 738 /* coarse grid */ 739 KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first; 740 Mat Lmat = Aarr[(level=pc_gamg->Nlevels-1)]; lidx = 0; 741 742 ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr); 743 ierr = KSPSetOperators(smoother, Lmat, Lmat);CHKERRQ(ierr); 744 if (!pc_gamg->use_parallel_coarse_grid_solver) { 745 ierr = KSPSetNormType(smoother, KSP_NORM_NONE);CHKERRQ(ierr); 746 ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr); 747 ierr = PCSetType(subpc, PCBJACOBI);CHKERRQ(ierr); 748 ierr = PCSetUp(subpc);CHKERRQ(ierr); 749 ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr); 750 if (ii != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D is not one",ii); 751 ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr); 752 ierr = PCSetType(pc2, PCLU);CHKERRQ(ierr); 753 ierr = PCFactorSetShiftType(pc2,MAT_SHIFT_INBLOCKS);CHKERRQ(ierr); 754 ierr = KSPSetTolerances(k2[0],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);CHKERRQ(ierr); 755 ierr = KSPSetType(k2[0], KSPPREONLY);CHKERRQ(ierr); 756 /* This flag gets reset by PCBJacobiGetSubKSP(), but our BJacobi really does the same algorithm everywhere (and in 757 * fact, all but one process will have zero dofs), so we reset the flag to avoid having PCView_BJacobi attempt to 758 * view every subdomain as though they were different. */ 759 ((PC_BJacobi*)subpc->data)->same_local_solves = PETSC_TRUE; 760 } 761 } 762 763 /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 764 ierr = PetscObjectOptionsBegin((PetscObject)pc);CHKERRQ(ierr); 765 ierr = PCSetFromOptions_MG(PetscOptionsObject,pc);CHKERRQ(ierr); 766 ierr = PetscOptionsEnd();CHKERRQ(ierr); 767 ierr = PCMGSetGalerkin(pc,PC_MG_GALERKIN_EXTERNAL);CHKERRQ(ierr); 768 769 /* setup cheby eigen estimates from SA */ 770 if (pc_gamg->use_sa_esteig==1) { 771 for (lidx = 1, level = pc_gamg->Nlevels-2; level >= 0 ; lidx++, level--) { 772 KSP smoother; 773 PetscBool ischeb; 774 775 savesetfromoptions[level] = NULL; 776 ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr); 777 ierr = PetscObjectTypeCompare((PetscObject)smoother,KSPCHEBYSHEV,&ischeb);CHKERRQ(ierr); 778 if (ischeb) { 779 KSP_Chebyshev *cheb = (KSP_Chebyshev*)smoother->data; 780 781 ierr = KSPSetFromOptions(smoother);CHKERRQ(ierr); /* let command line emax override using SA's eigenvalues */ 782 if (mg->max_eigen_DinvA[level] > 0 && cheb->emax == 0.) { 783 PC subpc; 784 PetscBool isjac; 785 ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr); 786 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCJACOBI,&isjac);CHKERRQ(ierr); 787 if (isjac && pc_gamg->use_sa_esteig==1) { 788 PetscReal emax,emin; 789 790 emin = mg->min_eigen_DinvA[level]; 791 emax = mg->max_eigen_DinvA[level]; 792 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); 793 cheb->emin_computed = emin; 794 cheb->emax_computed = emax; 795 ierr = KSPChebyshevSetEigenvalues(smoother, cheb->tform[2]*emin + cheb->tform[3]*emax, cheb->tform[0]*emin + cheb->tform[1]*emax);CHKERRQ(ierr); 796 797 /* We have set the eigenvalues and consumed the transformation values 798 prevent from flagging the recomputation of the eigenvalues again in PCSetUp_MG 799 below when setfromoptions will be called again */ 800 savesetfromoptions[level] = smoother->ops->setfromoptions; 801 smoother->ops->setfromoptions = NULL; 802 } 803 } 804 } 805 } 806 } 807 808 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 809 810 /* restore Chebyshev smoother for next calls */ 811 if (pc_gamg->use_sa_esteig==1) { 812 for (lidx = 1, level = pc_gamg->Nlevels-2; level >= 0 ; lidx++, level--) { 813 if (savesetfromoptions[level]) { 814 KSP smoother; 815 ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr); 816 smoother->ops->setfromoptions = savesetfromoptions[level]; 817 } 818 } 819 } 820 821 /* clean up */ 822 for (level=1; level<pc_gamg->Nlevels; level++) { 823 ierr = MatDestroy(&Parr[level]);CHKERRQ(ierr); 824 ierr = MatDestroy(&Aarr[level]);CHKERRQ(ierr); 825 } 826 } else { 827 KSP smoother; 828 829 ierr = PetscInfo(pc,"One level solver used (system is seen as DD). Using default solver.\n");CHKERRQ(ierr); 830 ierr = PCMGGetSmoother(pc, 0, &smoother);CHKERRQ(ierr); 831 ierr = KSPSetOperators(smoother, Aarr[0], Aarr[0]);CHKERRQ(ierr); 832 ierr = KSPSetType(smoother, KSPPREONLY);CHKERRQ(ierr); 833 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 834 } 835 PetscFunctionReturn(0); 836 } 837 838 /* ------------------------------------------------------------------------- */ 839 /* 840 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 841 that was created with PCCreate_GAMG(). 842 843 Input Parameter: 844 . pc - the preconditioner context 845 846 Application Interface Routine: PCDestroy() 847 */ 848 PetscErrorCode PCDestroy_GAMG(PC pc) 849 { 850 PetscErrorCode ierr; 851 PC_MG *mg = (PC_MG*)pc->data; 852 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 853 854 PetscFunctionBegin; 855 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 856 if (pc_gamg->ops->destroy) { 857 ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr); 858 } 859 ierr = PetscFree(pc_gamg->ops);CHKERRQ(ierr); 860 ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr); 861 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 862 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 863 PetscFunctionReturn(0); 864 } 865 866 /*@ 867 PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction. 868 869 Logically Collective on PC 870 871 Input Parameters: 872 + pc - the preconditioner context 873 - n - the number of equations 874 875 876 Options Database Key: 877 . -pc_gamg_process_eq_limit <limit> 878 879 Notes: 880 GAMG will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process 881 that has degrees of freedom 882 883 Level: intermediate 884 885 .seealso: PCGAMGSetCoarseEqLim(), PCGAMGSetRankReductionFactors() 886 @*/ 887 PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n) 888 { 889 PetscErrorCode ierr; 890 891 PetscFunctionBegin; 892 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 893 ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 894 PetscFunctionReturn(0); 895 } 896 897 static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n) 898 { 899 PC_MG *mg = (PC_MG*)pc->data; 900 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 901 902 PetscFunctionBegin; 903 if (n>0) pc_gamg->min_eq_proc = n; 904 PetscFunctionReturn(0); 905 } 906 907 /*@ 908 PCGAMGSetCoarseEqLim - Set maximum number of equations on coarsest grid. 909 910 Collective on PC 911 912 Input Parameters: 913 + pc - the preconditioner context 914 - n - maximum number of equations to aim for 915 916 Options Database Key: 917 . -pc_gamg_coarse_eq_limit <limit> 918 919 Notes: For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid 920 has less than 1000 unknowns. 921 922 Level: intermediate 923 924 .seealso: PCGAMGSetProcEqLim(), PCGAMGSetRankReductionFactors() 925 @*/ 926 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n) 927 { 928 PetscErrorCode ierr; 929 930 PetscFunctionBegin; 931 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 932 ierr = PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 933 PetscFunctionReturn(0); 934 } 935 936 static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n) 937 { 938 PC_MG *mg = (PC_MG*)pc->data; 939 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 940 941 PetscFunctionBegin; 942 if (n>0) pc_gamg->coarse_eq_limit = n; 943 PetscFunctionReturn(0); 944 } 945 946 /*@ 947 PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids 948 949 Collective on PC 950 951 Input Parameters: 952 + pc - the preconditioner context 953 - n - PETSC_TRUE or PETSC_FALSE 954 955 Options Database Key: 956 . -pc_gamg_repartition <true,false> 957 958 Notes: 959 this will generally improve the loading balancing of the work on each level 960 961 Level: intermediate 962 963 .seealso: () 964 @*/ 965 PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n) 966 { 967 PetscErrorCode ierr; 968 969 PetscFunctionBegin; 970 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 971 ierr = PetscTryMethod(pc,"PCGAMGSetRepartition_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 972 PetscFunctionReturn(0); 973 } 974 975 static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n) 976 { 977 PC_MG *mg = (PC_MG*)pc->data; 978 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 979 980 PetscFunctionBegin; 981 pc_gamg->repart = n; 982 PetscFunctionReturn(0); 983 } 984 985 /*@ 986 PCGAMGSetEstEigKSPMaxIt - Set number of KSP iterations in eigen estimator (for Cheby) 987 988 Collective on PC 989 990 Input Parameters: 991 + pc - the preconditioner context 992 - n - number of its 993 994 Options Database Key: 995 . -pc_gamg_esteig_ksp_max_it <its> 996 997 Notes: 998 999 Level: intermediate 1000 1001 .seealso: () 1002 @*/ 1003 PetscErrorCode PCGAMGSetEstEigKSPMaxIt(PC pc, PetscInt n) 1004 { 1005 PetscErrorCode ierr; 1006 1007 PetscFunctionBegin; 1008 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1009 ierr = PetscTryMethod(pc,"PCGAMGSetEstEigKSPMaxIt_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 1010 PetscFunctionReturn(0); 1011 } 1012 1013 static PetscErrorCode PCGAMGSetEstEigKSPMaxIt_GAMG(PC pc, PetscInt n) 1014 { 1015 PC_MG *mg = (PC_MG*)pc->data; 1016 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1017 1018 PetscFunctionBegin; 1019 pc_gamg->esteig_max_it = n; 1020 PetscFunctionReturn(0); 1021 } 1022 1023 /*@ 1024 PCGAMGSetUseSAEstEig - Use eigen estimate from smoothed aggregation for Cheby smoother 1025 1026 Collective on PC 1027 1028 Input Parameters: 1029 + pc - the preconditioner context 1030 - n - number of its 1031 1032 Options Database Key: 1033 . -pc_gamg_use_sa_esteig <true,false> 1034 1035 Notes: 1036 1037 Level: intermediate 1038 1039 .seealso: () 1040 @*/ 1041 PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool n) 1042 { 1043 PetscErrorCode ierr; 1044 1045 PetscFunctionBegin; 1046 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1047 ierr = PetscTryMethod(pc,"PCGAMGSetUseSAEstEig_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 1048 PetscFunctionReturn(0); 1049 } 1050 1051 static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool n) 1052 { 1053 PC_MG *mg = (PC_MG*)pc->data; 1054 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1055 1056 PetscFunctionBegin; 1057 pc_gamg->use_sa_esteig = n ? 1 : 0; 1058 PetscFunctionReturn(0); 1059 } 1060 1061 /*@C 1062 PCGAMGSetEstEigKSPType - Set type of KSP in eigen estimator (for Cheby) 1063 1064 Collective on PC 1065 1066 Input Parameters: 1067 + pc - the preconditioner context 1068 - t - ksp type 1069 1070 Options Database Key: 1071 . -pc_gamg_esteig_ksp_type <type> 1072 1073 Notes: 1074 1075 Level: intermediate 1076 1077 .seealso: () 1078 @*/ 1079 PetscErrorCode PCGAMGSetEstEigKSPType(PC pc, char t[]) 1080 { 1081 PetscErrorCode ierr; 1082 1083 PetscFunctionBegin; 1084 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1085 ierr = PetscTryMethod(pc,"PCGAMGSetEstEigKSPType_C",(PC,char[]),(pc,t));CHKERRQ(ierr); 1086 PetscFunctionReturn(0); 1087 } 1088 1089 static PetscErrorCode PCGAMGSetEstEigKSPType_GAMG(PC pc, char t[]) 1090 { 1091 PetscErrorCode ierr; 1092 PC_MG *mg = (PC_MG*)pc->data; 1093 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1094 1095 PetscFunctionBegin; 1096 ierr = PetscStrcpy(pc_gamg->esteig_type,t);CHKERRQ(ierr); 1097 PetscFunctionReturn(0); 1098 } 1099 1100 /*@ 1101 PCGAMGSetEigenvalues - Set eigenvalues 1102 1103 Collective on PC 1104 1105 Input Parameters: 1106 + pc - the preconditioner context 1107 - emax - max eigenvalue 1108 . emin - min eigenvalue 1109 1110 Options Database Key: 1111 . -pc_gamg_eigenvalues 1112 1113 Level: intermediate 1114 1115 .seealso: PCGAMGSetEstEigKSPMaxIt(), PCGAMGSetUseSAEstEig(), PCGAMGSetEstEigKSPType() 1116 @*/ 1117 PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax,PetscReal emin) 1118 { 1119 PetscErrorCode ierr; 1120 1121 PetscFunctionBegin; 1122 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1123 ierr = PetscTryMethod(pc,"PCGAMGSetEigenvalues_C",(PC,PetscReal,PetscReal),(pc,emax,emin));CHKERRQ(ierr); 1124 PetscFunctionReturn(0); 1125 } 1126 static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc,PetscReal emax,PetscReal emin) 1127 { 1128 PC_MG *mg = (PC_MG*)pc->data; 1129 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1130 1131 PetscFunctionBegin; 1132 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); 1133 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); 1134 pc_gamg->emax = emax; 1135 pc_gamg->emin = emin; 1136 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /*@ 1141 PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding algebraic multigrid preconditioner 1142 1143 Collective on PC 1144 1145 Input Parameters: 1146 + pc - the preconditioner context 1147 - n - PETSC_TRUE or PETSC_FALSE 1148 1149 Options Database Key: 1150 . -pc_gamg_reuse_interpolation <true,false> 1151 1152 Level: intermediate 1153 1154 Notes: 1155 this may negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows 1156 rebuilding the preconditioner quicker. 1157 1158 .seealso: () 1159 @*/ 1160 PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n) 1161 { 1162 PetscErrorCode ierr; 1163 1164 PetscFunctionBegin; 1165 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1166 ierr = PetscTryMethod(pc,"PCGAMGSetReuseInterpolation_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr); 1167 PetscFunctionReturn(0); 1168 } 1169 1170 static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n) 1171 { 1172 PC_MG *mg = (PC_MG*)pc->data; 1173 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1174 1175 PetscFunctionBegin; 1176 pc_gamg->reuse_prol = n; 1177 PetscFunctionReturn(0); 1178 } 1179 1180 /*@ 1181 PCGAMGASMSetUseAggs - Have the PCGAMG smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner. 1182 1183 Collective on PC 1184 1185 Input Parameters: 1186 + pc - the preconditioner context 1187 - flg - PETSC_TRUE to use aggregates, PETSC_FALSE to not 1188 1189 Options Database Key: 1190 . -pc_gamg_asm_use_agg 1191 1192 Level: intermediate 1193 1194 .seealso: () 1195 @*/ 1196 PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg) 1197 { 1198 PetscErrorCode ierr; 1199 1200 PetscFunctionBegin; 1201 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1202 ierr = PetscTryMethod(pc,"PCGAMGASMSetUseAggs_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 1203 PetscFunctionReturn(0); 1204 } 1205 1206 static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg) 1207 { 1208 PC_MG *mg = (PC_MG*)pc->data; 1209 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1210 1211 PetscFunctionBegin; 1212 pc_gamg->use_aggs_in_asm = flg; 1213 PetscFunctionReturn(0); 1214 } 1215 1216 /*@ 1217 PCGAMGSetUseParallelCoarseGridSolve - allow a parallel coarse grid solver 1218 1219 Collective on PC 1220 1221 Input Parameters: 1222 + pc - the preconditioner context 1223 - flg - PETSC_TRUE to not force coarse grid onto one processor 1224 1225 Options Database Key: 1226 . -pc_gamg_use_parallel_coarse_grid_solver 1227 1228 Level: intermediate 1229 1230 .seealso: PCGAMGSetCoarseGridLayoutType(), PCGAMGSetCpuPinCoarseGrids() 1231 @*/ 1232 PetscErrorCode PCGAMGSetUseParallelCoarseGridSolve(PC pc, PetscBool flg) 1233 { 1234 PetscErrorCode ierr; 1235 1236 PetscFunctionBegin; 1237 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1238 ierr = PetscTryMethod(pc,"PCGAMGSetUseParallelCoarseGridSolve_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 1239 PetscFunctionReturn(0); 1240 } 1241 1242 static PetscErrorCode PCGAMGSetUseParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg) 1243 { 1244 PC_MG *mg = (PC_MG*)pc->data; 1245 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1246 1247 PetscFunctionBegin; 1248 pc_gamg->use_parallel_coarse_grid_solver = flg; 1249 PetscFunctionReturn(0); 1250 } 1251 1252 /*@ 1253 PCGAMGSetCpuPinCoarseGrids - pin reduced grids to CPU 1254 1255 Collective on PC 1256 1257 Input Parameters: 1258 + pc - the preconditioner context 1259 - flg - PETSC_TRUE to pin coarse grids to CPU 1260 1261 Options Database Key: 1262 . -pc_gamg_cpu_pin_coarse_grids 1263 1264 Level: intermediate 1265 1266 .seealso: PCGAMGSetCoarseGridLayoutType(), PCGAMGSetUseParallelCoarseGridSolve() 1267 @*/ 1268 PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg) 1269 { 1270 PetscErrorCode ierr; 1271 1272 PetscFunctionBegin; 1273 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1274 ierr = PetscTryMethod(pc,"PCGAMGSetCpuPinCoarseGrids_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 1275 PetscFunctionReturn(0); 1276 } 1277 1278 static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg) 1279 { 1280 PC_MG *mg = (PC_MG*)pc->data; 1281 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1282 1283 PetscFunctionBegin; 1284 pc_gamg->cpu_pin_coarse_grids = flg; 1285 PetscFunctionReturn(0); 1286 } 1287 1288 /*@ 1289 PCGAMGSetCoarseGridLayoutType - place reduce grids on processors with natural order (compact type) 1290 1291 Collective on PC 1292 1293 Input Parameters: 1294 + pc - the preconditioner context 1295 - flg - Layout type 1296 1297 Options Database Key: 1298 . -pc_gamg_coarse_grid_layout_type 1299 1300 Level: intermediate 1301 1302 .seealso: PCGAMGSetUseParallelCoarseGridSolve(), PCGAMGSetCpuPinCoarseGrids() 1303 @*/ 1304 PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg) 1305 { 1306 PetscErrorCode ierr; 1307 1308 PetscFunctionBegin; 1309 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1310 ierr = PetscTryMethod(pc,"PCGAMGSetCoarseGridLayoutType_C",(PC,PCGAMGLayoutType),(pc,flg));CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 1314 static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg) 1315 { 1316 PC_MG *mg = (PC_MG*)pc->data; 1317 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1318 1319 PetscFunctionBegin; 1320 pc_gamg->layout_type = flg; 1321 PetscFunctionReturn(0); 1322 } 1323 1324 /*@ 1325 PCGAMGSetNlevels - Sets the maximum number of levels PCGAMG will use 1326 1327 Not collective on PC 1328 1329 Input Parameters: 1330 + pc - the preconditioner 1331 - n - the maximum number of levels to use 1332 1333 Options Database Key: 1334 . -pc_mg_levels 1335 1336 Level: intermediate 1337 1338 .seealso: () 1339 @*/ 1340 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n) 1341 { 1342 PetscErrorCode ierr; 1343 1344 PetscFunctionBegin; 1345 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1346 ierr = PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr); 1347 PetscFunctionReturn(0); 1348 } 1349 1350 static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n) 1351 { 1352 PC_MG *mg = (PC_MG*)pc->data; 1353 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1354 1355 PetscFunctionBegin; 1356 pc_gamg->Nlevels = n; 1357 PetscFunctionReturn(0); 1358 } 1359 1360 /*@ 1361 PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph 1362 1363 Not collective on PC 1364 1365 Input Parameters: 1366 + pc - the preconditioner context 1367 . 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 1368 - n - number of threshold values provided in array 1369 1370 Options Database Key: 1371 . -pc_gamg_threshold <threshold> 1372 1373 Notes: 1374 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. 1375 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. 1376 1377 If n is less than the total number of coarsenings (see PCGAMGSetNlevels()), then threshold scaling (see PCGAMGSetThresholdScale()) is used for each successive coarsening. 1378 In this case, PCGAMGSetThresholdScale() must be called before PCGAMGSetThreshold(). 1379 If n is greater than the total number of levels, the excess entries in threshold will not be used. 1380 1381 Level: intermediate 1382 1383 .seealso: PCGAMGFilterGraph(), PCGAMGSetSquareGraph() 1384 @*/ 1385 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n) 1386 { 1387 PetscErrorCode ierr; 1388 1389 PetscFunctionBegin; 1390 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1391 if (n) PetscValidRealPointer(v,2); 1392 ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal[],PetscInt),(pc,v,n));CHKERRQ(ierr); 1393 PetscFunctionReturn(0); 1394 } 1395 1396 static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n) 1397 { 1398 PC_MG *mg = (PC_MG*)pc->data; 1399 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1400 PetscInt i; 1401 PetscFunctionBegin; 1402 for (i=0; i<PetscMin(n,PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i]; 1403 for (; i<PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i-1]*pc_gamg->threshold_scale; 1404 PetscFunctionReturn(0); 1405 } 1406 1407 /*@ 1408 PCGAMGSetRankReductionFactors - Set manual schedual for process reduction on coarse grids 1409 1410 Collective on PC 1411 1412 Input Parameters: 1413 + pc - the preconditioner context 1414 . v - array of reduction factors. 0 for fist value forces a reduction to one process/device on first level in Cuda 1415 - n - number of values provided in array 1416 1417 Options Database Key: 1418 . -pc_gamg_rank_reduction_factors <factors> 1419 1420 Level: intermediate 1421 1422 .seealso: PCGAMGSetProcEqLim(), PCGAMGSetCoarseEqLim() 1423 @*/ 1424 PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n) 1425 { 1426 PetscErrorCode ierr; 1427 1428 PetscFunctionBegin; 1429 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1430 if (n) PetscValidIntPointer(v,2); 1431 ierr = PetscTryMethod(pc,"PCGAMGSetRankReductionFactors_C",(PC,PetscInt[],PetscInt),(pc,v,n));CHKERRQ(ierr); 1432 PetscFunctionReturn(0); 1433 } 1434 1435 static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n) 1436 { 1437 PC_MG *mg = (PC_MG*)pc->data; 1438 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1439 PetscInt i; 1440 PetscFunctionBegin; 1441 for (i=0; i<PetscMin(n,PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i]; 1442 for (; i<PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */ 1443 PetscFunctionReturn(0); 1444 } 1445 1446 /*@ 1447 PCGAMGSetThresholdScale - Relative threshold reduction at each level 1448 1449 Not collective on PC 1450 1451 Input Parameters: 1452 + pc - the preconditioner context 1453 - scale - the threshold value reduction, ussually < 1.0 1454 1455 Options Database Key: 1456 . -pc_gamg_threshold_scale <v> 1457 1458 Notes: 1459 The initial threshold (for an arbitrary number of levels starting from the finest) can be set with PCGAMGSetThreshold(). 1460 This scaling is used for each subsequent coarsening, but must be called before PCGAMGSetThreshold(). 1461 1462 Level: advanced 1463 1464 .seealso: PCGAMGSetThreshold() 1465 @*/ 1466 PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v) 1467 { 1468 PetscErrorCode ierr; 1469 1470 PetscFunctionBegin; 1471 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1472 ierr = PetscTryMethod(pc,"PCGAMGSetThresholdScale_C",(PC,PetscReal),(pc,v));CHKERRQ(ierr); 1473 PetscFunctionReturn(0); 1474 } 1475 1476 static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v) 1477 { 1478 PC_MG *mg = (PC_MG*)pc->data; 1479 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1480 PetscFunctionBegin; 1481 pc_gamg->threshold_scale = v; 1482 PetscFunctionReturn(0); 1483 } 1484 1485 /*@C 1486 PCGAMGSetType - Set solution method 1487 1488 Collective on PC 1489 1490 Input Parameters: 1491 + pc - the preconditioner context 1492 - type - PCGAMGAGG, PCGAMGGEO, or PCGAMGCLASSICAL 1493 1494 Options Database Key: 1495 . -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply 1496 1497 Level: intermediate 1498 1499 .seealso: PCGAMGGetType(), PCGAMG, PCGAMGType 1500 @*/ 1501 PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type) 1502 { 1503 PetscErrorCode ierr; 1504 1505 PetscFunctionBegin; 1506 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1507 ierr = PetscTryMethod(pc,"PCGAMGSetType_C",(PC,PCGAMGType),(pc,type));CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 1511 /*@C 1512 PCGAMGGetType - Get solution method 1513 1514 Collective on PC 1515 1516 Input Parameter: 1517 . pc - the preconditioner context 1518 1519 Output Parameter: 1520 . type - the type of algorithm used 1521 1522 Level: intermediate 1523 1524 .seealso: PCGAMGSetType(), PCGAMGType 1525 @*/ 1526 PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type) 1527 { 1528 PetscErrorCode ierr; 1529 1530 PetscFunctionBegin; 1531 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1532 ierr = PetscUseMethod(pc,"PCGAMGGetType_C",(PC,PCGAMGType*),(pc,type));CHKERRQ(ierr); 1533 PetscFunctionReturn(0); 1534 } 1535 1536 static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type) 1537 { 1538 PC_MG *mg = (PC_MG*)pc->data; 1539 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1540 1541 PetscFunctionBegin; 1542 *type = pc_gamg->type; 1543 PetscFunctionReturn(0); 1544 } 1545 1546 static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type) 1547 { 1548 PetscErrorCode ierr,(*r)(PC); 1549 PC_MG *mg = (PC_MG*)pc->data; 1550 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1551 1552 PetscFunctionBegin; 1553 pc_gamg->type = type; 1554 ierr = PetscFunctionListFind(GAMGList,type,&r);CHKERRQ(ierr); 1555 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown GAMG type %s given",type); 1556 if (pc_gamg->ops->destroy) { 1557 ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr); 1558 ierr = PetscMemzero(pc_gamg->ops,sizeof(struct _PCGAMGOps));CHKERRQ(ierr); 1559 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1560 /* cleaning up common data in pc_gamg - this should disapear someday */ 1561 pc_gamg->data_cell_cols = 0; 1562 pc_gamg->data_cell_rows = 0; 1563 pc_gamg->orig_data_cell_cols = 0; 1564 pc_gamg->orig_data_cell_rows = 0; 1565 ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr); 1566 pc_gamg->data_sz = 0; 1567 } 1568 ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr); 1569 ierr = PetscStrallocpy(type,&pc_gamg->gamg_type_name);CHKERRQ(ierr); 1570 ierr = (*r)(pc);CHKERRQ(ierr); 1571 PetscFunctionReturn(0); 1572 } 1573 1574 /* -------------------------------------------------------------------------- */ 1575 /* 1576 PCMGGetGridComplexity - compute coarse grid complexity of MG hierarchy 1577 1578 Input Parameter: 1579 . pc - the preconditioner context 1580 1581 Output Parameter: 1582 . gc - grid complexity = sum_i(nnz_i) / nnz_0 1583 1584 Level: advanced 1585 */ 1586 static PetscErrorCode PCMGGetGridComplexity(PC pc, PetscReal *gc) 1587 { 1588 PetscErrorCode ierr; 1589 PC_MG *mg = (PC_MG*)pc->data; 1590 PC_MG_Levels **mglevels = mg->levels; 1591 PetscInt lev; 1592 PetscLogDouble nnz0 = 0, sgc = 0; 1593 MatInfo info; 1594 1595 PetscFunctionBegin; 1596 if (!pc->setupcalled) { 1597 *gc = 0; 1598 PetscFunctionReturn(0); 1599 } 1600 if (!mg->nlevels) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MG has no levels"); 1601 for (lev=0; lev<mg->nlevels; lev++) { 1602 Mat dB; 1603 ierr = KSPGetOperators(mglevels[lev]->smoothd,NULL,&dB);CHKERRQ(ierr); 1604 ierr = MatGetInfo(dB,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); /* global reduction */ 1605 sgc += info.nz_used; 1606 if (lev==mg->nlevels-1) nnz0 = info.nz_used; 1607 } 1608 if (nnz0 > 0) *gc = (PetscReal)(sgc/nnz0); 1609 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number for grid points on finest level is not available"); 1610 PetscFunctionReturn(0); 1611 } 1612 1613 static PetscErrorCode PCView_GAMG(PC pc,PetscViewer viewer) 1614 { 1615 PetscErrorCode ierr,i; 1616 PC_MG *mg = (PC_MG*)pc->data; 1617 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1618 PetscReal gc=0; 1619 PetscFunctionBegin; 1620 ierr = PetscViewerASCIIPrintf(viewer," GAMG specific options\n");CHKERRQ(ierr); 1621 ierr = PetscViewerASCIIPrintf(viewer," Threshold for dropping small values in graph on each level =");CHKERRQ(ierr); 1622 for (i=0;i<mg->nlevels; i++) { 1623 ierr = PetscViewerASCIIPrintf(viewer," %g",(double)pc_gamg->threshold[i]);CHKERRQ(ierr); 1624 } 1625 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1626 ierr = PetscViewerASCIIPrintf(viewer," Threshold scaling factor for each level not specified = %g\n",(double)pc_gamg->threshold_scale);CHKERRQ(ierr); 1627 if (pc_gamg->use_aggs_in_asm) { 1628 ierr = PetscViewerASCIIPrintf(viewer," Using aggregates from coarsening process to define subdomains for PCASM\n");CHKERRQ(ierr); 1629 } 1630 if (pc_gamg->use_parallel_coarse_grid_solver) { 1631 ierr = PetscViewerASCIIPrintf(viewer," Using parallel coarse grid solver (all coarse grid equations not put on one process)\n");CHKERRQ(ierr); 1632 } 1633 #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) 1634 if (pc_gamg->cpu_pin_coarse_grids) { 1635 /* ierr = PetscViewerASCIIPrintf(viewer," Pinning coarse grids to the CPU)\n");CHKERRQ(ierr); */ 1636 } 1637 #endif 1638 /* if (pc_gamg->layout_type==PCGAMG_LAYOUT_COMPACT) { */ 1639 /* ierr = PetscViewerASCIIPrintf(viewer," Put reduced grids on processes in natural order (ie, 0,1,2...)\n");CHKERRQ(ierr); */ 1640 /* } else { */ 1641 /* ierr = PetscViewerASCIIPrintf(viewer," Put reduced grids on whole machine (ie, 0,1*f,2*f...,np-f)\n");CHKERRQ(ierr); */ 1642 /* } */ 1643 if (pc_gamg->ops->view) { 1644 ierr = (*pc_gamg->ops->view)(pc,viewer);CHKERRQ(ierr); 1645 } 1646 ierr = PCMGGetGridComplexity(pc,&gc);CHKERRQ(ierr); 1647 ierr = PetscViewerASCIIPrintf(viewer," Complexity: grid = %g\n",gc);CHKERRQ(ierr); 1648 PetscFunctionReturn(0); 1649 } 1650 1651 PetscErrorCode PCSetFromOptions_GAMG(PetscOptionItems *PetscOptionsObject,PC pc) 1652 { 1653 PetscErrorCode ierr; 1654 PC_MG *mg = (PC_MG*)pc->data; 1655 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 1656 PetscBool flag,f2; 1657 MPI_Comm comm; 1658 char prefix[256],tname[32]; 1659 PetscInt i,n; 1660 const char *pcpre; 1661 static const char *LayoutTypes[] = {"compact","spread","PCGAMGLayoutType","PC_GAMG_LAYOUT",NULL}; 1662 PetscFunctionBegin; 1663 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1664 ierr = PetscOptionsHead(PetscOptionsObject,"GAMG options");CHKERRQ(ierr); 1665 ierr = PetscOptionsFList("-pc_gamg_type","Type of AMG method","PCGAMGSetType",GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag);CHKERRQ(ierr); 1666 if (flag) { 1667 ierr = PCGAMGSetType(pc,tname);CHKERRQ(ierr); 1668 } 1669 ierr = PetscOptionsFList("-pc_gamg_esteig_ksp_type","Krylov method for eigen estimator","PCGAMGSetEstEigKSPType",KSPList,pc_gamg->esteig_type,tname,sizeof(tname),&flag);CHKERRQ(ierr); 1670 if (flag) { 1671 ierr = PCGAMGSetEstEigKSPType(pc,tname);CHKERRQ(ierr); 1672 } 1673 ierr = PetscOptionsBool("-pc_gamg_repartition","Repartion coarse grids","PCGAMGSetRepartition",pc_gamg->repart,&pc_gamg->repart,NULL);CHKERRQ(ierr); 1674 f2 = PETSC_TRUE; 1675 ierr = PetscOptionsBool("-pc_gamg_use_sa_esteig","Use eigen estimate from Smoothed aggregation for smoother","PCGAMGSetUseSAEstEig",f2,&f2,&flag);CHKERRQ(ierr); 1676 if (flag) pc_gamg->use_sa_esteig = f2 ? 1 : 0; 1677 ierr = PetscOptionsBool("-pc_gamg_reuse_interpolation","Reuse prolongation operator","PCGAMGReuseInterpolation",pc_gamg->reuse_prol,&pc_gamg->reuse_prol,NULL);CHKERRQ(ierr); 1678 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); 1679 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); 1680 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); 1681 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); 1682 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); 1683 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); 1684 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); 1685 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); 1686 n = PETSC_MG_MAXLEVELS; 1687 ierr = PetscOptionsRealArray("-pc_gamg_threshold","Relative threshold to use for dropping edges in aggregation graph","PCGAMGSetThreshold",pc_gamg->threshold,&n,&flag);CHKERRQ(ierr); 1688 if (!flag || n < PETSC_MG_MAXLEVELS) { 1689 if (!flag) n = 1; 1690 i = n; 1691 do {pc_gamg->threshold[i] = pc_gamg->threshold[i-1]*pc_gamg->threshold_scale;} while (++i<PETSC_MG_MAXLEVELS); 1692 } 1693 n = PETSC_MG_MAXLEVELS; 1694 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); 1695 if (!flag) i = 0; 1696 else i = n; 1697 do {pc_gamg->level_reduction_factors[i] = -1;} while (++i<PETSC_MG_MAXLEVELS); 1698 ierr = PetscOptionsInt("-pc_mg_levels","Set number of MG levels","PCGAMGSetNlevels",pc_gamg->Nlevels,&pc_gamg->Nlevels,NULL);CHKERRQ(ierr); 1699 { 1700 PetscReal eminmax[2] = {0., 0.}; 1701 n = 2; 1702 ierr = PetscOptionsRealArray("-pc_gamg_eigenvalues","extreme eigenvalues for smoothed aggregation","PCGAMGSetEigenvalues",eminmax,&n,&flag);CHKERRQ(ierr); 1703 if (flag) { 1704 if (n != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues"); 1705 ierr = PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0]);CHKERRQ(ierr); 1706 } 1707 } 1708 /* set options for subtype */ 1709 if (pc_gamg->ops->setfromoptions) {ierr = (*pc_gamg->ops->setfromoptions)(PetscOptionsObject,pc);CHKERRQ(ierr);} 1710 1711 ierr = PCGetOptionsPrefix(pc, &pcpre);CHKERRQ(ierr); 1712 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%spc_gamg_",pcpre ? pcpre : "");CHKERRQ(ierr); 1713 ierr = PetscOptionsTail();CHKERRQ(ierr); 1714 PetscFunctionReturn(0); 1715 } 1716 1717 /* -------------------------------------------------------------------------- */ 1718 /*MC 1719 PCGAMG - Geometric algebraic multigrid (AMG) preconditioner 1720 1721 Options Database Keys: 1722 + -pc_gamg_type <type> - one of agg, geo, or classical 1723 . -pc_gamg_repartition <true,default=false> - repartition the degrees of freedom accross the coarse grids as they are determined 1724 . -pc_gamg_reuse_interpolation <true,default=false> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations 1725 . -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 1726 . -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> 1727 equations on each process that has degrees of freedom 1728 . -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for. 1729 . -pc_gamg_threshold[] <thresh,default=0> - Before aggregating the graph GAMG will remove small values from the graph on each level 1730 - -pc_gamg_threshold_scale <scale,default=1> - Scaling of threshold on each coarser grid if not specified 1731 1732 Options Database Keys for default Aggregation: 1733 + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation 1734 . -pc_gamg_sym_graph <true,default=false> - symmetrize the graph before computing the aggregation 1735 - -pc_gamg_square_graph <n,default=1> - number of levels to square the graph before aggregating it 1736 1737 Multigrid options: 1738 + -pc_mg_cycles <v> - v or w, see PCMGSetCycleType() 1739 . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see PCMGSetDistinctSmoothUp() 1740 . -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade 1741 - -pc_mg_levels <levels> - Number of levels of multigrid to use. 1742 1743 1744 Notes: 1745 In order to obtain good performance for PCGAMG for vector valued problems you must 1746 Call MatSetBlockSize() to indicate the number of degrees of freedom per grid point 1747 Call MatSetNearNullSpace() (or PCSetCoordinates() if solving the equations of elasticity) to indicate the near null space of the operator 1748 See the Users Manual Chapter 4 for more details 1749 1750 Level: intermediate 1751 1752 .seealso: PCCreate(), PCSetType(), MatSetBlockSize(), PCMGType, PCSetCoordinates(), MatSetNearNullSpace(), PCGAMGSetType(), PCGAMGAGG, PCGAMGGEO, PCGAMGCLASSICAL, PCGAMGSetProcEqLim(), 1753 PCGAMGSetCoarseEqLim(), PCGAMGSetRepartition(), PCGAMGRegister(), PCGAMGSetReuseInterpolation(), PCGAMGASMSetUseAggs(), PCGAMGSetUseParallelCoarseGridSolve(), PCGAMGSetNlevels(), PCGAMGSetThreshold(), PCGAMGGetType(), PCGAMGSetReuseInterpolation(), PCGAMGSetUseSAEstEig(), PCGAMGSetEstEigKSPMaxIt(), PCGAMGSetEstEigKSPType() 1754 M*/ 1755 1756 PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc) 1757 { 1758 PetscErrorCode ierr,i; 1759 PC_GAMG *pc_gamg; 1760 PC_MG *mg; 1761 1762 PetscFunctionBegin; 1763 /* register AMG type */ 1764 ierr = PCGAMGInitializePackage();CHKERRQ(ierr); 1765 1766 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 1767 ierr = PCSetType(pc, PCMG);CHKERRQ(ierr); 1768 ierr = PetscObjectChangeTypeName((PetscObject)pc, PCGAMG);CHKERRQ(ierr); 1769 1770 /* create a supporting struct and attach it to pc */ 1771 ierr = PetscNewLog(pc,&pc_gamg);CHKERRQ(ierr); 1772 ierr = PCMGSetGalerkin(pc,PC_MG_GALERKIN_EXTERNAL);CHKERRQ(ierr); 1773 mg = (PC_MG*)pc->data; 1774 mg->innerctx = pc_gamg; 1775 1776 ierr = PetscNewLog(pc,&pc_gamg->ops);CHKERRQ(ierr); 1777 1778 /* these should be in subctx but repartitioning needs simple arrays */ 1779 pc_gamg->data_sz = 0; 1780 pc_gamg->data = NULL; 1781 1782 /* overwrite the pointers of PCMG by the functions of base class PCGAMG */ 1783 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 1784 pc->ops->setup = PCSetUp_GAMG; 1785 pc->ops->reset = PCReset_GAMG; 1786 pc->ops->destroy = PCDestroy_GAMG; 1787 mg->view = PCView_GAMG; 1788 1789 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG);CHKERRQ(ierr); 1790 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG);CHKERRQ(ierr); 1791 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetProcEqLim_C",PCGAMGSetProcEqLim_GAMG);CHKERRQ(ierr); 1792 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseEqLim_C",PCGAMGSetCoarseEqLim_GAMG);CHKERRQ(ierr); 1793 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRepartition_C",PCGAMGSetRepartition_GAMG);CHKERRQ(ierr); 1794 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetEstEigKSPType_C",PCGAMGSetEstEigKSPType_GAMG);CHKERRQ(ierr); 1795 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetEstEigKSPMaxIt_C",PCGAMGSetEstEigKSPMaxIt_GAMG);CHKERRQ(ierr); 1796 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetEigenvalues_C",PCGAMGSetEigenvalues_GAMG);CHKERRQ(ierr); 1797 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetUseSAEstEig_C",PCGAMGSetUseSAEstEig_GAMG);CHKERRQ(ierr); 1798 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetReuseInterpolation_C",PCGAMGSetReuseInterpolation_GAMG);CHKERRQ(ierr); 1799 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGASMSetUseAggs_C",PCGAMGASMSetUseAggs_GAMG);CHKERRQ(ierr); 1800 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetUseParallelCoarseGridSolve_C",PCGAMGSetUseParallelCoarseGridSolve_GAMG);CHKERRQ(ierr); 1801 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCpuPinCoarseGrids_C",PCGAMGSetCpuPinCoarseGrids_GAMG);CHKERRQ(ierr); 1802 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseGridLayoutType_C",PCGAMGSetCoarseGridLayoutType_GAMG);CHKERRQ(ierr); 1803 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThreshold_C",PCGAMGSetThreshold_GAMG);CHKERRQ(ierr); 1804 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRankReductionFactors_C",PCGAMGSetRankReductionFactors_GAMG);CHKERRQ(ierr); 1805 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThresholdScale_C",PCGAMGSetThresholdScale_GAMG);CHKERRQ(ierr); 1806 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetType_C",PCGAMGSetType_GAMG);CHKERRQ(ierr); 1807 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGGetType_C",PCGAMGGetType_GAMG);CHKERRQ(ierr); 1808 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNlevels_C",PCGAMGSetNlevels_GAMG);CHKERRQ(ierr); 1809 pc_gamg->repart = PETSC_FALSE; 1810 pc_gamg->reuse_prol = PETSC_FALSE; 1811 pc_gamg->use_aggs_in_asm = PETSC_FALSE; 1812 pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE; 1813 pc_gamg->cpu_pin_coarse_grids = PETSC_FALSE; 1814 pc_gamg->layout_type = PCGAMG_LAYOUT_SPREAD; 1815 pc_gamg->min_eq_proc = 50; 1816 pc_gamg->coarse_eq_limit = 50; 1817 for (i=0;i<PETSC_MG_MAXLEVELS;i++) pc_gamg->threshold[i] = 0.; 1818 pc_gamg->threshold_scale = 1.; 1819 pc_gamg->Nlevels = PETSC_MG_MAXLEVELS; 1820 pc_gamg->current_level = 0; /* don't need to init really */ 1821 ierr = PetscStrcpy(pc_gamg->esteig_type,NULL);CHKERRQ(ierr); 1822 pc_gamg->esteig_max_it = 10; 1823 pc_gamg->use_sa_esteig = -1; 1824 pc_gamg->emin = 0; 1825 pc_gamg->emax = 0; 1826 1827 pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG; 1828 1829 /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */ 1830 ierr = PCGAMGSetType(pc,PCGAMGAGG);CHKERRQ(ierr); 1831 PetscFunctionReturn(0); 1832 } 1833 1834 /*@C 1835 PCGAMGInitializePackage - This function initializes everything in the PCGAMG package. It is called 1836 from PCInitializePackage(). 1837 1838 Level: developer 1839 1840 .seealso: PetscInitialize() 1841 @*/ 1842 PetscErrorCode PCGAMGInitializePackage(void) 1843 { 1844 PetscErrorCode ierr; 1845 PetscInt l; 1846 1847 PetscFunctionBegin; 1848 if (PCGAMGPackageInitialized) PetscFunctionReturn(0); 1849 PCGAMGPackageInitialized = PETSC_TRUE; 1850 ierr = PetscFunctionListAdd(&GAMGList,PCGAMGGEO,PCCreateGAMG_GEO);CHKERRQ(ierr); 1851 ierr = PetscFunctionListAdd(&GAMGList,PCGAMGAGG,PCCreateGAMG_AGG);CHKERRQ(ierr); 1852 ierr = PetscFunctionListAdd(&GAMGList,PCGAMGCLASSICAL,PCCreateGAMG_Classical);CHKERRQ(ierr); 1853 ierr = PetscRegisterFinalize(PCGAMGFinalizePackage);CHKERRQ(ierr); 1854 1855 /* general events */ 1856 ierr = PetscLogEventRegister("PCGAMGGraph_AGG", 0, &PC_GAMGGraph_AGG);CHKERRQ(ierr); 1857 ierr = PetscLogEventRegister("PCGAMGGraph_GEO", PC_CLASSID, &PC_GAMGGraph_GEO);CHKERRQ(ierr); 1858 ierr = PetscLogEventRegister("PCGAMGCoarse_AGG", PC_CLASSID, &PC_GAMGCoarsen_AGG);CHKERRQ(ierr); 1859 ierr = PetscLogEventRegister("PCGAMGCoarse_GEO", PC_CLASSID, &PC_GAMGCoarsen_GEO);CHKERRQ(ierr); 1860 ierr = PetscLogEventRegister("PCGAMGProl_AGG", PC_CLASSID, &PC_GAMGProlongator_AGG);CHKERRQ(ierr); 1861 ierr = PetscLogEventRegister("PCGAMGProl_GEO", PC_CLASSID, &PC_GAMGProlongator_GEO);CHKERRQ(ierr); 1862 ierr = PetscLogEventRegister("PCGAMGPOpt_AGG", PC_CLASSID, &PC_GAMGOptProlongator_AGG);CHKERRQ(ierr); 1863 1864 ierr = PetscLogEventRegister("GAMG: createProl", PC_CLASSID, &petsc_gamg_setup_events[SET1]);CHKERRQ(ierr); 1865 ierr = PetscLogEventRegister(" Graph", PC_CLASSID, &petsc_gamg_setup_events[GRAPH]);CHKERRQ(ierr); 1866 /* PetscLogEventRegister(" G.Mat", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_MAT]); */ 1867 /* PetscLogEventRegister(" G.Filter", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_FILTER]); */ 1868 /* PetscLogEventRegister(" G.Square", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_SQR]); */ 1869 ierr = PetscLogEventRegister(" MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[SET4]);CHKERRQ(ierr); 1870 ierr = PetscLogEventRegister(" geo: growSupp", PC_CLASSID, &petsc_gamg_setup_events[SET5]);CHKERRQ(ierr); 1871 ierr = PetscLogEventRegister(" geo: triangle", PC_CLASSID, &petsc_gamg_setup_events[SET6]);CHKERRQ(ierr); 1872 ierr = PetscLogEventRegister(" search-set", PC_CLASSID, &petsc_gamg_setup_events[FIND_V]);CHKERRQ(ierr); 1873 ierr = PetscLogEventRegister(" SA: col data", PC_CLASSID, &petsc_gamg_setup_events[SET7]);CHKERRQ(ierr); 1874 ierr = PetscLogEventRegister(" SA: frmProl0", PC_CLASSID, &petsc_gamg_setup_events[SET8]);CHKERRQ(ierr); 1875 ierr = PetscLogEventRegister(" SA: smooth", PC_CLASSID, &petsc_gamg_setup_events[SET9]);CHKERRQ(ierr); 1876 ierr = PetscLogEventRegister("GAMG: partLevel", PC_CLASSID, &petsc_gamg_setup_events[SET2]);CHKERRQ(ierr); 1877 ierr = PetscLogEventRegister(" repartition", PC_CLASSID, &petsc_gamg_setup_events[SET12]);CHKERRQ(ierr); 1878 ierr = PetscLogEventRegister(" Invert-Sort", PC_CLASSID, &petsc_gamg_setup_events[SET13]);CHKERRQ(ierr); 1879 ierr = PetscLogEventRegister(" Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14]);CHKERRQ(ierr); 1880 ierr = PetscLogEventRegister(" Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15]);CHKERRQ(ierr); 1881 for (l=0;l<PETSC_MG_MAXLEVELS;l++) { 1882 char ename[32]; 1883 1884 ierr = PetscSNPrintf(ename,sizeof(ename),"PCGAMG Squ l%02d",l);CHKERRQ(ierr); 1885 ierr = PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0]);CHKERRQ(ierr); 1886 ierr = PetscSNPrintf(ename,sizeof(ename),"PCGAMG Gal l%02d",l);CHKERRQ(ierr); 1887 ierr = PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1]);CHKERRQ(ierr); 1888 ierr = PetscSNPrintf(ename,sizeof(ename),"PCGAMG Opt l%02d",l);CHKERRQ(ierr); 1889 ierr = PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2]);CHKERRQ(ierr); 1890 } 1891 /* PetscLogEventRegister(" PL move data", PC_CLASSID, &petsc_gamg_setup_events[SET13]); */ 1892 /* PetscLogEventRegister("GAMG: fix", PC_CLASSID, &petsc_gamg_setup_events[SET10]); */ 1893 /* PetscLogEventRegister("GAMG: set levels", PC_CLASSID, &petsc_gamg_setup_events[SET11]); */ 1894 /* create timer stages */ 1895 #if defined(GAMG_STAGES) 1896 { 1897 char str[32]; 1898 PetscInt lidx; 1899 sprintf(str,"MG Level %d (finest)",0); 1900 ierr = PetscLogStageRegister(str, &gamg_stages[0]);CHKERRQ(ierr); 1901 for (lidx=1; lidx<9; lidx++) { 1902 sprintf(str,"MG Level %d",(int)lidx); 1903 ierr = PetscLogStageRegister(str, &gamg_stages[lidx]);CHKERRQ(ierr); 1904 } 1905 } 1906 #endif 1907 PetscFunctionReturn(0); 1908 } 1909 1910 /*@C 1911 PCGAMGFinalizePackage - This function frees everything from the PCGAMG package. It is 1912 called from PetscFinalize() automatically. 1913 1914 Level: developer 1915 1916 .seealso: PetscFinalize() 1917 @*/ 1918 PetscErrorCode PCGAMGFinalizePackage(void) 1919 { 1920 PetscErrorCode ierr; 1921 1922 PetscFunctionBegin; 1923 PCGAMGPackageInitialized = PETSC_FALSE; 1924 ierr = PetscFunctionListDestroy(&GAMGList);CHKERRQ(ierr); 1925 PetscFunctionReturn(0); 1926 } 1927 1928 /*@C 1929 PCGAMGRegister - Register a PCGAMG implementation. 1930 1931 Input Parameters: 1932 + type - string that will be used as the name of the GAMG type. 1933 - create - function for creating the gamg context. 1934 1935 Level: advanced 1936 1937 .seealso: PCGAMGType, PCGAMG, PCGAMGSetType() 1938 @*/ 1939 PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC)) 1940 { 1941 PetscErrorCode ierr; 1942 1943 PetscFunctionBegin; 1944 ierr = PCGAMGInitializePackage();CHKERRQ(ierr); 1945 ierr = PetscFunctionListAdd(&GAMGList,type,create);CHKERRQ(ierr); 1946 PetscFunctionReturn(0); 1947 } 1948 1949