# how-to-optimize-gemm: row-major matrix multiplication, one kernel at a time

> This repository walks through optimizing a row-major matrix multiply by hand for aarch64, armv7, x86, CUDA, CUDA bf16, int8 and Vulkan, keeping the old and the new kernel side by side and plotting the results. The CUDA version is claimed to beat cuBLAS on a 3080, and the int8 work was merged into ncnn.

**tpoisonooo/how-to-optimize-gemm** — row-major matmul optimization

- Repository: https://github.com/tpoisonooo/how-to-optimize-gemm
- Website: https://zhuanlan.zhihu.com/p/65436463
- Stars: 759 · Forks: 94
- Language: C++
- License: GPL-3.0
- Published: 2026-10-21 · Updated: 2026-10-21 · Language: en
- Canonical page: https://hysenlabs.com/projects/tpoisonooo-how-to-optimize-gemm

## Every backend builds one old and one new kernel

how-to-optimize-gemm is a row-major matrix multiplication optimization tutorial, and the structure of the repository is what makes it work as a teaching tool. The top level holds one directory per backend, aarch64, aarch64-int8, armv7, cuda, cuda-bf16, cuda-int4, vulkan and x86, plus an OpenBLAS-0.2.20 directory, an mperf submodule, a MatrixMul.pdf and an images folder for the result plots.

The support table in the source README covers aarch64, aarch64-int8, armv7, cuda, cuda-bf16 and vulkan, with the x86 entry tracked as a separate line, and the tutorial table maps each backend to a written article. The aarch64 entries have two, a GEMM introduction and a GEMM caching article, while the armv7, cuda, cuda-int4 and vulkan entries each point at one article, and aarch64-int8 and cuda-bf16 have no linked tutorial at all.

The mechanism is the same in every backend directory. The makefile has an OLD variable and a NEW variable, each naming a kernel implementation, and on a first run you set both to the same implementation so you have a working baseline:

```bash
$ cd aarch64
$ cat makefile
OLD    := MMult_4x4_10
NEW   := MMult_4x4_10
..
```

Then a single target compiles and runs the kernel named by NEW and copies the result to output_new.m:

```bash
$ make run
$ cat output_new.m
```

The idea is that you change one kernel, rebuild, and compare. The OLD variable is where the previous version is kept, so the difference between two runs of the same benchmark is the difference between two kernels rather than the difference between two machines.

## Plotting the numbers instead of reading them

Reading raw numbers off a terminal is not intuitive, so the repository includes a plotting step. You install the two requirements, numpy and matplotlib, and run the script:

```bash
$ python3 -m pip install -r ../requirements.txt
$ python3 plot.py
```

The requirements file is genuinely minimal, containing just numpy and matplotlib, which tells you the plotter reads the .m result files that the makefile produced and nothing else. The generated images are referenced in the source README, including a comparison of an int8 implementation against the rk3399 single core fp32 peak and a comparison of the CUDA implementation against cuBLAS.

That last comparison is the one the repository leads with. The README states that the cuda version exceeds the speed of NVIDIA cuBLAS, and the image caption says the green line is the project's own implementation, the blue line is cuBLAS, and that both are on a 3080 without tensor cores.

Two caveats appear in the same section. You need to install the CUDA driver and nvcc yourself, and you need a CPU OpenBLAS build to serve as the baseline that verifies numerical correctness. The second point is the one that matters for anyone reproducing the result: a fast kernel that computes the wrong thing is worthless, and the OpenBLAS baseline is what tells you the kernel is right before the speed numbers mean anything.

On the ARM side the setup is lighter. The default architecture variable is native, so on an ARM Linux machine, whether a Raspberry Pi, an rk3399 board or an AWS ARM server, you can compile and run directly with:

```bash
$ cd armv8 && make run
```

The project also records an update note saying that in August 2023 the aarch64 backend gained cmake and mperf, with mperf enabled through -DMPERF_ENABLE=ON, because otherwise the performance work has no basis.

## The int8 variant and its boundary handling

The aarch64-int8 directory exists because readers of the GEMM introduction asked how to write an int8 gemm. The repository points at a separate project, chgemm, described as a usable int8 gemm library, and describes how it differs from the tutorial code.

There are four differences. The first is that edge cases are handled, unlike the tutorial code which only considers sizes that are multiples of four. The second is the measured throughput: the int8 path reaches 18.6 gflops, compared with a theoretical fp32 limit of 14.3 and a figure of roughly 12 to 14 gflops for gemmlowp. The third is the quantisation scheme, which is symmetric, so the input value range must fall inside [-127, +127] and the value -128 cannot appear. The fourth is that the library includes a small example of how to integrate it into Android Studio.

The performance figure is the interesting part of the comparison, and the plot it comes from is described in the README. The blue line is the chgemm implementation and the orange line is the rk3399 single core fp32 peak.

The claim that the chgemm work has been merged into the ncnn INT8 convolution implementation is worth noting for a practical reason. It is a third party project of the same author rather than part of this repository, so if you want the production int8 kernel rather than the tutorial one, ncnn is where it ended up.

## What the x86 and bf16 versions do and do not do

