1: static char help[] = "A placeholder for testing device matrices with over 2 billion nonzeros\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat          A, B;
  8:   Vec          x1, x2, y1, y2;
  9:   PetscInt     m = 1 << 6;
 10:   PetscInt     n = (1 << 6) + 2, nnz; // or m = 1 << 15, n = (1 << 16) + 2 to get > 2 billion nonzeros
 11:   PetscInt    *i, *j;
 12:   PetscScalar *a;
 13:   PetscReal    r1, r2;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 18:   // Create a 'dense' SEQAIJ matrix for simplicity, and also to save the row pointer memory cost
 19:   nnz = m * n;
 20:   PetscCall(PetscMalloc3(m + 1, &i, nnz, &j, nnz, &a));

 22:   i[0] = 0;
 23:   for (PetscInt k = 0; k < m; k++) {
 24:     i[k + 1] = i[k] + n;
 25:     for (PetscInt l = 0; l < n; l++) {
 26:       j[i[k] + l] = l;
 27:       a[i[k] + l] = (PetscScalar)(i[k] + l);
 28:     }
 29:   }
 30:   PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, m, n, i, j, a, &A));

 32:   PetscCall(MatCreateVecs(A, &x1, &y1));
 33:   PetscCall(VecSetRandom(x1, NULL));
 34:   PetscCall(MatMult(A, x1, y1));

 36:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
 37:   PetscCall(MatSetFromOptions(B));
 38:   PetscCall(MatDestroy(&A));
 39:   PetscCall(PetscFree3(i, j, a));

 41:   PetscCall(MatCreateVecs(B, &x2, &y2));
 42:   PetscCall(VecCopy(x1, x2));
 43:   PetscCall(MatMult(B, x2, y2));

 45:   PetscCall(VecNorm(y1, NORM_INFINITY, &r1));
 46:   PetscCall(VecAXPY(y2, -1.0, y1));
 47:   PetscCall(VecNorm(y2, NORM_INFINITY, &r2));
 48:   r2 /= r1;

 50:   PetscCheck(r2 < PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMult wrong with indices beyond 32-bit: relative error %g", (double)r2);

 52:   PetscCall(MatDestroy(&B));
 53:   PetscCall(VecDestroy(&x1));
 54:   PetscCall(VecDestroy(&x2));
 55:   PetscCall(VecDestroy(&y1));
 56:   PetscCall(VecDestroy(&y2));
 57:   PetscCall(PetscFinalize());
 58:   return 0;
 59: }

 61: /*TEST
 62:   testset:
 63:     nsize: 1
 64:     output_file: output/empty.out

 66:     test:
 67:       requires: kokkos_kernels single defined(PETSC_USE_64BIT_INDICES)
 68:       suffix: kokkos
 69:       args: -mat_type aijkokkos

 71:     test:
 72:       requires: cuda single defined(PETSC_USE_64BIT_INDICES)
 73:       suffix: cuda
 74:       args: -mat_type aijcusparse

 76:     test:
 77:       requires: hip single defined(PETSC_USE_64BIT_INDICES)
 78:       suffix: hip
 79:       args: -mat_type aijhipsparse
 80: TEST*/