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