Actual source code: ex82.c

  1: #include <petscksp.h>

  3: static char help[] = "Solves a linear system using PCHPDDM and MATHTOOL.\n\n";

  5: static PetscErrorCode GenEntries(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, void *ctx)
  6: {
  7:   PetscInt  d, j, k;
  8:   PetscReal diff = 0.0, *coords = (PetscReal *)(ctx);

 10: #if !PetscDefined(HAVE_OPENMP)
 12: #endif
 13:   for (j = 0; j < M; j++) {
 14:     for (k = 0; k < N; k++) {
 15:       diff = 0.0;
 16:       for (d = 0; d < sdim; d++) diff += (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]) * (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]);
 17:       ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff));
 18:     }
 19:   }
 20: #if !PetscDefined(HAVE_OPENMP)
 21:   return 0;
 22: #else
 23:   return 0;
 24: #endif
 25: }

 27: int main(int argc, char **argv)
 28: {
 29:   KSP            ksp;
 30:   PC             pc;
 31:   Vec            b, x;
 32:   Mat            A;
 33:   PetscInt       m = 100, dim = 3, M, begin = 0, n = 0, overlap = 1;
 34:   PetscMPIInt    size;
 35:   PetscReal     *coords, *gcoords;
 36:   MatHtoolKernel kernel = GenEntries;
 37:   PetscBool      flg, sym = PETSC_FALSE;
 38:   PetscRandom    rdm;

 41:   PetscInitialize(&argc, &argv, (char *)NULL, help);
 42:   PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL);
 43:   PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL);
 44:   PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
 45:   PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL);
 46:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 47:   M = size * m;
 48:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 49:   PetscMalloc1(m * dim, &coords);
 50:   PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
 51:   PetscRandomGetValuesReal(rdm, m * dim, coords);
 52:   PetscCalloc1(M * dim, &gcoords);
 53:   MPI_Exscan(&m, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD);
 54:   PetscArraycpy(gcoords + begin * dim, coords, m * dim);
 55:   MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);
 56:   MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A);
 57:   MatSetOption(A, MAT_SYMMETRIC, sym);
 58:   MatSetFromOptions(A);
 59:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 60:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 61:   MatCreateVecs(A, &b, &x);
 62:   VecSetRandom(b, rdm);
 63:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 64:   KSPSetOperators(ksp, A, A);
 65:   KSPSetFromOptions(ksp);
 66:   KSPGetPC(ksp, &pc);
 67:   PetscObjectTypeCompare((PetscObject)pc, PCHPDDM, &flg);
 68:   if (flg) {
 69: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
 70:     Mat aux;
 71:     IS  is;
 72:     MatGetOwnershipRange(A, &begin, &n);
 73:     n -= begin;
 74:     ISCreateStride(PETSC_COMM_SELF, n, begin, 1, &is);
 75:     MatIncreaseOverlap(A, 1, &is, overlap);
 76:     ISGetLocalSize(is, &n);
 77:     MatCreateDense(PETSC_COMM_SELF, n, n, n, n, NULL, &aux);
 78:     MatSetOption(aux, MAT_SYMMETRIC, sym);
 79:     MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY);
 80:     MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY);
 81:     MatShift(aux, 1.0); /* just the local identity matrix, not very meaningful numerically, but just testing that the necessary plumbing is there */
 82:     PCHPDDMSetAuxiliaryMat(pc, is, aux, NULL, NULL);
 83:     ISDestroy(&is);
 84:     MatDestroy(&aux);
 85: #endif
 86:   }
 87:   KSPSolve(ksp, b, x);
 88:   KSPDestroy(&ksp);
 89:   PetscRandomDestroy(&rdm);
 90:   VecDestroy(&b);
 91:   VecDestroy(&x);
 92:   MatDestroy(&A);
 93:   PetscFree(gcoords);
 94:   PetscFree(coords);
 95:   PetscFinalize();
 96:   return 0;
 97: }

 99: /*TEST

101:    build:
102:       requires: htool hpddm

104:    test:
105:       requires: htool hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
106:       nsize: 4
107:       # different numbers of iterations depending on PetscScalar type
108:       filter: sed -e "s/symmetry: S/symmetry: N/g" -e "/number of dense/d" -e "s/Linear solve converged due to CONVERGED_RTOL iterations 13/Linear solve converged due to CONVERGED_RTOL iterations 18/g"
109:       args: -ksp_view -ksp_converged_reason -mat_htool_epsilon 1e-2 -m_local 200 -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 1 -pc_hpddm_coarse_pc_type lu -pc_hpddm_levels_1_eps_gen_non_hermitian -symmetric {{false true}shared output} -overlap 2
110:       output_file: output/ex82_1.out

112: TEST*/