Actual source code: ex44.c
1: static char help[] = "Solves a tridiagonal linear system. Designed to compare SOR for different Mat impls.\n\n";
3: #include <petscksp.h>
5: int main(int argc, char **args)
6: {
7: KSP ksp; /* linear solver context */
8: Mat A; /* linear system matrix */
9: Vec x, b; /* approx solution, RHS */
10: PetscInt Ii, Istart, Iend;
11: PetscScalar v[3] = {-1. / 2., 1., -1. / 2.};
12: PetscInt j[3];
13: PetscInt k = 15;
14: PetscInt M, m = 420;
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-k", &k, NULL));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
21: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
22: PetscCall(KSPSetFromOptions(ksp));
23: PetscCall(KSPGetOperators(ksp, &A, NULL));
25: PetscCall(MatSetSizes(A, m, m, PETSC_DETERMINE, PETSC_DETERMINE));
26: PetscCall(MatSetFromOptions(A));
27: PetscCall(MatSetUp(A));
28: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
29: PetscCall(MatGetSize(A, &M, NULL));
30: for (Ii = Istart; Ii < Iend; Ii++) {
31: j[0] = Ii - k;
32: j[1] = Ii;
33: j[2] = (Ii + k) < M ? (Ii + k) : -1;
34: PetscCall(MatSetValues(A, 1, &Ii, 3, j, v, INSERT_VALUES));
35: }
36: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
37: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
38: PetscCall(MatCreateVecs(A, &x, &b));
40: PetscCall(VecSetFromOptions(b));
41: PetscCall(VecSet(b, 1.0));
42: PetscCall(VecSetFromOptions(x));
43: PetscCall(VecSet(x, 2.0));
45: PetscCall(KSPSolve(ksp, b, x));
47: PetscCall(VecDestroy(&b));
48: PetscCall(VecDestroy(&x));
49: PetscCall(KSPDestroy(&ksp));
51: PetscCall(PetscFinalize());
52: return 0;
53: }