Actual source code: util.c
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>
6: #include <petsc/private/kspimpl.h>
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
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
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
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;
32: PetscFunctionBegin;
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));
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 = (PetscScalar)data_in[dir * nloc + kk]; /* col oriented */
49: datas[dir * nnodes + kk] = PetscRealPart(crd); // get local part now
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: }
67: PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab)
68: {
69: PetscInt kk;
71: PetscFunctionBegin;
72: a_tab->size = a_size;
73: PetscCall(PetscMalloc2(a_size, &a_tab->table, a_size, &a_tab->data));
74: for (kk = 0; kk < a_size; kk++) a_tab->table[kk] = -1;
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab)
79: {
80: PetscFunctionBegin;
81: PetscCall(PetscFree2(a_tab->table, a_tab->data));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data)
86: {
87: PetscInt kk, idx;
89: PetscFunctionBegin;
90: PetscCheck(a_key >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Negative key %" PetscInt_FMT ".", a_key);
91: for (kk = 0, idx = GAMG_HASH(a_key); kk < a_tab->size; kk++, idx = (idx == (a_tab->size - 1)) ? 0 : idx + 1) {
92: if (a_tab->table[idx] == a_key) {
93: /* exists */
94: a_tab->data[idx] = a_data;
95: break;
96: } else if (a_tab->table[idx] == -1) {
97: /* add */
98: a_tab->table[idx] = a_key;
99: a_tab->data[idx] = a_data;
100: break;
101: }
102: }
103: if (kk == a_tab->size) {
104: /* this is not to efficient, waiting until completely full */
105: PetscInt oldsize = a_tab->size, new_size = 2 * a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data;
107: a_tab->size = new_size;
108: PetscCall(PetscMalloc2(a_tab->size, &a_tab->table, a_tab->size, &a_tab->data));
109: for (kk = 0; kk < a_tab->size; kk++) a_tab->table[kk] = -1;
110: for (kk = 0; kk < oldsize; kk++) {
111: if (oldtable[kk] != -1) PetscCall(PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]));
112: }
113: PetscCall(PetscFree2(oldtable, olddata));
114: PetscCall(PCGAMGHashTableAdd(a_tab, a_key, a_data));
115: }
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }