The same subroutine everywhere
Open a modern AI model, a physics engine, or a ray tracer, and underneath the surface you will usually find the same small, viciously optimized routine: a general matrix multiply, or GEMM. It is the multiplication of one grid of numbers by another, scaled and accumulated, repeated billions of times a second.
That routine is not glamorous. It does not look like a character, a light bounce, or a neural network layer. But it is the workhorse. If you care about speed at the bottom of the stack, you eventually have to care about how a CPU or a GPU multiplies matrices.
This post is a high-level tour. No proprietary code, no internal kernels, just the engineering shape of the problem and the patterns that are universal across every fast implementation.
"The best way to understand a fast computer is to ask what it is moving, not just what it is computing."
A two-minute primer
GEMM is the BLAS name for a matrix-matrix multiply with scaling. In plain terms, it computes:
Matrix A is M × K, B is K × N, and C is M × N. The scalars alpha and beta let you scale and accumulate, which is more useful than it sounds: it means one GEMM can do a linear layer, an attention projection, a rigid-body transform update, or a batch of ray-tracing instance transforms without extra passes.
GEMV is the vector cousin: y = alpha × A × x + beta × y. One matrix, one vector. It is a Level-2 BLAS operation (matrix-vector), while GEMM is Level-3 (matrix-matrix). The difference matters because GEMM reuses each loaded number many times; GEMV loads the whole matrix but only uses each entry once, so it is almost always memory-bound.
In linear algebra, BLAS is the Basic Linear Algebra Subprograms library family (Netlib BLAS, OpenBLAS, MKL) that defines GEMM and GEMV. In ray tracing, BLAS is a Bottom-Level Acceleration Structure — a tree of triangles. They share an acronym and nothing else. We will return to the ray-tracing meaning later.
Why the world runs on matrix multiply
A neural-network forward pass is, structurally, a sequence of GEMMs separated by element-wise activations. A transformer layer can spend more than 99% of its floating-point operations inside matrix multiplications. CNNs and physics simulators are no different: the heavy lifting reduces to multiplying large, structured grids of numbers.
The reason GEMM is so dominant in performance engineering is arithmetic intensity — the ratio of math done to bytes moved. For an M × N × K GEMM, the total work is about 2 × M × N × K floating-point operations, while the data moved is only the three matrices. For a 1024×1024×1024 multiply in FP32, that is roughly 170 operations for every byte loaded. Few kernels can match that.
High arithmetic intensity means the processor has enough work to hide memory latency. The data arrives, gets used hundreds of times, and only then is evicted. That is the property that lets GEMM saturate compute units while most other code is stuck waiting for DRAM.
The memory mountain
Every optimized kernel begins by admitting a hardware fact: memory is a pyramid. Registers are fast and tiny. Caches are larger and slower. DRAM is enormous and, from the processor's perspective, glacial.
A naïve matrix multiply keeps reading the same rows and columns from DRAM over and over. For a large matrix, that is catastrophic. The cache cannot hold the working set, so the processor spends most of its time waiting. The solution is to stop treating the matrix as one object and start treating it as a grid of tiles.
Tiling: how to multiply a matrix that does not fit in cache
Imagine a 4096×4096 matrix. It is far too large for any cache. But a 64×64 sub-matrix is tiny. The core idea of cache tiling is to compute one small block of the output at a time, loading only the A and B tiles that contribute to that block, doing all the work on them, and writing the result back once.
The classic three-level blocking strategy, shaped by Goto and van de Geijn and used in libraries like BLIS, OpenBLAS, and oneDNN, looks something like this:
Packing is the unsung hero. Real matrices are rarely stored in the exact order the micro-kernel wants. By copying a panel into a contiguous buffer, the inner loop becomes a simple sequential read, which the hardware prefetcher and cache can exploit.
The microkernel: where the work actually happens
The micro-kernel is the innermost loop — where the actual work happens. It is written in vector intrinsics or assembly, and its only job is to multiply a tiny tile of A by a tiny tile of B and accumulate into a tiny tile of C that lives in vector registers.
A common AVX2 micro-kernel computes an 8×8 tile of output. It keeps 64 floating-point accumulators in registers for the entire K-loop, loads one vector from B and broadcasts one vector from A, performs a fused multiply-add, and repeats. The 64 output values are written back to memory exactly once, at the very end.
On a CPU, the micro-kernel is driven by SIMD vector instructions: AVX2, AVX-512, NEON, or SVE. On a GPU, the same idea scales into a hierarchy of threads. One thread block (CUDA CTA / thread array, or workgroup in OpenCL/Vulkan compute) computes a larger tile cooperatively, using fast shared memory as a staging area. Inside the block, warps (NVIDIA, 32 threads) or wavefronts (AMD, 64 threads) compute sub-tiles. The smallest unit may be an MMA instruction feeding tensor cores on modern NVIDIA or AMD GPUs.
From physics to ray tracing
So far this has all been abstract. But the same tile-and-transform mathematics shows up everywhere in a game engine.
A rigid body updates its world-space vertices with a rotation matrix R(t) and a translation vector x(t): p(t) = R(t) × p₀ + x(t). In homogeneous coordinates this is a single 4×4 matrix multiply. The inertia tensor updates the same way. Every vertex, every bone, every rigid prop is a small matrix-vector or matrix-matrix operation batched across many instances.
Cloth and finite-element (FEM) simulations assemble a large sparse stiffness matrix and solve a linear system each frame. The assembly itself is a scatter of small element matrices into a global matrix, and the solver is dominated by sparse matrix-vector multiplications and preconditioned conjugate gradient iterations. The data structures are sparse, but the arithmetic is still built from the same multiply-accumulate patterns.
Rigid bodies
Transforms, inertia updates, and vertex skinning reduce to matrix-matrix and matrix-vector multiplies on 3×3 and 4×4 blocks.
Cloth / mass-spring
Implicit integration assembles a stiffness matrix and solves a sparse linear system; the inner loops are matrix-vector products.
FEM deformation
Element stiffness matrices are scattered into a global system. The solve is repeated matrix-vector work accelerated by the same tiled kernels.
Then there is ray tracing. In a ray-tracing pipeline, geometry lives inside a Bottom-Level Acceleration Structure (BLAS). Instances of that geometry are placed in the scene through a Top-Level Acceleration Structure (TLAS). Each TLAS instance carries a transformation matrix that scales, rotates, and translates the BLAS into world space.
When a ray is cast, the shader may transform the ray into an object's local space, or the bounding hierarchy may transform the geometry bounds into world space. Either way, the same 4×4 matrix math used for rigid-body physics is now being used to decide which triangle a ray will hit. The linear algebra does not care whether it is moving a vertex, solving a spring, or tracing a photon.
Again, the ray-tracing BLAS is a Bottom-Level Acceleration Structure, a tree of triangles. The linear-algebra BLAS is the Basic Linear Algebra Subprograms library family. They have nothing to do with each other except a three-letter collision. Context is the only disambiguator.
The hidden connection
What ties these domains together is not the matrix format. It is the realization that the world, at the bottom, is a set of coordinates being transformed. A neural-network weight matrix transforms an input vector into an output vector. A physics transform matrix moves a body from local space to world space. A TLAS instance matrix does the same for an entire mesh. A micro-kernel does not care which one it is; it only cares about the tile size, the data layout, and the instruction set.
"The matrix is not the application. The matrix is the shape the application takes when it has been flattened into arithmetic."
This is why high-performance GEMM matters. It is the compression layer between the idea — a ray, a neuron, a vertex, a spring — and the hardware. Optimize the matrix multiply, and you optimize everything that sits above it.
Further reading and related work
We have written about several pieces of this pipeline before, from hand-tuned AVX kernels to Vulkan compute to the instanced tile rendering that feeds the same spatial-grid thinking into a 2D engine.