Tested on more GPUs. The biggest improvement found so far over the standard matrix multiplication from the cuBLAS library is +43% when multiplying matrices of size 12288² on an A100 GPU.

  size    tflops_cublas  tflops_my  diff      gpu
  12288²  51.4           56.3       +9%       h100
  8192²   50.5           56.1       +11%      h100
  4096²   43.8           53.9       +23%      h100
  12288²  18.9           27.0       +43%      a100
  8192²   19.0           26.3       +38%      a100
  4096²   17.5           19.8       +13%      a100
  12288²  28.8           34.5       +20%      3090ti

Your `sum` array is only 64 elements but you're indexing with indices out of bounds, which is UB and the compiler knows it at compile time so it's skipping a bunch of work.

E.g. consider the line:

  sum[y2*8 + x2] += ...

In the final loop iteration when y2=15 and x2=7, the index is 127.

That's the trick. This is intended. The point is that the compiler does not notice oob access in the first stage, but notices it in the later stages, and compiles the code to a correctly working kernel. The result is correct, as checked by the function verify_matrix().

Sorry I'm a bit lost here, could you explain the reasoning behind this and why it works?

I want to have 512 threads per block, each thread calculating simultaneously 128 values. That's 65536 values per block. I can't accumulate each of these values in registers, because the GPU has the limit of max 65536 registers per block, and some additional registers are needed in the kernel. But if I find a way to trick the first stages of the compiler that it has sufficient amount of free registers, then sometimes, like in the case of this kernel, the later stages of the compiler are sufficiently smart to give me what I want: 512 threads per block, each calculating 128 values.

I hate to say it but that simply doesn't work: you can't write out of bounds to trick the compiler, it'll just ignore your out of bounds work.

You can look at the generated sass on godbolt: https://cuda.godbolt.org/z/19excTxM3

Note that there are 1024 FFMA instructions in the loop but you would expect 16*8*BK = 2048. This would suggest half the operations are skipped, which lines up with the half of writes that are out of bounds being omitted.

After the compute loop when you're calculating the final result and storing it, you can see that the FFMAs referencing out of bounds indices write QNAN instead of any real results.

Is it possible that the NANs are what are messing with your tests? Those are notoriously hard to deal with correctly, but you should assert that the result doesn't have any NANs whatsoever.

You are right. The function verify_matrix() from the original SGEMM_CUDA repository did not check for NANs. I deleted the repository. It was the 13th CUDA kernel I wrote in my life, and the whole endeavor teached me a lot. I appreciate the feedback.

Glad it was a learning experience for you and I apologize if I came off argumentative at all! I was mainly so incredulous because this is my day job haha, so I have a bit more experience than most in the area.

It definitely sucks to be led astray and have time wasted by a bug inherited from the original repo though, sorry to hear that :/