1c4762a1bSJed Brown static char help[] = "Tests PC and KSP on a tridiagonal matrix. Note that most\n\
2c4762a1bSJed Brown users should employ the KSP interface instead of using PC directly.\n\n";
3c4762a1bSJed Brown
4c4762a1bSJed Brown #include <petscksp.h>
5c4762a1bSJed Brown
main(int argc,char ** args)6d71ae5a4SJacob Faibussowitsch int main(int argc, char **args)
7d71ae5a4SJacob Faibussowitsch {
8c4762a1bSJed Brown Mat mat; /* matrix */
9c4762a1bSJed Brown Vec b, ustar, u; /* vectors (RHS, exact solution, approx solution) */
10c4762a1bSJed Brown PC pc; /* PC context */
11c4762a1bSJed Brown KSP ksp; /* KSP context */
12c4762a1bSJed Brown PetscInt n = 10, i, its, col[3];
13c4762a1bSJed Brown PetscScalar value[3];
14c4762a1bSJed Brown PCType pcname;
15c4762a1bSJed Brown KSPType kspname;
16c4762a1bSJed Brown PetscReal norm, tol = 1000. * PETSC_MACHINE_EPSILON;
17c4762a1bSJed Brown
18327415f7SBarry Smith PetscFunctionBeginUser;
19*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &args, NULL, help));
20c4762a1bSJed Brown /* Create and initialize vectors */
219566063dSJacob Faibussowitsch PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &b));
229566063dSJacob Faibussowitsch PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &ustar));
239566063dSJacob Faibussowitsch PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &u));
249566063dSJacob Faibussowitsch PetscCall(VecSet(ustar, 1.0));
259566063dSJacob Faibussowitsch PetscCall(VecSet(u, 0.0));
26c4762a1bSJed Brown
27c4762a1bSJed Brown /* Create and assemble matrix */
289566063dSJacob Faibussowitsch PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &mat));
299371c9d4SSatish Balay value[0] = -1.0;
309371c9d4SSatish Balay value[1] = 2.0;
319371c9d4SSatish Balay value[2] = -1.0;
32c4762a1bSJed Brown for (i = 1; i < n - 1; i++) {
339371c9d4SSatish Balay col[0] = i - 1;
349371c9d4SSatish Balay col[1] = i;
359371c9d4SSatish Balay col[2] = i + 1;
369566063dSJacob Faibussowitsch PetscCall(MatSetValues(mat, 1, &i, 3, col, value, INSERT_VALUES));
37c4762a1bSJed Brown }
389371c9d4SSatish Balay i = n - 1;
399371c9d4SSatish Balay col[0] = n - 2;
409371c9d4SSatish Balay col[1] = n - 1;
419566063dSJacob Faibussowitsch PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES));
429371c9d4SSatish Balay i = 0;
439371c9d4SSatish Balay col[0] = 0;
449371c9d4SSatish Balay col[1] = 1;
459371c9d4SSatish Balay value[0] = 2.0;
469371c9d4SSatish Balay value[1] = -1.0;
479566063dSJacob Faibussowitsch PetscCall(MatSetValues(mat, 1, &i, 2, col, value, INSERT_VALUES));
489566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
499566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
50c4762a1bSJed Brown
51c4762a1bSJed Brown /* Compute right-hand-side vector */
529566063dSJacob Faibussowitsch PetscCall(MatMult(mat, ustar, b));
53c4762a1bSJed Brown
54c4762a1bSJed Brown /* Create PC context and set up data structures */
559566063dSJacob Faibussowitsch PetscCall(PCCreate(PETSC_COMM_WORLD, &pc));
569566063dSJacob Faibussowitsch PetscCall(PCSetType(pc, PCNONE));
579566063dSJacob Faibussowitsch PetscCall(PCSetFromOptions(pc));
589566063dSJacob Faibussowitsch PetscCall(PCSetOperators(pc, mat, mat));
599566063dSJacob Faibussowitsch PetscCall(PCSetUp(pc));
60c4762a1bSJed Brown
61c4762a1bSJed Brown /* Create KSP context and set up data structures */
629566063dSJacob Faibussowitsch PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
639566063dSJacob Faibussowitsch PetscCall(KSPSetType(ksp, KSPRICHARDSON));
649566063dSJacob Faibussowitsch PetscCall(KSPSetFromOptions(ksp));
659566063dSJacob Faibussowitsch PetscCall(PCSetOperators(pc, mat, mat));
669566063dSJacob Faibussowitsch PetscCall(KSPSetPC(ksp, pc));
679566063dSJacob Faibussowitsch PetscCall(KSPSetUp(ksp));
68c4762a1bSJed Brown
69c4762a1bSJed Brown /* Solve the problem */
709566063dSJacob Faibussowitsch PetscCall(KSPGetType(ksp, &kspname));
719566063dSJacob Faibussowitsch PetscCall(PCGetType(pc, &pcname));
729566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "Running %s with %s preconditioning\n", kspname, pcname));
739566063dSJacob Faibussowitsch PetscCall(KSPSolve(ksp, b, u));
749566063dSJacob Faibussowitsch PetscCall(VecAXPY(u, -1.0, ustar));
759566063dSJacob Faibussowitsch PetscCall(VecNorm(u, NORM_2, &norm));
769566063dSJacob Faibussowitsch PetscCall(KSPGetIterationNumber(ksp, &its));
7748a46eb9SPierre Jolivet if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "2 norm of error %g Number of iterations %" PetscInt_FMT "\n", (double)norm, its));
78c4762a1bSJed Brown
79c4762a1bSJed Brown /* Free data structures */
809566063dSJacob Faibussowitsch PetscCall(KSPDestroy(&ksp));
819566063dSJacob Faibussowitsch PetscCall(VecDestroy(&u));
829566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ustar));
839566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b));
849566063dSJacob Faibussowitsch PetscCall(MatDestroy(&mat));
859566063dSJacob Faibussowitsch PetscCall(PCDestroy(&pc));
86c4762a1bSJed Brown
879566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
88b122ec5aSJacob Faibussowitsch return 0;
89c4762a1bSJed Brown }
90c4762a1bSJed Brown
91c4762a1bSJed Brown /*TEST
92c4762a1bSJed Brown
93c4762a1bSJed Brown test:
94c4762a1bSJed Brown args: -ksp_type cg -ksp_monitor_short
95c4762a1bSJed Brown
96c4762a1bSJed Brown TEST*/
97