The GEMM Built Cleanly. The Benchmark Was Still 24× Off.

Silent wrong answers, a broken tanh, and a build directory that broke CI for a reason no error message named.

There’s a gap between knowing that _mm256_fmadd_ps exists and knowing how to use it to build something fast, numerically correct, and honest about what it actually achieves. IntrinsicML was an attempt to close that gap: a hand-vectorized C++ microkernel library implementing the SIMD primitives, GEMM, activations, and reduced-precision kernels , that sit underneath every ML inference stack.

This covers two sessions of work spanning roughly three months. The first produced a working library with AVX2/AVX-512 GEMM, activation kernels, Python bindings, and a benchmark against OpenBLAS. The second fixed a CI pipeline with nine failing checks and implemented the reduced-precision kernels the first session had left as “not done.”

A note on process: I built this with Claude (Anthropic’s AI assistant) as a development partner across both sessions. I’ll be direct about what that looked like, including where Claude was wrong and I caught it. I’m not going to pretend the code emerged from a single human’s head, and I’m not going to oversell what AI-assisted development actually is.

What We Built

IntrinsicML is a C++17 library of hand-vectorized AVX2/AVX-512 microkernels for ML primitives, with Python/NumPy bindings via pybind11. By the end of session two, it includes:

• A GEMM implementation following the Goto/BLIS 5-loop structure with panel packing and register blocking, 8×8 AVX2 and 8×16 AVX-512 micro-kernels, plus an auto-dispatching small-matrix path that skips packing entirely when the working set fits in L2 cache
• FP16 GEMM using F16C conversion instructions (FP16 in, FP32 accumulate, FP32 out) and BF16 GEMM using only AVX2 integer instructions; no special ISA flag required
• GeLU, ReLU, SiLU, Softmax, and LayerNorm, all AVX2-vectorized
• A statistical benchmarking harness with 95% confidence intervals, JSON output, and a CI regression gate
• Python bindings with type stubs and per-call ISA forcing for benchmarking
• Four CI workflows and a working macOS Apple Silicon build path

The performance story: at N=256 and above, the packed GEMM runs at 55–59% of single-threaded OpenBLAS. At N=64, it originally ran at 3.7%, and that gap is the story of the first lesson below.

Lesson 1: Silent Wrong Answers Are the Worst Kind of Bug

The most dangerous bug in the project produced no crash, no error, and no suspicious benchmark number. It just silently computed the wrong thing.

The original gemm_micro_6x32_avx512 stub declared a 32-column output tile but only ever loaded and accumulated one __m512 register (16 floats) per row. The other 16 columns of every tile were never touched. With beta=0 (overwrite mode), those columns were zero from the initial memset , which happened to be the correct answer for the early test sizes, none of which reliably hit the full-block code path.

The discipline that caught it: before wiring any new kernel into the dispatch path, write a standalone unit test that (1) picks a dimension that’s an exact multiple of the tile width, so no edge blocks provide coincidental cover; (2) checks the “lo half” and “hi half” of every output tile independently against a scalar reference; and (3) poisons the initial C values with 1e30f instead of zero, so any unwritten column produces an obviously wrong answer. The original bug would have sailed through a lo-half-only check.

The general version of this: a test that lets an initial value survive unnoticed isn’t testing your kernel; it’s documenting what you hoped would happen. In SIMD code specifically, “the full-block path writes everything it claims to write” is almost never checked explicitly. It should be, with adversarial initial values, every time.

Lesson 2: The Tanh That Wasn’t Tanh

The GeLU implementation we inherited used a rational polynomial approximation of tanh with a maximum absolute error of 1.018 , which is to say, it wasn’t approximating tanh at all. It had the right shape near zero and nothing else.

The error only showed up for |x| > 2, where a meaningful fraction of GeLU's actual behavior lives. The existing tests passed because they compared the fast path against the scalar path, and the scalar path used the same broken polynomial. Two wrong implementations agreeing with each other looked exactly like correctness.

The fix used Cody–Waite range reduction into exp(2y), then tanh(y) = (exp(2y) − 1) / (exp(2y) + 1) , the same approach XNNPACK, Eigen, Intel oneDNN, and Sleef all independently converge on. Max absolute error dropped to under 2×10⁻⁷. The lesson wasn't "derive better coefficients"; it was "check what production libraries do before inventing your own." Cody and Waite published this in 1980; we used it.

The deeper trap here: when your reference implementation and your fast implementation share a bug, your tests will tell you they agree. Testing numerical code means cross-validating against something independent of both; for GeLU, that’s PyTorch’s F.gelu(x, approximate='tanh'), not your own scalar fallback.

Lesson 3: A Passing Tolerance Can Still Be a Broken Metric

A test with tolerance 1e-4 failed in CI with max_rel_err = 1.46e-4, using the naive metric |got − ref| / (|ref| + 1e-7).

