feat(backend-native-cpu): native FP16 matmul kernel by michalharakal · Pull Request #896 · SKaiNET-developers/SKaiNET

michalharakal · 2026-07-29T06:36:29Z

Fixes #887.

The issue's root cause is wrong

#887 concludes "the difference is the per-element decode, not the kernel structure or the dispatch." It is the dispatch.

Timing the two Panama kernels directly, bypassing all dispatch (ffn_up 8B, 4096x11008):

	m=1	m=16
`PanamaVectorBf16MatmulKernel`	143 ms	1796 ms
`PanamaVectorFp16MatmulKernel`	164 ms	2062 ms

Within ~15% of each other. But measured through the real dispatch, FP16 came out at 2060 ms — matching its Panama kernel exactly — while BF16 came out at 182 ms, 10x faster than its own Panama kernel. BF16 was never running the JVM kernel.

NativeKernelProvider (priority 100, FFM) implements matmulBf16() but had no matmulFp16() override, so it returned the interface default null and FP16 cascaded to PanamaVectorKernelProvider at priority 50. One narrow-float format was served by a native kernel and the other was not. That asymmetry is the entire 2-18x.

The fix

skainet_fp16_matmul takes the same caller contract and strides as skainet_bf16_matmul, differing in two places.

The dequant. BF16's conversion is free — it is the high half of an FP32. Binary16 needs exponent rebiasing and gradual underflow, so the conversion folds both special cases in with arithmetic masks and stays branch-free, keeping the inner loop a straight-line sequence the vectorizer can widen.

Deliberately no _Float16 and no F16C intrinsics: the x86_64 build carries no -march flag, so F16C cannot be assumed, and _Float16 without it lowers to libgcc helper calls that are slower than the bit math and block vectorization outright. AArch64 does build with +fp16, but a second path would double the surface to test for a handful of integer ops. Runtime ISA dispatch is where that belongs if it ever pays for itself.

The iteration order, which depends on m. i-p-j walks the whole of B once per row of A, re-decoding every element each time. BF16 can afford that at one shift per element; this kernel cannot — with i-p-j it was still 1.55-2.06x slower than FP32 at m=16 even running natively. So at m > 1 j is tiled and each B row is decoded once per tile into a 512-float stack buffer, then multiplied into all m rows of C: decodes drop from m*k*n to k*n and B is read once in total instead of m times, with no allocation entering the kernel. At m == 1 there is nothing to amortize and tiling only trades sequential streaming for a column-block walk, so the straight pass is kept. Accumulation into any given C element is p ascending on both paths, so both are bit-identical to the original i-p-j formulation, not merely close.

Two smaller pieces ride along:

The JVM PanamaVectorFp16MatmulKernel now widens through Float.float16ToFloat (a JDK 20+ HotSpot intrinsic) and does the widening a vector at a time instead of a lane at a time. This does not affect the numbers below — everything there goes native — but it is the path for consumers without the FFM library: non-JVM targets, sandboxes, JDKs without FFM.
Fp16Codec.decode loses its subnormal renormalization loop. Binary16 subnormals are mant * 2^-24 with mant in [1, 1023], both factors exact in FP32 and the product a normal FP32, so one multiply lets the hardware renormalize. This is what ScalarFp16MatmulKernel uses on targets with no intrinsic.
KernelProvider.supports gains its "Float16" arm, which was absent while every other matmul dtype was present — the same omission as the missing accessor, one layer up.

Measured

i7-9750H (AVX2), OpenJDK 21.0.11, median ms per call, NarrowFloatMatmulBenchmark in SKaiNET-transformers:

shape          batch     fp32   fp16 before   fp16 after
q_proj  1B         1     3.16         18.42         3.01
q_proj  1B        16    16.15        298.91        10.59
q_proj  8B         1    40.10         73.82        26.28
q_proj  8B        16    99.97       1182.08        63.76
ffn_up  8B         1   108.53        199.49        71.14
ffn_up  8B        16   270.27       3200.78       161.18
ffn_down 8B        1   106.48        199.13        70.25
ffn_down 8B       16   254.50       3191.62       158.59

FP16 is now 1.5-1.7x faster than FP32 everywhere except 2048x2048 at batch 1, where the two are a wash. Against the old numbers that is 18-28x at batch 16.

The last commit is why the batch-1 numbers are what they are. Tiling j is what buys the batch-16 win, but at m == 1 there is nothing to amortize and it cost 15% (ffn_up 82.1 vs 71.1 ms), so the kernel branches on m and keeps the straight pass for a single row.

The same amortization applies to BF16 and is not in this PR — it lands separately in #897, so BF16's measured behaviour changes with its own numbers rather than silently inside an FP16 fix.

Tests

Native parity covers every shape the BF16 parity test covers, plus two that random weights never reach: a weight set of subnormals, zeros and the format extremes, and an exhaustive sweep multiplying all 65536 bit patterns by 1.0 into a zeroed accumulator, which pins the C conversion against Fp16Codec over the whole domain rather than trusting sampled shapes to have hit a subnormal.
A provider test asserts NativeKernelProvider.matmulFp16() hands out the native kernel — the exact regression this PR exists to prevent, since the failure mode was a missing kernel that looked like a slow one.
Codec equivalence is asserted two ways: in commonTest the pre-change loop-based decode is kept verbatim as the oracle and must match on every non-NaN pattern, so it runs on every target; in a JVM-only test the codec must be bit-identical to Float.float16ToFloat on all 65536 inputs.
The Panama kernel gets the same exhaustive sweep, chunked at 999 columns so each chunk exercises both the vectorized body and the scalar tail.

Green on jvm and linuxX64; apiCheck clean.

NaN behaviour (deliberate change)

The codec and the hardware conversion disagreed on exactly 1022 inputs — every signaling NaN. Zero non-NaN differences. The hardware quiets sNaN; the old Fp16Codec.decode reproduced the signaling bit verbatim.

Fp16Codec.decode now quiets too. Fp16Codec.encode already forced NaN quiet, so decode preserving a signaling NaN was internally inconsistent; leaving it would have made the JVM and every other target disagree on those patterns; and it loses nothing real, since the first arithmetic use quiets the value anyway. The quiet bit is applied only for a non-zero mantissa, or ±Inf would decode as NaN.

Inside the kernels a signaling NaN stays signaling — the value is immediately multiplied and the FMA quiets it, so it is not observable. The kernel sweeps assert bit equality on non-NaN patterns and NaN-ness on the rest, which is exactly that contract.

Platform caveat

I could only build and run the native library for linux-x86_64. The new source compiles wherever the existing kernels do — it is plain C11 with no intrinsics — but the AArch64 NEON build and the MSVC build are unexercised locally and want CI confirmation.

github-actions · 2026-07-29T06:39:34Z

📖 Documentation Preview

The documentation has been built successfully for this PR.

Generated Files:

Operator documentation: docs/modules/operators/_generated_/
JSON schema output: operators.json

Artifacts:

Download the documentation-preview-896 artifact to view the complete documentation locally.

This comment will be updated automatically when the PR is updated.

github-actions · 2026-07-29T07:05:21Z

📖 Documentation Preview

The documentation has been built successfully for this PR.

Generated Files:

Operator documentation: docs/modules/operators/_generated_/
JSON schema output: operators.json

Artifacts:

Download the documentation-preview-896 artifact to view the complete documentation locally.

This comment will be updated automatically when the PR is updated.

github-actions · 2026-07-29T13:45:14Z

📖 Documentation Preview

The documentation has been built successfully for this PR.

Generated Files:

Operator documentation: docs/modules/operators/_generated_/
JSON schema output: operators.json

Artifacts:

Download the documentation-preview-896 artifact to view the complete documentation locally.

This comment will be updated automatically when the PR is updated.

