How to Optimize a CUDA Matmul Kernel for CuBLAS-Like Performance: A Worklog | Hacker News Reader

How to Optimize a CUDA Matmul Kernel for CuBLAS-Like Performance: A Worklog

siboehm.com

130 pointstodsacerdoti16 comments

This is a great post for people who are new to optimizing GPU code.

It is interesting to see that the author got this far without interchanging the innermost loop over k to the outermost loop, as is done in CUTLASS (https://github.com/NVIDIA/cutlass).

As you can see in this blog post the code ends up with a lot of compile-time constants (e.g. BLOCKSIZE, BM, BN, BK, TM, TN) one way to optimize this code further is to use an auto-tuner to find the optimal value for all of these parameters for your GPU and problem size, for example Kernel Tuner (https://github.com/KernelTuner/kernel_tuner)

Kernel Tuner is great! Remember going to a tutorial at SC21. Would highly recommend the tutorials they used to get familiar as well (https://github.com/KernelTuner/kernel_tuner_tutorial)

I did pretty much the same in HLSL just a few days ago. Here’s one idea to try.

You don’t actually need 3 tiles in groupshared memory.

With some care, it’s possible to go with 2 tiles there. One tile is loaded from the first source matrices. Another tile is accumulator for the output matrix tile. And then in the main function of the compute shader load values for the second source matrix from the global memory, and update all affected elements of the output tile with these mad() instructions.

Shader model 5.0 limits amount of group shared memory to 32kb, and that streaming trick allows to push to the limit, with 64x64 tiles. I think I got about 15% performance win from that trick on some of the GPUs I’m testing my software.

He mentions having bank conflicts in shared memory, there is a simple trick to avoid that. You allocate an extra column for your shared buffers:

    __shared__ float As[(CHUNKSIZE+1) * CHUNKSIZE];

If CHUNKSIZE is a power of 2 this will eliminate bank conflicts.

Author here: Seems like a good trick! Though won't this affect shared memory alignment and make me loose those LDS.128 instructions? Or do these not require alignment? There's so little good docs on SASS.

In general I'm still confused about whether vectorized load instructions (LDS.128) necessarily lead to bank conflicts or not. My impression was that consecutive 32b floats get mapped to different banks, so to avoid conflicts I'd want the warp to load 32*32b consecutive elements at each step.

Hmm, I think you might have to adjust the padding to be 128 bits then:

    __shared__ float As[(CHUNKSIZE+4) \* CHUNKSIZE]

Ultimately it's down to trial and error, like always with GPGPU.

I am curious about doing the same kind of thing for compute shaders. I'm aware of Kompute.cc (which is Vulkan based) but haven't looked at their GEMM kernels, and also of wonnx for WebGPU ([1] is their GEMM code).

I'm also curious whether warp shuffle operations might be useful to reduce some of the shared memory traffic.

[1]: https://github.com/webonnx/wonnx/blob/master/wonnx/templates...

I'll be benchmarking the WONNX GEMM in the near future - would be curious to see how far a GEMM on WebGPU can be taken.

Bram wrote some interesting exploration here: https://jott.live/markdown/webgpu_safari

Those diagrams are pretty good, wonder how did the author make them

In the conclusion, the author thanks https://excalidraw.com/.

At the end of the post, he links to excalidraw[0]

[0] https://excalidraw.com/

I implemented the same for a class last year, goodness, wish I had these visualizations for when I was banging my head against the wall to max out the roofline.

Wow, what class?

https://sites.google.com/eng.ucsd.edu/cse260fall2022/home

Albeit, much of the boilerplate was in place. Students were tasked with maximizing bandwidth via the optimization techniques mentioned in the article and other literature.

Really well written blog, explaining all the concepts introduced in details.

>> Profiling our naive kernel

I would highly recommend trying out NVIDIA NSight Systems to help profiling: it gives a lot of information while optimizing for CUDA (together with info on CPU, NVTX markers, Vulkan,OpenGL etc).

And then there was Nervana Systems's maxas effort that, in Maxwell days, exceeded cuBLAS and was edging theoretical FLOPs despite the penalty paid for address calculations which on that architecture compete with single precision FLOPS.

Reply on news.ycombinator.com