cuBLAS, in depth: column-major GEMMs from row-major frameworks, handles and workspace, cuBLASLt heuristics and epilogues, FP8 rules and reproducibility

The cuBLAS family

cuBLAS ships several APIs. The cuBLAS API is the classic BLAS interface with Level 1, 2 and 3 routines, plus extensions such as cublasGemmEx, cublasGemmStridedBatchedEx and cublasGemmGroupedBatchedEx for many differently sized GEMMs in one call. The cuBLASLt API handles only matmul but exposes far more control: opaque descriptors for the operation and each matrix, a heuristic query that returns ranked algorithms, fused epilogues such as bias, ReLU and GELU, and the scaling machinery that FP8 and FP4 need. The cuBLASXt API spreads host-resident BLAS calls across several GPUs and is rarely used in training. Deep learning frameworks use cuBLASLt for most work today; PyTorch exposes torch.backends.cuda.preferred_blas_library() to choose between backends.

Column-major, and how row-major frameworks use it

cuBLAS follows Fortran BLAS and treats every matrix as column-major: element (i, j) lives at ptr[i + j * ld], where the leading dimension ld is the stride between columns. PyTorch stores tensors row-major. The two are reconciled without copying by one identity: a row-major M x N matrix occupies exactly the same bytes as a column-major N x M matrix, its transpose. So to compute a row-major C = A B, ask cuBLAS for the column-major product C^T = B^T A^T, passing B first and swapping M and N.

A Linear layer computes y = x W^T + b, with x of shape [tokens, in] and W of shape [out, in], both row-major. In column-major terms the call is m = out, n = tokens, k = in, with A = W transposed and B = x not transposed. That is the "TN" layout, which matters again for FP8. In code:

// y[T, out] = x[T, in] * W[out, in]^T, all row-major BF16, FP32 accumulation.
// Column-major view: y^T[out, T] = W[in, out]^T * x^T[in, T]
const float alpha = 1.0f, beta = 0.0f;
cublasGemmEx(handle,
             CUBLAS_OP_T, CUBLAS_OP_N,
             out, T, in,                       // m, n, k
             &alpha,
             W, CUDA_R_16BF, in,               // A, lda = in (row stride of W)
             x, CUDA_R_16BF, in,               // B, ldb = in
             &beta,
             y, CUDA_R_16BF, out,              // C, ldc = out
             CUBLAS_COMPUTE_32F,
             CUBLAS_GEMM_DEFAULT);

Get the leading dimensions wrong and cuBLAS either returns CUBLAS_STATUS_INVALID_VALUE or, worse, silently reads a strided view of the wrong elements. Unit-test every new GEMM wrapper against a reference on non-square shapes, where a swapped m and n cannot hide.

Handles, streams and workspace

A handle carries the library context: the stream, the math mode, the pointer mode and the workspace. Creating one is expensive, so create one per thread and device and reuse it. Workspace is scratch memory some kernels need, for example to reduce split-K partial sums. The documentation recommends 32 MiB on Hopper and newer and 4 MiB on other architectures; a user workspace set with cublasSetWorkspace must be aligned to 256 bytes, and a size of at least 16 KiB avoids CUBLAS_STATUS_ALLOC_FAILED. Two details bite in practice. First, cublasSetStream unconditionally resets the workspace to the default pool, so set the stream first and the workspace second. Second, a user workspace shared between threads is unsafe because one call may run several kernels that rely on the workspace staying intact; give each stream its own handle and workspace.

For cuBLASLt the workspace is passed per call, and the size you advertise through CUBLASLT_MATMUL_PREF_MAX_WORKSPACE_BYTES limits which algorithms the heuristic may return. Advertising zero workspace quietly removes the fastest split-K kernels for small-M shapes.

Compute types, TF32 and emulation

Two settings decide the arithmetic: the data types of the matrices and the compute type, which sets accumulation precision and what the library may substitute. CUBLAS_COMPUTE_32F with BF16 or FP16 inputs is the training default: low-precision inputs, FP32 accumulation. With FP32 inputs, CUBLAS_COMPUTE_32F_FAST_TF32 lets tensor cores round inputs to TF32, and CUBLAS_COMPUTE_32F_FAST_16BF allows down-conversion to BF16. Recent releases add emulation: CUBLAS_COMPUTE_32F_EMULATED_16BFX9 computes FP32 GEMMs from BF16 tensor-core products, and a fixed-point scheme emulates FP64. Both can also be enabled by math mode or environment variables such as CUBLAS_EMULATE_SINGLE_PRECISION, which means a deployment can change numerics without a code change; pin these variables in your launch environment.

Frameworks map their switches onto these. In PyTorch, torch.set_float32_matmul_precision("high") or torch.backends.cuda.matmul.allow_tf32 = True lets FP32 matmuls use TF32. That is usually right for training and wrong for numerical code that expects true FP32. Mixed Precision Training covers loss scaling and master weights around these choices.

cuBLASLt: descriptors, heuristics and epilogues

cuBLASLt separates describing a problem from choosing an algorithm. You build an operation descriptor, a layout for each matrix and a preference object, ask the heuristic for ranked candidates, then execute. The heuristic costs host time, tens of microseconds according to the documentation, so query once per shape and cache the result. The example fuses bias and GELU into the epilogue, which saves a full read and write of the output compared with separate kernels:

cublasLtMatmulDesc_t op;
cublasLtMatmulDescCreate(&op, CUBLAS_COMPUTE_32F, CUDA_R_32F);
cublasOperation_t ta = CUBLAS_OP_T, tb = CUBLAS_OP_N;
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_TRANSA, &ta, sizeof(ta));
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_TRANSB, &tb, sizeof(tb));
cublasLtEpilogue_t epi = CUBLASLT_EPILOGUE_GELU_BIAS;
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_EPILOGUE, &epi, sizeof(epi));
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_BIAS_POINTER, &bias, sizeof(bias));

cublasLtMatrixLayout_t La, Lb, Ld;               // column-major shapes as stored
cublasLtMatrixLayoutCreate(&La, CUDA_R_16BF, in, out, in);   // W
cublasLtMatrixLayoutCreate(&Lb, CUDA_R_16BF, in, T, in);     // x
cublasLtMatrixLayoutCreate(&Ld, CUDA_R_16BF, out, T, out);   // y

cublasLtMatmulPreference_t pref;
cublasLtMatmulPreferenceCreate(&pref);
size_t ws_bytes = 32u << 20;
cublasLtMatmulPreferenceSetAttribute(pref, CUBLASLT_MATMUL_PREF_MAX_WORKSPACE_BYTES,
                                     &ws_bytes, sizeof(ws_bytes));
cublasLtMatmulHeuristicResult_t res[8]; int n = 0;
cublasLtMatmulAlgoGetHeuristic(lt, op, La, Lb, Ld, Ld, pref, 8, res, &n);
if (n == 0) { /* unsupported combination: fall back or fix alignment */ }

cublasLtMatmul(lt, op, &alpha, W, La, x, Lb, &beta, y, Ld, y, Ld,
               &res[0].algo, workspace, ws_bytes, stream);

The bias vector has one entry per row of D, which here is one per output feature, exactly what a Linear layer needs. Other epilogues produce the auxiliary output the backward pass needs, such as CUBLASLT_EPILOGUE_GELU_AUX_BIAS paired with CUBLASLT_EPILOGUE_DGELU_BGRAD in the backward GEMM. If you overlap GEMMs with communication kernels, CUBLASLT_MATMUL_DESC_SM_COUNT_TARGET tells the heuristic how many SMs to plan for, so it does not pick a kernel sized for a GPU it will not get.

FP8 and block scaling

FP8 has a narrow range, so every FP8 operand carries a scale. The documentation lists the requirements for tensor- or block-scaled FP8 kernels: pointers and dimensions must support 16-byte alignment, the compute type must be CUBLAS_COMPUTE_32F, the scale type must be CUDA_R_32F, and on compute capability 8.9, 9.0 and 12.x A must be transposed and B not, the TN layout a Linear layer already uses. Scales are passed with CUBLASLT_MATMUL_DESC_A_SCALE_POINTER and its B, C and D siblings, and CUBLASLT_MATMUL_DESC_AMAX_D_POINTER returns the output's absolute maximum so the next step's scale can be computed. The scale mode selects the granularity:

Which modes a GPU supports differs by architecture; check the scaling-mode support table for yours, and treat a zero-result heuristic query as the signal that the combination is unsupported.

Worked example: one MLP layer, prefill versus decode

