Add FP8 muladd with fast_acc argument; add i8/u8 mma tests

Author: AntonOresten

README.md

@@ -197,7 +197,23 @@ Tile IR operations.
 | Operation | Description |
 |-----------|-------------|
 | `a * b` | Matrix multiplication: `a @ b` |
-| `muladd(a, b, acc)` | Matrix multiply-accumulate: `a * b + acc` |
+| `muladd(a, b, acc; fast_acc=false)` | Matrix multiply-accumulate: `a * b + acc` |
+| `ct.muladd_scaled(a, a_scale, b, b_scale, acc)` | Block-scaled multiply-accumulate |
+
+Each operation follows `Base.:*` / `Base.muladd`'s shape rules, with the addition of allowing trailing batch dimensions.
+
+`fast_acc=true` enables fast accumulation for FP8 inputs, and has an effect only on Hopper (sm_90; silently ignored on other
+architectures), and requires Tile IR v13.3+.
+
+`ct.muladd_scaled` multiplies each operand by a low-precision block scale before the matmul:
+each scale element covers a contiguous block of `B = K ÷ K_s` elements along the K dimension.
+Requires Blackwell. The supported operand/scale/accumulator dtypes and block sizes are:
+
+| Input (`a`/`b`) | Scale | Acc/Output | B |
+|-----------------|-------|------------|--------|
+| `Float8_E4M3FN`, `Float8_E5M2` | `Float8_E8M0FNU` | `Float32` | 32 |
+| `Float4_E2M1FN` | `Float8_E8M0FNU` | `Float32` | 16, 32 |
+| `Float4_E2M1FN` | `Float8_E4M3FN` | `Float32` | 16 |
 
 ### Higher-Order Functions
 | Operation | Description |

ext/DLFP8TypesExt.jl

@@ -18,6 +18,10 @@ end
 ct.mma_allowed_acc_dtypes(::Type{Float8_E4M3FN}) = (Float16, Float32)
 ct.mma_allowed_acc_dtypes(::Type{Float8_E5M2})   = (Float16, Float32)
 
+# `fast_acc` (lower-precision MMA accumulation) is an FP8-only throughput hint.
+ct.mma_supports_fast_acc(::Type{Float8_E4M3FN}) = true
+ct.mma_supports_fast_acc(::Type{Float8_E5M2})   = true
+
 # Float ↔ FP8 scalar constructor overlays (for map/convert dispatch)
 const FP8Types = (Float8_E4M3FN, Float8_E5M2)
 const StandardFloats = (Float16, ct.BFloat16, Float32, ct.TFloat32, Float64)

ext/MicrofloatsExt.jl

@@ -32,6 +32,10 @@ ct.ftof_rounding_mode(::Type{Float8_E8M0FNU}) = ct.RoundingMode.Zero
 ct.mma_allowed_acc_dtypes(::Type{Float8_E4M3FN}) = (Float16, Float32)
 ct.mma_allowed_acc_dtypes(::Type{Float8_E5M2})   = (Float16, Float32)
 
+# `fast_acc` (lower-precision MMA accumulation) is an FP8-only throughput hint.
+ct.mma_supports_fast_acc(::Type{Float8_E4M3FN}) = true
+ct.mma_supports_fast_acc(::Type{Float8_E5M2})   = true
+
 # Float ↔ microfloat scalar constructor overlays (for map/convert dispatch).
 # Mirrors DLFP8TypesExt: route to `Intrinsics.ftof` so kernel-side conversions
 # lower to the FToFOp Tile IR intrinsic instead of Microfloats' Float32-fallback

src/compiler/intrinsics/core.jl

@@ -415,6 +415,17 @@ mma_allowed_acc_dtypes(::Type{TFloat32}) = (Float32,)
 mma_allowed_acc_dtypes(::Type{Float64})  = (Float64,)
 mma_allowed_acc_dtypes(@nospecialize(::Type)) = nothing
 
+"""
+    mma_supports_fast_acc(input_T) -> Bool
+
+Whether `mma`'s `fast_acc` hint (lower-precision accumulation for throughput)
+is valid for operand element type `input_T`. Only the FP8 dtypes qualify;
+extensions (Microfloats/DLFP8Types) register them by overloading this, like
+[`mma_allowed_acc_dtypes`](@ref). Mirrors cuTile Python's
+`use_fast_acc is only supported for fp8 input dtypes` check.
+"""
+mma_supports_fast_acc(@nospecialize(::Type)) = false
+
 # First (preferred) accumulator dtype for a given input dtype. Used by
 # matmul to pick `acc = zeros(first_allowed_acc(T), …)` so that the input
 # dtype constraint and tileiras's acc-dtype constraint stay consistent
