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*/
```