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: PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);
32: MatGetOwnershipRange(Gmat, &my0, &Iend);
33: nloc = Iend - my0;
34: VecGetLocalSize(mpimat->lvec, &num_ghosts);
35: nnodes = num_ghosts + nloc;
36: *a_stride = nnodes;
37: MatCreateVecs(Gmat, &tmp_crds, NULL);
39: PetscMalloc1(data_sz * nnodes, &datas);
40: for (dir = 0; dir < data_sz; dir++) {
41: /* set local, and global */
42: for (kk = 0; kk < nloc; kk++) {
43: PetscInt gid = my0 + kk;
44: PetscScalar crd = (PetscScalar)data_in[dir * nloc + kk]; /* col oriented */
45: datas[dir * nnodes + kk] = PetscRealPart(crd); // get local part now
47: VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);
48: }
49: VecAssemblyBegin(tmp_crds);
50: VecAssemblyEnd(tmp_crds);
51: /* scatter / gather ghost data and add to end of output data */
52: VecScatterBegin(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD);
53: VecScatterEnd(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD);
54: VecGetArray(mpimat->lvec, &data_arr);
55: for (kk = nloc, jj = 0; jj < num_ghosts; kk++, jj++) datas[dir * nnodes + kk] = PetscRealPart(data_arr[jj]);
56: VecRestoreArray(mpimat->lvec, &data_arr);
57: }
58: VecDestroy(&tmp_crds);
59: *a_data_out = datas;
60: return 0;
61: }
63: PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab)
64: {
65: PetscInt kk;
67: a_tab->size = a_size;
68: PetscMalloc2(a_size, &a_tab->table, a_size, &a_tab->data);
69: for (kk = 0; kk < a_size; kk++) a_tab->table[kk] = -1;
70: return 0;
71: }
73: PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab)
74: {
75: PetscFree2(a_tab->table, a_tab->data);
76: return 0;
77: }
79: PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data)
80: {
81: PetscInt kk, idx;
84: for (kk = 0, idx = GAMG_HASH(a_key); kk < a_tab->size; kk++, idx = (idx == (a_tab->size - 1)) ? 0 : idx + 1) {
85: if (a_tab->table[idx] == a_key) {
86: /* exists */
87: a_tab->data[idx] = a_data;
88: break;
89: } else if (a_tab->table[idx] == -1) {
90: /* add */
91: a_tab->table[idx] = a_key;
92: a_tab->data[idx] = a_data;
93: break;
94: }
95: }
96: if (kk == a_tab->size) {
97: /* this is not to efficient, waiting until completely full */
98: PetscInt oldsize = a_tab->size, new_size = 2 * a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data;
100: a_tab->size = new_size;
101: PetscMalloc2(a_tab->size, &a_tab->table, a_tab->size, &a_tab->data);
102: for (kk = 0; kk < a_tab->size; kk++) a_tab->table[kk] = -1;
103: for (kk = 0; kk < oldsize; kk++) {
104: if (oldtable[kk] != -1) PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]);
105: }
106: PetscFree2(oldtable, olddata);
107: PCGAMGHashTableAdd(a_tab, a_key, a_data);
108: }
109: return 0;
110: }