PanamaVectorFp16MatmulKernel ran at a flat ~0.5 GFLOP/s regardless of shape or batch, making FP16 KEEP_NATIVE 2-18x slower than the FP32 SGEMM it replaces, while the structurally identical BF16 kernel was 1.5-2.1x faster. Both fill a scratch lane array scalar-wise before the vector FMA; the difference is what that fill costs. BF16 widens with three integer ops, FP16 called Fp16Codec.decode, whose subnormal arm renormalizes in a data-dependent do/while. Flat throughput while work scales 16x with batch is the signature of a scalar operation dominating the inner loop. Call Float.float16ToFloat from the kernel. It is a HotSpot intrinsic since JDK 20 that lowers to vcvtph2ps where F16C is available; the kernel is JVM-only and its provider already gates on JDK 21+, so it is always reachable here. Make the codec straight-line too, for the targets that have no intrinsic and still go through ScalarFp16MatmulKernel. The loop is unnecessary: binary16 subnormals are mant * 2^-24 with mant in [1, 1023], both factors exact in FP32 and the product a normal FP32, so one multiply lets the hardware renormalize. mant == 0 falls out of the same expression as +-0. Quiet NaN on decode, which is a deliberate behaviour change. The two implementations otherwise agree bit-for-bit on all 65536 inputs, but the hardware conversion quiets signaling NaNs and the old decode reproduced them verbatim, so 1022 patterns would have differed between the JVM and every other target. Quieting costs nothing on the hot path, matches encode -- which already never emits a signaling binary16 NaN -- and loses nothing real, since the first arithmetic use quiets the value anyway. The quiet bit is applied only for a non-zero mantissa, or +-Inf would decode as NaN. Tests sweep the whole 16-bit domain rather than sampling. In commonTest the old loop-based decode is kept verbatim as the oracle and must match on every non-NaN pattern, with the NaN change pinned separately: still NaN, quiet bit set, payload and sign preserved, and exactly the 1022 signaling patterns differing from the old result. A JVM-only test asserts the codec is bit-identical to Float.float16ToFloat across all 65536 inputs, which is what licenses the kernel to substitute one for the other. A kernel sweep multiplies every pattern by 1.0 into a zeroed accumulator -- exact under FMA -- so the two decode paths are compared directly, chunked at 999 columns so each chunk exercises both the vectorized body and the scalar tail. Green on jvm and linuxX64. API is unchanged. Refs #887.

Calling the JDK intrinsic per element bought about 1.6x but left FP16 1.2-11x slower than FP32, so the decode was still the bottleneck rather than the conversion arithmetic. The BF16 kernel fills its scratch buffer scalar-wise too, and it is 7-11x faster, which points at the scalar fill loop itself: BF16's body is a shift the JIT can autovectorize, while a float16ToFloat call in the loop body cannot be lifted into the SIMD domain the same way, so the fill stays element-at-a-time. Fill the scratch buffer with the raw 16-bit patterns instead of decoded floats, load it as an IntVector, and widen a whole vector at a time: shift the sign-free pattern left by 13 to land binary16's fields in FP32 positions, rebias the exponent, then apply the two special cases under vector masks rather than branches. Inf/NaN takes a second rebias that saturates the exponent; zero and subnormals are bumped one exponent step and have 2^-14 subtracted, which makes the FPU renormalize them. All branch-free, so every lane costs the same. The Vector API offers nothing better on JDK 21: it has no half-float species or conversion, and ShortVector.fromByteArray is gone, leaving only fromMemorySegment, which would pull java.lang.foreign -- preview on 21 -- into a kernel that must run without --enable-preview. Filling an IntArray keeps the load portable and costs the same integer work the BF16 fill already does. A signaling NaN stays signaling on this path, where the codec quiets it. That is unobservable: every lane feeds the FMA, and the FMA quiets it. The exhaustive kernel sweep pins exactly that contract -- bit equality with the codec on all 63490 non-NaN patterns, NaN-ness on the rest -- and passes, which also validates the branch-free algorithm across the whole domain rather than on samples. Refs #887.

This is what #887 actually is. The issue reads the FP16/BF16 gap as a slow decode, but the two Panama kernels are within about 15% of each other head to head (143 vs 164 ms for ffn_up 8B at m=1). Measured through the real dispatch, FP16 matched its Panama kernel exactly while BF16 came out 10x faster than its own -- because NativeKernelProvider carried matmulBf16 but not matmulFp16, so BF16 resolved to the FFM kernel at priority 100 and FP16 cascaded to Panama at 50. One format was served natively and the other was not; that is the whole 2-18x. Add the missing side. skainet_fp16_matmul takes the same caller contract and strides as skainet_bf16_matmul and differs in two places. The dequant. BF16 gets its conversion for free as the high half of an FP32; binary16 needs rebiasing and gradual underflow, so the conversion folds both special cases in with arithmetic masks and stays branch-free, keeping the inner loop a straight-line sequence the vectorizer can widen. No _Float16 and no F16C intrinsics: the x86_64 build carries no -march flag, so F16C cannot be assumed, and _Float16 without it lowers to libgcc helper calls that are slower than the bit math and block vectorization outright. AArch64 does build with +fp16, but a second path would double the surface to test for a handful of integer ops; runtime ISA dispatch is where that belongs if it ever pays. The iteration order. i-p-j re-decodes every B element once per row of A, which BF16 can afford at one shift per element and this kernel cannot -- it left FP16 still 1.55-2.06x slower than FP32 at m=16 even natively. So j is tiled and each B row is decoded once per tile into a 512-float stack buffer, then multiplied into all m rows of C. Decodes drop from m*k*n to k*n with B traffic unchanged, and no allocation enters the kernel. Accumulation into any given C element is still p ascending, so the result is bit-identical to the i-p-j formulation, not merely close. Also fill in the "Float16" arm of KernelProvider.supports, which was absent while every other matmul dtype was present -- the same omission as the missing accessor, one layer up. Measured on i7-9750H / OpenJDK 21, median ms, against the FP32 SGEMM: shape batch fp32 fp16 before fp16 after q_proj 1B 16 15.59 298.91 10.31 q_proj 8B 1 38.91 73.82 30.00 q_proj 8B 16 94.79 1182.08 57.72 ffn_up 8B 1 107.37 199.49 82.07 ffn_up 8B 16 258.96 3200.78 155.33 ffn_down 8B 16 254.13 3191.62 157.11 FP16 is now 1.27-1.67x faster than FP32 everywhere except 2048x2048 at batch 1, where the two are a wash. At m=16 it also passes BF16, which still pays the per-row decode -- the same amortization would help there, deliberately left out so BF16's measured behaviour changes in its own commit with its own numbers. Parity covers the shapes the BF16 test covers, plus two the random ones never reach: a weight set of subnormals, zeros and the format extremes, and an exhaustive sweep multiplying all 65536 patterns by 1.0 into a zeroed accumulator, which pins the C conversion against Fp16Codec across the whole domain rather than trusting sampled shapes to have hit a subnormal. A signaling NaN stays signaling on both native and Panama paths; the multiply quiets it, so only NaN-ness is asserted there. Refs #887.

…matmul Tiling j amortizes the decode across rows of A, which is why it is there, but at m == 1 there is nothing to amortize -- every B element is used exactly once either way -- and it trades sequential row streaming for a column-block walk. Measured 15% slower at m == 1 on ffn_up 8B (71 ms i-p-j against 82 ms tiled). That is the decode step of inference, so it is the wrong place to lose 15%. Branch on m and keep plain i-p-j for the single-row case. Accumulation stays p ascending on both paths, so they remain bit-identical to each other, which the new cross-path test asserts on raw bits rather than within a tolerance. Two coverage gaps went with it. The exhaustive decode sweep uses m == 1, so after this change it no longer touches the tiled loop at all -- it is now run a second time with a zero-weighted second row. And every parity shape was either m == 1 or n <= 256, so the tiled path only ever ran as one full tile and the tile-boundary arithmetic was never exercised; n = 1100 adds two full tiles plus a 76-column remainder. Refs #887.

The common codec test counted how many patterns the new decode changes relative to the old one and asserted 1022. That fails on Kotlin/JS with actual 0, and the codec is not at fault: JS quiets a signaling NaN itself whenever a Float crosses float32/double, so on that target the old implementation and the new one produce identical bits and there is nothing to count. Verified by reproducing jsTest locally, and by checking the round trip directly -- 7f802000 comes back as 7fc02000 while payload and sign survive untouched. So the count moves to the JVM test, which is the only place no platform sits in between, and it now also asserts that every changed pattern *is* a signaling NaN rather than only counting them -- a change outside that set would be a regression, not the intended quieting. What stays in commonTest is the part that is genuinely portable and is the actual contract: a NaN pattern decodes to a NaN, quiet, with payload and sign preserved. That passes on jvm, js, wasmJs and linuxX64. Refs #887.

