Eigen vs. BLAS/LAPACK libraries, AMD Zen 5 limited to AVX2 (September 2026)

The cells of the AVX-512 page measured again with Eigen compiled for AVX2 and FMA and every library held to its AVX2 kernels, on the same machine: what the comparison looks like on a CPU without AVX-512.

Conditions

The limits took effect: GEMM in double at n = 1024 ran at 121 and 105 GFLOP/s in OpenBLAS and AOCL without them and at 75 and 75 GFLOP/s with them, against 74 GFLOP/s for Eigen; oneMKL ran at 68 GFLOP/s either way.

How to read the numbers. Every rate is nominal GFLOP/s: the textbook operation count of the operation, stated in each section below, divided by the median wall time of one call. Both arms of a comparison are built into the same binary and run on the same freshly randomized operands, and the reference library’s result is validated against Eigen’s before it is timed. The library rates are from the session of 2026-09-14 and 2026-09-15. Eigen was re-measured on 2026-09-16, after changes to its Schur decomposition and divide-and-conquer SVD reached master: the same three binaries, rebuilt at that commit, ran their Eigen cells only, one pass each. The Eigen line on each chart is the median of the three passes and the shaded band is their spread, the run-to-run noise on this host: within 3 % for 83 % of the cells and above 10 % for 5 %. A ratio divides that median by the library’s rate from the earlier session. Eigen’s cells in that session bound the drift between sessions: every operation but GEEV and GESDD repeats its rates within 1 % at the median cell, so ratio differences of a few percent are not significant. A few level-3 cells run at one of two rates, up to 2.5 times apart, depending on the binary rather than the session: SYRK and GEMM at n = 192, and in float at n = 256. Each rate is stable within a binary; the earlier session’s three binaries were split between the two, and the three rebuilt ones all run at the same one. This matches the per-call allocation of Eigen’s packing buffers described in libeigen/eigen#3201, measured there with the same AVX2 build on a Zen 4 core: at these sizes some processes page-fault the buffers in on every call, depending on their earlier allocations, and run at the lower rate. For the eigenvalue and singular value decompositions the count is the one of the classical algorithm, so a divide-and-conquer or blocked implementation is credited with work it does not do; the ratios between arms are unaffected.

The dotted vertical lines on the charts mark the sizes at which the operands of the call, inputs and outputs together, fill the core’s L1 data cache (48 KiB), its L2 cache (1 MiB) and the L3 cache as the guest sees it (32 MiB).

AXPY

y += alpha * x against ?axpy, alpha = 0.75; (2n) flops. Lengths run from 2 to (2^{24}) (128 MiB per double vector), through every cache level into DRAM.

Short vectors. Up to about n = 32 the AXPY and DOT rates measure the cost of a call rather than the kernel. Eigen’s expression compiles inline into the calling loop, while a call into a shared library pays a cost of its own before it reads an element. That cost is fixed per call while the cost per element is similar in every arm, so the ratios at the small end are ratios of fixed costs. The DOT cells in double at n = 2 show it directly: 0.50 ns per call for Eigen against 3.6 ns for OpenBLAS, 4.0 ns for oneMKL and 8.3 ns for AOCL. Eigen still computes the whole operation on every call; nothing is hoisted out of the timed loop.

AXPY’s lead is smaller because every call adds into the y written by the call before it, and reading that memory back puts a floor under every arm: up to n = 32, Eigen’s AXPY holds at 2.6 to 3.0 ns per call where n is a multiple of the vector width, 4 in double and 8 in float, and at 4.0 ns where it is not and the remaining elements are written and reloaded through masked moves, while the same kernel cycling through 64 vectors instead of one took 0.8 to 1.4 ns in a separate test on this machine.

DOT

x.dot(y) against ?dot; (2n) flops, the same lengths. Real scalars only: the Fortran convention for returning a complex function value is not fixed by the platform ABI and differs between library builds, so ?dotc cannot be declared portably.

Up to about n = 32 the rates measure the cost of a call; see AXPY.

