Actual source code: ex48.c
1: static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
3: /*
4: Test example that demonstrates how MINRES can produce a dp of zero
5: but still converge.
7: Provided by: Mark Filipiak <mjf@staffmail.ed.ac.uk>
8: */
9: #include <petscksp.h>
11: int main(int argc, char **args)
12: {
13: Vec x, b, u; /* approx solution, RHS, exact solution */
14: Mat A; /* linear system matrix */
15: KSP ksp; /* linear solver context */
16: PC pc; /* preconditioner context */
17: PetscReal norm;
18: PetscInt i, n = 2, col[3], its;
19: PetscMPIInt size;
20: PetscScalar one = 1.0, value[3], shift = 0.0;
21: PetscBool nonzeroguess = PETSC_FALSE;
23: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, NULL, help));
25: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
27: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
28: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, NULL));
29: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-shift", &shift, 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));
84: PetscCall(MatShift(A, -shift));
86: /*
87: Set constant right-hand-side vector.
88: */
89: PetscCall(VecSet(b, one));
90: /*
91: Solution = RHS for the matrix [[2 -1] [-1 2]] and RHS [1 1]
92: */
93: PetscCall(VecSet(u, one));
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create the linear solver and set various options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: /*
99: Create linear solver context
100: */
101: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
103: /*
104: Set operators. Here the matrix that defines the linear system
105: also serves as the preconditioning matrix.
106: */
107: PetscCall(KSPSetOperators(ksp, A, A));
109: /*
110: Set linear solver defaults for this problem (optional).
111: - By extracting the KSP and PC contexts from the KSP context,
112: we can then directly call any KSP and PC routines to set
113: various options.
114: - The following four statements are optional; all of these
115: parameters could alternatively be specified at runtime via
116: KSPSetFromOptions();
117: */
118: PetscCall(KSPGetPC(ksp, &pc));
119: PetscCall(PCSetType(pc, PCJACOBI));
120: PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
122: /*
123: Set runtime options, e.g.,
124: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
125: These options will override those specified above as long as
126: KSPSetFromOptions() is called _after_ any other customization
127: routines.
128: */
129: PetscCall(KSPSetFromOptions(ksp));
131: if (nonzeroguess) {
132: PetscScalar p = .5;
133: PetscCall(VecSet(x, p));
134: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
135: }
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Solve the linear system
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: /*
141: Solve linear system
142: */
143: PetscCall(KSPSolve(ksp, b, x));
145: /*
146: View solver info; we could instead use the option -ksp_view to
147: print this info to the screen at the conclusion of KSPSolve().
148: */
149: PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Check solution and clean up
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: /*
155: Check the error
156: */
157: PetscCall(VecAXPY(x, -1.0, u));
158: PetscCall(VecNorm(x, NORM_2, &norm));
159: PetscCall(KSPGetIterationNumber(ksp, &its));
160: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
162: /*
163: Free work space. All PETSc objects should be destroyed when they
164: are no longer needed.
165: */
166: PetscCall(VecDestroy(&x));
167: PetscCall(VecDestroy(&u));
168: PetscCall(VecDestroy(&b));
169: PetscCall(MatDestroy(&A));
170: PetscCall(KSPDestroy(&ksp));
172: /*
173: Always call PetscFinalize() before exiting a program. This routine
174: - finalizes the PETSc libraries as well as MPI
175: - provides summary and diagnostic information if certain runtime
176: options are chosen (e.g., -log_view).
177: */
178: PetscCall(PetscFinalize());
179: return 0;
180: }
182: /*TEST
184: test:
185: args: -ksp_monitor_short -ksp_converged_reason -ksp_type minres -pc_type jacobi -ksp_error_if_not_converged
187: test:
188: suffix: 2
189: args: -ksp_monitor_short -ksp_converged_reason -ksp_type minres -pc_type jacobi -ksp_view_eigenvalues -ksp_view_singularvalues -shift 1.5 -ksp_rtol 0 -ksp_atol 0 -n 5 -ksp_max_it 5
191: TEST*/