1 static char help[] = "Tests for bugs in A->offloadmask consistency for GPU matrices\n\n"; 2 3 #include <petscmat.h> 4 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