Actual source code: ex2.c
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";
5: #include <petscksp.h>
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;
19: PetscInitialize(&argc,&args,(char*)0,help);
20: /* Create and initialize vectors */
21: VecCreateSeq(PETSC_COMM_SELF,n,&b);
22: VecCreateSeq(PETSC_COMM_SELF,n,&ustar);
23: VecCreateSeq(PETSC_COMM_SELF,n,&u);
24: VecSet(ustar,1.0);
25: VecSet(u,0.0);
27: /* Create and assemble matrix */
28: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&mat);
29: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
30: for (i=1; i<n-1; i++) {
31: col[0] = i-1; col[1] = i; col[2] = i+1;
32: MatSetValues(mat,1,&i,3,col,value,INSERT_VALUES);
33: }
34: i = n - 1; col[0] = n - 2; col[1] = n - 1;
35: MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
36: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
37: MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
38: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
39: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
41: /* Compute right-hand-side vector */
42: MatMult(mat,ustar,b);
44: /* Create PC context and set up data structures */
45: PCCreate(PETSC_COMM_WORLD,&pc);
46: PCSetType(pc,PCNONE);
47: PCSetFromOptions(pc);
48: PCSetOperators(pc,mat,mat);
49: PCSetUp(pc);
51: /* Create KSP context and set up data structures */
52: KSPCreate(PETSC_COMM_WORLD,&ksp);
53: KSPSetType(ksp,KSPRICHARDSON);
54: KSPSetFromOptions(ksp);
55: PCSetOperators(pc,mat,mat);
56: KSPSetPC(ksp,pc);
57: KSPSetUp(ksp);
59: /* Solve the problem */
60: KSPGetType(ksp,&kspname);
61: PCGetType(pc,&pcname);
62: PetscPrintf(PETSC_COMM_SELF,"Running %s with %s preconditioning\n",kspname,pcname);
63: KSPSolve(ksp,b,u);
64: VecAXPY(u,-1.0,ustar);
65: VecNorm(u,NORM_2,&norm);
66: KSPGetIterationNumber(ksp,&its);
67: if (norm > tol) {
68: PetscPrintf(PETSC_COMM_SELF,"2 norm of error %g Number of iterations %D\n",(double)norm,its);
69: }
71: /* Free data structures */
72: KSPDestroy(&ksp);
73: VecDestroy(&u);
74: VecDestroy(&ustar);
75: VecDestroy(&b);
76: MatDestroy(&mat);
77: PCDestroy(&pc);
79: PetscFinalize();
80: return 0;
81: }
83: /*TEST
85: test:
86: args: -ksp_type cg -ksp_monitor_short
88: TEST*/