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