The failing input was x = -3.65, where GeLU(-3.65) ≈ -3.4×10⁻⁴ , deep in the asymptotic tail. The absolute error was 1.3×10⁻⁷, comfortably inside the Cody–Waite error budget. But dividing that tiny absolute error by an even tinier reference value inflated the relative error past the threshold. The kernel was correct. The metric was lying.

The fix is what np.allclose does internally: a combined bound, |diff| ≤ atol + rtol × |ref|. When |ref| is near zero, atol dominates and the bound stops exploding. We found the same anti-pattern ten separate times across GEMM, GeLU, SiLU, and LayerNorm tests; anywhere the function legitimately crosses zero, naive relative error is a flakiness bomb waiting for an unlucky random seed.

Lesson 4: Reporting the Number You Measured, Not the One You Expected

We built solid benchmarking infrastructure: RDTSC with LFENCE serialization, wall-clock timing with 95% confidence intervals, CPU pinning, a CI regression gate with committed JSON baselines. But the sharpest lesson wasn’t about the measurement mechanism. It was about what we did with an inconvenient result.

The first draft of the README claimed “~90% of OpenBLAS at N=64.” The actual measurement was 3.7% , a 24× gap between the claim and reality. Not malice, just an ungrounded guess based on the general intuition that “cache-resident GEMM should be fast.”

Source: Image by the author.

The mechanism, once we measured instead of assumed: panel packing costs O(MK + KN) regardless of output size. At N=64, that overhead swamps the FMA work it's supposed to accelerate; the packing buffer is roughly 128× the size of the actual data being multiplied. Doubling SIMD width via AVX-512 barely moved the needle (3.7% → 3.9%), which confirmed the regime wasn't compute-bound at all; it was packing-overhead-bound. The real fix, implemented in session two, was a direct unpacked kernel path for small matrices , which brought N=64 up to a healthy fraction of OpenBLAS without touching the packed path at all.

A performance claim unsupported by measurement isn’t marketing; it’s misinformation. Anyone who built on top of “90% of OpenBLAS” would have shipped something 24× slower than they planned for.

Lesson 5: A Wheel That Installs Cleanly Can Still Be Broken

We added cmake --install rules so downstream C++ consumers could use find_package(simd_kernels REQUIRED). The Python wheel broke as a direct result.

scikit-build-core builds wheels by running cmake --install and packaging whatever lands in the install tree. Before we added install rules, the wheel happened to contain the compiled .so by default. After we added rules for headers, cmake config, and benchmark executables, the install tree contained all of those, but the .so itself had no explicit install() rule, so it silently disappeared from the tree. Every pip install succeeded. Every import simd_kernels failed with ModuleNotFoundError.

The fix needed two things simultaneously: an explicit install(TARGETS ... COMPONENT PythonModule) for the .so, and pyproject.toml telling scikit-build-core to package only that component. "cmake --install exits zero" and "cmake --install produces an importable package" are different guarantees, and only one of them was being tested.

Lesson 6: Committed Build Artifacts Cause Bugs That Don’t Look Like Their Cause

This one produced nine simultaneous CI failures in session two, and none of the error messages pointed at the real problem.

A gitignored build/ directory had, at some point, been force-committed. Inside it, build/_deps/doctest-src , a FetchContent-cloned copy of the doctest repo , looked to GitHub's tooling like a submodule gitlink with no matching .gitmodules entry. actions/checkout --submodules recursive failed immediately with a fatal: No url found for submodule path error, before any build step ran. Separately, build/CMakeCache.txt had an absolute path baked in from the original dev machine; CMake refuses to reuse a cache generated from a different source directory, so every fresh cmake -S . -B build step failed with a path mismatch on the CI runner's different filesystem layout.

Both failures looked unrelated to their actual cause; one read as a source-control problem, the other as a build-configuration problem. The fix was one command: git rm -r --cached build/.

The broader point: a repository should hold inputs to your build, not its outputs. Generated files carry invisible state, paths, timestamps, and checksums , that was true on the machine that produced them and false everywhere else. If CI starts failing in ways that don’t connect to any recent code change, git ls-files --error-unmatch build/ will tell you in one line whether this is your problem.

Lesson 7: An ISA Override Is a Concurrency Design Decision, Not a Config Flag

The roadmap called for an isa= parameter to force sgemm() onto a specific kernel for benchmarking , compare AVX2 and AVX-512 on one machine without recompiling.

The obvious implementation is a global atomic: std::atomic<IsaMode> g_isa_mode. It's wrong. "Global" is the wrong scope for a per-call override; setting it on one thread creates a window where a different thread's unrelated call also picks up the forced ISA.

The correct scope is thread_local, with an RAII guard (save, set, compute, restore-on-destruction) so the override resets even if the call throws. One subtlety sits on top: detected_isa() , the hardware capability check , is computed once via std::call_once and cached, because hardware capability doesn't change at runtime. It intentionally ignores the per-call override, because "what can this CPU do" and "what should this specific call use" are different questions that happen to share a similar name.

Lesson 8: The Cleanest Code Comes From Understanding the Format, Not the Instructions

