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