1 static char help[] = "Solves a tridiagonal linear system. Designed to compare SOR for different Mat impls.\n\n";
2
3 #include <petscksp.h>
4
main(int argc,char ** args)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;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-k", &k, NULL));
19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
20
21 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
22 PetscCall(KSPSetFromOptions(ksp));
23 PetscCall(KSPGetOperators(ksp, &A, NULL));
24
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));
39
40 PetscCall(VecSetFromOptions(b));
41 PetscCall(VecSet(b, 1.0));
42 PetscCall(VecSetFromOptions(x));
43 PetscCall(VecSet(x, 2.0));
44
45 PetscCall(KSPSolve(ksp, b, x));
46
47 PetscCall(VecDestroy(&b));
48 PetscCall(VecDestroy(&x));
49 PetscCall(KSPDestroy(&ksp));
50
51 PetscCall(PetscFinalize());
52 return 0;
53 }
54