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 <petsc/private/kspimpl.h> 7 8 // PetscClangLinter pragma disable: -fdoc-sowing-chars 9 /* 10 PCGAMGGetDataWithGhosts - Get array of local + ghost data with local data 11 hacks into Mat MPIAIJ so this must have size > 1 12 13 Input Parameters: 14 + Gmat - MPIAIJ matrix for scatters 15 . data_sz - number of data terms per node (# cols in output) 16 - data_in - column-oriented local data of size nloc*data_sz 17 18 Output Parameters: 19 + a_stride - number of rows of output (locals+ghosts) 20 - a_data_out - output data with ghosts of size stride*data_sz 21 22 */ 23 PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat, PetscInt data_sz, PetscReal data_in[], PetscInt *a_stride, PetscReal **a_data_out) 24 { 25 Vec tmp_crds; 26 Mat_MPIAIJ *mpimat; 27 PetscInt nnodes, num_ghosts, dir, kk, jj, my0, Iend, nloc; 28 PetscScalar *data_arr; 29 PetscReal *datas; 30 PetscBool isMPIAIJ; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(Gmat, MAT_CLASSID, 1); 34 mpimat = (Mat_MPIAIJ *)Gmat->data; 35 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ)); 36 PetscCall(MatGetOwnershipRange(Gmat, &my0, &Iend)); 37 nloc = Iend - my0; 38 PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts)); 39 nnodes = num_ghosts + nloc; 40 *a_stride = nnodes; 41 PetscCall(MatCreateVecs(Gmat, &tmp_crds, NULL)); 42 43 PetscCall(PetscMalloc1(data_sz * nnodes, &datas)); 44 for (dir = 0; dir < data_sz; dir++) { 45 /* set local, and global */ 46 for (kk = 0; kk < nloc; kk++) { 47 PetscInt gid = my0 + kk; 48 PetscScalar crd = data_in[dir * nloc + kk]; /* col oriented */ 49 datas[dir * nnodes + kk] = PetscRealPart(crd); // get local part now 50 51 PetscCall(VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES)); 52 } 53 PetscCall(VecAssemblyBegin(tmp_crds)); 54 PetscCall(VecAssemblyEnd(tmp_crds)); 55 /* scatter / gather ghost data and add to end of output data */ 56 PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 57 PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 58 PetscCall(VecGetArray(mpimat->lvec, &data_arr)); 59 for (kk = nloc, jj = 0; jj < num_ghosts; kk++, jj++) datas[dir * nnodes + kk] = PetscRealPart(data_arr[jj]); 60 PetscCall(VecRestoreArray(mpimat->lvec, &data_arr)); 61 } 62 PetscCall(VecDestroy(&tmp_crds)); 63 *a_data_out = datas; 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66