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, PetscCtx ctx)
6: {
7: PetscInt d, j, k;
8: PetscReal diff = 0.0, *coords = (PetscReal *)(ctx);
10: #if !PetscDefined(HAVE_OPENMP)
11: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
22: #else
23: return PETSC_SUCCESS;
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, overlap = 1;
34: PetscMPIInt size;
35: PetscReal *coords, *gcoords;
36: MatHtoolKernelFn *kernel = GenEntries;
37: PetscBool flg, sym = PETSC_FALSE;
38: PetscRandom rdm;
40: PetscFunctionBeginUser;
41: PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help));
42: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL));
43: PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL));
44: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
45: PetscCall(PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL));
46: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
47: M = size * m;
48: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
49: PetscCall(PetscMalloc1(m * dim, &coords));
50: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
51: PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords));
52: PetscCall(PetscCalloc1(M * dim, &gcoords));
53: PetscCallMPI(MPI_Exscan(&m, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
54: PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim));
55: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
56: PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A));
57: PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym));
58: PetscCall(MatSetFromOptions(A));
59: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
60: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
61: PetscCall(MatCreateVecs(A, &b, &x));
62: PetscCall(VecSetRandom(b, rdm));
63: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
64: PetscCall(KSPSetOperators(ksp, A, A));
65: PetscCall(KSPSetFromOptions(ksp));
66: PetscCall(KSPGetPC(ksp, &pc));
67: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCHPDDM, &flg));
68: if (flg) {
69: #if PetscDefined(HAVE_HPDDM) && PetscDefined(HAVE_DYNAMIC_LIBRARIES) && PetscDefined(USE_SHARED_LIBRARIES)
70: Mat aux;
71: IS is;
72: PetscInt n;
74: PetscCall(MatGetOwnershipRange(A, &begin, &n));
75: n -= begin;
76: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, begin, 1, &is));
77: PetscCall(MatIncreaseOverlap(A, 1, &is, overlap));
78: PetscCall(ISGetLocalSize(is, &n));
79: PetscCall(MatCreateDense(PETSC_COMM_SELF, n, n, n, n, NULL, &aux));
80: PetscCall(MatSetOption(aux, MAT_SYMMETRIC, sym));
81: PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatShift(aux, 1.0)); /* just the local identity matrix, not very meaningful numerically, but just testing that the necessary plumbing is there */
84: PetscCall(PCHPDDMSetAuxiliaryMat(pc, is, aux, NULL, NULL));
85: PetscCall(ISDestroy(&is));
86: PetscCall(MatDestroy(&aux));
87: #endif
88: }
89: PetscCall(KSPSolve(ksp, b, x));
90: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &flg));
91: if (flg) {
92: Vec y;
93: PetscReal norm;
94: PetscScalar alpha = 2.0;
96: PetscCall(VecDuplicate(x, &y));
97: PetscCall(MatScale(A, alpha));
98: PetscCall(KSPSolve(ksp, b, y));
99: PetscCall(VecAXPY(x, -alpha, y));
100: PetscCall(VecNorm(x, NORM_INFINITY, &norm));
101: PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Scaled difference of both solves has nonzero norm, ||x - alpha y||_inf = %g > %g", (double)norm, (double)PETSC_SMALL);
102: PetscCall(VecDestroy(&y));
103: }
104: PetscCall(KSPDestroy(&ksp));
105: PetscCall(PetscRandomDestroy(&rdm));
106: PetscCall(VecDestroy(&b));
107: PetscCall(VecDestroy(&x));
108: PetscCall(MatDestroy(&A));
109: PetscCall(PetscFree(gcoords));
110: PetscCall(PetscFree(coords));
111: PetscCall(PetscFinalize());
112: return 0;
113: }
115: /*TEST
117: build:
118: requires: htool
120: test:
121: requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
122: nsize: 4
123: # different numbers of iterations depending on PetscScalar type
124: 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" -e "s/(1,0) : type=transpose/(1,0) : type=seqdense/g"
125: 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
126: output_file: output/ex82_1.out
128: test:
129: suffix: bjacobi
130: nsize: 4
131: args: -ksp_max_it 20 -mat_htool_epsilon 1e-2 -m_local 200 -ksp_error_if_not_converged
132: output_file: output/empty.out
134: test:
135: suffix: asm
136: nsize: 4
137: args: -ksp_max_it 20 -mat_htool_epsilon 1e-2 -m_local 200 -pc_type asm -pc_asm_sub_mat_type nest -sub_pc_type {{lu cholesky}shared output} -ksp_error_if_not_converged
138: output_file: output/empty.out
140: TEST*/