What is cuBLAS?

cuBLAS (CUDA Basic Linear Algebra Subroutines) is NVIDIA's proprietary, closed-source implementation of the classic Basic Linear Algebra Subprograms (BLAS) standard for GPUs. Rather than hand-writing and tuning routines like matrix multiplication yourself, you link against cuBLAS and call its functions from host code. The library ships a large catalog of kernels, each tuned for a specific data type (FP32, FP16, and so on), matrix size, and Streaming Multiprocessor architecture, and dispatches to whichever kernel actually drives a chip's Tensor Cores hardest for that shape, via internal heuristics NVIDIA has not published. That makes cuBLAS the default foundation under most high-performance numerical code on NVIDIA GPUs, including the core matrix operations inside PyTorch and other deep learning frameworks, usually alongside more specialized libraries like cuDNN for convolutions and CUTLASS where a fixed cuBLAS routine does not fit.

The column-major trap

The single most common mistake when calling cuBLAS is memory layout. BLAS was originally a Fortran standard, and Fortran stores two-dimensional arrays in column-major order: values from the same column sit next to each other in memory. C, C++, and Python (via NumPy's default) all use row-major order instead, so a matrix your program has stored one way is silently reinterpreted the other way if you hand it to cuBLAS without accounting for the difference. cuBLAS also needs the leading dimension of each matrix argument, the memory stride between the start of one column and the next, which is just the row count for a full matrix but becomes the parent matrix's row count when you are passing a sub-matrix out of a larger allocation.

The practical fix is a linear algebra identity rather than a physical reordering of your data: if C = A @ B, then C^T = B^T @ A^T. A matrix stored in row-major order is bit-for-bit identical in memory to its own transpose stored in column-major order. So if you hand cuBLAS your row-major A and B but swap their order (and their M/N dimensions) in the function call, it computes C^T in what it thinks is column-major layout, and that same block of memory read back as row-major is exactly the C you wanted, with no data movement required. A minimal single-precision GEMM wrapper using this trick looks like:

#include <cublas_v2.h>

// C = alpha * A @ B + beta * C, with A, B, C all row-major
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) {
  // cuBLAS sees A (M x K row-major) as A^T (K x M col-major), ld = K
  // cuBLAS sees B (K x N row-major) as B^T (N x K col-major), ld = N
  // cuBLAS sees C (M x N row-major) as C^T (N x M col-major), ld = N
  // swap A/B and M/N to compute C^T instead of C
  cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
              N, M, K,
              alpha, B, N,
              A, K,
              beta, C, N);
}

CUBLAS_OP_N tells cuBLAS to use each matrix exactly as passed, with no extra transpose applied. This row/column swap is the standard way most C++ and Python codebases interoperate with cuBLAS without ever materializing a transposed copy of anything.

To use cuBLAS from CUDA C++, link the library (-lcublas at compile time) and include its cublas_v2.h header; NVIDIA's own cuBLAS documentation is the reference for the full function catalog. For workloads outside cuBLAS's fixed routine set, or where you need to fuse a GEMM with custom pre- or post-processing, CUTLASS is the usual next step, and you can watch kernel-level utilization on either with nvidia-smi while you benchmark. Aquanode's marketplace has H100 and H200 instances with cuBLAS and the rest of the CUDA toolkit preinstalled if you want to run these comparisons yourself.