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