Hands-On BLAS Demo

4096³ GEMM on NVIDIA + Apple — what the numbers look like

The cleanest way to feel BLAS's value is to run the same logical operation on both ecosystems and read the throughput numbers side by side.

Level	RTX 4090 (FP32)	M3 Ultra (FP32)	M3 Ultra (FP16)
BLAS-1 (AXPY)	20 GF/s	5 GF/s	—
BLAS-2 (GEMV)	0.44 GF/s	279 GF/s (GPU)	1,136 GF/s
BLAS-3 (GEMM)	2.32 TF/s (non-Tensor-Core)	13.89 TF/s	18.15 TF/s

Three things to read from this table:

BLAS-1 is awful on GPU. You see double-digit GF/s on hardware capable of TF/s — bandwidth-bound, no reuse, no chance. cuBLAS detects this and runs it on the host or skips Tensor Cores.
The CUDA FP32 GEMM here uses the non-tensor path for a fair comparison. Enable Tensor Core FP16 via cublasGemmEx with CUBLAS_GEMM_DEFAULT_TENSOR_OP and the 4090 jumps past 30 TF/s.
M3 Ultra GEMM crushes because the matrix coprocessor is fed by huge unified memory bandwidth. Tested at 4096³.

Intensity trends are platform-agnostic: BLAS-1 → bandwidth-bound, BLAS-3 → compute-bound. The absolute numbers vary; the shape of the picture doesn't.

Code

Benchmark harness — same shape, two backends, fair compare·python

import time
import numpy as np
import torch

N = 4096
A_cpu = np.random.randn(N, N).astype(np.float32)
B_cpu = np.random.randn(N, N).astype(np.float32)

def bench(label, A, B, runs=10, sync=lambda: None):
    sync(); start = time.perf_counter()
    for _ in range(runs):
        C = A @ B
        sync()
    elapsed = (time.perf_counter() - start) / runs
    flops = 2 * N**3
    print(f'{label}: {elapsed*1000:.2f} ms, {flops/elapsed/1e12:.2f} TFLOP/s')

# numpy (CPU via OpenBLAS / MKL)
bench('numpy CPU', A_cpu, B_cpu)

# PyTorch on best available accelerator
if torch.cuda.is_available():
    A = torch.from_numpy(A_cpu).cuda()
    B = torch.from_numpy(B_cpu).cuda()
    bench('torch CUDA FP32', A, B, sync=torch.cuda.synchronize)
    bench('torch CUDA FP16', A.half(), B.half(), sync=torch.cuda.synchronize)
elif torch.backends.mps.is_available():
    A = torch.from_numpy(A_cpu).to('mps')
    B = torch.from_numpy(B_cpu).to('mps')
    bench('torch MPS FP32', A, B, sync=torch.mps.synchronize)
    bench('torch MPS FP16', A.half(), B.half(), sync=torch.mps.synchronize)

Exercise

Run the benchmark above on whichever box you have. Capture results in a small CSV (label, time_ms, tflops). The number you get is your machine's BLAS reference. Every hand-rolled kernel you'll write in the rest of this quest gets compared to it. If you have access to both NVIDIA and Apple boxes, run on both — the asymmetry between the two ecosystems is part of the lesson.

4096³ GEMM on NVIDIA + Apple — what the numbers look like

Code

External links

Exercise

Progress

Comments 0