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