blas-microbenchmark measures the routines of the BLAS. This page says what those are, enough to read a result: what a routine computes, what its letters mean, and why some reach a hundred GFLOP/s where others stay at a few.

An interface, many libraries

The BLAS (Basic Linear Algebra Subprograms) is a standard set of vector and matrix operations, fixed in the 1970s and 1980s: add a multiple of one vector to another, multiply a matrix by a vector, multiply two matrices, solve a triangular system. Every numerical code that does linear algebra calls them, directly or through LAPACK, NumPy, Julia or a solver, and does not care who wrote them.

So many libraries implement the same interface, each tuned its own way: OpenBLAS, BLIS and AMD's AOCL-BLAS built on it, Intel's MKL, NVIDIA's NVPL, Arm's ArmPL, and Netlib's reference implementation, written for clarity rather than speed. The same program runs against any of them; how fast it runs depends on which. That is what these benchmarks measure (The BLAS libraries).

The benchmarks call the BLAS through CBLAS, its C interface: cblas_dgemm(...), with the storage order of the matrices as an argument (row-major, as C stores them, or column-major, as Fortran does; --layout).

Three levels

The routines come in three levels, by what they work on. The level decides what limits them: the memory, or the arithmetic.

Level	Works on	Example	FLOPs a byte (double)	Limited by
1	vectors	daxpy: y ← αx + y	2n for 24n: 1/12	memory bandwidth
2	a matrix and vectors	dgemv: y ← αAx + βy	2n² for 8n²: 1/4	memory bandwidth
3	matrices	dgemm: C ← αAB + βC	2n³ for 24n²: n/12	arithmetic, once n is large

A level 1 or 2 routine reads each number once and does one or two operations with it: it waits for the memory, and its GB/s says more than its GFLOP/s. A level 3 routine uses each number n times: a good library keeps the operands in the caches and the vector units busy, and approaches the processor's peak. That is why the report draws a roofline for levels 1 and 2, and the peak for level 3 (The peak).

bmb_daxpy -v 1024:16777216
bmb_dgemm -m 256:2048

Run both: daxpy's rate falls as the vectors leave each cache, dgemm's rises with the size.

Four precisions

A routine's first letter is its precision:

Letter	Precision	Bytes an element
s	single, 32-bit real	4
d	double, 64-bit real	8
c	single complex	8
z	double complex	16

One routine mixes them: sbgemm, a GEMM whose A and B are in bfloat16 (2 bytes, 8 bits of significand) and whose C, and the sums, are in single precision, as machine-learning codes call it.

• Single moves half the bytes of double and fits twice as many numbers in a vector register: it can run about twice as fast, at about 7 significant digits rather than 16.
• Complex numbers are pairs of reals. A complex multiply-add is four multiplies and four adds, eight FLOPs on two numbers, where a real one is two FLOPs on one: four times the arithmetic on twice the data. A complex level 3 routine is therefore more compute-bound than its real counterpart.
• Two letters name a routine that takes one precision and returns another: scasum sums a single complex vector into a single real, dznrm2 takes the norm of a double complex vector. i?amax returns an index: idamax.

The FLOPs and bytes the benchmarks count follow these conventions, so that a GFLOP/s is comparable across precisions and libraries (How a number is measured).

The routines

A routine is a precision and an operation, a family: dgemm is gemm in double. α and β are scalars, x and y vectors, A, B and C matrices; Aᵀ is A transposed, Aᴴ its conjugate transpose.

Level 1: vectors

Family	Computes
axpy	y ← αx + y
scal	x ← αx
copy, swap	y ← x; x ↔ y (no arithmetic: only GB/s)
dot	xᵀy (real); dotu xᵀy and dotc xᴴy (complex)
nrm2	‖x‖₂, the Euclidean norm, without overflow
asum	Σ|xᵢ| (for a complex x, |re| + |im|)
iamax	the index of the largest |xᵢ| (no arithmetic counted)
rot	applies a plane rotation to the pairs (xᵢ, yᵢ)
axpby	y ← αx + βy: an extension, which not every library has

