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;

 18:   PetscInitialize(&argc, &args, (char *)0, help);
 19:   PetscOptionsGetInt(NULL, NULL, "-k", &k, NULL);
 20:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);

 22:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 23:   KSPSetFromOptions(ksp);
 24:   KSPGetOperators(ksp, &A, NULL);

 26:   MatSetSizes(A, m, m, PETSC_DETERMINE, PETSC_DETERMINE);
 27:   MatSetFromOptions(A);
 28:   MatSetUp(A);
 29:   MatGetOwnershipRange(A, &Istart, &Iend);
 30:   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:     MatSetValues(A, 1, &Ii, 3, j, v, INSERT_VALUES);
 36:   }
 37:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 38:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 39:   MatCreateVecs(A, &x, &b);

 41:   VecSetFromOptions(b);
 42:   VecSet(b, 1.0);
 43:   VecSetFromOptions(x);
 44:   VecSet(x, 2.0);

 46:   KSPSolve(ksp, b, x);

 48:   VecDestroy(&b);
 49:   VecDestroy(&x);
 50:   KSPDestroy(&ksp);

 52:   PetscFinalize();
 53:   return 0;
 54: }