Actual source code: ex4.c
2: static char help[] = "Solves a linear system in parallel with KSP and HMG.\n\
3: Input parameters include:\n\
4: -view_exact_sol : write exact solution vector to stdout\n\
5: -m <mesh_x> : number of mesh points in x-direction\n\
6: -n <mesh_y> : number of mesh points in y-direction\n\
7: -bs : number of variables on each mesh vertex \n\n";
9: /*
10: Simple example is used to test PCHMG
11: */
12: #include <petscksp.h>
14: int main(int argc, char **args)
15: {
16: Vec x, b, u; /* approx solution, RHS, exact solution */
17: Mat A; /* linear system matrix */
18: KSP ksp; /* linear solver context */
19: PetscReal norm; /* norm of solution error */
20: PetscInt i, j, Ii, J, Istart, Iend, m = 8, n = 7, its, bs = 1, II, JJ, jj;
21: PetscBool flg, test = PETSC_FALSE, reuse = PETSC_FALSE, viewexpl = PETSC_FALSE;
22: PetscScalar v;
23: PC pc;
25: PetscFunctionBeginUser;
26: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
27: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
28: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
29: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
30: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_hmg_interface", &test, NULL));
31: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_reuse_interpolation", &reuse, NULL));
32: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_explicit_mat", &viewexpl, NULL));
34: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
35: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n * bs, m * n * bs));
36: PetscCall(MatSetBlockSize(A, bs));
37: PetscCall(MatSetFromOptions(A));
38: PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
39: PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
40: #if defined(PETSC_HAVE_HYPRE)
41: PetscCall(MatHYPRESetPreallocation(A, 5, NULL, 5, NULL));
42: #endif
44: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
46: for (Ii = Istart / bs; Ii < Iend / bs; Ii++) {
47: v = -1.0;
48: i = Ii / n;
49: j = Ii - i * n;
50: if (i > 0) {
51: J = Ii - n;
52: for (jj = 0; jj < bs; jj++) {
53: II = Ii * bs + jj;
54: JJ = J * bs + jj;
55: PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
56: }
57: }
58: if (i < m - 1) {
59: J = Ii + n;
60: for (jj = 0; jj < bs; jj++) {
61: II = Ii * bs + jj;
62: JJ = J * bs + jj;
63: PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
64: }
65: }
66: if (j > 0) {
67: J = Ii - 1;
68: for (jj = 0; jj < bs; jj++) {
69: II = Ii * bs + jj;
70: JJ = J * bs + jj;
71: PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
72: }
73: }
74: if (j < n - 1) {
75: J = Ii + 1;
76: for (jj = 0; jj < bs; jj++) {
77: II = Ii * bs + jj;
78: JJ = J * bs + jj;
79: PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
80: }
81: }
82: v = 4.0;
83: for (jj = 0; jj < bs; jj++) {
84: II = Ii * bs + jj;
85: PetscCall(MatSetValues(A, 1, &II, 1, &II, &v, ADD_VALUES));
86: }
87: }
89: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
90: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
91: if (viewexpl) {
92: Mat E;
93: PetscCall(MatComputeOperator(A, MATAIJ, &E));
94: PetscCall(MatView(E, NULL));
95: PetscCall(MatDestroy(&E));
96: }
98: PetscCall(MatCreateVecs(A, &u, NULL));
99: PetscCall(VecSetFromOptions(u));
100: PetscCall(VecDuplicate(u, &b));
101: PetscCall(VecDuplicate(b, &x));
103: PetscCall(VecSet(u, 1.0));
104: PetscCall(MatMult(A, u, b));
106: flg = PETSC_FALSE;
107: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL));
108: if (flg) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
110: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
111: PetscCall(KSPSetOperators(ksp, A, A));
112: PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT));
114: if (test) {
115: PetscCall(KSPGetPC(ksp, &pc));
116: PetscCall(PCSetType(pc, PCHMG));
117: PetscCall(PCHMGSetInnerPCType(pc, PCGAMG));
118: PetscCall(PCHMGSetReuseInterpolation(pc, PETSC_TRUE));
119: PetscCall(PCHMGSetUseSubspaceCoarsening(pc, PETSC_TRUE));
120: PetscCall(PCHMGUseMatMAIJ(pc, PETSC_FALSE));
121: PetscCall(PCHMGSetCoarseningComponent(pc, 0));
122: }
124: PetscCall(KSPSetFromOptions(ksp));
125: PetscCall(KSPSolve(ksp, b, x));
127: if (reuse) {
128: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
129: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
130: PetscCall(KSPSolve(ksp, b, x));
131: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
132: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
133: PetscCall(KSPSolve(ksp, b, x));
134: /* Make sparsity pattern different and reuse interpolation */
135: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
136: PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
137: PetscCall(MatGetSize(A, &m, NULL));
138: n = 0;
139: v = 0;
140: m--;
141: /* Connect the last element to the first element */
142: PetscCall(MatSetValue(A, m, n, v, ADD_VALUES));
143: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
144: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
145: PetscCall(KSPSolve(ksp, b, x));
146: }
148: PetscCall(VecAXPY(x, -1.0, u));
149: PetscCall(VecNorm(x, NORM_2, &norm));
150: PetscCall(KSPGetIterationNumber(ksp, &its));
152: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));
154: PetscCall(KSPDestroy(&ksp));
155: PetscCall(VecDestroy(&u));
156: PetscCall(VecDestroy(&x));
157: PetscCall(VecDestroy(&b));
158: PetscCall(MatDestroy(&A));
160: PetscCall(PetscFinalize());
161: return 0;
162: }
164: /*TEST
166: build:
167: requires: !complex !single
169: test:
170: suffix: hypre
171: nsize: 2
172: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
173: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre
175: test:
176: suffix: hypre_bs4
177: nsize: 2
178: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
179: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1
181: test:
182: suffix: hypre_asm
183: nsize: 2
184: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
185: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_3_pc_type asm
187: test:
188: suffix: hypre_fieldsplit
189: nsize: 2
190: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
191: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -mg_levels_4_pc_type fieldsplit
193: test:
194: suffix: gamg
195: nsize: 2
196: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg
198: test:
199: suffix: gamg_bs4
200: nsize: 2
201: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1
203: test:
204: suffix: gamg_asm
205: nsize: 2
206: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_1_pc_type asm
208: test:
209: suffix: gamg_fieldsplit
210: nsize: 2
211: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -mg_levels_1_pc_type fieldsplit
213: test:
214: suffix: interface
215: nsize: 2
216: args: -ksp_monitor -ksp_rtol 1e-6 -test_hmg_interface 1 -bs 4
218: test:
219: suffix: reuse
220: nsize: 2
221: args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_reuse_interpolation 1 -test_reuse_interpolation 1 -hmg_inner_pc_type gamg
223: test:
224: suffix: component
225: nsize: 2
226: args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_coarsening_component 2 -pc_hmg_use_subspace_coarsening 1 -bs 4 -hmg_inner_pc_type gamg
228: testset:
229: output_file: output/ex4_expl.out
230: nsize: {{1 2}}
231: filter: grep -v " MPI process" | grep -v " type:" | grep -v "Mat Object"
232: args: -ksp_converged_reason -view_explicit_mat -pc_type none -ksp_type {{cg gmres}}
233: test:
234: suffix: expl_aij
235: args: -mat_type aij
236: test:
237: suffix: expl_hypre
238: requires: hypre
239: args: -mat_type hypre
241: test:
242: suffix: hypre_device
243: nsize: {{1 2}}
244: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
245: args: -mat_type hypre -ksp_converged_reason -pc_type hypre -m 13 -n 17
247: test:
248: suffix: hypre_device_cusparse
249: output_file: output/ex4_hypre_device.out
250: nsize: {{1 2}}
251: requires: hypre cuda defined(PETSC_HAVE_HYPRE_DEVICE)
252: args: -mat_type {{aij aijcusparse}} -vec_type cuda -ksp_converged_reason -pc_type hypre -m 13 -n 17
254: TEST*/