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
6: #include <stdio.h>
10: int main(int argc,char **args)
11: {
12: Mat mat; /* matrix */
13: Vec b,ustar,u; /* vectors (RHS, exact solution, approx solution) */
14: PC pc; /* PC context */
15: KSP ksp; /* KSP context */
17: PetscInt n = 10,i,its,col[3];
18: PetscScalar value[3];
19: const PCType pcname;
20: const KSPType kspname;
21: PetscReal norm;
23: PetscInitialize(&argc,&args,(char *)0,help);
25: /* Create and initialize vectors */
26: VecCreateSeq(PETSC_COMM_SELF,n,&b);
27: VecCreateSeq(PETSC_COMM_SELF,n,&ustar);
28: VecCreateSeq(PETSC_COMM_SELF,n,&u);
29: VecSet(ustar,1.0);
30: VecSet(u,0.0);
32: /* Create and assemble matrix */
33: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&mat);
34: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
35: for (i=1; i<n-1; i++) {
36: col[0] = i-1; col[1] = i; col[2] = i+1;
37: MatSetValues(mat,1,&i,3,col,value,INSERT_VALUES);
38: }
39: i = n - 1; col[0] = n - 2; col[1] = n - 1;
40: MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
41: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
42: MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
43: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
44: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
46: /* Compute right-hand-side vector */
47: MatMult(mat,ustar,b);
49: /* Create PC context and set up data structures */
50: PCCreate(PETSC_COMM_WORLD,&pc);
51: PCSetType(pc,PCNONE);
52: PCSetFromOptions(pc);
53: PCSetOperators(pc,mat,mat,DIFFERENT_NONZERO_PATTERN);
54: PCSetUp(pc);
56: /* Create KSP context and set up data structures */
57: KSPCreate(PETSC_COMM_WORLD,&ksp);
58: KSPSetType(ksp,KSPRICHARDSON);
59: KSPSetFromOptions(ksp);
60: PCSetOperators(pc,mat,mat,DIFFERENT_NONZERO_PATTERN);
61: KSPSetPC(ksp,pc);
62: KSPSetUp(ksp);
64: /* Solve the problem */
65: KSPGetType(ksp,&kspname);
66: PCGetType(pc,&pcname);
67: PetscPrintf(PETSC_COMM_SELF,"Running %s with %s preconditioning\n",kspname,pcname);
68: KSPSolve(ksp,b,u);
69: VecAXPY(u,-1.0,ustar);
70: VecNorm(u,NORM_2,&norm);
71: KSPGetIterationNumber(ksp,&its);
72: PetscPrintf(PETSC_COMM_SELF,"2 norm of error %A Number of iterations %D\n",norm,its);
74: /* Free data structures */
75: KSPDestroy(ksp);
76: VecDestroy(u);
77: VecDestroy(ustar);
78: VecDestroy(b);
79: MatDestroy(mat);
80: PCDestroy(pc);
82: PetscFinalize();
83: return 0;
84: }
85: