CUDA-Kernel for a Dense-Sparse matrix multiplication

Question

i have been working on a very big project for some time. Within this project, I wrote my own CUDA-kernels to do various operations. One of them is to perform a sparse affine transformation to a list of sparse inputs.

Basically my input is a list of sparse vectors which are always either 1 or 0. I know for a fact that I can have at most 32 ones in a single vector.

v1 = [0, 0, 0, 1, 0, 0, 1, ...]
v2 = [1, 0, 0, 0, 0, 0, 1, ...]
...

My idea now was to wrap all these vectors into a sparse format like:

3 2 4 2 ...
-----------
1 2 3 2 
5 4 4 4
9 . 5 . 
. . 8 .

Its basically a matrix. The first row coresponds to the amount of non-zero entries. The values below are the indices of the non-zero entries.

Now when performing the matrix-vector multiplication, all I have to do for each output element is look at the input matrix, get the weights at the given indices and add them up.

So far so good. I wrote the following kernel:

__global__ void sparse_affine_kernel(
    const float*        __restrict__ mat,
    const unsigned int* __restrict__ inp_col_indices,
    const unsigned int               inp_col_max_entries,
    const float*        __restrict__ bia,
          float*        __restrict__ res,
    const unsigned int               m,
    const unsigned int               n,
    const unsigned int               lda,
    const unsigned int               ldc){

    // clang-format on
    // compute which output value we are looking at
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    // skip out of bounds
    if (col >= n || row >= m)
        return;

    // get the offset at which we look into our sparse input
    int offset = col * (inp_col_max_entries + 1);
    // check how many values we are going to read
    int count = inp_col_indices[offset];

    // track the sum
    float sum = bia[row];

    // start at offset + 1 (offset contains the amount of values to read)
    for (int i = offset + 1; i < offset + 1 + count; i++) {

        // get the sparse index (set row of the input)
        auto b_row = inp_col_indices[i];
        // get the corresponding weight
        auto wgt = mat[MATRIX_INDEX(lda, row, b_row)];

        sum += wgt;
    }
    res[MATRIX_INDEX(ldc, row, col)] = sum;
};

Now the code should be somewhat straight forward. Id like to know the following things:

1. Do you see any concrete way of improving this somewhat straight forward operation?
2. Is there anything directly related to CUDA which I could use to improve the performance of this code? Maybe using shared memory? I tried using some shared-memory some time ago and simply remembered that I wasnt able to improve the performance of the code.

I am very happy for a review and optimization-ideas for my code :)

Greetings Finn

Answer 1

Thank you for offering this for review.

I understand you're primarily interested in performance. But I confess I found the code a little on the opaque side and not quite ready to invite lots of folks to collaborate on it.

Your introductory paragraphs, outside of the code artifact, were very clear and helpful.

---

I am reading the signature. It could be improved. Starting with a URL pointing to data structure documentation, similar to your opening paragraphs. The code artifact should be self-describing.

The terse nomenclature of the signature leaves me with several questions. I do not know what "bia", "lda", & "ldc" mean. For example, when reading them, should I mentally pronounce it "bias"? linear discriminant array? linked data column? IDK. Google offered no relevant abbreviation expansions.

---

Please spell the 5th argument result. Abbreviating its single use helps no one. It was clear enough, but identifiers in a public API have a higher documentation burden than locals. Consider adhering to the convention where input args tend to appear near beginning of signature and outputs near the end.

---

The 1st two lines of code use lovely identifiers and are wonderfully clear, thank you.

---

I am skeptical about the out-of-bounds return. Maybe it is conventional and the right thing to do. In other languages I would expect an exception to be raised. Here, I don't see so much as an errno or error counter being affected. We consult a pair of block globals and a thread global. I am concerned that a subset of threads will win, higher threads will lose, and we've just offered the app developer the gift of a silent Heisenbug.

In particular, from a DbC perspective, it does not appear to me that "caller was incorrect" if we're out of bounds. So responsibility is still with the library routine to fulfill the contract. That might be "set an error flag or side-effect the matrix", but that's not what we see implemented.

---

Consider eliding the track the sum comment, as it doesn't add anything beyond what the well-chosen identifier is telling us.

I found b_row slightly puzzling. Maybe it could be bit_row? But it seems to be used where a column might be expected.

Rather than wgt, please just call it weight, and then we probably don't need the comment to explain it.

---

I imagine that MATRIX_INDEX is a macro with a few adds and multiplies (or shifts), but you did not include it. I was hoping it would help me to better understand the lda / ldc distinction.

---

Overall?

This is simple enough code, but it's not a code base I would want to assign or accept maintenance tasks for, not yet.