1 static char help[] = "Illustrates use of PCGASM.\n\
2 The Generalized Additive Schwarz Method for solving a linear system in parallel with KSP. The\n\
3 code indicates the procedure for setting user-defined subdomains.\n\
4 See section 'ex62' below for command-line options.\n\
5 Without -user_set_subdomains, the general PCGASM options are meaningful:\n\
6 -pc_gasm_total_subdomains\n\
7 -pc_gasm_print_subdomains\n\
8 \n";
9
10 /*
11 Note: This example focuses on setting the subdomains for the GASM
12 preconditioner for a problem on a 2D rectangular grid. See ex1.c
13 and ex2.c for more detailed comments on the basic usage of KSP
14 (including working with matrices and vectors).
15
16 The GASM preconditioner is fully parallel. The user-space routine
17 CreateSubdomains2D that computes the domain decomposition is also parallel
18 and attempts to generate both subdomains straddling processors and multiple
19 domains per processor.
20
21 This matrix in this linear system arises from the discretized Laplacian,
22 and thus is not very interesting in terms of experimenting with variants
23 of the GASM preconditioner.
24 */
25
26 /*
27 Include "petscksp.h" so that we can use KSP solvers. Note that this file
28 automatically includes:
29 petscsys.h - base PETSc routines petscvec.h - vectors
30 petscmat.h - matrices
31 petscis.h - index sets petscksp.h - Krylov subspace methods
32 petscviewer.h - viewers petscpc.h - preconditioners
33 */
34 #include <petscksp.h>
35
36 PetscErrorCode AssembleMatrix(Mat, PetscInt m, PetscInt n);
37
main(int argc,char ** args)38 int main(int argc, char **args)
39 {
40 Vec x, b, u; /* approx solution, RHS, exact solution */
41 Mat A; /* linear system matrix */
42 KSP ksp; /* linear solver context */
43 PC pc; /* PC context */
44 IS *inneris, *outeris; /* array of index sets that define the subdomains */
45 PetscInt overlap; /* width of subdomain overlap */
46 PetscInt Nsub; /* number of subdomains */
47 PetscInt m, n; /* mesh dimensions in x- and y- directions */
48 PetscInt M, N; /* number of subdomains in x- and y- directions */
49 PetscMPIInt size;
50 PetscBool flg = PETSC_FALSE;
51 PetscBool user_set_subdomains = PETSC_FALSE;
52 PetscReal one, e;
53
54 PetscFunctionBeginUser;
55 PetscCall(PetscInitialize(&argc, &args, NULL, help));
56 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
57 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex62", "PCGASM");
58 m = 15;
59 PetscCall(PetscOptionsInt("-M", "Number of mesh points in the x-direction", "PCGASMCreateSubdomains2D", m, &m, NULL));
60 n = 17;
61 PetscCall(PetscOptionsInt("-N", "Number of mesh points in the y-direction", "PCGASMCreateSubdomains2D", n, &n, NULL));
62 user_set_subdomains = PETSC_FALSE;
63 PetscCall(PetscOptionsBool("-user_set_subdomains", "Use the user-specified 2D tiling of mesh by subdomains", "PCGASMCreateSubdomains2D", user_set_subdomains, &user_set_subdomains, NULL));
64 M = 2;
65 PetscCall(PetscOptionsInt("-Mdomains", "Number of subdomain tiles in the x-direction", "PCGASMSetSubdomains2D", M, &M, NULL));
66 N = 1;
67 PetscCall(PetscOptionsInt("-Ndomains", "Number of subdomain tiles in the y-direction", "PCGASMSetSubdomains2D", N, &N, NULL));
68 overlap = 1;
69 PetscCall(PetscOptionsInt("-overlap", "Size of tile overlap.", "PCGASMSetSubdomains2D", overlap, &overlap, NULL));
70 PetscOptionsEnd();
71
72 /* -------------------------------------------------------------------
73 Compute the matrix and right-hand-side vector that define
74 the linear system, Ax = b.
75 ------------------------------------------------------------------- */
76
77 /*
78 Assemble the matrix for the five point stencil, YET AGAIN
79 */
80 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
81 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
82 PetscCall(MatSetFromOptions(A));
83 PetscCall(MatSetUp(A));
84 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
85 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
86 PetscCall(AssembleMatrix(A, m, n));
87
88 /*
89 Create and set vectors
90 */
91 PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
92 PetscCall(VecSetSizes(b, PETSC_DECIDE, m * n));
93 PetscCall(VecSetFromOptions(b));
94 PetscCall(VecDuplicate(b, &u));
95 PetscCall(VecDuplicate(b, &x));
96 one = 1.0;
97 PetscCall(VecSet(u, one));
98 PetscCall(MatMult(A, u, b));
99
100 /*
101 Create linear solver context
102 */
103 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
104
105 /*
106 Set operators. Here the matrix that defines the linear system
107 also serves as the matrix from which the preconditioner is constructed.
108 */
109 PetscCall(KSPSetOperators(ksp, A, A));
110
111 /*
112 Set the default preconditioner for this program to be GASM
113 */
114 PetscCall(KSPGetPC(ksp, &pc));
115 PetscCall(PCSetType(pc, PCGASM));
116
117 /* -------------------------------------------------------------------
118 Define the problem decomposition
119 ------------------------------------------------------------------- */
120
121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122 Basic method, should be sufficient for the needs of many users.
123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124
125 Set the overlap, using the default PETSc decomposition via
126 PCGASMSetOverlap(pc,overlap);
127 Could instead use the option -pc_gasm_overlap <ovl>
128
129 Set the total number of blocks via -pc_gasm_blocks <blks>
130 Note: The GASM default is to use 1 block per processor. To
131 experiment on a single processor with various overlaps, you
132 must specify use of multiple blocks!
133 */
134
135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136 More advanced method, setting user-defined subdomains
137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138
139 Firstly, create index sets that define the subdomains. The utility
140 routine PCGASMCreateSubdomains2D() is a simple example, which partitions
141 the 2D grid into MxN subdomains with an optional overlap.
142 More generally, the user should write a custom routine for a particular
143 problem geometry.
144
145 Then call PCGASMSetLocalSubdomains() with resulting index sets
146 to set the subdomains for the GASM preconditioner.
147 */
148
149 if (user_set_subdomains) { /* user-control version */
150 PetscCall(PCGASMCreateSubdomains2D(pc, m, n, M, N, 1, overlap, &Nsub, &inneris, &outeris));
151 PetscCall(PCGASMSetSubdomains(pc, Nsub, inneris, outeris));
152 PetscCall(PCGASMDestroySubdomains(Nsub, &inneris, &outeris));
153 flg = PETSC_FALSE;
154 PetscCall(PetscOptionsGetBool(NULL, NULL, "-subdomain_view", &flg, NULL));
155 if (flg) {
156 PetscInt i;
157 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Nmesh points: %" PetscInt_FMT " x %" PetscInt_FMT "; subdomain partition: %" PetscInt_FMT " x %" PetscInt_FMT "; overlap: %" PetscInt_FMT "; Nsub: %" PetscInt_FMT "\n", m, n, M, N, overlap, Nsub));
158 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Outer IS:\n"));
159 for (i = 0; i < Nsub; i++) {
160 PetscCall(PetscPrintf(PETSC_COMM_SELF, " outer IS[%" PetscInt_FMT "]\n", i));
161 PetscCall(ISView(outeris[i], PETSC_VIEWER_STDOUT_SELF));
162 }
163 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Inner IS:\n"));
164 for (i = 0; i < Nsub; i++) {
165 PetscCall(PetscPrintf(PETSC_COMM_SELF, " inner IS[%" PetscInt_FMT "]\n", i));
166 PetscCall(ISView(inneris[i], PETSC_VIEWER_STDOUT_SELF));
167 }
168 }
169 } else { /* basic setup */
170 PetscCall(KSPSetFromOptions(ksp));
171 }
172
173 /* -------------------------------------------------------------------
174 Set the linear solvers for the subblocks
175 ------------------------------------------------------------------- */
176
177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178 Basic method, should be sufficient for the needs of most users.
179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180
181 By default, the GASM preconditioner uses the same solver on each
182 block of the problem. To set the same solver options on all blocks,
183 use the prefix -sub before the usual PC and KSP options, e.g.,
184 -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4
185
186 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187 Advanced method, setting different solvers for various blocks.
188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189
190 Note that each block's KSP context is completely independent of
191 the others, and the full range of uniprocessor KSP options is
192 available for each block.
193
194 - Use PCGASMGetSubKSP() to extract the array of KSP contexts for
195 the local blocks.
196 - See ex7.c for a simple example of setting different linear solvers
197 for the individual blocks for the block Jacobi method (which is
198 equivalent to the GASM method with zero overlap).
199 */
200
201 flg = PETSC_FALSE;
202 PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_set_subdomain_solvers", &flg, NULL));
203 if (flg) {
204 KSP *subksp; /* array of KSP contexts for local subblocks */
205 PetscInt i, nlocal, first; /* number of local subblocks, first local subblock */
206 PC subpc; /* PC context for subblock */
207 PetscBool isasm;
208
209 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "User explicitly sets subdomain solvers.\n"));
210
211 /*
212 Set runtime options
213 */
214 PetscCall(KSPSetFromOptions(ksp));
215
216 /*
217 Flag an error if PCTYPE is changed from the runtime options
218 */
219 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &isasm));
220 PetscCheck(isasm, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Cannot Change the PCTYPE when manually changing the subdomain solver settings");
221
222 /*
223 Call KSPSetUp() to set the block Jacobi data structures (including
224 creation of an internal KSP context for each block).
225
226 Note: KSPSetUp() MUST be called before PCGASMGetSubKSP().
227 */
228 PetscCall(KSPSetUp(ksp));
229
230 /*
231 Extract the array of KSP contexts for the local blocks
232 */
233 PetscCall(PCGASMGetSubKSP(pc, &nlocal, &first, &subksp));
234
235 /*
236 Loop over the local blocks, setting various KSP options
237 for each block.
238 */
239 for (i = 0; i < nlocal; i++) {
240 PetscCall(KSPGetPC(subksp[i], &subpc));
241 PetscCall(PCSetType(subpc, PCILU));
242 PetscCall(KSPSetType(subksp[i], KSPGMRES));
243 PetscCall(KSPSetTolerances(subksp[i], 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
244 }
245 } else {
246 /*
247 Set runtime options
248 */
249 PetscCall(KSPSetFromOptions(ksp));
250 }
251
252 /* -------------------------------------------------------------------
253 Solve the linear system
254 ------------------------------------------------------------------- */
255
256 PetscCall(KSPSolve(ksp, b, x));
257
258 /* -------------------------------------------------------------------
259 Assemble the matrix again to test repeated setup and solves.
260 ------------------------------------------------------------------- */
261
262 PetscCall(AssembleMatrix(A, m, n));
263 PetscCall(KSPSolve(ksp, b, x));
264
265 /* -------------------------------------------------------------------
266 Compare result to the exact solution
267 ------------------------------------------------------------------- */
268 PetscCall(VecAXPY(x, -1.0, u));
269 PetscCall(VecNorm(x, NORM_INFINITY, &e));
270
271 flg = PETSC_FALSE;
272 PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_error", &flg, NULL));
273 if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e));
274
275 /*
276 Free work space. All PETSc objects should be destroyed when they
277 are no longer needed.
278 */
279
280 PetscCall(KSPDestroy(&ksp));
281 PetscCall(VecDestroy(&u));
282 PetscCall(VecDestroy(&x));
283 PetscCall(VecDestroy(&b));
284 PetscCall(MatDestroy(&A));
285 PetscCall(PetscFinalize());
286 return 0;
287 }
288
AssembleMatrix(Mat A,PetscInt m,PetscInt n)289 PetscErrorCode AssembleMatrix(Mat A, PetscInt m, PetscInt n)
290 {
291 PetscInt i, j, Ii, J, Istart, Iend;
292 PetscScalar v;
293
294 PetscFunctionBegin;
295 PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
296 for (Ii = Istart; Ii < Iend; Ii++) {
297 v = -1.0;
298 i = Ii / n;
299 j = Ii - i * n;
300 if (i > 0) {
301 J = Ii - n;
302 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
303 }
304 if (i < m - 1) {
305 J = Ii + n;
306 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
307 }
308 if (j > 0) {
309 J = Ii - 1;
310 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
311 }
312 if (j < n - 1) {
313 J = Ii + 1;
314 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
315 }
316 v = 4.0;
317 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
318 }
319 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
320 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
321 PetscFunctionReturn(PETSC_SUCCESS);
322 }
323
324 /*TEST
325
326 test:
327 suffix: 2D_1
328 args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains
329
330 test:
331 suffix: 2D_2
332 nsize: 2
333 args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains
334
335 test:
336 suffix: 2D_3
337 nsize: 3
338 args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains
339
340 test:
341 suffix: hp
342 nsize: 4
343 requires: superlu_dist
344 args: -M 7 -N 9 -pc_gasm_overlap 1 -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist -ksp_monitor -print_error -pc_gasm_total_subdomains 2 -pc_gasm_use_hierachical_partitioning 1
345 output_file: output/empty.out
346 TODO: bug, triggers New nonzero at (0,15) caused a malloc in MatCreateSubMatrices_MPIAIJ_SingleIS_Local
347
348 test:
349 suffix: superlu_dist_1
350 requires: superlu_dist
351 args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 1 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist
352
353 test:
354 suffix: superlu_dist_2
355 nsize: 2
356 requires: superlu_dist
357 args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 1 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist
358
359 test:
360 suffix: superlu_dist_3
361 nsize: 3
362 requires: superlu_dist
363 args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 2 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist
364
365 test:
366 suffix: superlu_dist_4
367 nsize: 4
368 requires: superlu_dist
369 args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 2 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist
370
371 test:
372 suffix: gasm_view
373 nsize: 8
374 args: -pc_type gasm -ksp_view -Mdomains 2 -Ndomains 2 -user_set_subdomains
375
376 TEST*/
377