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