Actual source code: ex64.c
petsc-3.9.4 2018-09-11
1: /*
2: * ex64.c
3: *
4: * Created on: Aug 10, 2015
5: * Author: Fande Kong <fdkong.jd@gmail.com>
6: */
8: static char help[] = "Illustrates use of the preconditioner GASM.\n \
9: using hierarchical partitioning and MatIncreaseOverlapSplit \
10: -pc_gasm_total_subdomains\n \
11: -pc_gasm_print_subdomains\n \n";
13: /*
14: Note: This example focuses on setting the subdomains for the GASM
15: preconditioner for a problem on a 2D rectangular grid. See ex1.c
16: and ex2.c for more detailed comments on the basic usage of KSP
17: (including working with matrices and vectors).
19: The GASM preconditioner is fully parallel. The user-space routine
20: CreateSubdomains2D that computes the domain decomposition is also parallel
21: and attempts to generate both subdomains straddling processors and multiple
22: domains per processor.
25: This matrix in this linear system arises from the discretized Laplacian,
26: and thus is not very interesting in terms of experimenting with variants
27: of the GASM preconditioner.
28: */
30: /*T
31: Concepts: KSP^Additive Schwarz Method (GASM) with user-defined subdomains
32: Processors: n
33: T*/
35: /*
36: Include "petscksp.h" so that we can use KSP solvers. Note that this file
37: automatically includes:
38: petscsys.h - base PETSc routines petscvec.h - vectors
39: petscmat.h - matrices
40: petscis.h - index sets petscksp.h - Krylov subspace methods
41: petscviewer.h - viewers petscpc.h - preconditioners
42: */
43: #include <petscksp.h>
44: #include <petscmat.h>
47: int main(int argc,char **args)
48: {
49: Vec x,b,u; /* approx solution, RHS, exact solution */
50: Mat A,perA; /* linear system matrix */
51: KSP ksp; /* linear solver context */
52: PC pc; /* PC context */
53: PetscInt overlap; /* width of subdomain overlap */
54: PetscInt m,n; /* mesh dimensions in x- and y- directions */
55: PetscInt M,N; /* number of subdomains in x- and y- directions */
56: PetscInt i,j,Ii,J,Istart,Iend;
58: PetscMPIInt size;
59: PetscBool flg;
60: PetscBool user_set_subdomains;
61: PetscScalar v, one = 1.0;
62: MatPartitioning part;
63: IS coarseparts,fineparts;
64: IS is,isn,isrows;
65: MPI_Comm comm;
66: PetscReal e;
68: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
69: comm = PETSC_COMM_WORLD;
70: MPI_Comm_size(comm,&size);
71: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex62","PC");
72: m = 15;
73: PetscOptionsInt("-M", "Number of mesh points in the x-direction","PCGASMCreateSubdomains2D",m,&m,NULL);
74: n = 17;
75: PetscOptionsInt("-N","Number of mesh points in the y-direction","PCGASMCreateSubdomains2D",n,&n,NULL);
76: user_set_subdomains = PETSC_FALSE;
77: PetscOptionsBool("-user_set_subdomains","Use the user-specified 2D tiling of mesh by subdomains","PCGASMCreateSubdomains2D",user_set_subdomains,&user_set_subdomains,NULL);
78: M = 2;
79: PetscOptionsInt("-Mdomains","Number of subdomain tiles in the x-direction","PCGASMSetSubdomains2D",M,&M,NULL);
80: N = 1;
81: PetscOptionsInt("-Ndomains","Number of subdomain tiles in the y-direction","PCGASMSetSubdomains2D",N,&N,NULL);
82: overlap = 1;
83: PetscOptionsInt("-overlap","Size of tile overlap.","PCGASMSetSubdomains2D",overlap,&overlap,NULL);
84: PetscOptionsEnd();
86: /* -------------------------------------------------------------------
87: Compute the matrix and right-hand-side vector that define
88: the linear system, Ax = b.
89: ------------------------------------------------------------------- */
91: /*
92: Assemble the matrix for the five point stencil, YET AGAIN
93: */
94: MatCreate(comm,&A);
95: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
96: MatSetFromOptions(A);
97: MatSetUp(A);
98: MatGetOwnershipRange(A,&Istart,&Iend);
99: for (Ii=Istart; Ii<Iend; Ii++) {
100: v = -1.0; i = Ii/n; j = Ii - i*n;
101: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
102: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
103: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
104: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
105: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
106: }
107: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
108: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
110: /*
111: Partition the graph of the matrix
112: */
113: MatPartitioningCreate(comm,&part);
114: MatPartitioningSetAdjacency(part,A);
115: MatPartitioningSetType(part,MATPARTITIONINGHIERARCH);
116: MatPartitioningHierarchicalSetNcoarseparts(part,2);
117: MatPartitioningHierarchicalSetNfineparts(part,2);
118: MatPartitioningSetFromOptions(part);
119: /* get new processor owner number of each vertex */
120: MatPartitioningApply(part,&is);
121: /* get coarse parts according to which we rearrange the matrix */
122: MatPartitioningHierarchicalGetCoarseparts(part,&coarseparts);
123: /* fine parts are used with coarse parts */
124: MatPartitioningHierarchicalGetFineparts(part,&fineparts);
125: /* get new global number of each old global number */
126: ISPartitioningToNumbering(is,&isn);
127: ISBuildTwoSided(is,NULL,&isrows);
128: MatCreateSubMatrix(A,isrows,isrows,MAT_INITIAL_MATRIX,&perA);
129: MatCreateVecs(perA,&b,NULL);
130: VecSetFromOptions(b);
131: VecDuplicate(b,&u);
132: VecDuplicate(b,&x);
133: VecSet(u,one);
134: MatMult(perA,u,b);
135: KSPCreate(comm,&ksp);
136: /*
137: Set operators. Here the matrix that defines the linear system
138: also serves as the preconditioning matrix.
139: */
140: KSPSetOperators(ksp,perA,perA);
142: /*
143: Set the default preconditioner for this program to be GASM
144: */
145: KSPGetPC(ksp,&pc);
146: PCSetType(pc,PCGASM);
147: KSPSetFromOptions(ksp);
148: /* -------------------------------------------------------------------
149: Solve the linear system
150: ------------------------------------------------------------------- */
152: KSPSolve(ksp,b,x);
153: /* -------------------------------------------------------------------
154: Compare result to the exact solution
155: ------------------------------------------------------------------- */
156: VecAXPY(x,-1.0,u);
157: VecNorm(x,NORM_INFINITY, &e);
159: flg = PETSC_FALSE;
160: PetscOptionsGetBool(NULL,NULL,"-print_error",&flg,NULL);
161: if (flg) {
162: PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e);
163: }
164: /*
165: Free work space. All PETSc objects should be destroyed when they
166: are no longer needed.
167: */
168: KSPDestroy(&ksp);
169: VecDestroy(&u);
170: VecDestroy(&x);
171: VecDestroy(&b);
172: MatDestroy(&A);
173: MatDestroy(&perA);
174: ISDestroy(&is);
175: ISDestroy(&coarseparts);
176: ISDestroy(&fineparts);
177: ISDestroy(&isrows);
178: ISDestroy(&isn);
179: MatPartitioningDestroy(&part);
180: PetscFinalize();
181: return ierr;
182: }
185: /*TEST
187: test:
188: nsize: 4
189: requires: superlu_dist parmetis
190: args: -ksp_monitor -pc_gasm_overlap 1 -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist -pc_gasm_total_subdomains 2
192: TEST*/