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: /*
9: PCGAMGGetDataWithGhosts - Get array of local + ghost data with local data
10: hacks into Mat MPIAIJ so this must have size > 1
12: Input Parameter:
13: . Gmat - MPIAIJ matrix for scatters
14: . data_sz - number of data terms per node (# cols in output)
15: . data_in[nloc*data_sz] - column oriented local data
17: Output Parameter:
18: . a_stride - number of rows of output (locals+ghosts)
19: . a_data_out[stride*data_sz] - output data with ghosts
21: */
22: PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat, PetscInt data_sz, PetscReal data_in[], PetscInt *a_stride, PetscReal **a_data_out)
23: {
24: Vec tmp_crds;
25: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Gmat->data;
26: PetscInt nnodes, num_ghosts, dir, kk, jj, my0, Iend, nloc;
27: PetscScalar *data_arr;
28: PetscReal *datas;
29: PetscBool isMPIAIJ;
31: PetscFunctionBegin;
32: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ));
33: PetscCall(MatGetOwnershipRange(Gmat, &my0, &Iend));
34: nloc = Iend - my0;
35: PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
36: nnodes = num_ghosts + nloc;
37: *a_stride = nnodes;
38: PetscCall(MatCreateVecs(Gmat, &tmp_crds, NULL));
40: PetscCall(PetscMalloc1(data_sz * nnodes, &datas));
41: for (dir = 0; dir < data_sz; dir++) {
42: /* set local, and global */
43: for (kk = 0; kk < nloc; kk++) {
44: PetscInt gid = my0 + kk;
45: PetscScalar crd = (PetscScalar)data_in[dir * nloc + kk]; /* col oriented */
46: datas[dir * nnodes + kk] = PetscRealPart(crd); // get local part now
48: PetscCall(VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES));
49: }
50: PetscCall(VecAssemblyBegin(tmp_crds));
51: PetscCall(VecAssemblyEnd(tmp_crds));
52: /* scatter / gather ghost data and add to end of output data */
53: PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
54: PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
55: PetscCall(VecGetArray(mpimat->lvec, &data_arr));
56: for (kk = nloc, jj = 0; jj < num_ghosts; kk++, jj++) datas[dir * nnodes + kk] = PetscRealPart(data_arr[jj]);
57: PetscCall(VecRestoreArray(mpimat->lvec, &data_arr));
58: }
59: PetscCall(VecDestroy(&tmp_crds));
60: *a_data_out = datas;
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab)
65: {
66: PetscInt kk;
68: PetscFunctionBegin;
69: a_tab->size = a_size;
70: PetscCall(PetscMalloc2(a_size, &a_tab->table, a_size, &a_tab->data));
71: for (kk = 0; kk < a_size; kk++) a_tab->table[kk] = -1;
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab)
76: {
77: PetscFunctionBegin;
78: PetscCall(PetscFree2(a_tab->table, a_tab->data));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data)
83: {
84: PetscInt kk, idx;
86: PetscFunctionBegin;
87: PetscCheck(a_key >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Negative key %" PetscInt_FMT ".", a_key);
88: for (kk = 0, idx = GAMG_HASH(a_key); kk < a_tab->size; kk++, idx = (idx == (a_tab->size - 1)) ? 0 : idx + 1) {
89: if (a_tab->table[idx] == a_key) {
90: /* exists */
91: a_tab->data[idx] = a_data;
92: break;
93: } else if (a_tab->table[idx] == -1) {
94: /* add */
95: a_tab->table[idx] = a_key;
96: a_tab->data[idx] = a_data;
97: break;
98: }
99: }
100: if (kk == a_tab->size) {
101: /* this is not to efficient, waiting until completely full */
102: PetscInt oldsize = a_tab->size, new_size = 2 * a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data;
104: a_tab->size = new_size;
105: PetscCall(PetscMalloc2(a_tab->size, &a_tab->table, a_tab->size, &a_tab->data));
106: for (kk = 0; kk < a_tab->size; kk++) a_tab->table[kk] = -1;
107: for (kk = 0; kk < oldsize; kk++) {
108: if (oldtable[kk] != -1) PetscCall(PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]));
109: }
110: PetscCall(PetscFree2(oldtable, olddata));
111: PetscCall(PCGAMGHashTableAdd(a_tab, a_key, a_data));
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }