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*/