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