Actual source code: ex7.c
petsc-master 2021-01-19
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);CHKERRMPI(ierr);
51: MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
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);
78: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
80: /*
81: Create parallel vectors
82: */
83: MatCreateVecs(A,&u,&b);
84: VecDuplicate(u,&x);
86: /*
87: Set exact solution; then compute right-hand-side vector.
88: */
89: VecSet(u,one);
90: MatMult(A,u,b);
92: /*
93: Create linear solver context
94: */
95: KSPCreate(PETSC_COMM_WORLD,&ksp);
97: /*
98: Set operators. Here the matrix that defines the linear system
99: also serves as the preconditioning matrix.
100: */
101: KSPSetOperators(ksp,A,A);
103: /*
104: Set default preconditioner for this program to be block Jacobi.
105: This choice can be overridden at runtime with the option
106: -pc_type <type>
108: IMPORTANT NOTE: Since the inners solves below are constructed to use
109: iterative methods (such as KSPGMRES) the outer Krylov method should
110: be set to use KSPFGMRES since it is the only Krylov method (plus KSPFCG)
111: that allows the preconditioners to be nonlinear (that is have iterative methods
112: inside them). The reason these examples work is because the number of
113: iterations on the inner solves is left at the default (which is 10,000)
114: and the tolerance on the inner solves is set to be a tight value of around 10^-6.
115: */
116: KSPGetPC(ksp,&pc);
117: PCSetType(pc,PCBJACOBI);
120: /* -------------------------------------------------------------------
121: Define the problem decomposition
122: ------------------------------------------------------------------- */
124: /*
125: Call PCBJacobiSetTotalBlocks() to set individually the size of
126: each block in the preconditioner. This could also be done with
127: the runtime option
128: -pc_bjacobi_blocks <blocks>
129: Also, see the command PCBJacobiSetLocalBlocks() to set the
130: local blocks.
132: Note: The default decomposition is 1 block per processor.
133: */
134: PetscMalloc1(m,&blks);
135: for (i=0; i<m; i++) blks[i] = n;
136: PCBJacobiSetTotalBlocks(pc,m,blks);
137: PetscFree(blks);
139: /*
140: Set runtime options
141: */
142: KSPSetFromOptions(ksp);
145: /* -------------------------------------------------------------------
146: Set the linear solvers for the subblocks
147: ------------------------------------------------------------------- */
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Basic method, should be sufficient for the needs of most users.
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: By default, the block Jacobi method uses the same solver on each
154: block of the problem. To set the same solver options on all blocks,
155: use the prefix -sub before the usual PC and KSP options, e.g.,
156: -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4
157: */
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Advanced method, setting different solvers for various blocks.
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Note that each block's KSP context is completely independent of
164: the others, and the full range of uniprocessor KSP options is
165: available for each block. The following section of code is intended
166: to be a simple illustration of setting different linear solvers for
167: the individual blocks. These choices are obviously not recommended
168: for solving this particular problem.
169: */
170: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
171: if (isbjacobi) {
172: /*
173: Call KSPSetUp() to set the block Jacobi data structures (including
174: creation of an internal KSP context for each block).
176: Note: KSPSetUp() MUST be called before PCBJacobiGetSubKSP().
177: */
178: KSPSetUp(ksp);
180: /*
181: Extract the array of KSP contexts for the local blocks
182: */
183: PCBJacobiGetSubKSP(pc,&nlocal,&first,&subksp);
185: /*
186: Loop over the local blocks, setting various KSP options
187: for each block.
188: */
189: for (i=0; i<nlocal; i++) {
190: KSPGetPC(subksp[i],&subpc);
191: if (!rank) {
192: if (i%2) {
193: PCSetType(subpc,PCILU);
194: } else {
195: PCSetType(subpc,PCNONE);
196: KSPSetType(subksp[i],KSPBCGS);
197: KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
198: }
199: } else {
200: PCSetType(subpc,PCJACOBI);
201: KSPSetType(subksp[i],KSPGMRES);
202: KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
203: }
204: }
205: }
207: /* -------------------------------------------------------------------
208: Solve the linear system
209: ------------------------------------------------------------------- */
211: /*
212: Solve the linear system
213: */
214: KSPSolve(ksp,b,x);
216: /* -------------------------------------------------------------------
217: Check solution and clean up
218: ------------------------------------------------------------------- */
220: /*
221: Check the error
222: */
223: VecAXPY(x,none,u);
224: VecNorm(x,NORM_2,&norm);
225: KSPGetIterationNumber(ksp,&its);
226: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
228: /*
229: Free work space. All PETSc objects should be destroyed when they
230: are no longer needed.
231: */
232: KSPDestroy(&ksp);
233: VecDestroy(&u); VecDestroy(&x);
234: VecDestroy(&b); MatDestroy(&A);
235: PetscFinalize();
236: return ierr;
237: }
240: /*TEST
242: test:
243: suffix: 1
244: nsize: 2
245: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
247: test:
248: suffix: 2
249: nsize: 2
250: args: -ksp_view
252: test:
253: suffix: viennacl
254: requires: viennacl
255: args: -ksp_monitor_short -mat_type aijviennacl
256: output_file: output/ex7_mpiaijcusparse.out
258: test:
259: suffix: viennacl_2
260: nsize: 2
261: requires: viennacl
262: args: -ksp_monitor_short -mat_type aijviennacl
263: output_file: output/ex7_mpiaijcusparse_2.out
265: test:
266: suffix: mpiaijcusparse
267: requires: cuda
268: args: -ksp_monitor_short -mat_type aijcusparse
270: test:
271: suffix: mpiaijcusparse_2
272: nsize: 2
273: requires: cuda
274: args: -ksp_monitor_short -mat_type aijcusparse
276: test:
277: suffix: mpiaijcusparse_simple
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_simple_2
283: nsize: 2
284: requires: cuda
285: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -sub_ksp_type preonly -sub_pc_type ilu
287: test:
288: suffix: mpiaijcusparse_3
289: requires: cuda
290: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
292: test:
293: suffix: mpiaijcusparse_4
294: nsize: 2
295: requires: cuda
296: args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
298: testset:
299: args: -ksp_monitor_short -pc_type gamg -ksp_view
300: test:
301: suffix: gamg_cuda
302: nsize: {{1 2}separate output}
303: requires: cuda
304: args: -mat_type aijcusparse
305: # triggers cusparse MatTransposeMat operation when squaring the graph
306: args: -pc_gamg_sym_graph 0 -pc_gamg_threshold -1 -pc_gamg_square_graph 1
307: test:
308: suffix: gamg_kokkos
309: nsize: {{1 2}separate output}
310: requires: kokkos_kernels
311: args: -mat_type aijkokkos
313: TEST*/