feat: shared narrow-float (BF16 + FP16) layer

skainet-backends/skainet-backend-api/src/commonMain/kotlin/sk/ainet/backend/api/kernel/Bf16MatmulKernel.kt

@@ -1,5 +1,8 @@
 package sk.ainet.backend.api.kernel
 
+import sk.ainet.lang.types.Bf16Codec
+import sk.ainet.lang.types.NarrowFloatCodec
+
 /**
  * FP32 input × BF16-packed weight matrix multiplication:
  *
@@ -43,7 +46,10 @@ package sk.ainet.backend.api.kernel
  * and byte-stride accessors on the same `ByteArray`. For a contiguous
  * `(k, n)` BF16 matrix `bByteStride == n * 2`.
  */
-public interface Bf16MatmulKernel {
+public interface Bf16MatmulKernel : NarrowFloatMatmulKernel {
+
+    override val codec: NarrowFloatCodec get() = Bf16Codec
+
     /**
      * @param a left operand `(m, k)`, row-major FP32, stride [aStride]
      *   floats per row.
@@ -64,7 +70,7 @@ public interface Bf16MatmulKernel {
      * @param k contraction dimension (cols of A == rows of B). `k == 0`
      *   zeros the `m × n` output block.
      */
-    public fun matmul(
+    override fun matmul(
         a: FloatArray, aOffset: Int, aStride: Int,
         b: ByteArray, bByteOffset: Int, bByteStride: Int,
         out: FloatArray, outOffset: Int, outStride: Int,

skainet-backends/skainet-backend-api/src/commonMain/kotlin/sk/ainet/backend/api/kernel/KernelProvider.kt

@@ -61,6 +61,12 @@ public interface KernelProvider {
      */
     public fun matmulBf16(): Bf16MatmulKernel? = null
 
+    /**
+     * IEEE binary16 weight matmul, or `null` when this provider carries no FP16 kernel.
+     * Callers cascade to a lower-priority provider, ultimately the scalar reference.
+     */
+    public fun matmulFp16(): Fp16MatmulKernel? = null
+
     /**
      * F32 × Q8_0 matmul kernel exposed by this provider, or `null` if
      * this provider does not specialize Q8_0. Same fall-through pattern.

skainet-backends/skainet-backend-api/src/commonMain/kotlin/sk/ainet/backend/api/kernel/NarrowFloatMatmulKernel.kt

@@ -0,0 +1,61 @@
+package sk.ainet.backend.api.kernel
+
+import sk.ainet.lang.types.Fp16Codec
+import sk.ainet.lang.types.NarrowFloatCodec
+
+/**
+ * Matmul over a **16-bit float** right operand: `C(m,n) = A(m,k) · B(k,n)` where `B` is packed at
+ * 2 bytes per element and `A`/`C` are FP32.
+ *
+ * The shape, offset, aliasing and "fully overwrite the `m × n` block" obligations are identical to
+ * `Fp32MatmulKernel`. Only the inner-loop decode of `B` differs, which is what [codec] selects — so
+ * BF16 and FP16 differ by a decode step and nothing else.
+ *
+ * **Accumulation is always FP32.** The narrow format is a storage width, never an accumulate width;
+ * every implementation widens each `B` element before the multiply and sums in `Float`. This
+ * matches what PyTorch, JAX and tensor cores do, and is not a shortcut to be "optimized" away.
+ *
+ * ## Stride convention
+ *
+ * `A` and `out` use **float** strides (counting FP32 elements). `B` uses **byte** offsets and
+ * strides, matching the byte-packed convention of the block-quantized kernels and avoiding the
+ * foot-gun of mixing element- and byte-strides on the same `ByteArray`. For a contiguous `(k, n)`
+ * matrix, `bByteStride == n * 2`.
+ */
+public interface NarrowFloatMatmulKernel {
+
+    /** The 16-bit format this kernel decodes `B` from. */
+    public val codec: NarrowFloatCodec
+
+    /**
+     * @param a left operand `(m, k)`, row-major FP32, [aStride] floats per row.
+     * @param aOffset element offset into [a] of the (0, 0) entry.
+     * @param aStride distance in **floats** between consecutive rows of [a].
+     * @param b right operand `(k, n)`, row-major, packed little-endian 16-bit floats.
+     * @param bByteOffset **byte** offset into [b] of the (0, 0) entry.
+     * @param bByteStride **byte** distance between consecutive rows of [b].
+     * @param out output `(m, n)`, row-major FP32, [outStride] floats per row.
+     * @param outOffset element offset into [out].
+     * @param outStride distance in **floats** between consecutive rows of [out].
+     * @param m rows of A and C.
+     * @param n columns of B and C.
+     * @param k contraction dimension. `k == 0` zeros the `m × n` output block.
+     */
+    public fun matmul(
+        a: FloatArray, aOffset: Int, aStride: Int,
+        b: ByteArray, bByteOffset: Int, bByteStride: Int,
+        out: FloatArray, outOffset: Int, outStride: Int,
+        m: Int, n: Int, k: Int,
+    )
+}
+
+/**
+ * Matmul over an **IEEE binary16** right operand.
+ *
+ * The FP16 counterpart to [Bf16MatmulKernel]. Note that FP16 cannot use the bit-shift decode that
+ * makes BF16 nearly free — binary16 needs exponent rebiasing — so vectorized FP16 implementations
+ * lag their BF16 equivalents. Correctness first; throughput is a separate concern.
+ */
+public interface Fp16MatmulKernel : NarrowFloatMatmulKernel {
+    override val codec: NarrowFloatCodec get() = Fp16Codec
+}

skainet-backends/skainet-backend-cpu/api/jvm/skainet-backend-cpu.api

@@ -43,6 +43,13 @@ public final class sk/ainet/context/DirectCpuExecutionContext$Companion {
 
 public final class sk/ainet/exec/kernel/PanamaVectorBf16MatmulKernel : sk/ainet/backend/api/kernel/Bf16MatmulKernel {
 	public static final field INSTANCE Lsk/ainet/exec/kernel/PanamaVectorBf16MatmulKernel;
+	public fun getCodec ()Lsk/ainet/lang/types/NarrowFloatCodec;
+	public fun matmul ([FII[BII[FIIIII)V
+}
+
+public final class sk/ainet/exec/kernel/PanamaVectorFp16MatmulKernel : sk/ainet/backend/api/kernel/Fp16MatmulKernel {
+	public static final field INSTANCE Lsk/ainet/exec/kernel/PanamaVectorFp16MatmulKernel;
+	public fun getCodec ()Lsk/ainet/lang/types/NarrowFloatCodec;
 	public fun matmul ([FII[BII[FIIIII)V
 }
 
@@ -52,6 +59,7 @@ public final class sk/ainet/exec/kernel/PanamaVectorKernelProvider : sk/ainet/ba
 	public fun getPriority ()I
 	public fun isAvailable ()Z
 	public fun matmulBf16 ()Lsk/ainet/backend/api/kernel/Bf16MatmulKernel;
+	public fun matmulFp16 ()Lsk/ainet/backend/api/kernel/Fp16MatmulKernel;
 	public fun matmulFp32 ()Lsk/ainet/backend/api/kernel/Fp32MatmulKernel;
 	public fun matmulQ4K ()Lsk/ainet/backend/api/kernel/Q4KMatmulKernel;
 	public fun matmulQ4_0 ()Lsk/ainet/backend/api/kernel/Q4_0MatmulKernel;
@@ -69,6 +77,7 @@ public final class sk/ainet/exec/kernel/PanamaVectorKernelProviderFactory : sk/a
 	public fun getPriority ()I
 	public fun isAvailable ()Z
 	public fun matmulBf16 ()Lsk/ainet/backend/api/kernel/Bf16MatmulKernel;
+	public fun matmulFp16 ()Lsk/ainet/backend/api/kernel/Fp16MatmulKernel;
 	public fun matmulFp32 ()Lsk/ainet/backend/api/kernel/Fp32MatmulKernel;
 	public fun matmulQ4K ()Lsk/ainet/backend/api/kernel/Q4KMatmulKernel;
 	public fun matmulQ4_0 ()Lsk/ainet/backend/api/kernel/Q4_0MatmulKernel;
@@ -122,6 +131,13 @@ public final class sk/ainet/exec/kernel/PanamaVectorQ8_0MatmulKernel : sk/ainet/
 
 public final class sk/ainet/exec/kernel/ScalarBf16MatmulKernel : sk/ainet/backend/api/kernel/Bf16MatmulKernel {
 	public static final field INSTANCE Lsk/ainet/exec/kernel/ScalarBf16MatmulKernel;
+	public fun getCodec ()Lsk/ainet/lang/types/NarrowFloatCodec;
+	public fun matmul ([FII[BII[FIIIII)V
+}
+
+public final class sk/ainet/exec/kernel/ScalarFp16MatmulKernel : sk/ainet/backend/api/kernel/Fp16MatmulKernel {
+	public static final field INSTANCE Lsk/ainet/exec/kernel/ScalarFp16MatmulKernel;
+	public fun getCodec ()Lsk/ainet/lang/types/NarrowFloatCodec;
 	public fun matmul ([FII[BII[FIIIII)V
 }
 
@@ -131,6 +147,7 @@ public final class sk/ainet/exec/kernel/ScalarKernelProvider : sk/ainet/backend/
 	public fun getPriority ()I
 	public fun isAvailable ()Z
 	public fun matmulBf16 ()Lsk/ainet/backend/api/kernel/Bf16MatmulKernel;
+	public fun matmulFp16 ()Lsk/ainet/backend/api/kernel/Fp16MatmulKernel;
 	public fun matmulFp32 ()Lsk/ainet/backend/api/kernel/Fp32MatmulKernel;
 	public fun matmulQ4K ()Lsk/ainet/backend/api/kernel/Q4KMatmulKernel;
 	public fun matmulQ4_0 ()Lsk/ainet/backend/api/kernel/Q4_0MatmulKernel;
@@ -148,6 +165,7 @@ public final class sk/ainet/exec/kernel/ScalarKernelProviderFactory : sk/ainet/b
 	public fun getPriority ()I
 	public fun isAvailable ()Z
 	public fun matmulBf16 ()Lsk/ainet/backend/api/kernel/Bf16MatmulKernel;
+	public fun matmulFp16 ()Lsk/ainet/backend/api/kernel/Fp16MatmulKernel;
 	public fun matmulFp32 ()Lsk/ainet/backend/api/kernel/Fp32MatmulKernel;
 	public fun matmulQ4K ()Lsk/ainet/backend/api/kernel/Q4KMatmulKernel;
 	public fun matmulQ4_0 ()Lsk/ainet/backend/api/kernel/Q4_0MatmulKernel;

skainet-backends/skainet-backend-cpu/src/commonMain/kotlin/sk/ainet/exec/kernel/ScalarFp16MatmulKernel.kt

@@ -0,0 +1,49 @@
+package sk.ainet.exec.kernel
+
+import sk.ainet.backend.api.kernel.Fp16MatmulKernel
+import sk.ainet.lang.types.Fp16Codec
+
+/**
+ * Scalar reference implementation of [Fp16MatmulKernel] — three nested loops, decode inline, no
+ * SIMD. Always available on every KMP target, and the correctness reference an accelerated FP16
+ * kernel must match within FP-ordering tolerance.
+ *
+ * Mirrors [ScalarBf16MatmulKernel]; the only difference is the decode, which for binary16 needs
+ * exponent rebiasing rather than BF16's bit shift.
+ */
+public object ScalarFp16MatmulKernel : Fp16MatmulKernel {
+    override fun matmul(
+        a: FloatArray, aOffset: Int, aStride: Int,
+        b: ByteArray, bByteOffset: Int, bByteStride: Int,
+        out: FloatArray, outOffset: Int, outStride: Int,
+        m: Int, n: Int, k: Int,
+    ) {
+        require(m >= 0 && n >= 0 && k >= 0) {
+            "ScalarFp16MatmulKernel: m, n, k must be non-negative; got m=$m n=$n k=$k"
+        }
+        if (m == 0 || n == 0) return
+        // k == 0 is an explicit "zero the output block" per the SPI.
+        if (k == 0) {
+            for (i in 0 until m) {
+                val rowOff = outOffset + i * outStride
+                for (j in 0 until n) out[rowOff + j] = 0f
+            }
+            return
+        }
+        for (i in 0 until m) {
+            val aRowOff = aOffset + i * aStride
+            val outRowOff = outOffset + i * outStride
+            for (j in 0 until n) {
+                // Accumulate in FP32 — the narrow format is storage only.
+                var sum = 0f
+                for (p in 0 until k) {
+                    val bByteIdx = bByteOffset + p * bByteStride + j * 2
+                    val lo = b[bByteIdx].toInt() and 0xFF
+                    val hi = b[bByteIdx + 1].toInt() and 0xFF
+                    sum += a[aRowOff + p] * Fp16Codec.decode((hi shl 8) or lo)
+                }
+                out[outRowOff + j] = sum
+            }
+        }
+    }
+}

skainet-backends/skainet-backend-cpu/src/commonMain/kotlin/sk/ainet/exec/kernel/ScalarKernelProvider.kt

@@ -1,6 +1,7 @@
 package sk.ainet.exec.kernel
 
 import sk.ainet.backend.api.kernel.Bf16MatmulKernel
+import sk.ainet.backend.api.kernel.Fp16MatmulKernel
 import sk.ainet.backend.api.kernel.Fp32MatmulKernel
 import sk.ainet.backend.api.kernel.KernelProvider
 import sk.ainet.backend.api.kernel.Q4KMatmulKernel
@@ -30,6 +31,8 @@ public object ScalarKernelProvider : KernelProvider {
     override fun isAvailable(): Boolean = true
     override fun matmulFp32(): Fp32MatmulKernel = ScalarMatmulKernel
     override fun matmulBf16(): Bf16MatmulKernel = ScalarBf16MatmulKernel
+
+    override fun matmulFp16(): Fp16MatmulKernel = ScalarFp16MatmulKernel
     override fun matmulQ8_0(): Q8_0MatmulKernel = ScalarQ8_0MatmulKernel
     override fun matmulQ4_0(): Q4_0MatmulKernel = ScalarQ4_0MatmulKernel
     override fun matmulQ4K(): Q4KMatmulKernel = ScalarQ4_KMatmulKernel

skainet-backends/skainet-backend-cpu/src/commonMain/kotlin/sk/ainet/exec/tensor/ops/DefaultCpuOps.kt

@@ -30,8 +30,11 @@ import sk.ainet.lang.tensor.data.Q8_0BlockTensorData
 import sk.ainet.lang.tensor.data.TensorData
 import sk.ainet.lang.tensor.data.TensorDataFactory
 import sk.ainet.lang.tensor.ops.UpsampleMode
+import sk.ainet.lang.types.BF16
+import sk.ainet.lang.types.Bf16Codec
 import sk.ainet.lang.types.FP16
 import sk.ainet.lang.types.FP32
+import sk.ainet.lang.types.Fp16Codec
 import sk.ainet.lang.types.Int32
 import sk.ainet.lang.types.Int8
 import kotlin.math.floor
@@ -2685,18 +2688,34 @@ public open class DefaultCpuOpsBase(protected val dataFactory: TensorDataFactory
 
         @Suppress("UNCHECKED_CAST")
         val outData = when (targetClass) {
-            FP32::class, FP16::class -> dataFactory.fromFloatArray<TTo, Float>(
+            FP32::class -> dataFactory.fromFloatArray<TTo, Float>(
                 tensor.shape,
                 targetClass,
                 copyTensorValuesAsFloatArray(tensor)
             ) as TensorData<TTo, V>
+            // Narrowing to 16 bits must actually ROUND. This previously shared the FP32 branch,
+            // so `convert(x, FP16)` only re-tagged the dtype and kept full FP32 precision — the
+            // resulting tensor claimed to be FP16 while holding values no FP16 can represent, with
+            // no rounding, no overflow-to-Inf past the format's range, and no NaN handling.
+            // Storage stays float-backed (matching DenseTensorDataFactory's narrow tags); it is
+            // the VALUES that now reflect the target format.
+            FP16::class, BF16::class -> {
+                val codec = if (targetClass == BF16::class) Bf16Codec else Fp16Codec
+                val src = copyTensorValuesAsFloatArray(tensor)
+                val rounded = FloatArray(src.size) { codec.decode(codec.encode(src[it])) }
+                dataFactory.fromFloatArray<TTo, Float>(
+                    tensor.shape,
+                    targetClass,
+                    rounded
+                ) as TensorData<TTo, V>
+            }
             Int32::class, Int8::class -> dataFactory.fromIntArray<TTo, Int>(
                 tensor.shape,
                 targetClass,
                 copyTensorValuesAsIntArray(tensor)
             ) as TensorData<TTo, V>
             else -> throw IllegalArgumentException(
-                "convert supports FP32, FP16, Int32, and Int8 targets, got ${targetType.name}"
+                "convert supports FP32, FP16, BF16, Int32, and Int8 targets, got ${targetType.name}"
             )
         }
         return CpuTensor(outData, this, targetClass, GradState(requiresGrad = tensor.requiresGrad))

skainet-backends/skainet-backend-cpu/src/jvmMain/kotlin/sk/ainet/exec/kernel/PanamaVectorFp16MatmulKernel.kt

@@ -0,0 +1,87 @@
+package sk.ainet.exec.kernel
+
+import jdk.incubator.vector.FloatVector
+import jdk.incubator.vector.VectorSpecies
+import sk.ainet.backend.api.kernel.Fp16MatmulKernel
+import sk.ainet.lang.types.Fp16Codec
+
+/**
+ * Panama Vector API implementation of [Fp16MatmulKernel] — the FP16 counterpart to
+ * [PanamaVectorBf16MatmulKernel], using the same i-p-j outer-product shape and FP32 FMA
+ * accumulation.
+ *
+ * **Why the decode stays scalar here.** BF16 can be widened lane-wise with an integer shift
+ * (`IntVector.lanewise(LSHL, 16).reinterpretAsFloats()`), which is why the BF16 kernel vectorizes
+ * its dequant. Binary16 needs exponent rebiasing and gradual-underflow handling, and mainstream
+ * JDKs expose no FP16 vector species, so each element is decoded scalar into a lane-width scratch
+ * buffer before the vectorized multiply-accumulate. The FMA over `n` is still fully vectorized —
+ * only the widening is not. Expect this kernel to trail the BF16 one; that is inherent to the
+ * format on this platform, not a defect in this implementation.
+ *
+ * Numerical parity vs [ScalarFp16MatmulKernel] is asserted by
+ * `PanamaVectorFp16MatmulKernelParityTest`.
+ */
+public object PanamaVectorFp16MatmulKernel : Fp16MatmulKernel {
+
+    private val floatSpecies: VectorSpecies<Float> = FloatVector.SPECIES_PREFERRED
+
+    override fun matmul(
+        a: FloatArray, aOffset: Int, aStride: Int,
+        b: ByteArray, bByteOffset: Int, bByteStride: Int,
+        out: FloatArray, outOffset: Int, outStride: Int,
+        m: Int, n: Int, k: Int,
+    ) {
+        require(m >= 0 && n >= 0 && k >= 0) {
+            "PanamaVectorFp16MatmulKernel: m, n, k must be non-negative; got m=$m n=$n k=$k"
+        }
+        if (m == 0 || n == 0) return
+        if (k == 0) {
+            for (i in 0 until m) {
+                val rowOff = outOffset + i * outStride
+                for (j in 0 until n) out[rowOff + j] = 0f
+            }
+            return
+        }
+
+        val laneCount = floatSpecies.length()
+        val bound = floatSpecies.loopBound(n)
+        val scratch = FloatArray(laneCount)
+
+        // Zero the output block first — the i-p-j outer product accumulates into it.
+        for (i in 0 until m) {
+            val rowOff = outOffset + i * outStride
+            for (j in 0 until n) out[rowOff + j] = 0f
+        }
+
+        for (i in 0 until m) {
+            val aRowOff = aOffset + i * aStride
+            val outRowOff = outOffset + i * outStride
+            for (p in 0 until k) {
+                val aIp = a[aRowOff + p]
+                val aBcast = FloatVector.broadcast(floatSpecies, aIp)
+                val bRowByteOff = bByteOffset + p * bByteStride
+                var j = 0
+                while (j < bound) {
+                    val byteBase = bRowByteOff + j * 2
+                    for (lane in 0 until laneCount) {
+                        val lo = b[byteBase + lane * 2].toInt() and 0xFF
+                        val hi = b[byteBase + lane * 2 + 1].toInt() and 0xFF
+                        scratch[lane] = Fp16Codec.decode((hi shl 8) or lo)
+                    }
+                    val bVec = FloatVector.fromArray(floatSpecies, scratch, 0)
+                    val outVec = FloatVector.fromArray(floatSpecies, out, outRowOff + j)
+                    aBcast.fma(bVec, outVec).intoArray(out, outRowOff + j)
+                    j += laneCount
+                }
+                // Tail (scalar): up to laneCount - 1 trailing columns.
+                while (j < n) {
+                    val bByteIdx = bRowByteOff + j * 2
+                    val lo = b[bByteIdx].toInt() and 0xFF
+                    val hi = b[bByteIdx + 1].toInt() and 0xFF
+                    out[outRowOff + j] += aIp * Fp16Codec.decode((hi shl 8) or lo)
+                    j++
+                }
+            }
+        }
+    }
+}

skainet-backends/skainet-backend-cpu/src/jvmMain/kotlin/sk/ainet/exec/kernel/PanamaVectorKernelProvider.kt

@@ -1,6 +1,7 @@
 package sk.ainet.exec.kernel
 
 import sk.ainet.backend.api.kernel.Bf16MatmulKernel
+import sk.ainet.backend.api.kernel.Fp16MatmulKernel
 import sk.ainet.backend.api.kernel.Fp32MatmulKernel
 import sk.ainet.backend.api.kernel.KernelProvider
 import sk.ainet.backend.api.kernel.Q4KMatmulKernel
@@ -51,6 +52,9 @@ public object PanamaVectorKernelProvider : KernelProvider {
     override fun matmulBf16(): Bf16MatmulKernel? =
         if (isAvailable()) PanamaVectorBf16MatmulKernel else null
 
+    override fun matmulFp16(): Fp16MatmulKernel? =
+        if (isAvailable()) PanamaVectorFp16MatmulKernel else null
+
     override fun matmulQ8_0(): Q8_0MatmulKernel? =
         if (isAvailable()) PanamaVectorQ8_0MatmulKernel else null