Closing the cuBLAS gap

Profiling against cuBLAS identifies whether the remaining gap is memory, occupancy, or instruction issue.

Compare against a library you did not write

cuBLAS is column-major, so a row-major C = A * B is a column-major C^T = B^T * A^T, and the operands go in the other order:

cublasHandle_t h; cublasCreate(&h);
const float alpha = 1.f, beta = 0.f;
cublasSgemm(h, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K,
            &alpha, dB, N, dA, K, &beta, dC, N);

Swap the operands the wrong way and you get a plausible matrix, not an error message. Check both results against the same CPU reference before timing either.

One harness, interleaved

Same buffers, same shape, same warmup, alternating cuBLAS and your kernel, median with a spread. A ratio smaller than the spread is not a gap; it is noise with a decimal point.

Three places the gap can live

Memory: DRAM throughput at peak and sectors per request above four. Occupancy: achieved warps against theoretical, and the limiter the profiler names. Issue: scheduler slots used, pipe utilisation, and the stall with the largest share. Profile one kernel per run and compare the same section side by side; the section that is saturated in cuBLAS and slack in yours is the diagnosis. The rows look like this, with your values:

gpu__dram_throughput.avg.pct_of_peak_sustained_elapsed   <value>
smsp__warps_active.avg.pct_of_peak_sustained_active     <value>
launch__registers_per_thread                            <value>

Mode names and TF32 behaviour differ across CUDA 12 and 13; confirm against the cuBLAS guide for your toolkit before quoting one.

Things that break

On Ampere and later, cuBLAS can satisfy an FP32 GEMM with TF32 tensor-core math, which is faster and less accurate than the FP32 path. The default depends on the library version and the math mode, so set the mode explicitly and record which one you compared. Creating the handle or the first lazy initialisation inside the timed region inflates the library’s time. And clocks drift between a cuBLAS block and a kernel block, which is the reason the runs alternate.

Sources