1: static const char help[] = "Tests MatDenseGetColumnVec() and friends on dense matrices created from a VecType, including device dense matrices whose vectors are VECKOKKOS\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat                A;
  8:   Vec                v;
  9:   char               vtype[64] = VECSTANDARD;
 10:   PetscInt           M = 9, N = 3, lda, rstart, rend, i, j;
 11:   PetscReal          norm;
 12:   PetscScalar        sum;
 13:   const PetscScalar *array;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 17:   PetscCall(PetscOptionsGetString(NULL, NULL, "-vec_type", vtype, sizeof(vtype), NULL));
 18:   /* A VECKOKKOS type gives a MATDENSECUDA/MATDENSEHIP with VECKOKKOS vectors when the Kokkos backend runs on that device */
 19:   PetscCall(MatCreateDenseFromVecType(PETSC_COMM_WORLD, vtype, PETSC_DECIDE, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, &A));
 20:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));

 22:   /* Write column j through its column vector, then scale it through a read/write column vector: A(:, j) = 2 (j + 1) */
 23:   for (j = 0; j < N; j++) {
 24:     PetscCall(MatDenseGetColumnVecWrite(A, j, &v));
 25:     PetscCall(VecSet(v, (PetscScalar)(j + 1)));
 26:     PetscCall(MatDenseRestoreColumnVecWrite(A, j, &v));
 27:   }
 28:   for (j = 0; j < N; j++) {
 29:     PetscCall(MatDenseGetColumnVec(A, j, &v));
 30:     PetscCall(VecScale(v, 2.0));
 31:     PetscCall(MatDenseRestoreColumnVec(A, j, &v));
 32:   }

 34:   /* Check through read-only column vectors, and independently through the matrix itself */
 35:   for (j = 0; j < N; j++) {
 36:     PetscCall(MatDenseGetColumnVecRead(A, j, &v));
 37:     PetscCall(VecSum(v, &sum));
 38:     PetscCheck(PetscAbsScalar(sum - 2.0 * (j + 1) * M) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Column %" PetscInt_FMT " read back through MatDenseGetColumnVecRead() sums to %g, expected %g", j, (double)PetscRealPart(sum), 2.0 * (j + 1) * M);
 39:     PetscCall(MatDenseRestoreColumnVecRead(A, j, &v));
 40:   }
 41:   PetscCall(MatNorm(A, NORM_INFINITY, &norm));
 42:   PetscCheck(PetscAbsReal(norm - N * (N + 1)) < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatNorm() gives %g, expected %g", (double)norm, (double)(N * (N + 1)));
 43:   PetscCall(MatDenseGetLDA(A, &lda));
 44:   PetscCall(MatDenseGetArrayRead(A, &array));
 45:   for (j = 0; j < N; j++) {
 46:     for (i = 0; i < rend - rstart; i++)
 47:       PetscCheck(array[i + j * lda] == 2.0 * (j + 1), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Entry (%" PetscInt_FMT ", %" PetscInt_FMT ") is %g, expected %g", rstart + i, j, (double)PetscRealPart(array[i + j * lda]), 2.0 * (j + 1));
 48:   }
 49:   PetscCall(MatDenseRestoreArrayRead(A, &array));

 51:   PetscCall(MatDestroy(&A));
 52:   PetscCall(PetscFinalize());
 53:   return 0;
 54: }

 56: /*TEST

 58:   testset:
 59:     output_file: output/empty.out
 60:     nsize: {{1 2}}
 61:     test:
 62:       suffix: standard
 63:       args: -vec_type standard
 64:     test:
 65:       suffix: cuda
 66:       requires: cuda
 67:       args: -vec_type cuda
 68:     test:
 69:       suffix: hip
 70:       requires: hip
 71:       args: -vec_type hip
 72:     # On a CUDA or HIP build this is a MATDENSECUDA or MATDENSEHIP with VECKOKKOS column vectors
 73:     test:
 74:       suffix: kokkos
 75:       requires: kokkos_kernels !sycl
 76:       args: -vec_type kokkos

 78: TEST*/