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