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)
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.
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.
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.
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.
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.
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).
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.
16 comments
[ 3.3 ms ] story [ 46.1 ms ] threadIt 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)
[0] https://excalidraw.com/
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.
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.
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.
>> 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).
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...
Bram wrote some interesting exploration here: https://jott.live/markdown/webgpu_safari