Actual source code: ex3.c
1: /*$Id: ex3.c,v 1.53 2001/08/07 03:03:43 balay Exp $*/
3: static char help[] = "Demonstrates the use of fast Richardson for SOR. And\n\
4: also tests the MatRelax() routines. Input parameters are:\n\
5: -n <n> : problem dimension\n\n";
7: #include petscksp.h
8: #include petscpc.h
12: int main(int argc,char **args)
13: {
14: Mat mat; /* matrix */
15: Vec b,ustar,u; /* vectors (RHS, exact solution, approx solution) */
16: PC pc; /* PC context */
17: KSP ksp; /* KSP context */
18: int ierr,n = 10,i,its,col[3];
19: PetscScalar value[3],one = 1.0,zero = 0.0;
20: KSPType kspname;
21: PCType pcname;
23: PetscInitialize(&argc,&args,(char *)0,help);
24: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
26: /* Create and initialize vectors */
27: VecCreateSeq(PETSC_COMM_SELF,n,&b);
28: VecCreateSeq(PETSC_COMM_SELF,n,&ustar);
29: VecCreateSeq(PETSC_COMM_SELF,n,&u);
30: VecSet(&one,ustar);
31: VecSet(&zero,u);
33: /* Create and assemble matrix */
34: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&mat);
35: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
36: for (i=1; i<n-1; i++) {
37: col[0] = i-1; col[1] = i; col[2] = i+1;
38: MatSetValues(mat,1,&i,3,col,value,INSERT_VALUES);
39: }
40: i = n - 1; col[0] = n - 2; col[1] = n - 1;
41: MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
42: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
43: MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);
44: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
45: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
47: /* Compute right-hand-side vector */
48: MatMult(mat,ustar,b);
50: /* Create PC context and set up data structures */
51: PCCreate(PETSC_COMM_WORLD,&pc);
52: PCSetType(pc,PCNONE);
53: PCSetFromOptions(pc);
54: PCSetOperators(pc,mat,mat,DIFFERENT_NONZERO_PATTERN);
55: PCSetVector(pc,u);
56: PCSetUp(pc);
58: /* Create KSP context and set up data structures */
59: KSPCreate(PETSC_COMM_WORLD,&ksp);
60: KSPSetType(ksp,KSPRICHARDSON);
61: KSPSetFromOptions(ksp);
62: KSPSetSolution(ksp,u);
63: KSPSetRhs(ksp,b);
64: PCSetOperators(pc,mat,mat,DIFFERENT_NONZERO_PATTERN);
65: KSPSetPC(ksp,pc);
66: KSPSetUp(ksp);
68: /* Solve the problem */
69: KSPGetType(ksp,&kspname);
70: PCGetType(pc,&pcname);
71: PetscPrintf(PETSC_COMM_SELF,"Running %s with %s preconditioning\n",kspname,pcname);
72: KSPSolve(ksp,&its);
73: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations %d\n",its);
75: /* Free data structures */
76: KSPDestroy(ksp);
77: VecDestroy(u);
78: VecDestroy(ustar);
79: VecDestroy(b);
80: MatDestroy(mat);
81: PCDestroy(pc);
82: PetscFinalize();
83: return 0;
84: }
85: