Actual source code: ex301.c
2: static char help[] = "Tests for bugs in A->offloadmask consistency for GPU matrices\n\n";
4: #include <petscmat.h>
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;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
16: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
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));
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++) PetscCall(MatSetValue(A, j, j, one, INSERT_VALUES));
38: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
39: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
40: break;
41: case 2:
42: for (j = rstart; j < rend; j++) PetscCall(MatSetValuesRow(A, j, &one));
43: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
44: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
45: default:
46: break;
47: }
49: /* Compute y <- A*x and verify that the difference between y and x is negligible, as it should be since A is the identity. */
50: PetscCall(MatMult(A, x, y));
51: PetscCall(VecAXPY(y, negativeone, x));
52: PetscCall(VecNorm(y, NORM_2, &norm));
53: 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));
55: PetscCall(MatDestroy(&A));
56: PetscCall(VecDestroy(&x));
57: PetscCall(VecDestroy(&y));
58: }
60: PetscCall(PetscFinalize());
61: return 0;
62: }
64: /*TEST
66: test:
67: suffix: aijviennacl_1
68: nsize: 1
69: args: -mat_type aijviennacl
70: requires: viennacl
72: test:
73: suffix: aijviennacl_2
74: nsize: 2
75: args: -mat_type aijviennacl
76: requires: viennacl
78: test:
79: suffix: aijcusparse_1
80: nsize: 1
81: args: -mat_type aijcusparse
82: requires: cuda
84: test:
85: suffix: aijcusparse_2
86: nsize: 2
87: args: -mat_type aijcusparse
88: requires: cuda
89: TEST*/