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: }