Actual source code: ex7.c
1: static char help[] = "Block Jacobi preconditioner for solving a linear system in parallel with KSP.\n\
2: The code indicates the\n\
3: procedures for setting the particular block sizes and for using different\n\
4: linear solvers on the individual blocks.\n\n";
6: /*
7: Note: This example focuses on ways to customize the block Jacobi
8: preconditioner. See ex1.c and ex2.c for more detailed comments on
9: the basic usage of KSP (including working with matrices and vectors).
11: Recall: The block Jacobi method is equivalent to the ASM preconditioner
12: with zero overlap.
13: */
15: /*
16: Include "petscksp.h" so that we can use KSP solvers. Note that this file
17: automatically includes:
18: petscsys.h - base PETSc routines petscvec.h - vectors
19: petscmat.h - matrices
20: petscis.h - index sets petscksp.h - Krylov subspace methods
21: petscviewer.h - viewers petscpc.h - preconditioners
22: */
23: #include <petscksp.h>
25: int main(int argc, char **args)
26: {
27: Vec x, b, u; /* approx solution, RHS, exact solution */
28: Mat A; /* linear system matrix */
29: KSP ksp; /* KSP context */
30: KSP *subksp; /* array of local KSP contexts on this processor */
31: PC pc; /* PC context */
32: PC subpc; /* PC context for subdomain */
33: PetscReal norm; /* norm of solution error */
34: PetscInt i, j, Ii, J, *blks, m = 4, n;
35: PetscMPIInt rank, size;
36: PetscInt its, nlocal, first, Istart, Iend;
37: PetscScalar v, one = 1.0, none = -1.0;
38: PetscBool isbjacobi;
41: PetscInitialize(&argc, &args, (char *)0, help);
42: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
43: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
44: MPI_Comm_size(PETSC_COMM_WORLD, &size);
45: n = m + 2;
47: /* -------------------------------------------------------------------
48: Compute the matrix and right-hand-side vector that define
49: the linear system, Ax = b.
50: ------------------------------------------------------------------- */
52: /*
53: Create and assemble parallel matrix
54: */
55: MatCreate(PETSC_COMM_WORLD, &A);
56: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n);
57: MatSetFromOptions(A);
58: MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
59: MatSeqAIJSetPreallocation(A, 5, NULL);
60: MatGetOwnershipRange(A, &Istart, &Iend);
61: for (Ii = Istart; Ii < Iend; Ii++) {
62: v = -1.0;
63: i = Ii / n;
64: j = Ii - i * n;
65: if (i > 0) {
66: J = Ii - n;
67: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
68: }
69: if (i < m - 1) {
70: J = Ii + n;
71: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
72: }
73: if (j > 0) {
74: J = Ii - 1;
75: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
76: }
77: if (j < n - 1) {
78: J = Ii + 1;
79: MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
80: }
81: v = 4.0;
82: MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
83: }
84: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
86: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
88: /*
89: Create parallel vectors
90: */
91: MatCreateVecs(A, &u, &b);
92: VecDuplicate(u, &x);
94: /*
95: Set exact solution; then compute right-hand-side vector.
96: */
97: VecSet(u, one);
98: MatMult(A, u, b);
100: /*
101: Create linear solver context
102: */
103: KSPCreate(PETSC_COMM_WORLD, &ksp);
105: /*
106: Set operators. Here the matrix that defines the linear system
107: also serves as the preconditioning matrix.
108: */
109: KSPSetOperators(ksp, A, A);
111: /*
112: Set default preconditioner for this program to be block Jacobi.
113: This choice can be overridden at runtime with the option
114: -pc_type <type>
116: IMPORTANT NOTE: Since the inners solves below are constructed to use
117: iterative methods (such as KSPGMRES) the outer Krylov method should
118: be set to use KSPFGMRES since it is the only Krylov method (plus KSPFCG)
119: that allows the preconditioners to be nonlinear (that is have iterative methods
120: inside them). The reason these examples work is because the number of
121: iterations on the inner solves is left at the default (which is 10,000)
122: and the tolerance on the inner solves is set to be a tight value of around 10^-6.
123: */
124: KSPGetPC(ksp, &pc);
125: PCSetType(pc, PCBJACOBI);
127: /* -------------------------------------------------------------------
128: Define the problem decomposition
129: ------------------------------------------------------------------- */
131: /*
132: Call PCBJacobiSetTotalBlocks() to set individually the size of
133: each block in the preconditioner. This could also be done with
134: the runtime option
135: -pc_bjacobi_blocks <blocks>
136: Also, see the command PCBJacobiSetLocalBlocks() to set the
137: local blocks.
139: Note: The default decomposition is 1 block per processor.
140: */
141: PetscMalloc1(m, &blks);
142: for (i = 0; i < m; i++) blks[i] = n;
143: PCBJacobiSetTotalBlocks(pc, m, blks);
144: PetscFree(blks);
146: /*
147: Set runtime options
148: */
149: KSPSetFromOptions(ksp);
151: /* -------------------------------------------------------------------
152: Set the linear solvers for the subblocks
153: ------------------------------------------------------------------- */
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Basic method, should be sufficient for the needs of most users.
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: By default, the block Jacobi method uses the same solver on each
160: block of the problem. To set the same solver options on all blocks,
161: use the prefix -sub before the usual PC and KSP options, e.g.,
162: -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4
163: */
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Advanced method, setting different solvers for various blocks.
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: Note that each block's KSP context is completely independent of
170: the others, and the full range of uniprocessor KSP options is
171: available for each block. The following section of code is intended
172: to be a simple illustration of setting different linear solvers for
173: the individual blocks. These choices are obviously not recommended
174: for solving this particular problem.
175: */
176: PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &isbjacobi);
177: if (isbjacobi) {
178: /*
179: Call KSPSetUp() to set the block Jacobi data structures (including
180: creation of an internal KSP context for each block).
182: Note: KSPSetUp() MUST be called before PCBJacobiGetSubKSP().
183: */
184: KSPSetUp(ksp);
186: /*
187: Extract the array of KSP contexts for the local blocks
188: */
189: PCBJacobiGetSubKSP(pc, &nlocal, &first, &subksp);
191: /*
192: Loop over the local blocks, setting various KSP options
193: for each block.
194: */
195: for (i = 0; i < nlocal; i++) {
196: KSPGetPC(subksp[i], &subpc);
197: if (rank == 0) {
198: if (i % 2) {
199: PCSetType(subpc, PCILU);
200: } else {
201: PCSetType(subpc, PCNONE);
202: KSPSetType(subksp[i], KSPBCGS);
203: KSPSetTolerances(subksp[i], 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
204: }
205: } else {
206: PCSetType(subpc, PCJACOBI);
207: KSPSetType(subksp[i], KSPGMRES);
208: KSPSetTolerances(subksp[i], 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
209: }
210: }
211: }
213: /* -------------------------------------------------------------------
214: Solve the linear system
215: ------------------------------------------------------------------- */
217: /*
218: Solve the linear system
219: */
220: KSPSolve(ksp, b, x);
222: /* -------------------------------------------------------------------
223: Check solution and clean up
224: ------------------------------------------------------------------- */
226: /*
227: Check the error
228: */
229: VecAXPY(x, none, u);
230: VecNorm(x, NORM_2, &norm);
231: KSPGetIterationNumber(ksp, &its);
232: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its);
234: /*
235: Free work space. All PETSc objects should be destroyed when they
236: are no longer needed.
237: */
238: KSPDestroy(&ksp);
239: VecDestroy(&u);
240: VecDestroy(&x);
241: VecDestroy(&b);
242: MatDestroy(&A);
243: PetscFinalize();
244: return 0;
245: }
247: /*TEST
249: test:
250: suffix: 1
251: nsize: 2
252: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
254: test:
255: suffix: 2
256: nsize: 2
257: args: -ksp_view ::ascii_info_detail
259: test:
260: suffix: viennacl
261: requires: viennacl
262: args: -ksp_monitor_short -mat_type aijviennacl
263: output_file: output/ex7_mpiaijcusparse.out
265: test:
266: suffix: viennacl_2
267: nsize: 2
268: requires: viennacl
269: args: -ksp_monitor_short -mat_type aijviennacl
270: output_file: output/ex7_mpiaijcusparse_2.out
272: test:
273: suffix: mpiaijcusparse
274: requires: cuda
275: args: -ksp_monitor_short -mat_type aijcusparse
277: test:
278: suffix: mpiaijcusparse_2
279: nsize: 2
280: requires: cuda
281: args: -ksp_monitor_short -mat_type aijcusparse
283: test:
284: suffix: mpiaijcusparse_simple
285: requires: cuda
286: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -sub_ksp_type preonly -sub_pc_type ilu
288: test:
289: suffix: mpiaijcusparse_simple_2
290: nsize: 2
291: requires: cuda
292: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -sub_ksp_type preonly -sub_pc_type ilu
294: test:
295: suffix: mpiaijcusparse_3
296: requires: cuda
297: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
299: test:
300: suffix: mpiaijcusparse_4
301: nsize: 2
302: requires: cuda
303: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
305: testset:
306: args: -ksp_monitor_short -pc_type gamg -ksp_view -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10
307: test:
308: suffix: gamg_cuda
309: nsize: {{1 2}separate output}
310: requires: cuda
311: args: -mat_type aijcusparse
312: args: -pc_gamg_threshold -1
313: test:
314: suffix: gamg_kokkos
315: nsize: {{1 2}separate output}
316: requires: kokkos_kernels
317: args: -mat_type aijkokkos
319: TEST*/