Fast multiplication of large matrices

Programs for multiplying large square matrices in single precision (FP32) on multicore CPUs, in accurate, safe, fast and large-matrix variants.

The ideas of what to implement are mine; the code was written by AI. There is no new theory here: everything used is publicly known.

Programs for Linux x86-64 (tar.gz, 891 KB), version of 12 October 2026 (adds the large mode)
SHA-256 95b075591559bee023da924ef4a1f0f47e3b941d89f4de5876267a81b535975f
Previous version (28 September 2026, four modes): fastmm-linux-x86_64_2026-09-28.tar.gz, SHA-256 0187b3e104861359732f24f25fc428f6d3d4be94dfc6a249146b00e7f6890a94
Licence: PolyForm Strict 1.0.0. Free to download and use for noncommercial purposes; no changes, redistribution or sale.

The variants

The program bin/fastmm_multiply has five modes. Times for two 30 720 × 30 720 matrices with entries uniform in [−1, 1] on one AMD EPYC 7713 (64 cores); error is the relative error against an FP64 reference, as a multiple of the error of the vendor library (tuned AOCL-BLAS 5.3.2).

mode	use	time	speed-up	error
`accurate`	error at or below the vendor library’s	14.6 s	1.09×	0.97×
`safe` (default)	faster; robust on inputs with large means	11.6 s	1.36×	2.3×
`fast`	fastest of the original modes; ordinary data only (see below)	11.1 s	1.43×	4.1×
`large` (new)	for big matrices (fastest from N = 61 440); robust like `safe`	11.5 s	1.39×	5.1×
`classical`	the vendor library, for comparison	15.9 s	1	1

The gain grows with the size. The same programs on larger matrices of the same kind, on the same machine:

N = 61 440

mode	time	speed-up	error	peak memory
`accurate`	100.5 s	1.29×	0.87×	63.6 GB
`safe`	83.4 s	1.56×	2.1×	47.3 GB
`fast`	80.1 s	1.62×	3.8×	61.4 GB
`large` (new)	73.3 s	1.77×	5.6×	46 GB
`classical` (vendor library)	129.9 s	1	1	42.4 GB

N = 100 000

mode	time	speed-up	error	peak memory
`accurate`	418.6 s	1.34×	0.74×	164 GB
`safe`	335.7 s	1.66×	1.9×	123 GB
`fast`	326.2 s	1.71×	3.3×	161 GB
`large` (new)	269 s	2.08×	5.9×	118 GB
`classical` (vendor library)	558.8 s	1	1	112 GB

N = 122 880

mode	time	speed-up	error	peak memory
`accurate`	743.8 s	1.44×	0.71×	245 GB
`safe`	628.8 s	1.70×	1.8×	185 GB
`fast`	611.0 s	1.75×	3.2×	242 GB
`large` (new)	513 s	2.09×	6.2×	176 GB
`classical` (vendor library)	1070 s	1	1	169 GB

Measured on 28 September 2026 with the programs offered here; the large rows on 12 October 2026 with the program offered here, on the same machine with the same inputs and protocol, against the 28 September times of the vendor library (the safe and fast times reproduced within 2 % on the same day). Below about N = 16 000 there is no gain.

large is the mode to use for big matrices when its error level (about five to six times the vendor library’s on ordinary data, 1.5× that of fast) is acceptable: at N = 30 720 it lies between safe and fast, from N = 61 440 on it is the fastest and its gain grows with the size, it needs the least working memory, and it is protected against inputs with large means and badly scaled rows and columns in the same way as safe. It needs at least 30 blocks per side; with the default block size, matrices below N = 30 720 are multiplied in safe mode (the program says so). Its results are not bitwise reproducible from run to run (the parallel order of additions varies); the error level is.

fast is not protected against cancellation. On inputs with large means and heavy cancellation (entries of A near 1, rows of B alternating near ±1) its relative error reached 1.4 %, against 5·10⁻⁵ for the vendor library and 3·10⁻⁸ for safe and accurate. Use safe or accurate for such data.

Requirements

Linux on x86-64 with AVX2 and FMA (AMD Zen 2 or newer, Intel Haswell or newer).
glibc 2.34 or newer (for example Ubuntu 22.04, Debian 12, RHEL 9 or newer).
OpenBLAS: the shared library libopenblas.so.0 from your distribution’s OpenBLAS package.
AMD’s AOCL-BLAS (the multithreaded library libblis-mt.so), loaded at run time.
numactl (recommended); python3 for the self-test.

Install and test

tar xzf fastmm-linux-x86_64.tar.gz
cd fastmm
export FASTMM_BLIS_LIB=/full/path/to/libblis-mt.so    # your AOCL-BLAS library
./SELFTEST.sh                                         # 10 correctness checks, a few minutes

Multiplying your own matrices

numactl --interleave=all bin/fastmm_multiply A.npy B.npy C.npy --mode large --check

Input: two square float32 matrices of the same size, each a NumPy .npy file (float32, C order) or a raw file of little-endian float32 values in row-major order.
Output: the product, in the format given by the extension of the third file name, and a report with the multiplication time, the peak memory and, with --check, the relative error against an FP64 reference on sampled rows and entries.
Options: --mode accurate|safe|fast|large|classical; --workers W (default 64, at most 64; set it to the number of physical cores on smaller machines); --leaf L (block size, default 1024).
Any size is accepted; sides are zero-padded internally, and below N = 8192 the classical route is used.
The report printed after the multiplication gives the multiplication time separately from reading and writing, the peak memory, the time of each internal phase, and, with --check, the error against an FP64 reference.
Memory: the three matrices plus up to 60 % working space.
Integer-valued inputs give exact results while all sums stay below 2²⁴.
On a machine with two sockets or more than 64 cores, confine the run to one socket (numactl --cpunodebind=0 --interleave=0).

Also in the package

bin/fastmm_multiply_modp: exact product of two integer matrices modulo a prime (a modest gain, 1.02–1.11×).

The README.md in the archive has the full options and more measurements. Tested on one AMD EPYC 7713 (Zen 3) under openSUSE Leap 16. The programs are built for AVX2; on processors where the vendor library uses AVX-512 the gains will be smaller. There is no GPU version.