Level 2: a matrix and vectors

Family	Computes
gemv	y ← αAx + βy, or with Aᵀ or Aᴴ (--trans)
gbmv	the same with A a band matrix: only K diagonals each side of the main one are stored (-M is K)
symv, hemv	y ← αAx + βy with A symmetric (real) or Hermitian (complex), one triangle stored (--uplo); spmv with that triangle packed, column after column, without the other half
ger	A ← αxyᵀ + A, a rank-1 update; geru and gerc (xyᴴ) in complex
syr, her	A ← αxxᵀ + A (xxᴴ), on one triangle; syr2, her2: αxyᵀ + αyxᵀ + A
trmv	x ← Ax with A triangular (--uplo, --trans, --diag: a unit diagonal, not stored)
trsv	solves Ax = b for x, in place, A triangular; tpsv with A packed, tbsv with A banded

Level 3: matrices

Family	Computes
gemm	C ← αAB + βC, each of A and B transposed or not (--trans NT, …). The routine most of linear algebra spends its time in. -K sets the inner dimension apart; --batch-count does many at once
symm, hemm	C ← αAB + βC with A symmetric or Hermitian, on the left or the right (--side)
syrk, herk	C ← αAAᵀ + βC (AAᴴ), on one triangle of C: a rank-k update
syr2k, her2k	C ← αABᵀ + αBAᵀ + βC, the rank-2k update
trmm	B ← αAB with A triangular, on either side
trsm	solves AX = αB for X, A triangular: the core of LU and Cholesky solves
gemmt	C ← αAB + βC on one triangle of C (--uplo), for a product known to be symmetric: an extension
sbgemm	gemm with A and B in bfloat16, C in single precision: an extension

The extensions are calls libraries add to the BLAS. Their benchmarks are built whatever the library; where it lacks the call, they measure nothing and say which call is missing.

Not every family exists in every precision: dot, ger, symv, syr, syr2 and spmv are real, their complex counterparts are dotu and dotc, geru and gerc, hemv, her, her2; and here rot, gbmv, tbsv and tpsv are measured in real precision only. The manual lists the 113 benchmarks by level, with what -m, -M and -K are for each, and Options the arguments each takes.

Matrix structures

BLAS has no option for the kind of matrix: the structure is the routine's. A general matrix goes to gemm or gemv; the others each have their own routines, which read only the part that matters:

Matrix	Routines	What is used
Triangular	trmv, trsv, trmm, trsm	one triangle (--uplo), its diagonal taken as ones with --diag U
Symmetric, Hermitian	symv, symm, syr, syr2, syrk, syr2k; hemv, hemm, her, her2, herk, her2k	one triangle of the symmetric or Hermitian matrix (--uplo), the other being its mirror, conjugated for a Hermitian one; for syr, syr2, her and her2 that matrix is A, which they update; for syrk, syr2k, herk and her2k it is C, the triangle written
Band	gbmv, tbsv	the diagonals near the main one, in band storage: K each side for gbmv, K on one side for tbsv, which takes --uplo and --diag too (-M is K, the band's half-width)
Packed	spmv, tpsv	one triangle (--uplo; --diag for tpsv), stored without the other's place

rot applies a plane (Givens) rotation to two vectors: the rotation is given by its cosine and sine, not as a matrix. Sparse matrices, mostly zeros stored as lists of their non-zeros, are not in the dense BLAS these libraries implement: they have routines and libraries of their own (the Sparse BLAS, MKL's sparse routines, cuSPARSE, PETSc), and no benchmark here takes one.

For most dense calls the values barely change the speed: the same operations run in the same order whatever the operands hold. A library may still take shortcuts on some values, as the Reference BLAS does on a zero element, and subnormals, NaN and infinities slow some processors down. The benchmarks fill every operand themselves, the same for every library, with none of those (The operands).