Actual source code: ex64.c
petsc-3.7.5 2017-01-01
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>
49: int main(int argc,char **args)
50: {
51: Vec x,b,u; /* approx solution, RHS, exact solution */
52: Mat A,perA; /* linear system matrix */
53: KSP ksp; /* linear solver context */
54: PC pc; /* PC context */
55: PetscInt overlap; /* width of subdomain overlap */
56: PetscInt m,n; /* mesh dimensions in x- and y- directions */
57: PetscInt M,N; /* number of subdomains in x- and y- directions */
58: PetscInt i,j,Ii,J,Istart,Iend;
60: PetscMPIInt size;
61: PetscBool flg;
62: PetscBool user_set_subdomains;
63: PetscScalar v, one = 1.0;
64: MatPartitioning part;
65: IS coarseparts,fineparts;
66: IS is,isn,isrows;
67: MPI_Comm comm;
68: PetscReal e;
70: PetscInitialize(&argc,&args,(char*)0,help);
71: comm = PETSC_COMM_WORLD;
72: MPI_Comm_size(comm,&size);
73: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex62","PC");
74: m = 15;
75: PetscOptionsInt("-M", "Number of mesh points in the x-direction","PCGASMCreateSubdomains2D",m,&m,NULL);
76: n = 17;
77: PetscOptionsInt("-N","Number of mesh points in the y-direction","PCGASMCreateSubdomains2D",n,&n,NULL);
78: user_set_subdomains = PETSC_FALSE;
79: PetscOptionsBool("-user_set_subdomains","Use the user-specified 2D tiling of mesh by subdomains","PCGASMCreateSubdomains2D",user_set_subdomains,&user_set_subdomains,NULL);
80: M = 2;
81: PetscOptionsInt("-Mdomains","Number of subdomain tiles in the x-direction","PCGASMSetSubdomains2D",M,&M,NULL);
82: N = 1;
83: PetscOptionsInt("-Ndomains","Number of subdomain tiles in the y-direction","PCGASMSetSubdomains2D",N,&N,NULL);
84: overlap = 1;
85: PetscOptionsInt("-overlap","Size of tile overlap.","PCGASMSetSubdomains2D",overlap,&overlap,NULL);
86: PetscOptionsEnd();
88: /* -------------------------------------------------------------------
89: Compute the matrix and right-hand-side vector that define
90: the linear system, Ax = b.
91: ------------------------------------------------------------------- */
93: /*
94: Assemble the matrix for the five point stencil, YET AGAIN
95: */
96: MatCreate(comm,&A);
97: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
98: MatSetFromOptions(A);
99: MatSetUp(A);
100: MatGetOwnershipRange(A,&Istart,&Iend);
101: for (Ii=Istart; Ii<Iend; Ii++) {
102: v = -1.0; i = Ii/n; j = Ii - i*n;
103: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
104: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
105: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
106: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
107: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
108: }
109: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
110: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
112: /*
113: Partition the graph of the matrix
114: */
115: MatPartitioningCreate(comm,&part);
116: MatPartitioningSetAdjacency(part,A);
117: MatPartitioningSetType(part,MATPARTITIONINGHIERARCH);
118: MatPartitioningHierarchicalSetNcoarseparts(part,2);
119: MatPartitioningHierarchicalSetNfineparts(part,2);
120: MatPartitioningSetFromOptions(part);
121: /* get new processor owner number of each vertex */
122: MatPartitioningApply(part,&is);
123: /* get coarse parts according to which we rearrange the matrix */
124: MatPartitioningHierarchicalGetCoarseparts(part,&coarseparts);
125: /* fine parts are used with coarse parts */
126: MatPartitioningHierarchicalGetFineparts(part,&fineparts);
127: /* get new global number of each old global number */
128: ISPartitioningToNumbering(is,&isn);
129: ISBuildTwoSided(is,NULL,&isrows);
130: MatGetSubMatrix(A,isrows,isrows,MAT_INITIAL_MATRIX,&perA);
131: MatCreateVecs(perA,&b,NULL);
132: VecSetFromOptions(b);
133: VecDuplicate(b,&u);
134: VecDuplicate(b,&x);
135: VecSet(u,one);
136: MatMult(perA,u,b);
137: KSPCreate(comm,&ksp);
138: /*
139: Set operators. Here the matrix that defines the linear system
140: also serves as the preconditioning matrix.
141: */
142: KSPSetOperators(ksp,perA,perA);
144: /*
145: Set the default preconditioner for this program to be GASM
146: */
147: KSPGetPC(ksp,&pc);
148: PCSetType(pc,PCGASM);
149: KSPSetFromOptions(ksp);
150: /* -------------------------------------------------------------------
151: Solve the linear system
152: ------------------------------------------------------------------- */
154: KSPSolve(ksp,b,x);
155: /* -------------------------------------------------------------------
156: Compare result to the exact solution
157: ------------------------------------------------------------------- */
158: VecAXPY(x,-1.0,u);
159: VecNorm(x,NORM_INFINITY, &e);
161: flg = PETSC_FALSE;
162: PetscOptionsGetBool(NULL,NULL,"-print_error",&flg,NULL);
163: if (flg) {
164: PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e);
165: }
166: /*
167: Free work space. All PETSc objects should be destroyed when they
168: are no longer needed.
169: */
170: KSPDestroy(&ksp);
171: VecDestroy(&u);
172: VecDestroy(&x);
173: VecDestroy(&b);
174: MatDestroy(&A);
175: MatDestroy(&perA);
176: ISDestroy(&is);
177: ISDestroy(&coarseparts);
178: ISDestroy(&fineparts);
179: ISDestroy(&isrows);
180: ISDestroy(&isn);
181: MatPartitioningDestroy(&part);
182: PetscFinalize();
183: return 0;
184: }