Eigen vs. BLAS/LAPACK libraries, AMD Zen 5 (October 2026)

Eigen’s dense kernels measured head-to-head against three BLAS/LAPACK libraries on one machine, one thread, the same binary driver and the same flop count for every arm: eight BLAS operations from level 1 to level 3 and nine LAPACK factorizations and solvers. The same cells with Eigen and every library limited to AVX2 are on the AVX2 page.

This page was first published with a measurement of 2026-09-10 to 2026-09-16 at Eigen master 15f227178 on Ubuntu 24.04. It now shows a full re-measurement of Eigen and of every library, taken after the host moved to Ubuntu 26.04 with GCC 15.2, at Eigen master 674e4a4e5, and adds five operations: TRMM, SYMM, PSTRF, SYTRF and GELS. Against September’s last Eigen measurement, on the cells both pages share, Eigen’s AVX-512 rates repeat within 4 % (geometric mean per operation and scalar) for AXPY, DOT, GEMM, GESDD and TRSM, and for SYRK in double; GEMV is 10 to 13 % slower, POTRF 6 to 8 %, SYEV 7 to 8 % and GETRF 4 to 7 % slower; GEEV is 14 to 21 % faster, GEQRF about 8 % faster and SYRK in float 9 % faster. The libraries changed with the host, so the ratios of the two versions of the page differ for that reason as well. The September version of the page and its result files are in the history of the website repository.

Conditions

GEMM at n = 1024 ran at 144 GFLOP/s in double and 300 in float in Eigen, against 122 and 287 for OpenBLAS, 93 and 284 for AOCL, and 68 and 135 for oneMKL.

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. A ratio divides Eigen’s rate by the library’s rate from the same binary and the same run. Each of the three binaries, one per library, also measured every Eigen cell, so Eigen was measured three times, minutes to hours apart: 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. The spread is within 3 % for 76 % of the cells and above 10 % for 5 %, mostly AXPY, DOT and GEMV cells and the smallest sizes. 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.

Each run covered one operation against one library and ran under a fleet-quiet lease that kept other work off the host. It started only after the host passed a gate (no compilers or containers, the guest’s busy-core count below 0.5 and its one-minute load below 0.9, the Windows host’s CPU load at most 8 %) and was discarded and repeated if the guest averaged 1.8 busy cores or more in any 3 s sample, or a container started. The Windows host’s load is invisible inside the guest.

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.60 ns per call for Eigen against 3.9 ns for OpenBLAS, 4.2 ns for oneMKL and 8.7 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: Eigen’s AXPY holds at 3.0 to 4.3 ns per call from n = 2 to 32; a September test on this machine that cycled the same kernel through 64 vectors instead of one measured 0.8 to 1.5 ns.

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.

TRMM

X.noalias() = A.triangularView<Lower>() * B against ?trmm with the m-by-m lower triangle on the left and n columns in B, no transposition, non-unit diagonal; (m(m+1)n) flops, TRSM’s count. ?trmm overwrites B, so the reference arm copies B into the destination and multiplies in place; the Eigen arm writes the product to the destination, which noalias() zero-fills before accumulating. Each arm makes the one pass over the destination its interface requires, and neither allocates inside the timed region.

The same shapes as TRSM: square, then the sweeps with a fixed triangle of order 1024 or 4096 and a growing number of columns, then the remaining shapes with 4096 columns against a small triangle.

SYMM

C.noalias() += A.selfadjointView<Lower>() * B against ?symm (?hemm for a complex scalar) with the m-by-m symmetric matrix stored in its lower triangle, on the left, and n columns in B; (2m^2n) flops, the count of the general product of the same shape. C accumulates across iterations and is never reset, as in GEMM.

The same shapes as TRSM.

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.

PSTRF

LDLT<Ref<MatrixXd>> in place against ?pstrf, lower triangle, symmetric positive definite operand; (n(n-1)(n-2)/3 + 2n(n-1)) flops for both, (n^3/3) to leading order. The two compute different factors of the same pivoted matrix, (P^T A P = L L^T) for LAPACK and (L D L^T) for Eigen, and pivot on the same entry, the largest diagonal element of the remaining Schur complement. Both arms copy the operand into a destination allocated once and factorize the copy, as GETRF does. LDLT::compute also evaluates the operand’s L1 norm for rcond(), which ?pstrf does not. Each result is checked by applying the stored factors to a random vector and comparing with (Ax).

