static char help[] = "Tests for bugs in A->offloadmask consistency for GPU matrices\n\n";

#include <petscmat.h>

int main(int argc, char **args)
{
  Mat         A;
  PetscInt    i, j, rstart, rend, m = 3;
  PetscScalar one = 1.0, zero = 0.0, negativeone = -1.0;
  PetscReal   norm;
  Vec         x, y;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &args, NULL, help));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));

  for (i = 0; i < 2; i++) {
    /* Create the matrix and set it to contain explicit zero entries on the diagonal. */
    PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * m, m * m));
    PetscCall(MatSetFromOptions(A));
    PetscCall(MatSetUp(A));
    PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
    PetscCall(MatCreateVecs(A, &x, &y));
    PetscCall(VecSet(x, one));
    PetscCall(VecSet(y, zero));
    PetscCall(MatDiagonalSet(A, y, INSERT_VALUES));

    /* Now set A to be the identity using various approaches.
     * Note that there may be other approaches that should be added here. */
    switch (i) {
    case 0:
      PetscCall(MatDiagonalSet(A, x, INSERT_VALUES));
      break;
    case 1:
      for (j = rstart; j < rend; j++) PetscCall(MatSetValue(A, j, j, one, INSERT_VALUES));
      PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
      PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
      break;
    case 2:
      for (j = rstart; j < rend; j++) PetscCall(MatSetValuesRow(A, j, &one));
      PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
      PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    default:
      break;
    }

    /* Compute y <- A*x and verify that the difference between y and x is negligible, as it should be since A is the identity. */
    PetscCall(MatMult(A, x, y));
    PetscCall(VecAXPY(y, negativeone, x));
    PetscCall(VecNorm(y, NORM_2, &norm));
    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));

    PetscCall(MatDestroy(&A));
    PetscCall(VecDestroy(&x));
    PetscCall(VecDestroy(&y));
  }

  PetscCall(PetscFinalize());
  return 0;
}

/*TEST
  testset:
    output_file: output/empty.out

    test:
      suffix: aijviennacl_1
      nsize: {{1 2}}
      args: -mat_type aijviennacl
      requires: viennacl

    test:
      suffix: aijcusparse_1
      nsize: {{1 2}}
      args: -mat_type aijcusparse
      requires: cuda

    test:
      suffix: kokkos
      nsize: {{1 2}}
      args: -mat_type aijkokkos
      requires: kokkos_kernels

TEST*/