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, NULL, 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
 64:   testset:
 65:     output_file: output/empty.out

 67:     test:
 68:       suffix: aijviennacl_1
 69:       nsize: {{1 2}}
 70:       args: -mat_type aijviennacl
 71:       requires: viennacl

 73:     test:
 74:       suffix: aijcusparse_1
 75:       nsize: {{1 2}}
 76:       args: -mat_type aijcusparse
 77:       requires: cuda

 79:     test:
 80:       suffix: kokkos
 81:       nsize: {{1 2}}
 82:       args: -mat_type aijkokkos
 83:       requires: kokkos_kernels

 85: TEST*/