Take the up-projection of a transformer MLP with in = 4,096 and out = 14,336 in BF16. During prefill of 8,192 tokens the GEMM does 2 x 8,192 x 14,336 x 4,096, about 0.96 TFLOP, and moves the two inputs and the output once: about 420 MB. That is roughly 2,300 FLOPs per byte, far above the ratio of peak compute to memory bandwidth, so it is compute-bound and its time is set by tensor-core throughput. At an achieved 700 TFLOP/s it takes about 1.4 ms.

During decode with a batch of 8 sequences the same layer does 2 x 8 x 14,336 x 4,096, about 0.94 GFLOP, but must still read all 117 MB of weights. Intensity falls to about 8 FLOPs per byte and the GEMM is memory-bound: at 3.35 TB/s of HBM bandwidth it takes about 35 microseconds no matter how fast the tensor cores are. Different shapes, different winning kernels: for decode, the heuristic's split-K choices and your workspace size matter, and batching more sequences is nearly free until intensity crosses the machine balance point. For mixture-of-experts layers, where every expert sees a different M, see Grouped GEMM for MoE.

Reproducibility and logging

cuBLAS guarantees bit-wise identical results from run to run only for the same toolkit version on GPUs of the same architecture and SM count, and only with a single stream. With several concurrent streams, give each its own handle and workspace and set CUBLAS_WORKSPACE_CONFIG to :16:8 (may cost performance) or :4096:8 (about 24 MiB more GPU memory). PyTorch's torch.use_deterministic_algorithms(True) requires that variable for this reason. Routines such as symv lose reproducibility if you enable their faster atomics path with cublasSetAtomicsMode. A newer, experimental option, CUBLASLT_MATMUL_DESC_BATCH_INVARIANCE_FLAGS, makes results independent of how a problem is split along M and N, which is what an inference server needs when batch composition changes; it is documented for compute capability 10.x only.

To see what the library actually did, turn on logging. CUBLASLT_LOG_LEVEL ranges from 0 (off) through 1 (errors), 2 (trace of kernel-launching calls), 3 (performance hints) and 4 (info, including heuristic status) to 5 (full API trace); CUBLASLT_LOG_FILE redirects it. The classic API uses CUBLAS_LOGINFO_DBG and CUBLAS_LOGDEST_DBG. Level 3 hints are the cheapest performance review you will ever get.

Failure modes

• CUBLAS_STATUS_NOT_SUPPORTED after a dtype change. Usually an FP8 rule: alignment, layout or compute type. Log at level 4 to see why the heuristic returned nothing.
• Slow GEMMs on odd shapes. Dimensions that are not multiples of 8 or 16 elements lose 16-byte alignment and the fast kernels. Pad vocabulary and hidden sizes.
• Workspace silently reset. Calling cublasSetStream after cublasSetWorkspace drops the user workspace.
• Corruption with shared workspace across threads. Interleaved kernels clobber each other's scratch. One handle and workspace per stream.
• Numerics changed by the environment. TF32 or emulation enabled by a framework default or variable. Record matmul precision settings with every run.
• Run-to-run drift. Multiple streams without a workspace config, or a different SM count after a MIG change. Pin both when you need determinism.

Trade-offs

cuBLAS and cuBLASLt give you NVIDIA-tuned kernels for every architecture on day one with no compile step, at the cost of a closed set of fusions: what the epilogue enum offers is what you get. When you need a fusion it lacks, CUTLASS or Triton is the next step; CUTLASS in depth compares the three. For convolution and attention graphs the analogous library is cuDNN.

What to do next

1. Write a GEMM wrapper test that checks non-square shapes against a reference, so layout and leading-dimension mistakes fail loudly.
2. Create one handle per stream; set the stream, then a 256-byte-aligned workspace of the recommended size.
3. Decide TF32 and emulation policy explicitly and record it in run metadata.
4. Cache cuBLASLt heuristic results per shape instead of querying every call.
5. Fuse bias and activation into epilogues where the enum offers them.
6. Run once with CUBLASLT_LOG_LEVEL=3 and act on the hints.
7. Classify your hot GEMMs as compute- or memory-bound with the intensity arithmetic above before tuning anything.
8. If you need reproducibility, set CUBLAS_WORKSPACE_CONFIG and pin the toolkit version and GPU configuration.