xref: /petsc/src/mat/tests/ex3.c (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
1 
2 static char help[] = "Tests relaxation for dense matrices.\n\n";
3 
4 #include <petscmat.h>
5 
6 int main(int argc,char **args)
7 {
8   Mat            C;
9   Vec            u,x,b,e;
10   PetscInt       i,n = 10,midx[3];
11   PetscErrorCode ierr;
12   PetscScalar    v[3];
13   PetscReal      omega = 1.0,norm;
14 
15   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
16   CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-omega",&omega,NULL));
17   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
18 
19   CHKERRQ(MatCreate(PETSC_COMM_SELF,&C));
20   CHKERRQ(MatSetSizes(C,n,n,n,n));
21   CHKERRQ(MatSetType(C,MATSEQDENSE));
22   CHKERRQ(MatSetUp(C));
23   CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&b));
24   CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&x));
25   CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&u));
26   CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&e));
27   CHKERRQ(VecSet(u,1.0));
28   CHKERRQ(VecSet(x,0.0));
29 
30   v[0] = -1.; v[1] = 2.; v[2] = -1.;
31   for (i=1; i<n-1; i++) {
32     midx[0] = i-1; midx[1] = i; midx[2] = i+1;
33     CHKERRQ(MatSetValues(C,1,&i,3,midx,v,INSERT_VALUES));
34   }
35   i    = 0; midx[0] = 0; midx[1] = 1;
36   v[0] = 2.0; v[1] = -1.;
37   CHKERRQ(MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES));
38   i    = n-1; midx[0] = n-2; midx[1] = n-1;
39   v[0] = -1.0; v[1] = 2.;
40   CHKERRQ(MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES));
41 
42   CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
43   CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
44 
45   CHKERRQ(MatMult(C,u,b));
46 
47   for (i=0; i<n; i++) {
48     CHKERRQ(MatSOR(C,b,omega,SOR_FORWARD_SWEEP,0.0,1,1,x));
49     CHKERRQ(VecWAXPY(e,-1.0,x,u));
50     CHKERRQ(VecNorm(e,NORM_2,&norm));
51     CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"2-norm of error %g\n",(double)norm));
52   }
53   CHKERRQ(MatDestroy(&C));
54   CHKERRQ(VecDestroy(&x));
55   CHKERRQ(VecDestroy(&b));
56   CHKERRQ(VecDestroy(&u));
57   CHKERRQ(VecDestroy(&e));
58   ierr = PetscFinalize();
59   return ierr;
60 }
61 
62 /*TEST
63 
64    test:
65 
66 TEST*/
67