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 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