@@ -425,7 +436,7 @@ first_allowed_acc_dtype(::Type{T}) where {T} =
     end
 
 """
-    Intrinsics.mma(a::Tile, b::Tile, acc::Tile) -> typeof(acc)
+    Intrinsics.mma(a::Tile, b::Tile, acc::Tile, fast_acc::Bool=false) -> typeof(acc)
 
 Matrix-multiply-accumulate computing `a*b + acc`. Dispatches at codegen
 based on element types:
@@ -436,15 +447,23 @@ based on element types:
 - `i8` `a`/`b` with `i32` `acc` lower to `cuda_tile.mmai`. Per-input
   signedness is derived from the Julia type (`Int8` → signed,
   `UInt8` → unsigned); `acc` and the result are always signed `i32`.
+
+`fast_acc` enables fast accumulation (trading accumulator precision for
+throughput). It is only valid for FP8 inputs (see
+[`mma_supports_fast_acc`](@ref)) and requires Tile IR v13.3+.
 """
-@intrinsic mma(a::Tile, b::Tile, acc::Tile)
-tfunc(𝕃, ::typeof(Intrinsics.mma), @nospecialize(a), @nospecialize(b), @nospecialize(acc)) = CC.widenconst(acc)
+@intrinsic mma(a::Tile, b::Tile, acc::Tile, fast_acc::Bool=false)
+tfunc(𝕃, ::typeof(Intrinsics.mma), @nospecialize(a), @nospecialize(b), @nospecialize(acc),
+      @nospecialize(rest...)) = CC.widenconst(acc)
 function emit_intrinsic!(ctx::CGCtx, ::typeof(Intrinsics.mma), args)
     cb = ctx.cb
 
     lhs = emit_value!(ctx, args[1])
     rhs = emit_value!(ctx, args[2])
     acc = emit_value!(ctx, args[3])
+    fast_acc = length(args) >= 4 ?
+        (@something get_constant(ctx, args[4]) throw(IRError("mma: fast_acc must be a compile-time constant"))) :
+        false
 
     (lhs === nothing || rhs === nothing || acc === nothing) && throw(IRError("Cannot resolve operands for mma()"))
 
@@ -467,9 +486,15 @@ function emit_intrinsic!(ctx::CGCtx, ::typeof(Intrinsics.mma), args)
         acc_elem in allowed_acc ||
             throw(IRError("mma: acc dtype $acc_elem is not allowed for input dtype " *
                           "$lhs_elem; tileiras requires acc ∈ $allowed_acc"))
-        encode_MmaFOp!(cb, acc.type_id, lhs.v, rhs.v, acc.v)
+        # fast_acc is an FP8-only throughput hint; reject it elsewhere with a
+        # clear error rather than emitting IR tileiras would reject.
+        fast_acc && !Base.invokelatest(mma_supports_fast_acc, lhs_elem) &&
+            throw(IRError("mma: fast_acc is only supported for fp8 input dtypes " *
+                          "(f8e4m3fn, f8e5m2), got $lhs_elem"))
+        encode_MmaFOp!(cb, acc.type_id, lhs.v, rhs.v, acc.v; fast_acc)
     elseif lhs_elem <: Union{Int8, UInt8} && rhs_elem <: Union{Int8, UInt8} &&
            acc_elem === Int32
