Actual source code: ex301.c
1: static char help[] = "Tests for bugs in A->offloadmask consistency for GPU matrices\n\n";
3: #include <petscmat.h>
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;
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
15: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
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));
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: }
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));
54: PetscCall(MatDestroy(&A));
55: PetscCall(VecDestroy(&x));
56: PetscCall(VecDestroy(&y));
57: }
59: PetscCall(PetscFinalize());
60: return 0;
61: }
63: /*TEST
65: test:
66: suffix: aijviennacl_1
67: nsize: 1
68: args: -mat_type aijviennacl
69: requires: viennacl
71: test:
72: suffix: aijviennacl_2
73: nsize: 2
74: args: -mat_type aijviennacl
75: requires: viennacl
77: test:
78: suffix: aijcusparse_1
79: nsize: 1
80: args: -mat_type aijcusparse
81: requires: cuda
83: test:
84: suffix: aijcusparse_2
85: nsize: 2
86: args: -mat_type aijcusparse
87: requires: cuda
88: TEST*/