Actual source code: ex11.c
1: static char help[] = "Solves a linear system in parallel with KSP.\n\n";
3: /*
4: Description: Solves a complex linear system in parallel with KSP.
6: The model problem:
7: Solve Helmholtz equation on the unit square: (0,1) x (0,1)
8: -delta u - sigma1*u + i*sigma2*u = f,
9: where delta = Laplace operator
10: Dirichlet b.c.'s on all sides
11: Use the 2-D, five-point finite difference stencil.
13: Compiling the code:
14: This code uses the complex numbers version of PETSc, so configure
15: must be run to enable this
16: */
18: /*
19: Include "petscksp.h" so that we can use KSP solvers. Note that this file
20: automatically includes:
21: petscsys.h - base PETSc routines petscvec.h - vectors
22: petscmat.h - matrices
23: petscis.h - index sets petscksp.h - Krylov subspace methods
24: petscviewer.h - viewers petscpc.h - preconditioners
25: */
26: #include <petscksp.h>
28: int main(int argc, char **args)
29: {
30: Vec x, b, u; /* approx solution, RHS, exact solution */
31: Mat A; /* linear system matrix */
32: KSP ksp; /* linear solver context */
33: PetscReal norm; /* norm of solution error */
34: PetscInt dim, i, j, Ii, J, Istart, Iend, n = 6, its, use_random;
35: PetscScalar v, none = -1.0, sigma2, pfive = 0.5, *xa;
36: PetscRandom rctx;
37: PetscReal h2, sigma1 = 100.0;
38: PetscBool flg = PETSC_FALSE;
40: PetscFunctionBeginUser;
41: PetscCall(PetscInitialize(&argc, &args, NULL, help));
42: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
43: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
44: dim = n * n;
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Compute the matrix and right-hand-side vector that define
48: the linear system, Ax = b.
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: /*
51: Create parallel matrix, specifying only its global dimensions.
52: When using MatCreate(), the matrix format can be specified at
53: runtime. Also, the parallel partitioning of the matrix is
54: determined by PETSc at runtime.
55: */
56: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
57: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
58: PetscCall(MatSetFromOptions(A));
59: PetscCall(MatSetUp(A));
61: /*
62: Currently, all PETSc parallel matrix formats are partitioned by
63: contiguous chunks of rows across the processors. Determine which
64: rows of the matrix are locally owned.
65: */
66: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
68: /*
69: Set matrix elements in parallel.
70: - Each processor needs to insert only elements that it owns
71: locally (but any non-local elements will be sent to the
72: appropriate processor during matrix assembly).
73: - Always specify global rows and columns of matrix entries.
74: */
76: PetscCall(PetscOptionsGetBool(NULL, NULL, "-norandom", &flg, NULL));
77: if (flg) use_random = 0;
78: else use_random = 1;
79: if (use_random) {
80: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
81: PetscCall(PetscRandomSetFromOptions(rctx));
82: PetscCall(PetscRandomSetInterval(rctx, 0.0, PETSC_i));
83: } else sigma2 = 10.0 * PETSC_i;
84: h2 = 1.0 / ((n + 1) * (n + 1));
85: for (Ii = Istart; Ii < Iend; Ii++) {
86: v = -1.0;
87: i = Ii / n;
88: j = Ii - i * n;
89: if (i > 0) {
90: J = Ii - n;
91: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
92: }
93: if (i < n - 1) {
94: J = Ii + n;
95: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
96: }
97: if (j > 0) {
98: J = Ii - 1;
99: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
100: }
101: if (j < n - 1) {
102: J = Ii + 1;
103: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
104: }
105: if (use_random) PetscCall(PetscRandomGetValue(rctx, &sigma2));
106: v = 4.0 - sigma1 * h2 + sigma2 * h2;
107: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
108: }
109: if (use_random) PetscCall(PetscRandomDestroy(&rctx));
111: /*
112: Assemble matrix, using the 2-step process:
113: MatAssemblyBegin(), MatAssemblyEnd()
114: Computations can be done while messages are in transition
115: by placing code between these two statements.
116: */
117: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
118: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
120: /*
121: Create parallel vectors.
122: - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
123: we specify only the vector's global
124: dimension; the parallel partitioning is determined at runtime.
125: - Note: We form 1 vector from scratch and then duplicate as needed.
126: */
127: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
128: PetscCall(VecSetSizes(u, PETSC_DECIDE, dim));
129: PetscCall(VecSetFromOptions(u));
130: PetscCall(VecDuplicate(u, &b));
131: PetscCall(VecDuplicate(b, &x));
133: /*
134: Set exact solution; then compute right-hand-side vector.
135: */
137: if (use_random) {
138: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
139: PetscCall(PetscRandomSetFromOptions(rctx));
140: PetscCall(VecSetRandom(u, rctx));
141: } else PetscCall(VecSet(u, pfive));
142: PetscCall(MatMult(A, u, b));
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Create the linear solver and set various options
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: /*
149: Create linear solver context
150: */
151: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
153: /*
154: Set operators. Here the matrix that defines the linear system
155: also serves as the matrix from which the preconditioner is constructed.
156: */
157: PetscCall(KSPSetOperators(ksp, A, A));
159: /*
160: Set runtime options, e.g.,
161: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
162: */
163: PetscCall(KSPSetFromOptions(ksp));
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Solve the linear system
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: PetscCall(KSPSolve(ksp, b, x));
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Check solution and clean up
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: /*
176: Print the first 3 entries of x; this demonstrates extraction of the
177: real and imaginary components of the complex vector, x.
178: */
179: flg = PETSC_FALSE;
180: PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_x3", &flg, NULL));
181: if (flg) {
182: PetscCall(VecGetArray(x, &xa));
183: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "The first three entries of x are:\n"));
184: for (i = 0; i < 3; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x[%" PetscInt_FMT "] = %g + %g i\n", i, (double)PetscRealPart(xa[i]), (double)PetscImaginaryPart(xa[i])));
185: PetscCall(VecRestoreArray(x, &xa));
186: }
188: /*
189: Check the error
190: */
191: PetscCall(VecAXPY(x, none, u));
192: PetscCall(VecNorm(x, NORM_2, &norm));
193: PetscCall(KSPGetIterationNumber(ksp, &its));
194: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));
196: /*
197: Free work space. All PETSc objects should be destroyed when they
198: are no longer needed.
199: */
200: PetscCall(KSPDestroy(&ksp));
201: if (use_random) PetscCall(PetscRandomDestroy(&rctx));
202: PetscCall(VecDestroy(&u));
203: PetscCall(VecDestroy(&x));
204: PetscCall(VecDestroy(&b));
205: PetscCall(MatDestroy(&A));
206: PetscCall(PetscFinalize());
207: return 0;
208: }
210: /*TEST
212: build:
213: requires: complex
215: test:
216: args: -n 6 -norandom -pc_type none -ksp_monitor -ksp_gmres_cgs_refinement_type refine_always
218: testset:
219: suffix: deflation
220: args: -norandom -pc_type deflation -ksp_monitor
221: requires: superlu_dist
223: test:
224: nsize: 6
226: test:
227: nsize: 3
228: args: -pc_deflation_compute_space {{db2 aggregation}}
230: test:
231: suffix: pc_deflation_init_only-0
232: nsize: 4
233: args: -ksp_type fgmres -pc_deflation_compute_space db4 -pc_deflation_compute_space_size 2 -pc_deflation_levels 2 -deflation_ksp_max_it 10
234: #TODO remove suffix and next test when this works
235: #args: -pc_deflation_init_only {{0 1}separate output}
236: args: -pc_deflation_init_only 0
238: test:
239: suffix: pc_deflation_init_only-1
240: nsize: 4
241: args: -ksp_type fgmres -pc_deflation_compute_space db4 -pc_deflation_compute_space_size 2 -pc_deflation_levels 2 -deflation_ksp_max_it 10
242: args: -pc_deflation_init_only 1
244: TEST*/