GEMV

y.noalias() += A * x against ?gemv, column-major, no transposition; (2mn) flops.

Square matrices, then the short-fat (n = 10000, m grows) and tall-skinny (m = 10000, n grows) sweeps.

GEMM

C.noalias() += A * B against ?gemm, column-major, no transposition; (2mnk) flops.

Square matrices first; then the two tall-skinny sweeps, where m stays at 4096 or 10000 while n = k grows; then the assorted non-square shapes that fit neither.

TRSM

A.triangularView<Lower>().solveInPlace(X) against ?trsm with the m-by-m lower triangle on the left and n right-hand sides, no transposition, non-unit diagonal; (m(m+1)n) flops. The solve is destructive, so the timed region includes one copy of B for both arms.

Square shapes (m = n), then the sweeps with a fixed triangle of order 1024 or 4096 and a growing number of right-hand sides, then the remaining shapes with 4096 right-hand sides against a small triangle.

SYRK

C.selfadjointView<Lower>().rankUpdate(A) against ?syrk, lower triangle, no transposition, with C n-by-n and A n-by-k; (n(n+1)k) flops.

Square shapes (n = k), then the sweeps: a 4096-by-4096 matrix updated by a growing number of columns, and the Gram matrix of a short, wide factor with 4096 columns.

POTRF

LLT<Ref<MatrixXd>> in place against ?potrf, lower triangle, symmetric positive definite operand; (n^3/3) flops. The timed region includes one copy of the operand for both arms, since the factorization is destructive.

GETRF

PartialPivLU<Ref<MatrixXd>> in place against ?getrf, square operand; (2n^3/3) flops. Both arms copy the operand into a destination allocated once and factorize the copy; neither allocates inside the timed region.

GEQRF

HouseholderQR<Ref<MatrixXd>> in place against ?geqrf, m ≥ n; (2mn^2 - 2n^3/3) flops, the same copy-then-factorize protocol as GETRF.

Square shapes, then the tall-skinny sweeps with m fixed at 4096, 10000 and 100000.

SYEV

SelfAdjointEigenSolver with eigenvectors against ?syev (jobz = 'V', lower triangle), the like-for-like pair: tridiagonalization followed by implicit QR iteration on both sides; (9n^3) flops, the count of that algorithm (Golub & Van Loan, section 8.3.3). Neither works in place: the solver copies the operand internally and the reference arm pays the same one copy. Sizes to 2048.

GESDD

BDCSVD<MatrixXd, ComputeThinU | ComputeThinV> against ?gesdd with jobz = 'S', divide and conquer on both sides, m ≥ n, the thin U and V computed; (14mn^2 + 8n^3) flops, the Golub-Reinsch count (Golub & Van Loan, table 8.6.1). Square sizes to 2048, then tall-skinny sweeps with m fixed at 4096 and 10000.

GEEV

EigenSolver with eigenvectors against ?geev (jobvr = 'V'), a dense random nonsymmetric operand: Hessenberg reduction, QR iteration to the Schur form and back-substitution for the eigenvectors on both sides; (25n^3) flops, the count of the Hessenberg QR algorithm with Schur vectors (Golub & Van Loan, section 7.5.6). Sizes to 1024, where a call of EigenSolver takes 2.9 s in double.

Summary

Geometric mean of the Eigen ÷ library rate ratio over the square shapes in each size band. Above 1, Eigen is faster. Only the vector operations reach the last band.

Caveats

• Linux runs under WSL2, so absolute rates are a lower bound on bare-metal Linux.
• The limits emulate an AVX2 CPU on Zen 5 silicon: the caches, memory bandwidth and branch predictor are those of this machine, not of an AVX2-only part.
• oneMKL as shipped does not use its AVX2 or AVX-512 kernels on AMD CPUs.

Full tables

GFLOP/s per shape and library, with the Eigen ÷ library ratio: Eigen’s median over its three passes of 2026-09-16 against the library’s rate from the earlier session.