SYTRF

BunchKaufman<Ref<MatrixXd>> in place against ?sytrf (?hetrf for a complex scalar), lower triangle, symmetric indefinite operand; (n(n-1)(n-2)/3 + 2n(n-1)) flops for both, (n^3/3) to leading order. Both arms copy the operand into a destination allocated once and factorize the copy, as GETRF does; pivot and workspace arrays are set up once per shape. Each result is checked by applying the stored factors, with the pivot blocks of D, to a random vector and comparing with (Ax).

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.

GELS

Full-rank least squares (\min \lVert AX - B \rVert) through Householder QR with A m-by-n, m ≥ n, and B m-by-nrhs: HouseholderQR<Ref<MatrixXd>> followed by qr.solve(B) against ?gels (trans = 'N'); (2mn^2 - 2n^3/3) flops for the factorization plus (\text{nrhs},(4mn - n^2)) for applying (Q^T) and solving with R. ?gels overwrites A and B, so the reference arm copies both into destinations allocated once. The Eigen arm factorizes a copy of A through the in-place HouseholderQR and writes qr.solve(B) to a preallocated X; solve() copies B into a temporary it allocates, which is the cost of the expression a user writes, so it is timed. The solution is checked through the normal equations on a well-conditioned operand.

Square shapes with one right-hand side, then the tall-skinny sweeps with m fixed at 4096, 10000 and 100000, then right-hand-side sweeps: 16 to 256 columns of B at m = n = 1024 and 16 to 1024 at 4096 by 256.

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.1 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.

Observations

• Where Eigen leads. In the middle of the size range (n = 24 to 1024) Eigen is ahead of all three libraries in nearly every BLAS operation and in the Cholesky-, LU- and QR-based factorizations; the exceptions in double are GEMM against OpenBLAS and AOCL at n = 24 to 128 (0.90 and 0.89), TRSM against AOCL (0.83 to 0.96) and GESDD (0.75 to 0.98). oneMKL runs its generic code path on this CPU: Eigen is 2.0 to 2.9 times ahead of it on GEMM, TRSM, SYMM and SYRK at n = 24 to 1024.
• The new factorizations. SYTRF (Bunch-Kaufman) is 1.1 to 2.5 times faster in double than OpenBLAS, oneMKL and AOCL at every size band above n = 16; PSTRF (pivoted Cholesky against Eigen’s LDLT) 1.1 to 2.3 times, and GELS 0.9 to 2.8 times.
• TRMM is the exception. For triangles of order 2048 and up Eigen’s rate falls to 62, 43 and 37 GFLOP/s in double at m = 2048, 3072 and 4096, while OpenBLAS holds 113 to 117: ratios of 0.54, 0.38 and 0.32, and 0.29 to 0.48 against AOCL. The cliff sits between m = 1536 (108 GFLOP/s, 0.95) and 2048, and shows in the sweeps too (0.54 at 1024 by 4096 and 4096 by 1024). For triangles of order 2 to 12 Eigen is also 0.3 to 0.5 times OpenBLAS. In between it is 1.0 to 1.2 times OpenBLAS.
• GEEV trails all three libraries at n = 192 to 1024 (0.58 to 0.81 in double), as it did in September; it leads at the small end.

Caveats

• Linux runs under WSL2, so absolute rates are a lower bound on bare-metal Linux.
• oneMKL as shipped does not use its AVX-512 kernels on AMD CPUs.
• Ubuntu 26.04’s GCC 15.2 targets x86-64-v3 by default, so even the Eigen build for the AVX2 page has that baseline; the -mavx512dq -mfma of this page’s build adds to it.
• The AOCL LAPACK arm is AOCL-libFLAME with the AOCL-BLIS it links; the OpenBLAS LAPACK arm is the reference LAPACK that OpenBLAS bundles, on top of OpenBLAS’s BLAS kernels. LAPACK ratios therefore compare Eigen’s factorizations with different libraries’ blocking as well as their kernels.
• PSTRF compares two algorithms that pivot on the same entry but produce different factors, (L L^T) and (L D L^T), and charges both the same flop count.

Full tables

GFLOP/s per shape and library, with the Eigen ÷ library ratio. The Eigen column is the median over its three passes; each ratio is formed within the library’s own binary and run.