Matrix Multiplication

1. Basic Syntax

Tensors are created and filled on the GPU (Rand, RandN, Full, Zeros). Host memory is touched only to read results back — for validation or inspection — never to upload input data. Execution happens at the OaContext::Scope boundary.

// Inputs generated ON THE GPU — no host upload.
OaMatrix a = OaFnMatrix::Rand(OaShape2D(M, K));
OaMatrix b = OaFnMatrix::Rand(OaShape2D(N, K));

OaMatrix c;
{
    OaContext::Scope scope(OaContext::GetDefault());
    c = OaFnMatrix::MatMul(a, b);  // C = A @ B^T, shape [M, N]
}                                  // scope exit: Execute() + Sync()

// Read result back to host (validation only).
std::vector<float> host(M * N);
OaFnMatrix::CopyToHost(c, host.data(), host.size() * sizeof(float));

The shape convention matches current OA GEMM: A is [M, K],B is stored as [N, K], result is [M, N]. The GEMM router, tuner, and cache are internal — the user never calls them directly.

2. Correctness Pass

Start with small deterministic shapes and compare every output element against a CPU reference. Use normalized error to avoid false failures near zero:

for (OaU32 m = 0; m < M; ++m) {
    for (OaU32 n = 0; n < N; ++n) {
        float sum = 0.0f;
        for (OaU32 k = 0; k < K; ++k) {
            sum += aHost[m * K + k] * bHost[n * K + k];
        }
        cRef[m * N + n] = sum;
    }
}
// Normalized error (handles near-zero output entries)
float norm_err = max_i|ref[i] - gpu[i]| / max(max_i|ref[i]|, 1e-6f);

Measured Correctness — RTX 5090 Laptop (BF16 CoopMat, tol = 3e-2)

All shapes land at norm_err ≈ 5e-3, consistent with BF16's 7-bit mantissa — well inside the 3e-2 tolerance. On FP32-only devices the router picks tiled/naive kernels and norm_err drops below 1e-4.

3. Performance Pass

After correctness, run timed shapes that exercise the router. Each shape runs 50 iterations after warmup. GFLOP/s = 2·M·N·K / wall_time.

Single-Dispatch Throughput — RTX 5090 Laptop (BF16 CoopMat, 50 iters)

Device theoretical peak: 64 TFLOPS (BF16 tensor cores). Throughput climbs with problem size as fixed per-dispatch overhead amortizes, peaking near 19.2 TFLOP/s at 2048³. The gemv-decode row is memory-bound — a single output row can't saturate tensor cores. This is the autoregressive single-token path.

4. Batch Dispatch

Recording multiple operations into a single OaContext scope before executing reduces CPU→GPU submission overhead. For small operations where GPU time ≈ submit latency, batching provides 4–8× speedup.

// Batch dispatch: 4 operations, single submit
OaVec<OaMatrix> results;
for (int i = 0; i < 4; ++i) {
    results.PushBack(OaFnMatrix::MatMul(aVec[i], bVec[i]));
}
ctx.Execute();  // All 4 matmuls execute in one submission
ctx.Sync();

5. Cross-Stack Comparison

All measurements are wall-clock time including full public API path (submit + sync overhead).

OA BF16 CoopMat achieves 19.2 TFLOP/s at 2048³ (30.0% of 64 TFLOP theoretical peak). The reference TF32 implementation achieves 38.7 TFLOP/s (60% of peak). There is room for improvement through tile tuning and shared memory optimization — the kernel path uses the same tensor core hardware. Intel iGPU runs the same binary at 259 GFLOP/s (FP32 fallback, 10.4% of 2.5 TFLOP theoretical peak).

6. Cross-Device Portability

Same binary, same source — runtime device selection:

cmake --build Build/Release --target TutorialCoreMatMulIntro -j
./Bin/Release/Tutorial/Core/TutorialCoreMatMulIntro

# Intel iGPU
OA_DEVICE=integrated ./Bin/Release/Tutorial/Core/TutorialCoreMatMulIntro

See also: detailed benchmark results with cross-stack comparison and batch-dispatch tables.