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;
23: PetscInitialize(&argc, &args, (char *)0, help);
24: MPI_Comm_size(PETSC_COMM_WORLD, &size);
26: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
27: 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: VecCreate(PETSC_COMM_WORLD, &x);
39: PetscObjectSetName((PetscObject)x, "Solution");
40: VecSetSizes(x, PETSC_DECIDE, n);
41: VecSetFromOptions(x);
42: VecDuplicate(x, &b);
43: 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: MatCreate(PETSC_COMM_WORLD, &A);
54: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
55: MatSetFromOptions(A);
56: 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: 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: 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: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
80: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
83: /*
84: jroman: added matrices
85: */
86: MatCreate(PETSC_COMM_WORLD, &B);
87: MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n);
88: MatSetFromOptions(B);
89: MatSetUp(B);
90: for (i = 0; i < n; i++) {
91: MatSetValue(B, i, i, value[1], INSERT_VALUES);
92: if (n - i + n / 3 < n) {
93: MatSetValue(B, n - i + n / 3, i, value[0], INSERT_VALUES);
94: MatSetValue(B, i, n - i + n / 3, value[0], INSERT_VALUES);
95: }
96: }
97: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
98: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
99: MatDuplicate(A, MAT_COPY_VALUES, &C);
100: MatAXPY(C, 2.0, B, DIFFERENT_NONZERO_PATTERN);
102: /*
103: Set exact solution; then compute right-hand-side vector.
104: */
105: VecSet(u, one);
106: MatMult(C, u, b);
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Create the linear solver and set various options
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: /*
112: Create linear solver context
113: */
114: 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: 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: KSPGetPC(ksp, &pc);
132: PCSetType(pc, PCJACOBI);
133: 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: KSPSetFromOptions(ksp);
144: if (nonzeroguess) {
145: PetscScalar p = .5;
146: VecSet(x, p);
147: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
148: }
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Solve the linear system
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: /*
154: Solve linear system
155: */
156: KSPSolve(ksp, b, x);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Check solution and clean up
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: /*
162: Check the error
163: */
164: VecAXPY(x, -1.0, u);
165: VecNorm(x, NORM_2, &norm);
166: KSPGetIterationNumber(ksp, &its);
167: 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: VecDestroy(&x);
174: VecDestroy(&u);
175: VecDestroy(&b);
176: MatDestroy(&A);
177: MatDestroy(&B);
178: MatDestroy(&C);
179: 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: 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*/