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