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;
49: PetscMPIInt size;
50: PetscBool flg;
51: PetscBool user_set_subdomains;
52: PetscScalar v, one = 1.0;
53: MatPartitioning part;
54: IS coarseparts, fineparts;
55: IS is, isn, isrows;
56: MPI_Comm comm;
57: 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;
93: i = Ii / n;
94: j = Ii - i * n;
95: if (i > 0) {
96: J = Ii - n;
97: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
98: }
99: if (i < m - 1) {
100: J = Ii + n;
101: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
102: }
103: if (j > 0) {
104: J = Ii - 1;
105: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
106: }
107: if (j < n - 1) {
108: J = Ii + 1;
109: MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
110: }
111: v = 4.0;
112: MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
113: }
114: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
117: /*
118: Partition the graph of the matrix
119: */
120: MatPartitioningCreate(comm, &part);
121: MatPartitioningSetAdjacency(part, A);
122: MatPartitioningSetType(part, MATPARTITIONINGHIERARCH);
123: MatPartitioningHierarchicalSetNcoarseparts(part, 2);
124: MatPartitioningHierarchicalSetNfineparts(part, 2);
125: MatPartitioningSetFromOptions(part);
126: /* get new processor owner number of each vertex */
127: MatPartitioningApply(part, &is);
128: /* get coarse parts according to which we rearrange the matrix */
129: MatPartitioningHierarchicalGetCoarseparts(part, &coarseparts);
130: /* fine parts are used with coarse parts */
131: MatPartitioningHierarchicalGetFineparts(part, &fineparts);
132: /* get new global number of each old global number */
133: ISPartitioningToNumbering(is, &isn);
134: ISBuildTwoSided(is, NULL, &isrows);
135: MatCreateSubMatrix(A, isrows, isrows, MAT_INITIAL_MATRIX, &perA);
136: MatCreateVecs(perA, &b, NULL);
137: VecSetFromOptions(b);
138: VecDuplicate(b, &u);
139: VecDuplicate(b, &x);
140: VecSet(u, one);
141: MatMult(perA, u, b);
142: KSPCreate(comm, &ksp);
143: /*
144: Set operators. Here the matrix that defines the linear system
145: also serves as the preconditioning matrix.
146: */
147: KSPSetOperators(ksp, perA, perA);
149: /*
150: Set the default preconditioner for this program to be GASM
151: */
152: KSPGetPC(ksp, &pc);
153: PCSetType(pc, PCGASM);
154: KSPSetFromOptions(ksp);
155: /* -------------------------------------------------------------------
156: Solve the linear system
157: ------------------------------------------------------------------- */
159: KSPSolve(ksp, b, x);
160: /* -------------------------------------------------------------------
161: Compare result to the exact solution
162: ------------------------------------------------------------------- */
163: VecAXPY(x, -1.0, u);
164: VecNorm(x, NORM_INFINITY, &e);
166: flg = PETSC_FALSE;
167: PetscOptionsGetBool(NULL, NULL, "-print_error", &flg, NULL);
168: if (flg) PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e);
169: /*
170: Free work space. All PETSc objects should be destroyed when they
171: are no longer needed.
172: */
173: KSPDestroy(&ksp);
174: VecDestroy(&u);
175: VecDestroy(&x);
176: VecDestroy(&b);
177: MatDestroy(&A);
178: MatDestroy(&perA);
179: ISDestroy(&is);
180: ISDestroy(&coarseparts);
181: ISDestroy(&fineparts);
182: ISDestroy(&isrows);
183: ISDestroy(&isn);
184: MatPartitioningDestroy(&part);
185: PetscFinalize();
186: return 0;
187: }
189: /*TEST
191: test:
192: nsize: 4
193: requires: superlu_dist parmetis
194: args: -ksp_monitor -pc_gasm_overlap 1 -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist -pc_gasm_total_subdomains 2
196: TEST*/