+        fast_acc && throw(IRError("mma: fast_acc is not supported for integer (mmai) inputs"))
         s_lhs = lhs_elem <: Signed ? Signedness.Signed : Signedness.Unsigned
         s_rhs = rhs_elem <: Signed ? Signedness.Signed : Signedness.Unsigned
         encode_MmaIOp!(cb, acc.type_id, lhs.v, rhs.v, acc.v;

src/language/operations.jl

@@ -1296,7 +1296,7 @@ end
 =============================================================================#
 
 """
-    muladd(a::Tile, b::Tile, acc::Tile) -> Tile
+    muladd(a::Tile, b::Tile, acc::Tile; fast_acc::Bool=false) -> Tile
 
 Matrix multiply-accumulate `a * b + acc` over tiles, lowering to
 `cuda_tile.mmaf` (float) or `cuda_tile.mmai` (`i8 × i8 → i32`).
@@ -1306,45 +1306,49 @@ promoted (vec-mat / mat-vec) and any trailing dimensions (≥3-D) are treated as
 broadcast batch dims, lifting `Base.muladd`'s shape rules to tiles. `acc`
 carries the result dtype, which must be one tileiras allows for the input dtype
 (f16/f32 for f16 and f8; f32 for bf16/tf32; f64 for f64; i32 for i8).
+
+`fast_acc` enables fast accumulation (lower accumulator precision for
+throughput); it is valid only for FP8 inputs and requires Tile IR v13.3+.
 """
-@inline function Base.muladd(a::Tile{T1, SA}, b::Tile{T2, SB}, acc::Tile{T3, SC}) where {T1, T2, T3, SA, SB, SC}
+@inline function Base.muladd(a::Tile{T1, SA}, b::Tile{T2, SB}, acc::Tile{T3, SC};
+                             fast_acc::Bool=false) where {T1, T2, T3, SA, SB, SC}
     # SA, SB, SC type parameters avoid ambiguity with the scalar `muladd`
     # methods during codegen.
-    _muladd(a, b, acc, Val(ndims(a)), Val(ndims(b)))
+    _muladd(a, b, acc, Val(ndims(a)), Val(ndims(b)), fast_acc)
 end
 
 # 2D × 2D: MmaFOp with swapped operands for row-major Tile IR
 # Julia (M,K)*(K,N) → TileIR (K,M)*(N,K) → mmaf(b,a,acc) → TileIR (N,M) → Julia (M,N)
-@inline function _muladd(a::Tile, b::Tile, acc::Tile, ::Val{2}, ::Val{2})
-    Intrinsics.mma(b, a, acc)
+@inline function _muladd(a::Tile, b::Tile, acc::Tile, ::Val{2}, ::Val{2}, fast_acc::Bool)
+    Intrinsics.mma(b, a, acc, fast_acc)
 end
 
 # Vec-mat (1D × 2D): reshape (M,) → (M, 1), MmaFOp, acc is already (M, N)
-@inline function _muladd(a::Tile, b::Tile, acc::Tile, ::Val{1}, ::Val{2})
+@inline function _muladd(a::Tile, b::Tile, acc::Tile, ::Val{1}, ::Val{2}, fast_acc::Bool)
     a2d = reshape(a, (size(a, 1), 1))
-    _muladd(a2d, b, acc, Val(2), Val(2))
+    _muladd(a2d, b, acc, Val(2), Val(2), fast_acc)
 end
 
 # Mat-vec (2D × 1D): reshape b (K,) → (K, 1), acc (M,) → (M, 1), MmaFOp, squeeze back
-@inline function _muladd(a::Tile, b::Tile, acc::Tile, ::Val{2}, ::Val{1})
+@inline function _muladd(a::Tile, b::Tile, acc::Tile, ::Val{2}, ::Val{1}, fast_acc::Bool)
     M, K = size(a, 1), size(b, 1)
     b2d = reshape(b, (K, 1))
     acc2d = reshape(acc, (M, 1))
-    result = _muladd(a, b2d, acc2d, Val(2), Val(2))
+    result = _muladd(a, b2d, acc2d, Val(2), Val(2), fast_acc)
     reshape(result, (M,))
 end
 
 # Vec-vec (1D × 1D): not supported
-@generated function _muladd(::Tile, ::Tile, ::Tile, ::Val{1}, ::Val{1})
+@generated function _muladd(::Tile, ::Tile, ::Tile, ::Val{1}, ::Val{1}, ::Bool)
     return :(throw(ArgumentError("Vector-vector multiply-accumulate is not supported.")))
 end
 
 # Batched mat-vec / vec-mat (≥3D × 1D or 1D × ≥3D): not supported, unsqueeze manually
-@generated function _muladd(::Tile, ::Tile, ::Tile, ::Val{1}, ::Val{NB}) where {NB}
+@generated function _muladd(::Tile, ::Tile, ::Tile, ::Val{1}, ::Val{NB}, ::Bool) where {NB}
     NB >= 3 || return :(throw(ArgumentError("unreachable")))
     return :(throw(ArgumentError("Batched vec-mat is not supported. Reshape the 1D operand to 2D first.")))
 end
-@generated function _muladd(::Tile, ::Tile, ::Tile, ::Val{NA}, ::Val{1}) where {NA}
+@generated function _muladd(::Tile, ::Tile, ::Tile, ::Val{NA}, ::Val{1}, ::Bool) where {NA}
     NA >= 3 || return :(throw(ArgumentError("unreachable")))
     return :(throw(ArgumentError("Batched mat-vec is not supported. Reshape the 1D operand to 2D first.")))
 end
@@ -1357,7 +1361,7 @@ end
 #   3. MmaFOp with swapped operands: mmaf(b, a, acc)
 #   4. Unflatten batch dims via reshape
 @generated function _muladd(a::Tile{T1, SA}, b::Tile{T2, SB}, acc::Tile{T3, SC},
-                            ::Val{NA}, ::Val{NB}) where {T1, T2, T3, SA, SB, SC, NA, NB}
+                            ::Val{NA}, ::Val{NB}, fast_acc::Bool) where {T1, T2, T3, SA, SB, SC, NA, NB}
     sa = Tuple(SA.parameters)
     sb = Tuple(SB.parameters)
 