github-actions · 2026-07-30T06:33:05Z

📖 Documentation Preview

The documentation has been built successfully for this PR.

Generated Files:

Operator documentation: docs/modules/operators/_generated_/
JSON schema output: operators.json

Artifacts:

Download the documentation-preview-896 artifact to view the complete documentation locally.

This comment will be updated automatically when the PR is updated.

i-p-j walks the whole of B once per row of A. For ffn_up 8B at m=16 that is 16 passes over 90 MiB, 1.4 GiB of traffic to do 1.4 GFLOP. Tile j instead and widen each B row once per tile into a small stack buffer, then multiply it into all m rows of C, so B is read once in total. Keep plain i-p-j at m == 1. There every B element is used exactly once either way, so tiling only trades sequential streaming for a column-block walk -- it cost the FP16 kernel 15% at m == 1, and m == 1 is the decode step of inference. Measured on i7-9750H / OpenJDK 21, median ms, fp32 column as the scale (the two runs differ by ~3% on the baseline): shape batch fp32 bf16 before bf16 after q_proj 1B 16 16.15 10.52 9.30 q_proj 8B 16 99.97 65.35 59.42 ffn_up 8B 16 270.27 178.94 145.41 ffn_down 8B 16 254.50 174.43 143.35 ffn_up 8B 1 108.53 56.78 56.94 9-19% at m=16 and unchanged at m=1. Worth noting that is far less than cutting memory traffic 16x would suggest: at 90 MiB the remaining B traffic is a few ms, so what is left is the FMA chain. This kernel is compute-bound at m=16, not bandwidth-bound, and the next real win there is a blocked microkernel or bfdot on ARMv8.6-A+, not more layout work. Accumulation into any given C element stays p ascending on both paths, so results are bit-identical to the previous formulation, not merely within tolerance. The new cross-path test asserts that on raw bits. Two coverage gaps closed while here: every existing parity shape was either m == 1 or n <= 256, so a tiled path would only ever have run as a single full tile with its boundary arithmetic never exercised -- n = 1100 adds two full tiles plus a 76-column remainder. Follows the same change to skainet_fp16_matmul in SKaiNET-developers#896.

michalharakal mentioned this pull request Jul 29, 2026

perf(backend-native-cpu): read B once per matmul in the BF16 kernel #897

Merged

michalharakal requested a review from aharakal July 30, 2026 04:23

michalharakal added 5 commits July 30, 2026 08:29

michalharakal force-pushed the fix/fp16-decode-intrinsic-887 branch from 0d94ba3 to 1b83d12 Compare July 30, 2026 06:30

aharakal approved these changes Jul 30, 2026

View reviewed changes

michalharakal merged commit 0237082 into develop Jul 30, 2026
21 checks passed

michalharakal deleted the fix/fp16-decode-intrinsic-887 branch July 30, 2026 07:43

michalharakal mentioned this pull request Jul 30, 2026

release: prepare 0.38.0 #904

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

feat(backend-native-cpu): native FP16 matmul kernel - #896

feat(backend-native-cpu): native FP16 matmul kernel#896
michalharakal merged 5 commits into
developfrom
fix/fp16-decode-intrinsic-887

michalharakal commented Jul 29, 2026 •

edited

Loading

Uh oh!

github-actions Bot commented Jul 29, 2026

Uh oh!

github-actions Bot commented Jul 29, 2026

Uh oh!

github-actions Bot commented Jul 29, 2026

Uh oh!

github-actions Bot commented Jul 30, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

michalharakal commented Jul 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

The issue's root cause is wrong

The fix

Measured

Tests

NaN behaviour (deliberate change)

Platform caveat

Uh oh!

github-actions Bot commented Jul 29, 2026

Uh oh!

github-actions Bot commented Jul 29, 2026

Uh oh!

github-actions Bot commented Jul 29, 2026

Uh oh!

github-actions Bot commented Jul 30, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

michalharakal commented Jul 29, 2026 •

edited

Loading