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