@@ -1384,7 +1388,7 @@ end
         b_3d = reshape(b_bc, $((K, N, B_flat)))
         acc_3d = reshape(acc_bc, $((M, N, B_flat)))
         # MmaFOp with swapped operands for row-major convention
-        result_3d = Intrinsics.mma(b_3d, a_3d, acc_3d)
+        result_3d = Intrinsics.mma(b_3d, a_3d, acc_3d, fast_acc)
         # Unflatten batch dims
         reshape(result_3d, $((M, N, batch_shape...)))
     end

test/device/integration.jl

@@ -33,3 +33,40 @@ using CUDA
 
     @test c_cpu ≈ c_ref
 end
+
+@testset "i8/u8 matmul (mmai)" begin
+    # i8/u8 × i8/u8 → i32 lowers to `cuda_tile.mmai`. The accumulator must be
+    # Int32 (the `*` operator would pick the input dtype as acc and fail), so we
+    # use `muladd` with an explicit i32 acc. Per-operand signedness is derived
+    # from the Julia element type, so we feed values that differ under signed vs
+    # unsigned interpretation (negatives, and magnitudes > 127) to confirm each
+    # operand is interpreted correctly. K = 16, |product| ≤ 255² · 16 ≈ 1.0e6,
+    # well within Int32, so the result is exact.
+    function mmai(a::ct.TileArray{T1,2}, b::ct.TileArray{T2,2}, c::ct.TileArray{Int32,2}) where {T1,T2}
+        ta = ct.load(a, (1, 1), (16, 16))
+        tb = ct.load(b, (1, 1), (16, 16))
+        tc = muladd(ta, tb, zeros(Int32, (16, 16)))
+        ct.store(c, (1, 1), tc)
+        return
+    end
+
+    M = K = N = 16
+    @testset "signed × signed" begin
+        a = rand(Int8(-128):Int8(127), M, K); b = rand(Int8(-128):Int8(127), K, N)
+        c = CUDA.zeros(Int32, M, N)
+        @cuda backend=cuTile blocks=1 mmai(CuArray(a), CuArray(b), c)
+        @test Array(c) == Int32.(a) * Int32.(b)
+    end
+    @testset "unsigned × unsigned" begin
+        a = rand(UInt8(0):UInt8(255), M, K); b = rand(UInt8(0):UInt8(255), K, N)
+        c = CUDA.zeros(Int32, M, N)
+        @cuda backend=cuTile blocks=1 mmai(CuArray(a), CuArray(b), c)
+        @test Array(c) == Int32.(a) * Int32.(b)
+    end
+    @testset "unsigned × signed" begin
+        a = rand(UInt8(0):UInt8(255), M, K); b = rand(Int8(-128):Int8(127), K, N)
+        c = CUDA.zeros(Int32, M, N)
+        @cuda backend=cuTile blocks=1 mmai(CuArray(a), CuArray(b), c)
+        @test Array(c) == Int32.(a) * Int32.(b)
+    end
+end

test/extensions/DLFP8Types.jl

@@ -45,14 +45,27 @@ end
     end
 end
 
