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