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