[oneMKL][spblas] Update sparse blas API

Rbiessy commented 7 months ago

Update the sparse blas API to better support other sparse backends.

This is a large change to the existing sparse API. Most notable changes are:

Introduce handle opaque types that store the USM pointer or SYCL buffer.
Introduce matrix properties and views.
Merge gemv, symv and trmv to spmv operation; gemm is renamed to spmm; trsv is renamed to spsv.
Drop support for gemvdot as there didn't seem to be enough use-cases nor support in other backends.
Add support for external workspace provided by the user which size is queried by *_buffer_size functions.

I have rendered the spec in html and verified it looks sane.

@gajanan-choudhary and @spencerpatty FYI.

andrewtbarker commented 6 months ago

Can you share your motivation for introducing the dense_matrix_handle_t type? The dense BLAS parts of the oneAPI spec do not use such a handle and I'm not sure it's helpful for standardization, but I am willing to be convinced.

gajanan-choudhary commented 6 months ago

Can you share your motivation for introducing the dense_matrix_handle_t type? The dense BLAS parts of the oneAPI spec do not use such a handle and I'm not sure it's helpful for standardization, but I am willing to be convinced.

@andrewtbarker, both cuSPARSE and rocSPARSE have introduced "generic APIs" in the last few years. These contain (likely lightweight) wrappers over dense matrices. For example, see section 6.4. Dense Matrix APIs of cuSPARSE documentation. Mainly, the wrapper stores nrows, ncols, layout, floating point type, and the matrix pointer.

While I agree that the dense BLAS domain does not have such wrappers over dense matrices, the APIs in the sparse domain look cleaner with these wrappers. Without these wrappers, implementations of the oneAPI specification that dispatch to cuSPARSE/rocSPARSE APIs in the backend must repeatedly create and destroy the dense matrix wrapper object any and every time their API is called.

The BLAS domain is fortunately well-standardized, so adding dense_matrix_handle_t at oneapi::mkl namespace level may not be appropriate, but the sparse BLAS domain isn't, and so we'd add such a handle under the oneapi::mkl::sparse namespace to indicate it's intention of use only with sparse BLAS APIs.

andrewtbarker commented 6 months ago

@gajanan-choudhary Thanks for your explanation, that makes a lot of sense.

My only remaining concern would be to make sure that the actual data in the dense_matrix_handle_t is transparent to users so they can use these objects in BLAS functions. For sparse_matrix_handle_t it is intentionally opaque because you do not want to expose implementation details (eg CSR versus COO) to users, but for the dense matrix I think it needs to be transparent for interoperability.

spencerpatty commented 6 months ago

@gajanan-choudhary Thanks for your explanation, that makes a lot of sense.

My only remaining concern would be to make sure that the actual data in the dense_matrix_handle_t is transparent to users so they can use these objects in BLAS functions. For sparse_matrix_handle_t it is intentionally opaque because you do not want to expose implementation details (eg CSR versus COO) to users, but for the dense matrix I think it needs to be transparent for interoperability.

We will want to have constructors that can take an std::mdspan or the raw data arrays/sizes/leading_dimensions/layouts for our dense_matrix_handle_t/dense_vector_handle_t and we will want to have an operator function like the following in the dense_matrix_handle_t class that returns an mdspan view:

operator std::mdspan<T, std::dextents<2>> const {return  {data, nrow, ncols}; }

using whatever our internal representation is to construct an mdspan from our internal matrix representation, since that is the direction of C++ for dense linear algebra (whenever we start supporting mdspans :) ) ... likewise an std::mdspan<T, std::dextents<1> > const {return {data, length}; } for the dense_vector_t as well ... That way our local object is lightweight and non-owning and usable in any function that takes in mdspans ...

of course it might be that we need to return a slightly more complicated mdspan with the layout and user defined strides, but that is the idea, that it is easily interpreted as an mdspan, so it can be immediately plugged into any other function that uses mdspans.

These really could be thought of as a dense_matrix_view_t or dense_vector_view_t which are lightweight and non-owning ... so maybe we want to reconsider our name "handle" (C style name) to switch over to "view" which is C++ style name indicating it is a wrapper object that is lightweight and non-owning...

Rbiessy commented 6 months ago

@gajanan-choudhary Thanks for your explanation, that makes a lot of sense. My only remaining concern would be to make sure that the actual data in the dense_matrix_handle_t is transparent to users so they can use these objects in BLAS functions. For sparse_matrix_handle_t it is intentionally opaque because you do not want to expose implementation details (eg CSR versus COO) to users, but for the dense matrix I think it needs to be transparent for interoperability.

We will want to have constructors that can take an std::mdspan or the raw data arrays/sizes/leading_dimensions/layouts for our dense_matrix_handle_t/dense_vector_handle_t and we will want to have an operator function like the following in the dense_matrix_handle_t class that returns an mdspan view:
operator std::mdspan<T, std::dextents<2>> const {return  {data, nrow, ncols}; }
using whatever our internal representation is to construct an mdspan from our internal matrix representation, since that is the direction of C++ for dense linear algebra (whenever we start supporting mdspans :) ) ... likewise an std::mdspan<T, std::dextents<1> > const {return {data, length}; } for the dense_vector_t as well ... That way our local object is lightweight and non-owning and usable in any function that takes in mdspans ...

of course it might be that we need to return a slightly more complicated mdspan with the layout and user defined strides, but that is the idea, that it is easily interpreted as an mdspan, so it can be immediately plugged into any other function that uses mdspans.

These really could be thought of as a dense_matrix_view_t or dense_vector_view_t which are lightweight and non-owning ... so maybe we want to reconsider our name "handle" (C style name) to switch over to "view" which is C++ style name indicating it is a wrapper object that is lightweight and non-owning...

This is interesting. I like the idea of making the dense handles compatible with mdspan. Given that there are some open questions on how this will work exactly and that the rest of oneMKL does not support mdspan yet, is it fair to add this feature in the future?

I wouldn't mind renaming handle to view. I agree that this would better fit C++ however these types are unfortunately closer to C-style with "init" and "destroy" functions. We would also need to rename the current matrix_view to something else. I'd say handle is fine since we document that it is not owning the data.