4 ms·
Fast Multidimensional Matrix Multiplication on CPU from Scratch (2022)
- ap4 2y agoRelated: 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 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
- ap4 2y agoTested 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.
- ladberg 2y agoYour `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.
- ap4 2y agoThat'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().
- ladberg 2y agoSorry I'm a bit lost here, could you explain the reasoning behind this and why it works?
- ap4 2y agoI 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.
- deleted 2y ago[deleted]
- ladberg 2y agoI 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 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.
- hedgehog 2y agoFor those interested in going deeper I think the classic reference in this area is the GotoBLAS paper: https://www.cs.utexas.edu/~pingali/CS378/2008sp/papers/gotoPaper.pdf https://www.cs.utexas.edu/~pingali/CS378/2008sp/papers/gotoP...
- gnufx 2y agohttps://github.com/flame/blis#citations https://github.com/flame/blis#citations is the current set of references for various aspects of the approach. The article, or something similar, was posted recently with discussion.
- Remnant44 2y agoI honestly didn't realize how performant the decades-old 2013 Haswell architecture is on vector workloads. 250GFLOP/core is no joke - He also cross-compared to an M1 Pro, that when not using the secret matrix coprocessor achieves effectively the same vector throughput, a decade later...
- Cold_Miserable 2y ago250 is nonsense. 2xFMA per cycle @ ~4.5Ghz = 32*4.5 = ~144 Gflops Beating cuBlas is unlikely. You probably made a mistake. Last I tested it, it was even better than MKL in efficiency.
- Remnant44 2y agoYes, the figure was 250GFLOP for 4 cores instead of per core, I misread. Still impressive but more reasonable
- adrian_b 2y agoThe floating-point FMA throughput per desktop CPU socket and per clock cycle has been doubled every few years in the sequence: Athlon 64 (2003) => Athlon 64 X2 (2005) => Core 2 (2006) => Nehalem (2008) => Sandy Bridge (2011) => Haswell (2013) => Coffee Lake Refresh (2018) => Ryzen 9 3950X (2019) => Ryzen 9 9950X (2024), going from 1 FP64 FMA/socket/clock cycle until 256 FP64 FMA/socket/clock cycle, with double numbers for FP32 FMA (1 FMA/s is counted as 2 Flop/s).
- kvemkon 2y agoI'd wish memory bandwidth could also be doubled so often on desktops. Instead of 256 (even more due to 2-3 times higher core frequency) only 14 times increase: DDR-400 6.4 GB/s => DDR5-5600 89.6 GB/s. The machine balance keeps falling even further. While flash memory became so fast in recent years, I haven't heard of any break-through technology prototypes to bring some significant progress into RAM. Let alone the RAM latency, which remained constant (+/- few ns) through all the years.
- kiririn 2y agoi7 6700 is skylake not haswell
- adrian_b 2y agoTrue, but it has the same floating-point throughput per socket and per clock cycle as Haswell. Until Haswell, for many years Intel had the target of doubling the FP throughput per desktop socket every 2 years, in order to avoid the risk of AMD catching up again with them. After launching Haswell, Intel has relaxed and they did not increase again the throughput for about a half of decade, until they had to start again to do that, under the menace of AMD Zen, when for lack of other better means they have increased every year the number of cores per socket from Kaby Lake (4 cores in 2016) to Comet Lake (10 cores in 2019).
- canjobear 2y agoQuasi-related: Do BLAS libraries ever actually implement Strassen's Algorithm?
- gyrovagueGeist 2y agoTypically no, although BLAS software engineers occasionally write HPC strassen's implementations and papers about them. In my opinion there's a few reasons why they're not common in BLAS libraries: - They're a bit less numerically stable, which used to matter more for common BLAS use cases. Less so these days. - The memory access patterns and algorithm parallelism make it much harder to reach as high a fraction of peak performance as standard GEMM. Matrix size restrictions for recursive mult algorithms is also an issue.
- gnufx 2y agoOne such paper is https://jianyuhuang.com/papers/sc16.pdf https://jianyuhuang.com/papers/sc16.pdf; I don't know if there's anything more recent/relevant.
- deleted 2y ago[deleted]
- deleted 2y ago[deleted]
- namibj 2y agoThe quoted assembly looks L1D$ bandwidth bound; on most common and vaguely recent architectures one has to use register tiling to saturate the FMA units, as a system unable to do more than one vector load and one vector store each per cycle can't ever fully saturate a single FMA unit on GEMM; for 2 FMA units even 2 vector loads and a vector store per cycle won't be enough without register tiling.