Actual source code: ex58.c
1: static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
3: /*
4: Modified from ex1.c for testing matrix operations when matrix structure is changed.
5: Contributed by Jose E. Roman, Feb. 2012.
6: */
7: #include <petscksp.h>
9: int main(int argc, char **args)
10: {
11: Vec x, b, u; /* approx solution, RHS, exact solution */
12: Mat A, B, C; /* linear system matrix */
13: KSP ksp; /* linear solver context */
14: PC pc; /* preconditioner context */
15: PetscReal norm; /* norm of solution error */
16: PetscInt i, n = 20, col[3], its;
17: PetscMPIInt size;
18: PetscScalar one = 1.0, value[3];
19: PetscBool nonzeroguess = PETSC_FALSE;
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
23: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
26: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonzero_guess", &nonzeroguess, NULL));
28: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: Compute the matrix and right-hand-side vector that define
30: the linear system, Ax = b.
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: /*
34: Create vectors. Note that we form 1 vector from scratch and
35: then duplicate as needed.
36: */
37: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
38: PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
39: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
40: PetscCall(VecSetFromOptions(x));
41: PetscCall(VecDuplicate(x, &b));
42: PetscCall(VecDuplicate(x, &u));
44: /*
45: Create matrix. When using MatCreate(), the matrix format can
46: be specified at runtime.
48: Performance tuning note: For problems of substantial size,
49: preallocation of matrix memory is crucial for attaining good
50: performance. See the matrix chapter of the users manual for details.
51: */
52: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
53: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
54: PetscCall(MatSetFromOptions(A));
55: PetscCall(MatSetUp(A));
57: /*
58: Assemble matrix
59: */
60: value[0] = -1.0;
61: value[1] = 2.0;
62: value[2] = -1.0;
63: for (i = 1; i < n - 1; i++) {
64: col[0] = i - 1;
65: col[1] = i;
66: col[2] = i + 1;
67: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
68: }
69: i = n - 1;
70: col[0] = n - 2;
71: col[1] = n - 1;
72: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
73: i = 0;
74: col[0] = 0;
75: col[1] = 1;
76: value[0] = 2.0;
77: value[1] = -1.0;
78: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
79: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
80: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
82: /*
83: jroman: added matrices
84: */
85: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
86: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n));
87: PetscCall(MatSetFromOptions(B));
88: PetscCall(MatSetUp(B));
89: for (i = 0; i < n; i++) {
90: PetscCall(MatSetValue(B, i, i, value[1], INSERT_VALUES));
91: if (n - i + n / 3 < n) {
92: PetscCall(MatSetValue(B, n - i + n / 3, i, value[0], INSERT_VALUES));
93: PetscCall(MatSetValue(B, i, n - i + n / 3, value[0], INSERT_VALUES));
94: }
95: }
96: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
97: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
98: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &C));
99: PetscCall(MatAXPY(C, 2.0, B, DIFFERENT_NONZERO_PATTERN));
101: /*
102: Set exact solution; then compute right-hand-side vector.
103: */
104: PetscCall(VecSet(u, one));
105: PetscCall(MatMult(C, u, b));
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create the linear solver and set various options
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Create linear solver context
112: */
113: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
115: /*
116: Set operators. Here the matrix that defines the linear system
117: also serves as the preconditioning matrix.
118: */
119: PetscCall(KSPSetOperators(ksp, C, C));
121: /*
122: Set linear solver defaults for this problem (optional).
123: - By extracting the KSP and PC contexts from the KSP context,
124: we can then directly call any KSP and PC routines to set
125: various options.
126: - The following four statements are optional; all of these
127: parameters could alternatively be specified at runtime via
128: KSPSetFromOptions();
129: */
130: PetscCall(KSPGetPC(ksp, &pc));
131: PetscCall(PCSetType(pc, PCJACOBI));
132: PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
134: /*
135: Set runtime options, e.g.,
136: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
137: These options will override those specified above as long as
138: KSPSetFromOptions() is called _after_ any other customization
139: routines.
140: */
141: PetscCall(KSPSetFromOptions(ksp));
143: if (nonzeroguess) {
144: PetscScalar p = .5;
145: PetscCall(VecSet(x, p));
146: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
147: }
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Solve the linear system
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: /*
153: Solve linear system
154: */
155: PetscCall(KSPSolve(ksp, b, x));
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Check solution and clean up
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: /*
161: Check the error
162: */
163: PetscCall(VecAXPY(x, -1.0, u));
164: PetscCall(VecNorm(x, NORM_2, &norm));
165: PetscCall(KSPGetIterationNumber(ksp, &its));
166: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
168: /*
169: Free work space. All PETSc objects should be destroyed when they
170: are no longer needed.
171: */
172: PetscCall(VecDestroy(&x));
173: PetscCall(VecDestroy(&u));
174: PetscCall(VecDestroy(&b));
175: PetscCall(MatDestroy(&A));
176: PetscCall(MatDestroy(&B));
177: PetscCall(MatDestroy(&C));
178: PetscCall(KSPDestroy(&ksp));
180: /*
181: Always call PetscFinalize() before exiting a program. This routine
182: - finalizes the PETSc libraries as well as MPI
183: - provides summary and diagnostic information if certain runtime
184: options are chosen (e.g., -log_view).
185: */
186: PetscCall(PetscFinalize());
187: return 0;
188: }
190: /*TEST
192: test:
193: args: -mat_type aij
194: output_file: output/ex58.out
196: test:
197: suffix: baij
198: args: -mat_type baij
199: output_file: output/ex58.out
201: test:
202: suffix: sbaij
203: args: -mat_type sbaij
204: output_file: output/ex58.out
206: TEST*/