cuBLAS & MPS Cheat Sheets

Two libraries, one API shape

cuBLAS and MPSMatrixMultiplication are different products from different companies but they implement essentially the same operation. The differences worth memorizing:

Aspect	cuBLAS (NVIDIA)	MPS (Apple)
Layout default	Column-major (Fortran)	Row-major
Default precision focus	FP32, with FP16/BF16 via cublasGemmEx	Half-precision-first; FP32 also supported
Lifecycle	Handle per thread, async by default	Command buffer per dispatch, sync at `waitUntilCompleted`
α/β scalars	Yes (full GEMM signature)	Yes
Algorithm picker	Heuristic via cublasLt; manual via algo IDs	Library picks; less direct control

Rules of the road for cuBLAS

Create one cublasHandle_t per thread (or per CUDA stream); reuse it across calls.
Async by default — call cudaStreamSynchronize if you need the result before the next CPU code.
Column-major math means most C/C++ users compute Cᵀ = Bᵀ · Aᵀ via CUBLAS_OP_T rather than transposing buffers.

Rules of the road for MPS

Row-major by default; same direction your numpy / Swift arrays naturally store.
Command buffer lifecycle is explicit: encode, commit, wait. The waiting step is non-optional.
Half-precision MPSMatrix is a first-class citizen — for inference, prefer FP16 unless numerical stability demands FP32.

Code

cuBLAS sgemm — column-major dance·cuda

#include <cublas_v2.h>

cublasHandle_t handle;
cublasCreate(&handle);

const float alpha = 1.f, beta = 0.f;
// Logical: C = A·B, where A is (m,k), B is (k,n), C is (m,n), all row-major.
// cuBLAS column-major: compute Cᵀ = Bᵀ · Aᵀ → swap operand order, use OP_T:
cublasSgemm(handle,
            CUBLAS_OP_T, CUBLAS_OP_T,    // op(A) = Aᵀ, op(B) = Bᵀ
            n, m, k,                     // result is (n × m) in column-major = (m × n) row-major
            &alpha,
            dB, n,                       // ldb = original cols of B (= n)
            dA, k,                       // lda = original cols of A (= k)
            &beta,
            dC, n);

cudaDeviceSynchronize();
cublasDestroy(handle);

MPSMatrixMultiplication — Swift, row-major, half precision·swift

import Metal
import MetalPerformanceShaders

let device = MTLCreateSystemDefaultDevice()!
let queue = device.makeCommandQueue()!

// Wrap your existing MTLBuffers as MPSMatrix descriptors.
let descA = MPSMatrixDescriptor(rows: m, columns: k,
    rowBytes: k * MemoryLayout<Float16>.stride, dataType: .float16)
let descB = MPSMatrixDescriptor(rows: k, columns: n,
    rowBytes: n * MemoryLayout<Float16>.stride, dataType: .float16)
let descC = MPSMatrixDescriptor(rows: m, columns: n,
    rowBytes: n * MemoryLayout<Float16>.stride, dataType: .float16)

let A = MPSMatrix(buffer: bufA, descriptor: descA)
let B = MPSMatrix(buffer: bufB, descriptor: descB)
let C = MPSMatrix(buffer: bufC, descriptor: descC)

let mm = MPSMatrixMultiplication(device: device,
    transposeLeft: false, transposeRight: false,
    resultRows: m, resultColumns: n,
    interiorColumns: k,
    alpha: 1.0, beta: 0.0)

let cb = queue.makeCommandBuffer()!
mm.encode(commandBuffer: cb, leftMatrix: A, rightMatrix: B, resultMatrix: C)
cb.commit()
cb.waitUntilCompleted()

Exercise

Pick whichever side you have access to and write a 'hello cuBLAS' or 'hello MPS' GEMM that multiplies two 1024×1024 random matrices. Compare the result to numpy.dot output element-wise within a 1e-5 tolerance. Time the call (excluding data transfer). On a 4090 you should see ~3–7 TF/s FP32, on an M3 Ultra ~13 TF/s FP32 / ~18 TF/s FP16. Note the numbers; you'll compare your hand-rolled tiled GEMM against these in Track 8.