Actual source code: ex62.c
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";
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).
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.
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: */
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>
36: PetscErrorCode AssembleMatrix(Mat, PetscInt m, PetscInt n);
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;
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();
72: /* -------------------------------------------------------------------
73: Compute the matrix and right-hand-side vector that define
74: the linear system, Ax = b.
75: ------------------------------------------------------------------- */
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));
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));
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 matrix from which the preconditioner is constructed.
108: */
109: PetscCall(KSPSetOperators(ksp, A, A));
111: /*
112: Set the default preconditioner for this program to be GASM
113: */
114: PetscCall(KSPGetPC(ksp, &pc));
115: PetscCall(PCSetType(pc, PCGASM));
117: /* -------------------------------------------------------------------
118: Define the problem decomposition
119: ------------------------------------------------------------------- */
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Basic method, should be sufficient for the needs of many users.
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Set the overlap, using the default PETSc decomposition via
126: PCGASMSetOverlap(pc,overlap);
127: Could instead use the option -pc_gasm_overlap <ovl>
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: */
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: More advanced method, setting user-defined subdomains
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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.
145: Then call PCGASMSetLocalSubdomains() with resulting index sets
146: to set the subdomains for the GASM preconditioner.
147: */
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: 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));
157: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Outer IS:\n"));
158: for (PetscInt i = 0; i < Nsub; i++) {
159: PetscCall(PetscPrintf(PETSC_COMM_SELF, " outer IS[%" PetscInt_FMT "]\n", i));
160: PetscCall(ISView(outeris[i], PETSC_VIEWER_STDOUT_SELF));
161: }
162: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Inner IS:\n"));
163: for (PetscInt i = 0; i < Nsub; i++) {
164: PetscCall(PetscPrintf(PETSC_COMM_SELF, " inner IS[%" PetscInt_FMT "]\n", i));
165: PetscCall(ISView(inneris[i], PETSC_VIEWER_STDOUT_SELF));
166: }
167: }
168: } else { /* basic setup */
169: PetscCall(KSPSetFromOptions(ksp));
170: }
172: /* -------------------------------------------------------------------
173: Set the linear solvers for the subblocks
174: ------------------------------------------------------------------- */
176: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: Basic method, should be sufficient for the needs of most users.
178: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: By default, the GASM preconditioner uses the same solver on each
181: block of the problem. To set the same solver options on all blocks,
182: use the prefix -sub before the usual PC and KSP options, e.g.,
183: -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Advanced method, setting different solvers for various blocks.
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Note that each block's KSP context is completely independent of
190: the others, and the full range of uniprocessor KSP options is
191: available for each block.
193: - Use PCGASMGetSubKSP() to extract the array of KSP contexts for
194: the local blocks.
195: - See ex7.c for a simple example of setting different linear solvers
196: for the individual blocks for the block Jacobi method (which is
197: equivalent to the GASM method with zero overlap).
198: */
200: flg = PETSC_FALSE;
201: PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_set_subdomain_solvers", &flg, NULL));
202: if (flg) {
203: KSP *subksp; /* array of KSP contexts for local subblocks */
204: PetscInt i, nlocal, first; /* number of local subblocks, first local subblock */
205: PC subpc; /* PC context for subblock */
206: PetscBool isasm;
208: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "User explicitly sets subdomain solvers.\n"));
210: /*
211: Set runtime options
212: */
213: PetscCall(KSPSetFromOptions(ksp));
215: /*
216: Flag an error if PCTYPE is changed from the runtime options
217: */
218: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &isasm));
219: PetscCheck(isasm, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Cannot Change the PCTYPE when manually changing the subdomain solver settings");
221: /*
222: Call KSPSetUp() to set the block Jacobi data structures (including
223: creation of an internal KSP context for each block).
225: Note: KSPSetUp() MUST be called before PCGASMGetSubKSP().
226: */
227: PetscCall(KSPSetUp(ksp));
229: /*
230: Extract the array of KSP contexts for the local blocks
231: */
232: PetscCall(PCGASMGetSubKSP(pc, &nlocal, &first, &subksp));
234: /*
235: Loop over the local blocks, setting various KSP options
236: for each block.
237: */
238: for (i = 0; i < nlocal; i++) {
239: PetscCall(KSPGetPC(subksp[i], &subpc));
240: PetscCall(PCSetType(subpc, PCILU));
241: PetscCall(KSPSetType(subksp[i], KSPGMRES));
242: PetscCall(KSPSetTolerances(subksp[i], 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
243: }
244: } else {
245: /*
246: Set runtime options
247: */
248: PetscCall(KSPSetFromOptions(ksp));
249: }
251: /* -------------------------------------------------------------------
252: Solve the linear system
253: ------------------------------------------------------------------- */
255: PetscCall(KSPSolve(ksp, b, x));
257: /* -------------------------------------------------------------------
258: Assemble the matrix again to test repeated setup and solves.
259: ------------------------------------------------------------------- */
261: PetscCall(AssembleMatrix(A, m, n));
262: PetscCall(KSPSolve(ksp, b, x));
264: /* -------------------------------------------------------------------
265: Compare result to the exact solution
266: ------------------------------------------------------------------- */
267: PetscCall(VecAXPY(x, -1.0, u));
268: PetscCall(VecNorm(x, NORM_INFINITY, &e));
270: flg = PETSC_FALSE;
271: PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_error", &flg, NULL));
272: if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e));
274: /*
275: Free work space. All PETSc objects should be destroyed when they
276: are no longer needed.
277: */
279: PetscCall(KSPDestroy(&ksp));
280: PetscCall(VecDestroy(&u));
281: PetscCall(VecDestroy(&x));
282: PetscCall(VecDestroy(&b));
283: PetscCall(MatDestroy(&A));
284: PetscCall(PetscFinalize());
285: return 0;
286: }
288: PetscErrorCode AssembleMatrix(Mat A, PetscInt m, PetscInt n)
289: {
290: PetscInt i, j, Ii, J, Istart, Iend;
291: PetscScalar v;
293: PetscFunctionBegin;
294: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
295: for (Ii = Istart; Ii < Iend; Ii++) {
296: v = -1.0;
297: i = Ii / n;
298: j = Ii - i * n;
299: if (i > 0) {
300: J = Ii - n;
301: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
302: }
303: if (i < m - 1) {
304: J = Ii + n;
305: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
306: }
307: if (j > 0) {
308: J = Ii - 1;
309: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
310: }
311: if (j < n - 1) {
312: J = Ii + 1;
313: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
314: }
315: v = 4.0;
316: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
317: }
318: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
319: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*TEST
325: test:
326: suffix: 2D_1
327: args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains
329: test:
330: suffix: 2D_2
331: nsize: 2
332: args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains
334: test:
335: suffix: 2D_3
336: nsize: 3
337: args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains
339: test:
340: suffix: hp
341: nsize: 4
342: requires: superlu_dist
343: 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
344: output_file: output/empty.out
345: TODO: bug, triggers New nonzero at (0,15) caused a malloc in MatCreateSubMatrices_MPIAIJ_SingleIS_Local
347: test:
348: suffix: superlu_dist_1
349: requires: superlu_dist
350: 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: test:
353: suffix: superlu_dist_2
354: nsize: 2
355: requires: superlu_dist
356: 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: test:
359: suffix: superlu_dist_3
360: nsize: 3
361: requires: superlu_dist
362: 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: test:
365: suffix: superlu_dist_4
366: nsize: 4
367: requires: superlu_dist
368: 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: test:
371: suffix: gasm_view
372: nsize: 8
373: args: -pc_type gasm -ksp_view -Mdomains 2 -Ndomains 2 -user_set_subdomains
375: TEST*/