The two reduced-precision kernels, implemented in session two after being marked “not done” in session one , turned out to be some of the cleanest code in the project, precisely because understanding the bit layout made the implementation nearly obvious.

FP16. IEEE 754 half-precision is 1 sign bit, 5 exponent bits, 10 mantissa bits. Converting 8 FP16 values to FP32 needs the F16C extension (present on every CPU with AVX2), via a single instruction, _mm256_cvtph_ps. Accumulation stays in FP32 throughout , FP16 accumulation loses roughly 1 ULP per step, invisible at K=4, obvious at K=128.

BF16. Brain Float 16 is 1 sign bit, 8 exponent bits, 7 mantissa bits , the same exponent width as FP32, which means identical dynamic range. The key fact: BF16 is literally the top 16 bits of FP32.

FP32 [31:0]:   [S:1][Exp:8][Mant:23]
BF16 [15:0]:   [S:1][Exp:8][Mant:7 ]   ← same sign + exponent, top 7 mantissa bits

Converting BF16 to FP32 is: zero-extend the 16 bits into 32, then shift left 16. Two AVX2 integer instructions for eight values, correct for normals, subnormals, infinity, and NaN , no F16C, no special compiler flag.

Source: Image by the author.

static inline __m256 bf16_to_fp32_x8(const uint16_t* ptr) noexcept {
    __m128i v16      = _mm_loadu_si128(reinterpret_cast<const __m128i*>(ptr));
    __m256i v32      = _mm256_cvtepu16_epi32(v16);    // zero-extend u16 → u32
    __m256i fp32bits = _mm256_slli_epi32(v32, 16);    // shift into FP32 upper half
    return  _mm256_castsi256_ps(fp32bits);             // reinterpret (free)
}

BF16 wasn’t an arbitrary design , Google Brain chose 8 exponent bits specifically to make this conversion trivial. Knowing why a format was designed a certain way is usually a shortcut to implementing it correctly.

Lesson 9: A New Platform Forces Implicit Assumptions Into the Open

The macOS Apple Silicon CI failure looked, in isolation, like a small annoyance: every AVX2 kernel source unconditionally #include <immintrin.h>, an x86-only header, so the arm64 build failed with dozens of "undeclared identifier _mm256_loadu_ps" errors.

The real fix was architectural: a SIMD_ML_ENABLE_X86 CMake option that auto-detects processor and compiler capability, then conditionally compiles kernel sources, flags, and test files, gated behind a SIMD_ML_X86_AVAILABLE preprocessor macro. The macOS CI job explicitly passes -DSIMD_ML_ENABLE_X86=OFF.

This is more code than the original unconditional version. It’s also honest about the codebase’s actual requirements , every x86-specific file now lives in an explicitly named build block, instead of being an implicit assumption that surfaces as a confusing compile error on the one platform where it’s wrong. Adding a new platform is a forcing function: it turns implicit assumptions into documented ones, which is worth doing even for platforms you never ship on.

Working With Claude, With No Memory Between Sessions

Session two started from zero; Claude had no memory of session one’s work. What made it productive was that session one had left behind a well-documented repository: a roadmap that named the small-matrix problem precisely, comments explaining why decisions were made, and CI logs that, read carefully, told most of the diagnostic story without needing to run code against the actual CI environment.

What Claude contributed: diagnosing all nine CI failures from logs alone, implementing the FP16/BF16 kernels and the ARM conditional-compilation structure, and writing the corresponding tests , all auditable line by line in the diffs. What stayed human judgment: deciding the BF16 path shouldn’t attempt an unpacked vdpbf16ps implementation (the interleaved-pair layout requirement makes it a net loss versus the AVX2 zero-extend path), deciding a FlashAttention microkernel was a separate session's worth of scope and shouldn't be started half-finished; and deciding which CI logs were worth sending as context versus which conclusions were already clear from the code.

The division of labor held across both sessions: Claude was strong at implementing well-specified problems and diagnosing issues given complete context. Deciding what to build, in what order, and with which tradeoffs stayed on the human side, and the project only worked because that side was actively reading the code, not approving output.

What Made This Project Worth Doing

The algorithms here aren’t secret; they’re in Goto’s 2008 paper and the Intel intrinsics reference. What’s less common is a library that shows its work: source comments explaining why at the point of each decision; benchmark numbers published, including the embarrassing one; and tests built specifically against the failure modes SIMD code actually has, wrong tile half, wrong edge-block stride, and relative error blowing up near zero.

The 55–59% of OpenBLAS at N≥256 is real. The N=64 gap was real, and now has a working fix. The Cody–Waite tanh matches Intel SVML’s error bound. None of these numbers are in the documentation because they flatter the project; they’re there because a reference implementation that hides its gaps is marketing, and one that documents them is actually useful to the next person who reads the code.

The codebase is at github.com/Mattral/SIMD-Microkernels-for-ML-Workloads. If you find a bug, open an issue. If you benchmark on hardware not represented in the docs, a PR with those numbers is genuinely useful , CONTRIBUTING.md has the checklist.