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;
26: PetscInitialize(&argc, &args, (char *)0, help);
27: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
28: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
29: PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
30: PetscOptionsGetBool(NULL, NULL, "-test_hmg_interface", &test, NULL);
31: PetscOptionsGetBool(NULL, NULL, "-test_reuse_interpolation", &reuse, NULL);
32: PetscOptionsGetBool(NULL, NULL, "-view_explicit_mat", &viewexpl, NULL);
34: MatCreate(PETSC_COMM_WORLD, &A);
35: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n * bs, m * n * bs);
36: MatSetBlockSize(A, bs);
37: MatSetFromOptions(A);
38: MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
39: MatSeqAIJSetPreallocation(A, 5, NULL);
40: #if defined(PETSC_HAVE_HYPRE)
41: MatHYPRESetPreallocation(A, 5, NULL, 5, NULL);
42: #endif
44: 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: 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: 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: 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: 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: MatSetValues(A, 1, &II, 1, &II, &v, ADD_VALUES);
86: }
87: }
89: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
91: if (viewexpl) {
92: Mat E;
93: MatComputeOperator(A, MATAIJ, &E);
94: MatView(E, NULL);
95: MatDestroy(&E);
96: }
98: MatCreateVecs(A, &u, NULL);
99: VecSetFromOptions(u);
100: VecDuplicate(u, &b);
101: VecDuplicate(b, &x);
103: VecSet(u, 1.0);
104: MatMult(A, u, b);
106: flg = PETSC_FALSE;
107: PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL);
108: if (flg) VecView(u, PETSC_VIEWER_STDOUT_WORLD);
110: KSPCreate(PETSC_COMM_WORLD, &ksp);
111: KSPSetOperators(ksp, A, A);
112: KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT);
114: if (test) {
115: KSPGetPC(ksp, &pc);
116: PCSetType(pc, PCHMG);
117: PCHMGSetInnerPCType(pc, PCGAMG);
118: PCHMGSetReuseInterpolation(pc, PETSC_TRUE);
119: PCHMGSetUseSubspaceCoarsening(pc, PETSC_TRUE);
120: PCHMGUseMatMAIJ(pc, PETSC_FALSE);
121: PCHMGSetCoarseningComponent(pc, 0);
122: }
124: KSPSetFromOptions(ksp);
125: KSPSolve(ksp, b, x);
127: if (reuse) {
128: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
130: KSPSolve(ksp, b, x);
131: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
133: KSPSolve(ksp, b, x);
134: /* Make sparsity pattern different and reuse interpolation */
135: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
136: MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
137: MatGetSize(A, &m, NULL);
138: n = 0;
139: v = 0;
140: m--;
141: /* Connect the last element to the first element */
142: MatSetValue(A, m, n, v, ADD_VALUES);
143: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
144: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
145: KSPSolve(ksp, b, x);
146: }
148: VecAXPY(x, -1.0, u);
149: VecNorm(x, NORM_2, &norm);
150: KSPGetIterationNumber(ksp, &its);
152: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its);
154: KSPDestroy(&ksp);
155: VecDestroy(&u);
156: VecDestroy(&x);
157: VecDestroy(&b);
158: MatDestroy(&A);
160: 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*/