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 /* 9 PCGAMGGetDataWithGhosts - Get array of local + ghost data with local data 10 hacks into Mat MPIAIJ so this must have size > 1 11 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 16 17 Output Parameter: 18 . a_stride - number of rows of output (locals+ghosts) 19 . a_data_out[stride*data_sz] - output data with ghosts 20 21 */ 22 PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat, PetscInt data_sz, PetscReal data_in[], PetscInt *a_stride, PetscReal **a_data_out) { 23 Vec tmp_crds; 24 Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)Gmat->data; 25 PetscInt nnodes, num_ghosts, dir, kk, jj, my0, Iend, nloc; 26 PetscScalar *data_arr; 27 PetscReal *datas; 28 PetscBool isMPIAIJ; 29 30 PetscFunctionBegin; 31 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ)); 32 PetscCall(MatGetOwnershipRange(Gmat, &my0, &Iend)); 33 nloc = Iend - my0; 34 PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts)); 35 nnodes = num_ghosts + nloc; 36 *a_stride = nnodes; 37 PetscCall(MatCreateVecs(Gmat, &tmp_crds, NULL)); 38 39 PetscCall(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 46 47 PetscCall(VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES)); 48 } 49 PetscCall(VecAssemblyBegin(tmp_crds)); 50 PetscCall(VecAssemblyEnd(tmp_crds)); 51 /* scatter / gather ghost data and add to end of output data */ 52 PetscCall(VecScatterBegin(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 53 PetscCall(VecScatterEnd(mpimat->Mvctx, tmp_crds, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 54 PetscCall(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 PetscCall(VecRestoreArray(mpimat->lvec, &data_arr)); 57 } 58 PetscCall(VecDestroy(&tmp_crds)); 59 *a_data_out = datas; 60 PetscFunctionReturn(0); 61 } 62 63 PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab) { 64 PetscInt kk; 65 66 PetscFunctionBegin; 67 a_tab->size = a_size; 68 PetscCall(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 PetscFunctionReturn(0); 71 } 72 73 PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab) { 74 PetscFunctionBegin; 75 PetscCall(PetscFree2(a_tab->table, a_tab->data)); 76 PetscFunctionReturn(0); 77 } 78 79 PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data) { 80 PetscInt kk, idx; 81 82 PetscFunctionBegin; 83 PetscCheck(a_key >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Negative key %" PetscInt_FMT ".", a_key); 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; 99 100 a_tab->size = new_size; 101 PetscCall(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) { PetscCall(PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk])); } 105 } 106 PetscCall(PetscFree2(oldtable, olddata)); 107 PetscCall(PCGAMGHashTableAdd(a_tab, a_key, a_data)); 108 } 109 PetscFunctionReturn(0); 110 } 111