The x86 implementation is derived from a project called flame, described in the README as the original implementation, and the differences are listed. The original was a column-major SSE version, while this repository's version is row-major. The current MMult_4x4_17.c reaches 70 percent of CPU peak, according to the README. Edge cases are not handled, and only cases where MNK are all multiples of four are considered, with the sub_kernel written only in the simplest assembly form, so practical use needs some adjustment. The plotting also dropped octave, because configuring an environment on an embedded device is troublesome, in favour of python.

The CUDA bf16 version has a different input contract. The input is not a standard GEMM but is aligned with nn.Linear, where the input shape is [M,K], the weight shape is [N,K] and the dtype is bf16.

The performance numbers there are honest about the gap. With random initialisation, MMult_cuda_bf16_5 at H800 4096 cubed is about 65 TFLOPS, which is roughly 12 times below nn.Linear at around 766 TFLOPS, against a theoretical peak of 1979 TFLOPS. The README's summary of that is that the current approach only gets them that far, and that from version six onwards you have to write it yourself.

The Vulkan backend depends on the kompute API wrapper for the build, with the build documentation in the vulkan directory of the repository. The result is that there was not time to reach peak performance, and the value of the backend is more as an introduction to learning compute shaders.

None of these three are finished products, and the README is consistent about that across all of them.

## The tooling the project points you at

The README closes with a short list of tools, and it reads like a reading list rather than a dependency list, since none of these are required to build the kernels.

megpeak is listed for measuring hardware limit performance, with support for arm, x86 and OpenCL. perf is pointed at as being in the basic Linux package, for system level performance analysis and disassembly. YHs_Sample is credited as an implementation by a well known author. mperf is described as a performance optimization guide, and it is the one that is actually vendored into the tree as a submodule alongside OpenBLAS.

The pattern across the whole repository is that measurement comes first. A kernel is written, benchmarked against a known peak, plotted, and only then compared with a reference implementation such as cuBLAS or OpenBLAS. The int8 section is the clearest example, where the comparison is drawn against the rk3399 single core fp32 peak rather than against another int8 library.

The licence is GPLv3, which is the practical constraint for anyone who wants to reuse a kernel. The repository is the tutorial, and the int8 library that came out of it is the separately licensed chgemm project. The primary language is C++ with assembly for the inner kernels, the default branch is master, the last push was on 2026-05-14, and there are no GitHub releases, so the tutorial articles on the linked site are the versioned history of the ideas while the repository is the current state of the code.

## Conclusion

Adopt this repository as a learning reference if you want to see how a matrix multiply kernel is written and improved step by step for a specific architecture, since every backend keeps the old and new implementations in one build and plots them together, and the int8 variant has already been merged into ncnn. Keep away if you need a production BLAS, since the code is tutorial grade: the x86 version does not handle edge cases, only sizes that are multiples of four, and the README says the same for the int8 comparison. Before you run the aarch64 path, note that the project expects an ARM Linux machine and that mperf has to be enabled with -DMPERF_ENABLE=ON, and check the last push date of 2026-05-14 given that the repository publishes no releases.

## FAQ

### How do I build and benchmark a kernel in how-to-optimize-gemm?

In the backend directory you want, edit the makefile so the OLD and NEW variables name the same implementation for a first run, then run make run. The makefile compiles and runs the kernel named by NEW and copies the result to output_new.m, and you plot the numbers with python3 plot.py after installing numpy and matplotlib.

### Does the CUDA version of how-to-optimize-gemm beat cuBLAS?

The README states that the cuda version exceeds the speed of NVIDIA cuBLAS, with a plot showing the project's own implementation in green and cuBLAS in blue, both measured on a 3080 without tensor cores. Running the comparison requires the CUDA driver and nvcc, plus a CPU OpenBLAS build to verify numerical correctness.

### What are the limits of the int8 implementation in how-to-optimize-gemm?

The int8 gemm is a separate project, chgemm, which handles edge cases that the tutorial code ignores, reaches 18.6 gflops against a 14.3 fp32 theoretical limit, and uses symmetric quantisation, so inputs must fall inside [-127, +127] and -128 cannot occur. That work has been merged into the ncnn INT8 convolution implementation.

### Which backends does how-to-optimize-gemm support?

The repository has directories for aarch64, aarch64-int8, armv7, cuda, cuda-bf16, cuda-int4, vulkan and x86, and the README notes small differences between them, such as the name chosen for NEW varying and the vulkan and int4 backends needing dependencies installed in advance.

## Sources

- [Issues](https://github.com/tpoisonooo/how-to-optimize-gemm/issues)
- [License: GPL-3.0](https://github.com/tpoisonooo/how-to-optimize-gemm/blob/master/LICENSE)
- [Project website](https://zhuanlan.zhihu.com/p/65436463)
- [README](https://github.com/tpoisonooo/how-to-optimize-gemm/blob/master/README.md)
- [tpoisonooo/how-to-optimize-gemm on GitHub](https://github.com/tpoisonooo/how-to-optimize-gemm)

---

Hysen Labs editorial analysis, written from the project's own repository and release notes. Cite the canonical page: https://hysenlabs.com/projects/tpoisonooo-how-to-optimize-gemm