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];
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
20: /* Create and assemble matrix. */
21: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
22: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N));
23: PetscCall(MatSetFromOptions(A));
24: PetscCall(MatSetUp(A));
25: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 2, col + 1, value + 1, INSERT_VALUES));
36: } else if (i == N - 1) {
37: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
38: } else {
39: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
40: }
41: }
42: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
43: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
44: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
46: /* Create vectors */
47: PetscCall(MatCreateVecs(A, &x, &b));
48: PetscCall(VecDuplicate(x, &u));
50: /* Set exact solution; then compute right-hand-side vector. */
51: PetscCall(VecSet(u, 1.0));
52: PetscCall(MatMult(A, u, b));
54: /* Create the linear solver and set various options. */
55: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
56: PetscCall(KSPGetPC(ksp, &pc));
57: PetscCall(PCSetType(pc, PCJACOBI));
58: PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
59: PetscCall(KSPSetOperators(ksp, A, A));
60: PetscCall(KSPSetFromOptions(ksp));
62: num_numfac = 1;
63: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 1, &i, &one, ADD_VALUES));
69: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
71: /* Update b */
72: PetscCall(MatMult(A, u, b));
74: /* Solve the linear system */
75: PetscCall(KSPSolve(ksp, b, x));
77: /* Check the solution and clean up */
78: PetscCall(VecAXPY(x, -1.0, u));
79: PetscCall(VecNorm(x, NORM_2, &norm));
80: PetscCall(KSPGetIterationNumber(ksp, &its));
81: if (norm > 100 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
82: }
84: /* Free work space. */
85: PetscCall(VecDestroy(&x));
86: PetscCall(VecDestroy(&u));
87: PetscCall(VecDestroy(&b));
88: PetscCall(MatDestroy(&A));
89: PetscCall(KSPDestroy(&ksp));
91: PetscCall(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*/