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 */
PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt * a_stride,PetscReal ** a_data_out)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