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