+# `fast_acc=true` is an FP8-only hint; it still lowers to `mmaf` (13.3+).
+@test @filecheck begin
+    @check_label "entry"
+    code_tiled(Tuple{ct.TileArray{Float8_E4M3FN,2,spec2d}, ct.TileArray{Float8_E4M3FN,2,spec2d},
+                     ct.TileArray{Float32,2,spec2d}}; bytecode_version=v"13.3") do a, b, c
+        ta = ct.load(a, (1, 1), (16, 16))
+        tb = ct.load(b, (1, 1), (16, 16))
+        @check "mmaf"
+        ct.store(c, (1, 1), muladd(ta, tb, zeros(Float32, (16, 16)); fast_acc=true))
+        return
+    end
 end
 
-# FP8 (e4m3/e5m2) conversions and matmul need Hopper (sm_90+).
-@testset "execution" begin
-if capability(device()) >= v"9"
+end
 
-# Round-trip Float32 → FP8 → Float32 on values exactly representable in
-# the target FP8 type — result must match input bit-for-bit.
+# Execution kernels are plain top-level functions, each defined next to the
+# test that exercises it. Kernels parametric on accumulator dtype must stay at
+# top level — defining them inside a testset scope boxes them into closures.
+
+# Round-trip Float32 → FP8 → Float32 on values exactly representable in the
+# target FP8 type — result must match input bit-for-bit.
 function rt_e4m3(a::ct.TileArray{Float32,1}, b::ct.TileArray{Float32,1})
     pid = ct.bid(1)
     tile = ct.load(a, pid, (16,))
@@ -65,19 +78,8 @@ function rt_e5m2(a::ct.TileArray{Float32,1}, b::ct.TileArray{Float32,1})
     ct.store(b, pid, convert(ct.Tile{Float32}, convert(ct.Tile{Float8_E5M2}, tile)))
     return
 end
-
-representable = Float32[0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 8.0,
-                        16.0, 32.0, 64.0, 128.0, 256.0, -1.0, -2.0, -0.5]
-let a = CuArray(representable), b = CUDA.zeros(Float32, length(representable))
-    @cuda backend=cuTile blocks=1 rt_e4m3(a, b)
-    @test Array(b) == representable
-    @cuda backend=cuTile blocks=1 rt_e5m2(a, b)
-    @test Array(b) == representable
-end
-
 # FMA in FP8: load Float32, convert to FP8, multiply-add in FP8, convert back.
-# Uses inputs whose products and sums also stay representable, so the result
-# is exact.
+# Inputs whose products and sums also stay representable, so the result is exact.
 function fma_e4m3(a::ct.TileArray{Float32,1}, b::ct.TileArray{Float32,1},
                   c::ct.TileArray{Float32,1}, d::ct.TileArray{Float32,1})
     pid = ct.bid(1)
