Actual source code: ex23.c
2: static char help[] = "Solves a tridiagonal linear system.\n\n";
4: /*
5: Include "petscksp.h" so that we can use KSP solvers. Note that this file
6: automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
12: Note: The corresponding uniprocessor example is ex1.c
13: */
14: #include <petscksp.h>
16: int main(int argc, char **args)
17: {
18: Vec x, b, u; /* approx solution, RHS, exact solution */
19: Mat A; /* linear system matrix */
20: KSP ksp; /* linear solver context */
21: PC pc; /* preconditioner context */
22: PetscReal norm, tol = 1000. * PETSC_MACHINE_EPSILON; /* norm of solution error */
23: PetscInt i, n = 10, col[3], its, rstart, rend, nlocal;
24: PetscScalar one = 1.0, value[3];
26: PetscFunctionBeginUser;
27: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
28: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
30: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: Compute the matrix and right-hand-side vector that define
32: the linear system, Ax = b.
33: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
35: /*
36: Create vectors. Note that we form 1 vector from scratch and
37: then duplicate as needed. For this simple case let PETSc decide how
38: many elements of the vector are stored on each processor. The second
39: argument to VecSetSizes() below causes PETSc to decide.
40: */
41: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
42: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
43: PetscCall(VecSetFromOptions(x));
44: PetscCall(VecDuplicate(x, &b));
45: PetscCall(VecDuplicate(x, &u));
47: /* Identify the starting and ending mesh points on each
48: processor for the interior part of the mesh. We let PETSc decide
49: above. */
51: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
52: PetscCall(VecGetLocalSize(x, &nlocal));
54: /*
55: Create matrix. When using MatCreate(), the matrix format can
56: be specified at runtime.
58: Performance tuning note: For problems of substantial size,
59: preallocation of matrix memory is crucial for attaining good
60: performance. See the matrix chapter of the users manual for details.
62: We pass in nlocal as the "local" size of the matrix to force it
63: to have the same parallel layout as the vector created above.
64: */
65: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
66: PetscCall(MatSetSizes(A, nlocal, nlocal, n, n));
67: PetscCall(MatSetFromOptions(A));
68: PetscCall(MatSetUp(A));
70: /*
71: Assemble matrix.
73: The linear system is distributed across the processors by
74: chunks of contiguous rows, which correspond to contiguous
75: sections of the mesh on which the problem is discretized.
76: For matrix assembly, each processor contributes entries for
77: the part that it owns locally.
78: */
80: if (!rstart) {
81: rstart = 1;
82: i = 0;
83: col[0] = 0;
84: col[1] = 1;
85: value[0] = 2.0;
86: value[1] = -1.0;
87: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
88: }
89: if (rend == n) {
90: rend = n - 1;
91: i = n - 1;
92: col[0] = n - 2;
93: col[1] = n - 1;
94: value[0] = -1.0;
95: value[1] = 2.0;
96: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
97: }
99: /* Set entries corresponding to the mesh interior */
100: value[0] = -1.0;
101: value[1] = 2.0;
102: value[2] = -1.0;
103: for (i = rstart; i < rend; i++) {
104: col[0] = i - 1;
105: col[1] = i;
106: col[2] = i + 1;
107: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
108: }
110: /* Assemble the matrix */
111: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
112: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
114: /*
115: Set exact solution; then compute right-hand-side vector.
116: */
117: PetscCall(VecSet(u, one));
118: PetscCall(MatMult(A, u, b));
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Create the linear solver and set various options
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: /*
124: Create linear solver context
125: */
126: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
128: /*
129: Set operators. Here the matrix that defines the linear system
130: also serves as the preconditioning matrix.
131: */
132: PetscCall(KSPSetOperators(ksp, A, A));
134: /*
135: Set linear solver defaults for this problem (optional).
136: - By extracting the KSP and PC contexts from the KSP context,
137: we can then directly call any KSP and PC routines to set
138: various options.
139: - The following four statements are optional; all of these
140: parameters could alternatively be specified at runtime via
141: KSPSetFromOptions();
142: */
143: PetscCall(KSPGetPC(ksp, &pc));
144: PetscCall(PCSetType(pc, PCJACOBI));
145: PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
147: /*
148: Set runtime options, e.g.,
149: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
150: These options will override those specified above as long as
151: KSPSetFromOptions() is called _after_ any other customization
152: routines.
153: */
154: PetscCall(KSPSetFromOptions(ksp));
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Solve the linear system
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: /*
160: Solve linear system
161: */
162: PetscCall(KSPSolve(ksp, b, x));
164: /*
165: View solver info; we could instead use the option -ksp_view to
166: print this info to the screen at the conclusion of KSPSolve().
167: */
168: PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Check solution and clean up
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: /*
174: Check the error
175: */
176: PetscCall(VecAXPY(x, -1.0, u));
177: PetscCall(VecNorm(x, NORM_2, &norm));
178: PetscCall(KSPGetIterationNumber(ksp, &its));
179: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
181: /*
182: Free work space. All PETSc objects should be destroyed when they
183: are no longer needed.
184: */
185: PetscCall(VecDestroy(&x));
186: PetscCall(VecDestroy(&u));
187: PetscCall(VecDestroy(&b));
188: PetscCall(MatDestroy(&A));
189: PetscCall(KSPDestroy(&ksp));
191: /*
192: Always call PetscFinalize() before exiting a program. This routine
193: - finalizes the PETSc libraries as well as MPI
194: - provides summary and diagnostic information if certain runtime
195: options are chosen (e.g., -log_view).
196: */
197: PetscCall(PetscFinalize());
198: return 0;
199: }
201: /*TEST
203: build:
204: requires: !complex !single
206: test:
207: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
209: test:
210: suffix: 2
211: nsize: 3
212: args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
214: test:
215: suffix: 3
216: nsize: 2
217: args: -ksp_monitor_short -ksp_rtol 1e-6 -ksp_type pipefgmres
219: TEST*/