Actual source code: ex6.c
2: static char help[] = "Solves a tridiagonal linear system with KSP. \n\
3: It illustrates how to do one symbolic factorization and multiple numeric factorizations using same matrix structure. \n\n";
5: #include <petscksp.h>
6: int main(int argc, char **args)
7: {
8: Vec x, b, u; /* approx solution, RHS, exact solution */
9: Mat A; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PC pc; /* preconditioner context */
12: PetscReal norm; /* norm of solution error */
13: PetscInt i, col[3], its, rstart, rend, N = 10, num_numfac;
14: PetscScalar value[3];
17: PetscInitialize(&argc, &args, (char *)0, help);
18: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
20: /* Create and assemble matrix. */
21: MatCreate(PETSC_COMM_WORLD, &A);
22: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N);
23: MatSetFromOptions(A);
24: MatSetUp(A);
25: MatGetOwnershipRange(A, &rstart, &rend);
27: value[0] = -1.0;
28: value[1] = 2.0;
29: value[2] = -1.0;
30: for (i = rstart; i < rend; i++) {
31: col[0] = i - 1;
32: col[1] = i;
33: col[2] = i + 1;
34: if (i == 0) {
35: MatSetValues(A, 1, &i, 2, col + 1, value + 1, INSERT_VALUES);
36: } else if (i == N - 1) {
37: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
38: } else {
39: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
40: }
41: }
42: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
44: MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
46: /* Create vectors */
47: MatCreateVecs(A, &x, &b);
48: VecDuplicate(x, &u);
50: /* Set exact solution; then compute right-hand-side vector. */
51: VecSet(u, 1.0);
52: MatMult(A, u, b);
54: /* Create the linear solver and set various options. */
55: KSPCreate(PETSC_COMM_WORLD, &ksp);
56: KSPGetPC(ksp, &pc);
57: PCSetType(pc, PCJACOBI);
58: KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
59: KSPSetOperators(ksp, A, A);
60: KSPSetFromOptions(ksp);
62: num_numfac = 1;
63: PetscOptionsGetInt(NULL, NULL, "-num_numfac", &num_numfac, NULL);
64: while (num_numfac--) {
65: /* An example on how to update matrix A for repeated numerical factorization and solve. */
66: PetscScalar one = 1.0;
67: PetscInt i = 0;
68: MatSetValues(A, 1, &i, 1, &i, &one, ADD_VALUES);
69: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
71: /* Update b */
72: MatMult(A, u, b);
74: /* Solve the linear system */
75: KSPSolve(ksp, b, x);
77: /* Check the solution and clean up */
78: VecAXPY(x, -1.0, u);
79: VecNorm(x, NORM_2, &norm);
80: KSPGetIterationNumber(ksp, &its);
81: if (norm > 100 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
82: }
84: /* Free work space. */
85: VecDestroy(&x);
86: VecDestroy(&u);
87: VecDestroy(&b);
88: MatDestroy(&A);
89: KSPDestroy(&ksp);
91: PetscFinalize();
92: return 0;
93: }
95: /*TEST
97: test:
98: args: -num_numfac 2 -pc_type lu
100: test:
101: suffix: 2
102: args: -num_numfac 2 -pc_type lu -pc_factor_mat_solver_type mumps
103: requires: mumps
105: test:
106: suffix: 3
107: nsize: 3
108: args: -num_numfac 2 -pc_type lu -pc_factor_mat_solver_type mumps
109: requires: mumps
111: TEST*/