xref: /petsc/src/mat/tests/ex301.c (revision ebead697dbf761eb322f829370bbe90b3bd93fa3)
1 
2 static char help[] = "Tests for bugs in A->offloadmask consistency for GPU matrices\n\n";
3 
4 #include <petscmat.h>
5 
6 int main(int argc,char **args)
7 {
8   Mat            A;
9   PetscInt       i,j,rstart,rend,m = 3;
10   PetscScalar    one = 1.0,zero = 0.0,negativeone = -1.0;
11   PetscReal      norm;
12   Vec            x,y;
13 
14   PetscFunctionBeginUser;
15   PetscCall(PetscInitialize(&argc,&args,(char*)0,help));
16   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
17 
18   for (i=0; i<2; i++) {
19     /* Create the matrix and set it to contain explicit zero entries on the diagonal. */
20     PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
21     PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*m,m*m));
22     PetscCall(MatSetFromOptions(A));
23     PetscCall(MatSetUp(A));
24     PetscCall(MatGetOwnershipRange(A,&rstart,&rend));
25     PetscCall(MatCreateVecs(A,&x,&y));
26     PetscCall(VecSet(x,one));
27     PetscCall(VecSet(y,zero));
28     PetscCall(MatDiagonalSet(A,y,INSERT_VALUES));
29 
30     /* Now set A to be the identity using various approaches.
31      * Note that there may be other approaches that should be added here. */
32     switch (i) {
33     case 0:
34       PetscCall(MatDiagonalSet(A,x,INSERT_VALUES));
35       break;
36     case 1:
37       for (j=rstart; j<rend; j++) {
38         PetscCall(MatSetValue(A,j,j,one,INSERT_VALUES));
39       }
40       PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
41       PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
42       break;
43     case 2:
44       for (j=rstart; j<rend; j++) {
45         PetscCall(MatSetValuesRow(A,j,&one));
46       }
47       PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
48       PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
49     default:
50       break;
51     }
52 
53     /* Compute y <- A*x and verify that the difference between y and x is negligible, as it should be since A is the identity. */
54     PetscCall(MatMult(A,x,y));
55     PetscCall(VecAXPY(y,negativeone,x));
56     PetscCall(VecNorm(y,NORM_2,&norm));
57     if (norm > PETSC_SQRT_MACHINE_EPSILON) {
58       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test %" PetscInt_FMT ": Norm of error is %g, but should be near 0.\n",i,(double)norm));
59     }
60 
61     PetscCall(MatDestroy(&A));
62     PetscCall(VecDestroy(&x));
63     PetscCall(VecDestroy(&y));
64   }
65 
66   PetscCall(PetscFinalize());
67   return 0;
68 }
69 
70 /*TEST
71 
72    test:
73       suffix: aijviennacl_1
74       nsize: 1
75       args: -mat_type aijviennacl
76       requires: viennacl
77 
78    test:
79       suffix: aijviennacl_2
80       nsize: 2
81       args: -mat_type aijviennacl
82       requires: viennacl
83 
84    test:
85       suffix: aijcusparse_1
86       nsize: 1
87       args: -mat_type aijcusparse
88       requires: cuda
89 
90    test:
91       suffix: aijcusparse_2
92       nsize: 2
93       args: -mat_type aijcusparse
94       requires: cuda
95 TEST*/
96