1 static char help[] = "Tests PC and KSP on a tridiagonal matrix. Note that most\n\ 2 users should employ the KSP interface instead of using PC directly.\n\n"; 3 4 #include <petscksp.h> 5 6 int main(int argc, char **args) 7 { 8 Mat mat; /* matrix */ 9 Vec b, ustar, u; /* vectors (RHS, exact solution, approx solution) */ 10 PC pc; /* PC context */ 11 KSP ksp; /* KSP context */ 12 PetscInt n = 10, i, its, col[3]; 13 PetscScalar value[3]; 14 PCType pcname; 15 KSPType kspname; 16 PetscReal norm, tol = 1000. * PETSC_MACHINE_EPSILON; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 20 /* Create and initialize vectors */ 21 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &b)); 22 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &ustar)); 23 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &u)); 24 PetscCall(VecSet(ustar, 1.0)); 25 PetscCall(VecSet(u, 0.0)); 26 27 /* Create and assemble matrix */ 28 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &mat)); 29 value[0] = -1.0; 30 value[1] = 2.0; 31 value[2] = -1.0; 32 for (i = 1; i < n - 1; i++) { 33 col[0] = i - 1; 34 col[1] = i; 35 col[2] = i + 1; 36 PetscCall(MatSetValues(mat, 1, &i, 3, col, value, INSERT_VALUES)); 37 } 38 i = n - 1; 39 col[0] = n - 2; 40 col[1] = n - 1; 41 PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES)); 42 i = 0; 43 col[0] = 0; 44 col[1] = 1; 45 value[0] = 2.0; 46 value[1] = -1.0; 47 PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES)); 48 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 49 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 50 51 /* Compute right-hand-side vector */ 52 PetscCall(MatMult(mat, ustar, b)); 53 54 /* Create PC context and set up data structures */ 55 PetscCall(PCCreate(PETSC_COMM_WORLD, &pc)); 56 PetscCall(PCSetType(pc, PCNONE)); 57 PetscCall(PCSetFromOptions(pc)); 58 PetscCall(PCSetOperators(pc, mat, mat)); 59 PetscCall(PCSetUp(pc)); 60 61 /* Create KSP context and set up data structures */ 62 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 63 PetscCall(KSPSetType(ksp, KSPRICHARDSON)); 64 PetscCall(KSPSetFromOptions(ksp)); 65 PetscCall(PCSetOperators(pc, mat, mat)); 66 PetscCall(KSPSetPC(ksp, pc)); 67 PetscCall(KSPSetUp(ksp)); 68 69 /* Solve the problem */ 70 PetscCall(KSPGetType(ksp, &kspname)); 71 PetscCall(PCGetType(pc, &pcname)); 72 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Running %s with %s preconditioning\n", kspname, pcname)); 73 PetscCall(KSPSolve(ksp, b, u)); 74 PetscCall(VecAXPY(u, -1.0, ustar)); 75 PetscCall(VecNorm(u, NORM_2, &norm)); 76 PetscCall(KSPGetIterationNumber(ksp, &its)); 77 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "2 norm of error %g Number of iterations %" PetscInt_FMT "\n", (double)norm, its)); 78 79 /* Free data structures */ 80 PetscCall(KSPDestroy(&ksp)); 81 PetscCall(VecDestroy(&u)); 82 PetscCall(VecDestroy(&ustar)); 83 PetscCall(VecDestroy(&b)); 84 PetscCall(MatDestroy(&mat)); 85 PetscCall(PCDestroy(&pc)); 86 87 PetscCall(PetscFinalize()); 88 return 0; 89 } 90 91 /*TEST 92 93 test: 94 args: -ksp_type cg -ksp_monitor_short 95 96 TEST*/ 97