Actual source code: ex2.c
1: static char help[] = "Solves a linear system in parallel with KSP.\n\
2: Input parameters include:\n\
3: -view_exact_sol : write exact solution vector to stdout\n\
4: -m <mesh_x> : number of mesh points in x-direction\n\
5: -n <mesh_y> : number of mesh points in y-direction\n\n";
7: /*
8: Include "petscksp.h" so that we can use KSP solvers.
9: */
10: #include <petscksp.h>
12: int main(int argc, char **args)
13: {
14: Vec x, b, u; /* approx solution, RHS, exact solution */
15: Mat A; /* linear system matrix */
16: KSP ksp; /* linear solver context */
17: PetscReal norm; /* norm of solution error */
18: PetscInt i, j, Ii, J, Istart, Iend, m = 8, n = 7, its;
19: PetscBool flg;
20: PetscScalar v;
22: PetscFunctionBeginUser;
23: PetscCall(PetscInitialize(&argc, &args, NULL, help));
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
26: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: Compute the matrix and right-hand-side vector that define
28: the linear system, Ax = b.
29: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
30: /*
31: Create parallel matrix, specifying only its global dimensions.
32: When using MatCreate(), the matrix format can be specified at
33: runtime. Also, the parallel partitioning of the matrix is
34: determined by PETSc at runtime.
36: Performance tuning note: For problems of substantial size,
37: preallocation of matrix memory is crucial for attaining good
38: performance. See the matrix chapter of the users manual for details.
39: */
40: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
41: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
42: PetscCall(MatSetFromOptions(A));
43: PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
44: PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
45: PetscCall(MatSeqSBAIJSetPreallocation(A, 1, 5, NULL));
46: PetscCall(MatMPISBAIJSetPreallocation(A, 1, 5, NULL, 5, NULL));
47: PetscCall(MatMPISELLSetPreallocation(A, 5, NULL, 5, NULL));
48: PetscCall(MatSeqSELLSetPreallocation(A, 5, NULL));
50: /*
51: Currently, all PETSc parallel matrix formats are partitioned by
52: contiguous chunks of rows across the processors. Determine which
53: rows of the matrix are locally owned.
54: */
55: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
57: /*
58: Set matrix elements for the 2-D, five-point stencil in parallel.
59: - Each processor needs to insert only elements that it owns
60: locally (but any non-local elements will be sent to the
61: appropriate processor during matrix assembly).
62: - Always specify global rows and columns of matrix entries.
64: Note: this uses the less common natural ordering that orders first
65: all the unknowns for x = h then for x = 2h etc; Hence you see J = Ii +- n
66: instead of J = I +- m as you might expect. The more standard ordering
67: would first do all variables for y = h, then y = 2h etc.
69: */
70: for (Ii = Istart; Ii < Iend; Ii++) {
71: v = -1.0;
72: i = Ii / n;
73: j = Ii - i * n;
74: if (i > 0) {
75: J = Ii - n;
76: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
77: }
78: if (i < m - 1) {
79: J = Ii + n;
80: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
81: }
82: if (j > 0) {
83: J = Ii - 1;
84: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
85: }
86: if (j < n - 1) {
87: J = Ii + 1;
88: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
89: }
90: v = 4.0;
91: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
92: }
94: /*
95: Assemble matrix, using the 2-step process:
96: MatAssemblyBegin(), MatAssemblyEnd()
97: Computations can be done while messages are in transition
98: by placing code between these two statements.
99: */
100: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
101: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
103: /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
104: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
106: /*
107: Create parallel vectors.
108: - We form 1 vector from scratch and then duplicate as needed.
109: - When using VecCreate(), VecSetSizes and VecSetFromOptions()
110: in this example, we specify only the
111: vector's global dimension; the parallel partitioning is determined
112: at runtime.
113: - When solving a linear system, the vectors and matrices MUST
114: be partitioned accordingly. PETSc automatically generates
115: appropriately partitioned matrices and vectors when MatCreate()
116: and VecCreate() are used with the same communicator.
117: - The user can alternatively specify the local vector and matrix
118: dimensions when more sophisticated partitioning is needed
119: (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
120: below).
121: */
122: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
123: PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n));
124: PetscCall(VecSetFromOptions(u));
125: PetscCall(VecDuplicate(u, &b));
126: PetscCall(VecDuplicate(b, &x));
128: /*
129: Set exact solution; then compute right-hand-side vector.
130: By default we use an exact solution of a vector with all
131: elements of 1.0;
132: */
133: PetscCall(VecSet(u, 1.0));
134: PetscCall(MatMult(A, u, b));
136: /*
137: View the exact solution vector if desired
138: */
139: flg = PETSC_FALSE;
140: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL));
141: if (flg) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Create the linear solver and set various options
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
148: /*
149: Set operators. Here the matrix that defines the linear system
150: also serves as the matrix from which the preconditioner is constructed.
151: */
152: PetscCall(KSPSetOperators(ksp, A, A));
154: /*
155: Set linear solver defaults for this problem (optional).
156: - By extracting the KSP and PC contexts from the KSP context,
157: we can then directly call any KSP and PC routines to set
158: various options.
159: - The following statement is optional; all of these
160: parameters could alternatively be specified at runtime via
161: KSPSetFromOptions(). All of these defaults can be
162: overridden at runtime, as indicated below.
163: */
164: PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), 1.e-50, PETSC_CURRENT, PETSC_CURRENT));
166: /*
167: Set runtime options, e.g.,
168: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
169: These options will override those specified above as long as
170: KSPSetFromOptions() is called _after_ any other customization
171: routines.
172: */
173: PetscCall(KSPSetFromOptions(ksp));
175: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: Solve the linear system
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: PetscCall(KSPSolve(ksp, b, x));
181: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: Check the solution and clean up
183: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184: PetscCall(VecAXPY(x, -1.0, u));
185: PetscCall(VecNorm(x, NORM_2, &norm));
186: PetscCall(KSPGetIterationNumber(ksp, &its));
188: /*
189: Print convergence information. PetscPrintf() produces a single
190: print statement from all processes that share a communicator.
191: An alternative is PetscFPrintf(), which prints to a file.
192: */
193: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));
195: /*
196: Free work space. All PETSc objects should be destroyed when they
197: are no longer needed.
198: */
199: PetscCall(KSPDestroy(&ksp));
200: PetscCall(VecDestroy(&u));
201: PetscCall(VecDestroy(&x));
202: PetscCall(VecDestroy(&b));
203: PetscCall(MatDestroy(&A));
205: /*
206: Always call PetscFinalize() before exiting a program. This routine
207: - finalizes the PETSc libraries as well as MPI
208: - provides summary and diagnostic information if certain runtime
209: options are chosen (e.g., -log_view).
210: */
211: PetscCall(PetscFinalize());
212: return 0;
213: }
215: /*TEST
217: build:
218: requires: !single
220: test:
221: suffix: chebyest_1
222: args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev -ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_monitor_short
224: test:
225: suffix: chebyest_2
226: args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev -ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_esteig_ksp_type cg -ksp_monitor_short
228: test:
229: args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always
231: test:
232: suffix: 2
233: nsize: 2
234: args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always
236: test:
237: suffix: 3
238: args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
240: test:
241: suffix: 4
242: args: -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
244: test:
245: suffix: 5
246: nsize: 2
247: args: -ksp_monitor_short -m 5 -n 5 -mat_view draw -ksp_gmres_cgs_refinement_type refine_always -nox
248: output_file: output/ex2_2.out
250: test:
251: suffix: bjacobi
252: nsize: 4
253: args: -pc_type bjacobi -pc_bjacobi_blocks 1 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres
255: test:
256: suffix: bjacobi_2
257: nsize: 4
258: args: -pc_type bjacobi -pc_bjacobi_blocks 2 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres -ksp_view
260: test:
261: suffix: bjacobi_3
262: nsize: 4
263: args: -pc_type bjacobi -pc_bjacobi_blocks 4 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres
265: test:
266: suffix: qmrcgs
267: args: -ksp_type qmrcgs -pc_type ilu
268: output_file: output/ex2_fbcgs.out
270: test:
271: suffix: qmrcgs_2
272: nsize: 3
273: args: -ksp_type qmrcgs -pc_type bjacobi
274: output_file: output/ex2_fbcgs_2.out
276: test:
277: suffix: fbcgs
278: args: -ksp_type fbcgs -pc_type ilu
280: test:
281: suffix: fbcgs_2
282: nsize: 3
283: args: -ksp_type fbcgsr -pc_type bjacobi
285: test:
286: suffix: groppcg
287: args: -ksp_monitor_short -ksp_type groppcg -m 9 -n 9
289: test:
290: suffix: mkl_pardiso_cholesky
291: requires: mkl_pardiso
292: args: -ksp_type preonly -pc_type cholesky -mat_type sbaij -pc_factor_mat_solver_type mkl_pardiso
294: test:
295: suffix: mkl_pardiso_lu
296: requires: mkl_pardiso
297: args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mkl_pardiso
299: test:
300: suffix: pipebcgs
301: args: -ksp_monitor_short -ksp_type pipebcgs -m 9 -n 9
303: test:
304: suffix: pipecg
305: args: -ksp_monitor_short -ksp_type pipecg -m 9 -n 9
307: test:
308: suffix: pipecgrr
309: args: -ksp_monitor_short -ksp_type pipecgrr -m 9 -n 9
311: test:
312: suffix: pipecr
313: args: -ksp_monitor_short -ksp_type pipecr -m 9 -n 9
315: test:
316: suffix: pipelcg
317: args: -ksp_monitor_short -ksp_type pipelcg -m 9 -n 9 -pc_type none -ksp_pipelcg_pipel 2 -ksp_pipelcg_lmax 2
318: filter: grep -v "sqrt breakdown in iteration"
320: test:
321: suffix: sell
322: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -m 9 -n 9 -mat_type sell
324: test:
325: requires: mumps
326: suffix: sell_mumps
327: args: -ksp_type preonly -m 9 -n 12 -mat_type sell -pc_type lu -pc_factor_mat_solver_type mumps -pc_factor_mat_ordering_type natural
329: test:
330: suffix: telescope
331: nsize: 4
332: args: -m 100 -n 100 -ksp_converged_reason -pc_type telescope -pc_telescope_reduction_factor 4 -telescope_pc_type bjacobi
334: test:
335: suffix: umfpack
336: requires: suitesparse
337: args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type umfpack
339: test:
340: suffix: spqr
341: requires: suitesparse
342: args: -ksp_type preonly -pc_type qr -pc_factor_mat_solver_type spqr
344: test:
345: suffix: pc_symmetric
346: args: -m 10 -n 9 -ksp_converged_reason -ksp_type gmres -ksp_pc_side symmetric -pc_type cholesky
348: test:
349: suffix: pipeprcg
350: args: -ksp_monitor_short -ksp_type pipeprcg -m 9 -n 9
352: test:
353: suffix: pipeprcg_rcw
354: args: -ksp_monitor_short -ksp_type pipeprcg -recompute_w false -m 9 -n 9
356: test:
357: suffix: pipecg2
358: requires: !defined(PETSC_HAVE_THREADSAFETY)
359: args: -ksp_monitor_short -ksp_type pipecg2 -m 9 -n 9 -ksp_norm_type {{preconditioned unpreconditioned natural}}
361: test:
362: suffix: pipecg2_2
363: requires: !defined(PETSC_HAVE_THREADSAFETY)
364: nsize: 4
365: args: -ksp_monitor_short -ksp_type pipecg2 -m 15 -n 9 -ksp_norm_type {{preconditioned unpreconditioned natural}}
367: test:
368: suffix: hpddm
369: nsize: 4
370: requires: hpddm !__float128 !__fp16
371: filter: sed -e "s/ iterations 9/ iterations 8/g"
372: args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{single double}shared output} -ksp_pc_side {{left right}shared output}
374: test:
375: suffix: hpddm___float128
376: output_file: output/ex2_hpddm.out
377: nsize: 4
378: requires: hpddm __float128
379: filter: sed -e "s/ iterations 9/ iterations 8/g"
380: args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{double quadruple}shared output} -ksp_pc_side {{left right}shared output}
382: test:
383: suffix: symmetric_pc
384: nsize: 1
385: args: -ksp_monitor -ksp_type gmres -pc_type bjacobi -sub_pc_type icc -ksp_pc_side symmetric
387: test:
388: suffix: symmetric_pc2
389: nsize: 1
390: args: -ksp_monitor -ksp_type gmres -pc_type bjacobi -sub_pc_type icc -ksp_pc_side symmetric -pc_bjacobi_blocks 2
392: test:
393: suffix: help
394: requires: !openblas !blis !mkl !hpddm !complex !kokkos_kernels !amgx !ml !spai !hypre !viennacl !parms !h2opus !metis !parmetis !superlu_dist !mkl_sparse_optimize !mkl_sparse !mkl_pardiso !mkl_cpardiso !cuda !hip defined(PETSC_USE_LOG) defined(PETSC_USE_INFO) cxx
395: nsize: 1
396: filter: grep -v -e "BLAS options" -e "blas_view" -e "^\-\-\-"
397: args: -ksp_monitor -help -petsc_ci_portable_error_output -error_output_stdout
399: test:
400: suffix: redundant
401: nsize: 2
402: args: -ksp_view -pc_type redundant -redundant_pc_type redundant -redundant_redundant_pc_type redundant -redundant_redundant_redundant_pc_type redundant -redundant_redundant_redundant_redundant_pc_type redundant
404: TEST*/