Eigen vs. BLAS/LAPACK libraries, Apple M4 with NEON (October 2026)

The cells of the Apple M4 SME page measured again on the same machine with Eigen compiled for NEON alone, the default for this platform, against OpenBLAS: what the comparison looks like without the SME unit.

Apple Accelerate and Arm Performance Libraries are not on this page: neither can be held to NEON through documented means. For the comparison with them, see the SME page.

This page was first published with a measurement of 2026-09-20 at Eigen master 46a2d0c56 on macOS 27.0. It now shows Eigen measured again at master 3cf994918 on macOS 27.0.1, the tree of the SME page’s October measurement. Against September, Eigen’s NEON rates repeat within 3 % (geometric mean per operation and scalar) except TRSM, 12 % faster, and SYEV in double and GEEV, 7 to 9 % slower. The GEMV ratios in double above n = 1000 fell because OpenBLAS’s dgemv measured faster in the October session than in September’s (2.6× at n = 1024). 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 60 GFLOP/s in double and 121 in float in Eigen, against 57 and 114 for OpenBLAS.

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. Eigen was measured three times, in the three binaries of the same build, those linked against OpenBLAS, Accelerate and Arm Performance Libraries, where only the Eigen cells ran. 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 95 % of the cells and above 10 % for 0.2 %. The widest spread depends on the binary rather than on noise: GETRF in float at n = 100 runs at 23 GFLOP/s in the binaries linked against OpenBLAS and Accelerate and at 36.5 in the third, while the ten repetitions within each binary span at most 4 % and the neighbouring sizes n = 96 and 128 agree across all three.

A ratio divides that median by OpenBLAS’s rate from the SME page’s session, 22 to 33 hours earlier for any one operation, since the library does not change with Eigen’s build. OpenBLAS’s own rates bound the drift between sessions on this machine: between the SME page’s September and October sessions, 13 days apart, they repeated within a geometric mean of 0.97 to 1.03 for every operation and scalar except dgemv (1.11).

Each run covered one operation in one binary. A run started only once the machine had been idle (other processes using less than 0.35 of a core over two 5 s windows, no keyboard or mouse input for 10 minutes) with no thermal or performance warning from pmset. It was first discarded and repeated if the whole machine averaged 1.35 busy cores or more during it, or reached 2.0 in any 5 s window, the limits of the SME page. Background services of the logged-in session (malware scans, search and photo indexing, mail) kept the machine near those limits through the night: 11 runs passed, and 14 failed four times each, although their discarded attempts agreed with the passed runs within 1 % (geometric mean over the cells of a run). The limits were then raised to 1.45 and 3.0. The fourth attempt of each of the 14 runs meets them (averages 1.27 to 1.40, peaks at most 2.76) and was kept; each agrees with the other two passes of its cells within 0.5 %. The remaining 11 runs were measured under the raised limits, 1 of them twice. Spotlight indexing was switched off for the session.

The dotted vertical lines on the charts mark the sizes at which the operands of the call, inputs and outputs together, fill a performance core’s L1 data cache (128 KiB) and the L2 cache that the four performance cores share (16 MiB); the M4 has no L3.

For the eigenvalue and singular value decompositions the operation 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.

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, so the ratios at the small end are ratios of fixed costs. The DOT cells in double at n = 2 show it directly: 0.49 ns per call for Eigen against 12.4 ns for OpenBLAS. Eigen still computes the whole operation on every call; nothing is hoisted out of the timed loop.

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, where a call of SelfAdjointEigenSolver takes 3.0 s in double.

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 3.7 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

• Accelerate and Arm Performance Libraries are left out because neither can be held to NEON through documented means; the SME page has the comparison with them.
• The OpenBLAS ratios are formed across two sessions on the same machine, 22 to 33 hours apart for any one operation; the drift bound above applies to them.
• macOS cannot bind a thread to a core, so a single-threaded rate is that of whichever core the scheduler chose; a busy foreground thread normally stays on a performance core.
• The machine is a desktop with a logged-in session and its background services. Runs that they disturbed beyond the limits above were repeated; smaller disturbances remain in the numbers. This session ran busier than the SME page’s (1.27 to 1.40 busy cores averaged over a kept run, against 1.10 to 1.13), and Eigen’s AXPY, DOT and GEMV came out 2 to 3 % below their September rates; whether that is the load or the code is not separated here.

Full tables

GFLOP/s per shape and library, with the Eigen ÷ library ratio: Eigen’s median over its three passes against the library’s rate.