@@ -87,6 +89,33 @@ function fma_e4m3(a::ct.TileArray{Float32,1}, b::ct.TileArray{Float32,1},
     ct.store(d, pid, convert(ct.Tile{Float32}, muladd.(ta, tb, tc)))
     return
 end
+# Non-scaled FP8 matmul with both allowed accumulator dtypes (f16 and f32).
+function mma_dl_fp8(A::ct.TileArray{Float8_E4M3FN,2}, B::ct.TileArray{Float8_E4M3FN,2},
+                    C::ct.TileArray{Tacc,2}, D::ct.TileArray{Float32,2}) where {Tacc<:Union{Float16,Float32}}
+    a = ct.load(A, (1, 1), (16, 16)); b = ct.load(B, (1, 1), (16, 16)); c = ct.load(C, (1, 1), (16, 16))
+    ct.store(D, (1, 1), convert(ct.Tile{Float32}, muladd(a, b, c)))
+    return
+end
+function mma_dl_fast(A::ct.TileArray{Float8_E4M3FN,2}, B::ct.TileArray{Float8_E4M3FN,2},
+                     C::ct.TileArray{Float32,2}, D::ct.TileArray{Float32,2})
+    a = ct.load(A, (1, 1), (16, 16)); b = ct.load(B, (1, 1), (16, 16)); c = ct.load(C, (1, 1), (16, 16))
+    ct.store(D, (1, 1), convert(ct.Tile{Float32}, muladd(a, b, c; fast_acc=true)))
+    return
+end
+
+# FP8 (e4m3/e5m2) conversions and matmul need Hopper (sm_90+).
+@testset "execution" begin
+if capability(device()) >= v"9"
+
+representable = Float32[0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 8.0,
+                        16.0, 32.0, 64.0, 128.0, 256.0, -1.0, -2.0, -0.5]
+let a = CuArray(representable), b = CUDA.zeros(Float32, length(representable))
+    @cuda backend=cuTile blocks=1 rt_e4m3(a, b)
+    @test Array(b) == representable
+    @cuda backend=cuTile blocks=1 rt_e5m2(a, b)
+    @test Array(b) == representable
+end
+
 let av = Float32[1.0, 2.0, 0.5, 4.0, 1.5, 2.0, -1.0, -0.5, 3.0, 0.5, 1.0, 2.0, -2.0, 1.0, 0.5, 4.0],
     bv = Float32[2.0, 1.0, 4.0, 0.5, 2.0, 3.0,  2.0,  4.0, 1.0, 2.0, 1.0, 0.5,  2.0, 1.0, 2.0, 1.0],
     cv = Float32[0.0, 1.0, 0.0, 0.0, 1.0, 1.0,  0.0,  0.0, 1.0, 0.0, 0.0, 1.0,  0.0, 0.0, 1.0, 0.0]
@@ -96,29 +125,30 @@ let av = Float32[1.0, 2.0, 0.5, 4.0, 1.5, 2.0, -1.0, -0.5, 3.0, 0.5, 1.0, 2.0, -
     @test Array(d) == av .* bv .+ cv
 end
 
-# Non-scaled FP8 matmul with both allowed accumulator dtypes (f16 and f32).
-function mma_dl_f32(A::ct.TileArray{Float8_E4M3FN,2}, B::ct.TileArray{Float8_E4M3FN,2},
-                    C::ct.TileArray{Float32,2}, D::ct.TileArray{Float32,2})
-    a = ct.load(A, (1, 1), (16, 16)); b = ct.load(B, (1, 1), (16, 16)); c = ct.load(C, (1, 1), (16, 16))
-    ct.store(D, (1, 1), convert(ct.Tile{Float32}, muladd(a, b, c)))
-    return
-end
-function mma_dl_f16(A::ct.TileArray{Float8_E4M3FN,2}, B::ct.TileArray{Float8_E4M3FN,2},
-                    C::ct.TileArray{Float16,2}, D::ct.TileArray{Float32,2})
-    a = ct.load(A, (1, 1), (16, 16)); b = ct.load(B, (1, 1), (16, 16)); c = ct.load(C, (1, 1), (16, 16))
-    ct.store(D, (1, 1), convert(ct.Tile{Float32}, muladd(a, b, c)))
-    return
-end
-@testset "mma → $Tacc acc" for (Tacc, kern) in ((Float32, mma_dl_f32), (Float16, mma_dl_f16))
+@testset "mma → $Tacc acc" for Tacc in (Float32, Float16)
     M = 16
     ah = Float8_E4M3FN.(Float32.(rand(0:2, M, M)) ./ 2)
     bh = Float8_E4M3FN.(Float32.(rand(0:2, M, M)) ./ 2)
     ch = Tacc.(Float32.(rand(0:2, M, M)))
     ref = Float32.(ah) * Float32.(bh) .+ Float32.(ch)
     D = CUDA.zeros(Float32, M, M)
-    @cuda backend=cuTile blocks=1 kern(CuArray(ah), CuArray(bh), CuArray(ch), D)
+    @cuda backend=cuTile blocks=1 mma_dl_fp8(CuArray(ah), CuArray(bh), CuArray(ch), D)
     @test Array(D) == ref
 end
 
+# fast_acc only has an effect on Hopper (sm_90); ignored elsewhere. So off
+# Hopper we assert the exact result (the flag must ride through without
+# perturbing the output); on Hopper we make no numeric claim.
+@testset "mma fast_acc (exact off Hopper)" begin
+    M = 16
+    ah = Float8_E4M3FN.(Float32.(rand(0:2, M, M)) ./ 2)
+    bh = Float8_E4M3FN.(Float32.(rand(0:2, M, M)) ./ 2)
+    ch = Float32.(rand(0:2, M, M))
+    ref = Float32.(ah) * Float32.(bh) .+ ch
+    D = CUDA.zeros(Float32, M, M)
+    @cuda backend=cuTile blocks=1 mma_dl_fast(CuArray(ah), CuArray(bh), CuArray(ch), D)
+    @test (Array(D) == ref) || (v"9" <= capability(device()) < v"10")
+end
+
 end
 end