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:   PetscFunctionBeginUser;
41:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
42:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
43:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
44:   PetscCallMPI(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:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
56:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
57:   PetscCall(MatSetFromOptions(A));
58:   PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
59:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
60:   PetscCall(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:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
68:     }
69:     if (i < m - 1) {
70:       J = Ii + n;
71:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
72:     }
73:     if (j > 0) {
74:       J = Ii - 1;
75:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
76:     }
77:     if (j < n - 1) {
78:       J = Ii + 1;
79:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
80:     }
81:     v = 4.0;
82:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
83:   }
84:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
85:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
86:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));

88:   /*
89:      Create parallel vectors
90:   */
91:   PetscCall(MatCreateVecs(A, &u, &b));
92:   PetscCall(VecDuplicate(u, &x));

94:   /*
95:      Set exact solution; then compute right-hand-side vector.
96:   */
97:   PetscCall(VecSet(u, one));
98:   PetscCall(MatMult(A, u, b));

100:   /*
101:      Create linear solver context
102:   */
103:   PetscCall(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:   PetscCall(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:   PetscCall(KSPGetPC(ksp, &pc));
125:   PetscCall(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:   PetscCall(PetscMalloc1(m, &blks));
142:   for (i = 0; i < m; i++) blks[i] = n;
143:   PetscCall(PCBJacobiSetTotalBlocks(pc, m, blks));
144:   PetscCall(PetscFree(blks));

146:   /*
147:     Set runtime options
148:   */
149:   PetscCall(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:   PetscCall(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:     PetscCall(KSPSetUp(ksp));

186:     /*
187:        Extract the array of KSP contexts for the local blocks
188:     */
189:     PetscCall(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:       PetscCall(KSPGetPC(subksp[i], &subpc));
197:       if (rank == 0) {
198:         if (i % 2) {
199:           PetscCall(PCSetType(subpc, PCILU));
200:         } else {
201:           PetscCall(PCSetType(subpc, PCNONE));
202:           PetscCall(KSPSetType(subksp[i], KSPBCGS));
203:           PetscCall(KSPSetTolerances(subksp[i], 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
204:         }
205:       } else {
206:         PetscCall(PCSetType(subpc, PCJACOBI));
207:         PetscCall(KSPSetType(subksp[i], KSPGMRES));
208:         PetscCall(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:   PetscCall(KSPSolve(ksp, b, x));

222:   /* -------------------------------------------------------------------
223:                       Check solution and clean up
224:      ------------------------------------------------------------------- */

226:   /*
227:      Check the error
228:   */
229:   PetscCall(VecAXPY(x, none, u));
230:   PetscCall(VecNorm(x, NORM_2, &norm));
231:   PetscCall(KSPGetIterationNumber(ksp, &its));
232:   PetscCall(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:   PetscCall(KSPDestroy(&ksp));
239:   PetscCall(VecDestroy(&u));
240:   PetscCall(VecDestroy(&x));
241:   PetscCall(VecDestroy(&b));
242:   PetscCall(MatDestroy(&A));
243:   PetscCall(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*/
```