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