diff --git a/main.c b/main.c index a3686cf..f2a2b1c 100644 --- a/main.c +++ b/main.c @@ -1,99 +1,94 @@ -static char help[] = "Solves a tridiagonal linear system with KSP.\n\n"; +static char help[] = "Solves a tridiagonal linear system.\n\n"; #include -#include int main(int argc, char **args) { - Vec x, b, u; /* approx solution, RHS, exact solution */ - Mat A; /* linear system matrix */ - KSP ksp; /* linear solver context */ - PC pc; /* preconditioner context */ - PetscReal norm; /* norm of solution error */ - PetscInt i, n = 10, col[3], its; - PetscMPIInt size; - PetscScalar value[3]; + Vec x, b, u; /* approx solution, RHS, exact solution */ + Mat A; /* linear system matrix */ + KSP ksp; /* linear solver context */ + PC pc; /* preconditioner context */ + PetscReal norm, tol = 1000. * PETSC_MACHINE_EPSILON; /* norm of solution error */ + PetscInt i, n = 10, col[3], its, rstart, rend, nlocal; + PetscScalar one = 1.0, value[3]; PetscFunctionBeginUser; PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); - PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); - PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); - PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); - PetscCall(VecCreate(PETSC_COMM_SELF, &x)); - PetscCall(PetscObjectSetName((PetscObject)x, "Solution")); + PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); PetscCall(VecSetSizes(x, PETSC_DECIDE, n)); PetscCall(VecSetFromOptions(x)); PetscCall(VecDuplicate(x, &b)); PetscCall(VecDuplicate(x, &u)); - PetscCall(MatCreate(PETSC_COMM_SELF, &A)); - PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)); + PetscCall(VecGetOwnershipRange(x, &rstart, &rend)); + PetscCall(VecGetLocalSize(x, &nlocal)); + + PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); + PetscCall(MatSetSizes(A, nlocal, nlocal, n, n)); PetscCall(MatSetFromOptions(A)); PetscCall(MatSetUp(A)); + if (!rstart) { + rstart = 1; + i = 0; + col[0] = 0; + col[1] = 1; + value[0] = 2.0; + value[1] = -1.0; + PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); + } + if (rend == n) { + rend = n - 1; + i = n - 1; + col[0] = n - 2; + col[1] = n - 1; + value[0] = -1.0; + value[1] = 2.0; + PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); + } + value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; - for (i = 1; i < n - 1; i++) { + for (i = rstart; i < rend; i++) { col[0] = i - 1; col[1] = i; col[2] = i + 1; PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); } - i = n - 1; - col[0] = n - 2; - col[1] = n - 1; - PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); - i = 0; - col[0] = 0; - col[1] = 1; - value[0] = 2.0; - value[1] = -1.0; - PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); + PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); - PetscCall(VecSet(u, 1.0)); + PetscCall(VecSet(u, one)); PetscCall(MatMult(A, u, b)); - PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp)); + PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); PetscCall(KSPSetOperators(ksp, A, A)); PetscCall(KSPGetPC(ksp, &pc)); PetscCall(PCSetType(pc, PCJACOBI)); - PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); + PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); PetscCall(KSPSetFromOptions(ksp)); PetscCall(KSPSolve(ksp, b, x)); + PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD)); + PetscCall(VecAXPY(x, -1.0, u)); PetscCall(VecNorm(x, NORM_2, &norm)); PetscCall(KSPGetIterationNumber(ksp, &its)); - PetscCall(PetscPrintf(PETSC_COMM_SELF, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its)); - - PetscCall(MatShift(A, 2.0)); - PetscCall(KSPSolve(ksp, b, x)); - - PetscCall(KSPDestroy(&ksp)); - - if (PCMPIServerActive) { - PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp)); - PetscCall(KSPSetOptionsPrefix(ksp, "prefix_test_")); - PetscCall(MatSetOptionsPrefix(A, "prefix_test_")); - PetscCall(KSPSetOperators(ksp, A, A)); - PetscCall(KSPSetFromOptions(ksp)); - PetscCall(KSPSolve(ksp, b, x)); - PetscCall(KSPDestroy(&ksp)); - } + if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its)); PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&u)); PetscCall(VecDestroy(&b)); PetscCall(MatDestroy(&A)); + PetscCall(KSPDestroy(&ksp)); PetscCall(PetscFinalize()); return 0;