Using device_matrix_data and faster SpMV

[Ben90001]

Hey,
I've been working on my performance testing since my last post and a few more questions came up:

1. device_matrix_data:
   In my previous post you suggested using the device_matrix_data to assemble the matrix on the GPU. I was unable to figure out a way to utilize it with significant speedup. As the access to the arrays of the device_matrix_data strucutre come with a performance warning the only way to use it at all I could come up with was by copying a matrix_data using .create_from_host() right before creating the matrix:

gko::matrix_data<double,gko::int64> mtx_data{gko::dim<2>(N)};
(...)
mtx_data.nonzeros.emplace_back(...); // inserting all elements
(...)
auto dmd = gko::device_matrix_data<double,gko::int64>::create_from_host(exec,mtx_data);
auto pA = mtx::create(exec);
pA->read(dmd);

Is there a better way to use device_matrix_data I missed or was it intended for different usecases?

2. SpMV overhead
   Compared to DUNE-ISTL Ginkgo shows significant overhead when comparing SpMV.
   Average time per SpMV of a 1x1 Matrix:
   DUNE: 1.6ns
   Ginkgo ref: 2600ns
   Ginkgo 1omp: 3200ns (1omp being the omp executor with to OMP_NUM_THREADS=1).
   The only explanation I could come up with is the runtime kernel and operator selection mentioned here (slide 2&6).
   I am using the follwowing lambda which gets executed a set amount of time (until it runs for at least 0.25 seconds):

auto experiment_SpMV = \
      [&pA,&x,&y](){
        pA->apply(x, y);  
      };

Is there a way to further optimize this? Maybe by converting to LinOp beforehand and save some overhead or maybe set a fixed amount of times to run the SpMV for better performance?

3. Sorting of matrix_data
   In the documentation it says to make sure the stored values are sorted in row major order.
   My results (SpMV as well as jacobi and ILU0 preconditioned CG solvers) are correct.
   The print out of the Matrix Market format however when using matrix_data is still unsorted.
   It doesn't really matter too much as my results are correct, but I was still wondering as to why that is because I would assume the .mtxdata would be (re)constructed from the matrix - is the matrix_data still stored together with the matrix itself somehow?

Any hint and help is highly appreciated
Thanks in advance
Ben

[upsj]

I'll try to answer your questions:

1. I'm not sure what you mean by "performance warning", but if you want to do matrix assembly on the device, you need to know beforehand the how many non-zero entries (may contain duplicates and be unordered, we have sum_duplicates and sort_row_major for that) you want to output precisely, and then write a GPU kernel to fill those storage entries with row/column positions and values. If you use matrix_data somewhere in this process, you will not get good performance out of it, since the computations still happen on the CPU. This is an advanced topic, since it requires you to write GPU kernels and know the storage requirements of your assembly process, e.g. when doing FEM, you need to know the number of elements, and figure out an indexing schemes for your vertices, edges, faces, ...
2. Ginkgo is not intended for tiny matrices - microseconds are not really something we would consider important. We can investigate if there is a way to reduce the overhead, but this is not unexpected. If you run a profiler, you will see a number of different things happening there - runtime dispatch, check for necessary data movement, logging infrastructure for high-level and low-level operations, ...
3. Your matrix_data always needs to be sorted in row-major order, i.e. the entries need to be grouped by row in an ascending order. The column indices may be unsorted (that's what the ILU and Jacobi preconditioners have the skip_sorting parameter for), and this order will preserved by most matrix formats (except probably Hybrid and Dense). The matrix_data itself is not stored anywhere in the matrix.

[MarcelKoch]

If you are interested in an example on how to use the gko::device_matrix_data to assemble on a GPU you can take a look at the example examples/kokkos-assembly. Specifically the generate_stencil_matrix function, which uses a Kokkos::parallel_for kernel in combination with device_matrix_data to assemble a simple 3pt stencil on the GPU.