cuBLAS

cuBLAS (CUDA Basic Linear Algebra Subroutines) is NVIDIA's high-performance implementation of the Basic Linear Algebra Subprograms (BLAS) standard. It is a proprietary library of highly optimized kernels for common linear algebra operations.

Instead of writing and tuning operations like matrix multiplication from scratch, you call cuBLAS functions from host code.

How it picks kernels#

cuBLAS contains many kernels, each fine-tuned for specific data types (e.g. FP32, FP16), matrix sizes and Streaming Multiprocessor (SM) architectures. At runtime, it uses internal heuristics (not publicly known) to select the fastest kernel and its optimal launch parameters.

This makes cuBLAS the foundation of most high-performance numerical computing on NVIDIA GPUs. Deep learning frameworks like PyTorch rely on it heavily for their core operations, alongside more specialized kernel libraries like cuDNN.

Matrix layout: the common pitfall#

The most common source of errors with cuBLAS is matrix data layout.

• Column-major order. For compatibility with the original BLAS standard (written in Fortran), cuBLAS expects matrices in column-major order. C, C++ and Python commonly use row-major.
• Leading dimension. Besides the operation's size (M, N, K), a BLAS function needs to know where each column starts in memory. The leading dimension (e.g. lda) is the stride between consecutive columns. For a whole allocated matrix, it's the number of rows. For a submatrix, it's the number of rows of the larger parent matrix.

Using row-major matrices without reordering#

For compute-heavy kernels like GEMM, you don't need to convert row-major matrices to column-major. Use the identity: if C = A × B, then Cᵀ = Bᵀ × Aᵀ.

Key idea A row-major matrix has exactly the same memory layout as its transpose stored column-major.

So pass your row-major A and B to cuBLAS in swapped order (and swap their dimensions). cuBLAS computes Cᵀ in column-major order, and that block of memory, read as row-major, is exactly the C you want:

#include <cublas_v2.h>

// performs single-precision C = alpha * A @ B + beta * C
// on row-major matrices using cublasSgemm
void sgemm_row_major(cublasHandle_t handle, int M, int N, int K,
                     const float *alpha,
                     const float *A, const float *B,
                     const float *beta,
                     float *C) {

  // A is M x K (row-major), cuBLAS sees it as A^T (K x M, column-major),
  //   the leading dimension of A^T is K
  // B is K x N (row-major), cuBLAS sees it as B^T (N x K, column-major),
  //   the leading dimension of B^T is N
  // C is M x N (row-major), cuBLAS sees it as C^T (N x M, column-major),
  //   the leading dimension of C^T is N

  // note the swapped A and B, and the swapped M and N
  cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
              N, M, K,
              alpha,
              B, N,  // leading dimension of B^T
              A, K,  // leading dimension of A^T
              beta,
              C, N); // leading dimension of C^T
}

The CUBLAS_OP_N flag tells the kernel to use the matrices as provided, with no additional transpose from its perspective. A complete, runnable version is on Godbolt.

Building with cuBLAS#

Functions are exposed via the cublas_v2.h header, and you must link the library (e.g. with -lcublas when compiling with nvcc).

See the official cuBLAS documentation.

Related#

• cuDNN: NVIDIA's deep learning kernel library
• CUTLASS: open-source templates for writing your own GEMM kernels
• Kernel: what cuBLAS functions launch