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