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