ap4

I have deleted it. The function verify_matrix() from the original SGEMM_CUDA repository did not check for NANs, and the kernel was returning NANs.

I have no way to delete this submission.

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.

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.

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().

For all doubters, I open-sourced it: https://github.com/arekpaterek/Faster_SGEMM_CUDA

I shared the code: https://github.com/arekpaterek/Faster_SGEMM_CUDA

I did more tests on various GPUs. The TFLOPS values for A100 exceed the maximum from the spec of the GPU, so perhaps TFLOPS is calculated in a different way in the spec.

  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
  16384²  28.8           34.9       +21%      3090ti
  12288²  28.8           34.5       +20%      3090ti
  8192²   29.3           33.3       +14%      3090ti
  4096²   27.9           26.7       -4%       3090ti

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

Edit: The values for A100 are fishy. They should not exceed 19.5 TFLOPS, which is the maximum from the spec. I have to take a closer look into that.

Also on other GPUS.

So far I have tested on:

  size    tflops_cublas  tflops_my  diff      gpu
  4096²   28.7-28.8      32.5       +13%      4070ts
  8192²   27.7-28.2      33.5       +19-21%   4070ts
  4096²   9.9-10.0       10.1-10.2  +1-2%     1080ti
  4096²   3.8-4.3        6.7        +56-76%   T4

Related: I created a CUDA kernel typically much faster than kernels from cuBLAS when multiplying large square float32 matrices. Tested mostly on a 4090 GPU so far.

Source code: https://github.com/arekpaterek/Faster_SGEMM_CUDA

  size    tflops_cublas  tflops_my  diff      gpu
  4096²   50.8-50.9      61.8       +21%      4090
  6144²   55.3           59.8       +8%       4090
  8192²   56.3-56.5      67.1       +19%      4090
  12288²  53.7           66.7       +24%      4090
  16384²  53.6           66.7       +24%      4090
  4096²   28.7-28.8      32.5       +13%      4070ts
  4096²   3.8-4.3        6.7        +56-76%   T4

Believe it or not.

On a 4090 gpu, average of 20 runs of SGEMM_CUDA:

  size    tflops_cublas  tflops_my  diff
  4096²   50.8-50.9      61.8       +21%
  8192²   56.3-56.4      67.1       +19%
  16384²  53.6           66.7       +24%

I guess the right thing to do now would be to hire a B2B salesman and figure out, which company needs it.

It is apples-to-apples. The code from the article runs a chosen kernel 50 times in a loop. In the case of CuBLAS it runs the function cublasGemmEx() 50 times.

I swapped one of the kernels in the code from the article to my kernel, and left only the multiplication of matrices of size 4096².

On average over 20 runs:

CuBLAS (./sgemm 0) has 50.9 TFLOPS.

My kernel has 61.8 TFLOPS, so it's actually +21% speedup in this benchmark.

How do I collect my paycheck?

Beating CuBLAS is easy. I wrote a 30-line kernel that multiplies two 4096² matrices 14% faster than CuBLAS on a 4090. The question is how to earn money on that.