oneMKL: transposed sparse gemv reads y when beta is zero #703

// oneMKL sparse::gemv on a CSR handle with op = trans/conjtrans reads y when beta == 0
#include <sycl/sycl.hpp>
#include <oneapi/mkl.hpp>
#include <mkl.h>
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <utility>

using oneapi::mkl::transpose;
namespace sparse = oneapi::mkl::sparse;

int main() {
    sycl::queue q{sycl::default_selector_v, sycl::property::queue::in_order()};
    char ver[256];
    mkl_get_version_string(ver, sizeof(ver));
    auto dev = q.get_device();
    printf("%s\n%s (driver %s)\n\n", ver, dev.get_info<sycl::info::device::name>().c_str(),
           dev.get_info<sycl::info::device::driver_version>().c_str());

    // 4x3 matrix in one-based CSR; row 2 and column 2 are empty
    //   [1 0 4]
    //   [0 0 0]
    //   [0 0 2]
    //   [3 0 0]
    const int64_t m = 4, n = 3, nnz = 4;
    int32_t *row_ptr = sycl::malloc_shared<int32_t>(m + 1, q);
    int32_t *col_ind = sycl::malloc_shared<int32_t>(nnz, q);
    float *val = sycl::malloc_shared<float>(nnz, q);
    const int32_t rp[] = {1, 3, 3, 4, 5}, ci[] = {1, 3, 3, 1};
    const float v[] = {1, 4, 2, 3};
    std::copy(rp, rp + m + 1, row_ptr);
    std::copy(ci, ci + nnz, col_ind);
    std::copy(v, v + nnz, val);

    sparse::matrix_handle_t A;
    sparse::init_matrix_handle(&A);
    sparse::set_csr_data(q, A, m, n, nnz, oneapi::mkl::index_base::one, row_ptr, col_ind, val).wait();

    // y = 1 * op(A) * x + 0 * y with x = ones and y prefilled; gemm uses 2 columns
    const float expect_N[] = {5, 0, 2, 3}, expect_T[] = {4, 0, 6};
    float *x = sycl::malloc_shared<float>(2 * m, q), *y = sycl::malloc_shared<float>(2 * m, q);
    const std::pair<transpose, const char *> ops[] = {
        {transpose::nontrans, "nontrans"}, {transpose::trans, "trans"}, {transpose::conjtrans, "conjtrans"}};
    for (int64_t ncols : {1, 2})
        for (auto [op, opname] : ops)
            for (float fill : {NAN, 42.0f}) {
                const bool N = op == transpose::nontrans;
                const int64_t nx = N ? n : m, ny = N ? m : n;
                std::fill(x, x + nx * ncols, 1.0f);
                std::fill(y, y + ny * ncols, fill);
                if (ncols == 1)
                    sparse::gemv(q, op, 1.0f, A, x, 0.0f, y).wait();
                else
                    sparse::gemm(q, oneapi::mkl::layout::col_major, op, transpose::nontrans, 1.0f, A,
                                 x, ncols, nx, 0.0f, y, ny).wait();
                bool ok = true;
                printf("%s %-9s y_in=%-3g y =", ncols == 1 ? "gemv" : "gemm", opname, fill);
                for (int64_t i = 0; i < ny * ncols; i++) {
                    printf(" %g", y[i]);
                    ok &= y[i] == (N ? expect_N : expect_T)[i % ny];
                }
                printf("  %s\n", ok ? "OK" : "WRONG");
            }

    sparse::release_matrix_handle(q, &A).wait();
    for (void *p : {(void *)row_ptr, (void *)col_ind, (void *)val, (void *)x, (void *)y})
        sycl::free(p, q);
}