[blas] Add bfloat16-output GEMM support

[zjin-lcf]

Summary

• add the specified bfloat16 A/B/C GEMM overload with float scalars for buffer and USM APIs
• implement BF16 storage with FP32 compute in cuBLAS and expose the existing rocBLAS implementation for both layouts and dispatch modes
• add deterministic transpose/layout coverage and version the expanded backend dispatch-table ABI

This is independent of #769. Both PRs expand the dispatch table, so whichever merges second will need its shared ABI-version changes rebased.

Test plan

• H100 (CUDA 13.2): cuBLAS CT and RT bfloat16-output buffer + USM tests, 4/4 each
• MI210 (HIP 7.1): rocBLAS CT and RT bfloat16-output buffer + USM tests, 4/4 each (CT repeated 3x without contention)
• H100 and MI210: broader ordinary GemmTestSuite + GemmUsmTestSuite, CT and RT
• HeCBench PR 323 blas-gemmEx-sycl oneMath path on H100 and MI210: BF16 output PASS at 128x128x128 (out-of-scope int8 calls locally guarded for validation)
• git diff --check

Toolchain: DPC++ clang 23.0.0git (fbf7d1fd2cbf), CUDA toolkit 12.5 / CUDA runtime 13.2, ROCm/HIP 7.1.

Made with Cursor

[melonakos]

Reviewed alongside #769, since the two are structurally the same change to different routines. The ABI handling here is done correctly, same as there — but the interaction between these two PRs is the thing that matters most, and it can fail silently.

Critical: this and #769 both claim SPEC_VERSION 2

Both PRs insert new function pointers into the middle of function_table_t, and both bump SPEC_VERSION from 1 to 2 plus the same 28 WRAPPER_VERSION files.

Mid-struct insertion is safe on its own, because the loader rejects a mismatched backend:

if (t->version != SPEC_VERSION)
    throw math::specification_mismatch();

But if both land at version 2, there are two different table layouts both asserting the same version. The check then passes while the offsets disagree, and a backend built from one tree loaded against the other calls the wrong function pointer — the exact corruption the guard exists to catch, now with no diagnostic whatsoever.

The textual merge conflict will show up in function_table.hpp and on the version line, so it should be visible — but it would be very easy to resolve that conflict and leave both PRs at 2, which is the failure case.

Whichever of these merges second needs SPEC_VERSION 3 and another pass over all 28 wrapper files. My suggestion: land #769 first (the dot surface is simpler), then rebase this one on top and take it to 3. I've flagged the same thing there. Please decide the order deliberately rather than letting merge timing decide it.

What's right here

The mechanics match #769 and I checked them the same way: SPEC_VERSION bumped, all 28 real backend wrapper files bumped, and the only two files mentioning WRAPPER_VERSION that you left alone are the backend-authoring templates carrying the <Wrapper version number> placeholder. Nothing missed.

As with #769, most of the 56-file diff is one-line wrapper regeneration rather than 56 files of logic — the actual review surface is function_table.hpp, blas.hxx, blas_loader.cpp, and the per-backend level3 implementations. Worth saying for whoever reviews this second: don't be put off by the file count.

The bfloat16 output question

For dot in #769 the decisive detail was that cublasDotEx runs with CUDA_R_32F as the execution type, so accumulation happens in fp32 and only the final store narrows. The same question applies here and matters more, because GEMM accumulates over k rather than over a single reduction: the accumulate type must stay fp32 with only the C output written as bfloat16. Please confirm that's what the cuBLAS and generic paths do, and add a short comment stating it — for a bfloat16-output GEMM that's the entire numerical contract, and it shouldn't require a reader to decode a positional argument to find it.

Also worth confirming consistency across the CPU backends you touched (armpl_level3.cxx, netlib_level3.cxx, openblas_level3.cxx, mkl_level3.cxx, generic_level3_bfloat16.cpp). If any of them accumulates in bfloat16 rather than float, the same call will give materially different answers depending on which backend is loaded, and that's the kind of discrepancy that takes a long time to track down from a user's bug report.

Same maintainer question as #769

This extends the public API in include/oneapi/math/blas.hxx with bfloat16-output GEMM overloads. Whether that's oneMath-spec'd or a vendor extension is a governance decision, not a review one.

@sknepper — same flag as on #769. If these overloads aren't in the spec, deciding now whether they live on the main BLAS surface or behind an extensions header is far cheaper than revisiting it after release.

Test coverage

You've extended gemm.cpp and gemm_usm.cpp plus reference_blas_templates.hpp (+30) and test_common.hpp (+14), which is the right set of places. Given the tolerance work you did in #761 for int8→float, I'd expect bfloat16 output to need similar thought — bfloat16 has 8 mantissa bits, so the output rounding dominates and a plain relative check against a float reference will be loose or flaky depending on the values. Did you reuse the abs_bound mechanism from #761 here, or derive a separate tolerance? If the latter, it's worth explaining why in a comment, since the two PRs are solving closely related problems and a future reader will wonder why they differ.

[ndingle-arm]

Sorry, this occurred to me after I left my previous comments: OpenBLAS also has BGEMM and SBGEMM, so you could apply the same change here as you did for ArmPL and support both this and lines 679-690 in OpenBLAS as well. I don't know if it has a similar ifdef to ArmPL to allow you to check if BF16 is enabled in a certain build though.