Actual source code: ex48.c
2: static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
4: /*
5: Test example that demonstrates how MINRES can produce a dp of zero
6: but still converge.
8: Provided by: Mark Filipiak <mjf@staffmail.ed.ac.uk>
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: PC pc; /* preconditioner context */
18: PetscReal norm;
19: PetscInt i, n = 2, col[3], its;
20: PetscMPIInt size;
21: PetscScalar one = 1.0, value[3];
22: PetscBool nonzeroguess = PETSC_FALSE;
24: PetscFunctionBeginUser;
25: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
26: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
27: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
28: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, NULL));
31: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: Compute the matrix and right-hand-side vector that define
33: the linear system, Ax = b.
34: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
36: /*
37: Create vectors. Note that we form 1 vector from scratch and
38: then duplicate as needed.
39: */
40: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
41: PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
42: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
43: PetscCall(VecSetFromOptions(x));
44: PetscCall(VecDuplicate(x, &b));
45: PetscCall(VecDuplicate(x, &u));
47: /*
48: Create matrix. When using MatCreate(), the matrix format can
49: be specified at runtime.
51: Performance tuning note: For problems of substantial size,
52: preallocation of matrix memory is crucial for attaining good
53: performance. See the matrix chapter of the users manual for details.
54: */
55: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
56: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
57: PetscCall(MatSetFromOptions(A));
58: PetscCall(MatSetUp(A));
60: /*
61: Assemble matrix
62: */
63: value[0] = -1.0;
64: value[1] = 2.0;
65: value[2] = -1.0;
66: for (i = 1; i < n - 1; i++) {
67: col[0] = i - 1;
68: col[1] = i;
69: col[2] = i + 1;
70: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
71: }
72: i = n - 1;
73: col[0] = n - 2;
74: col[1] = n - 1;
75: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
76: i = 0;
77: col[0] = 0;
78: col[1] = 1;
79: value[0] = 2.0;
80: value[1] = -1.0;
81: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
82: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
85: /*
86: Set constant right-hand-side vector.
87: */
88: PetscCall(VecSet(b, one));
89: /*
90: Solution = RHS for the matrix [[2 -1] [-1 2]] and RHS [1 1]
91: */
92: PetscCall(VecSet(u, one));
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create the linear solver and set various options
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: /*
98: Create linear solver context
99: */
100: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
102: /*
103: Set operators. Here the matrix that defines the linear system
104: also serves as the preconditioning matrix.
105: */
106: PetscCall(KSPSetOperators(ksp, A, A));
108: /*
109: Set linear solver defaults for this problem (optional).
110: - By extracting the KSP and PC contexts from the KSP context,
111: we can then directly call any KSP and PC routines to set
112: various options.
113: - The following four statements are optional; all of these
114: parameters could alternatively be specified at runtime via
115: KSPSetFromOptions();
116: */
117: PetscCall(KSPGetPC(ksp, &pc));
118: PetscCall(PCSetType(pc, PCJACOBI));
119: PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
121: /*
122: Set runtime options, e.g.,
123: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
124: These options will override those specified above as long as
125: KSPSetFromOptions() is called _after_ any other customization
126: routines.
127: */
128: PetscCall(KSPSetFromOptions(ksp));
130: if (nonzeroguess) {
131: PetscScalar p = .5;
132: PetscCall(VecSet(x, p));
133: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
134: }
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Solve the linear system
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: /*
140: Solve linear system
141: */
142: PetscCall(KSPSolve(ksp, b, x));
144: /*
145: View solver info; we could instead use the option -ksp_view to
146: print this info to the screen at the conclusion of KSPSolve().
147: */
148: PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Check solution and clean up
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: /*
154: Check the error
155: */
156: PetscCall(VecAXPY(x, -1.0, u));
157: PetscCall(VecNorm(x, NORM_2, &norm));
158: PetscCall(KSPGetIterationNumber(ksp, &its));
159: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
161: /*
162: Free work space. All PETSc objects should be destroyed when they
163: are no longer needed.
164: */
165: PetscCall(VecDestroy(&x));
166: PetscCall(VecDestroy(&u));
167: PetscCall(VecDestroy(&b));
168: PetscCall(MatDestroy(&A));
169: PetscCall(KSPDestroy(&ksp));
171: /*
172: Always call PetscFinalize() before exiting a program. This routine
173: - finalizes the PETSc libraries as well as MPI
174: - provides summary and diagnostic information if certain runtime
175: options are chosen (e.g., -log_view).
176: */
177: PetscCall(PetscFinalize());
178: return 0;
179: }
181: /*TEST
183: test:
184: args: -ksp_monitor_short -ksp_converged_reason -ksp_type minres -pc_type jacobi -ksp_error_if_not_converged
186: TEST*/