I will be using x86 naming throughout this post, but the concepts apply to other architectures too. Basic familiarity with SIMD and x86 assembly is assumed; uops.info and Intel’s intrinsics guide explain any instruction you don’t recognize.
A lot of SIMD code, particularly in the graphics programming space, use vec3/vec4 types all over the place. The issue is that this typically doesn’t scale beyond the SSE instruction set. In this post I want to share how you can “transpose” those vectors into lanes to scale your math operations beyond just those typical 4 component XMM registers, and have the efficiency scale up to 8 or 16 wide vectors (or maybe 32 wide vectors, if we ever get them).
To avoid the use of intrinsics, I’ll use Clang’s vector extensions (ext_vector_type) to represent the operations, since they’re easier to read and the logic can easily be reimplemented with intrinsics.
We’re going to look at a simple case on how to do N dot products at the same time (where N is the architecture’s SIMD width), and a typical CPU frustum culler as a more advanced example.
Frustum culling should ideally be done on a GPU, but that’s beyond this article
All measurements below have been made on a Ryzen 7 5800 (Zen 3) CPU, with Clang 23.1.0 targeting Windows x64. The compile flags are -std=c++20 -O3 -ffast-math, along with -march=x86-64-v3 (AVX2, FMA and BMI2) for the 8 lane code, and -march=x86-64-v2 (SSE4.2) for the 4 lane equivalents, for anyone who wants to follow along. Every number is the best of 12 runs, and varies by about 2% from run to run. All of the code is on Compiler Explorer, with one pane for each instruction set.
-ffast-mathlets the compiler assume no value is NaN or infinite. If your bounding boxes can have NaN or infinite positions or sizes, don’t use the flag: the culler’s comparisons may then give any answer for those boxes.
The basic example
Let’s take a look at what one dot product typically compiles to:
using vec3 = float __attribute__((ext_vector_type(3)));
using vec4 = float __attribute__((ext_vector_type(4)));
using vec8 = float __attribute__((ext_vector_type(8)));
float dot(vec3 a, vec3 b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
float dot(vec4 a, vec4 b) {
return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
}The assembly:
dot(float vector[3], float vector[3]):
vmovaps xmm0, xmmword ptr [rcx]
vmulps xmm0, xmm0, xmmword ptr [rdx]
vmovshdup xmm1, xmm0
vaddss xmm1, xmm0, xmm1
vshufpd xmm0, xmm0, xmm0, 1
vaddss xmm0, xmm1, xmm0
ret
dot(float vector[4], float vector[4]):
vmovaps xmm0, xmmword ptr [rcx]
vmulps xmm0, xmm0, xmmword ptr [rdx]
vshufpd xmm1, xmm0, xmm0, 1
vaddps xmm0, xmm0, xmm1
vmovshdup xmm1, xmm0
vaddss xmm0, xmm0, xmm1
retNot bad. And surprisingly, a vec3 dot costs the same as a vec4 dot. Regardless, the problem is that if we want to do multiple dot products, the cost scales linearly.
A tip for Windows: the default x64 calling convention passes any vector argument by pointer, which is where the
vmovapsand the memory operand above come from. Marking the function__vectorcall(or compiling with/Gvin clang-cl and MSVC) passes vectors in registers instead, as the System V ABI on Linux already does./Gvleaves member functions out, so a member function that takes vectors by value needs__vectorcalleither way. With it, the dot product compiles to:dotVectorcall(float vector[4], float vector[4]): vmulps xmm0, xmm1, xmm0 vshufpd xmm1, xmm0, xmm0, 1 vaddps xmm0, xmm0, xmm1 vmovshdup xmm1, xmm0 vaddss xmm0, xmm0, xmm1 ret
Let’s try transposing the vec4 one, so that each component of the vectors are stored together as one vector:
struct Vec4Lanes {
vec4 x;
vec4 y;
vec4 z;
vec4 w;
};
vec4 dot(const Vec4Lanes& a, const Vec4Lanes& b) {
return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
}The assembly:
dot(Vec4Lanes const&, Vec4Lanes const&):
vmovaps xmm0, xmmword ptr [rdx]
vmovaps xmm1, xmmword ptr [rdx + 16]
vmovaps xmm2, xmmword ptr [rdx + 32]
vmovaps xmm3, xmmword ptr [rdx + 48]
vmulps xmm0, xmm0, xmmword ptr [rcx]
vfmadd231ps xmm0, xmm1, xmmword ptr [rcx + 16]
vfmadd231ps xmm0, xmm2, xmmword ptr [rcx + 32]
vfmadd231ps xmm0, xmm3, xmmword ptr [rcx + 48]
retHardly any extra code generated, except we now generate 4 dot products at once, at the cost of a few loads, since the lanes are passed by reference.
What about doing 8 at the same time?
struct Vec8Lanes {
vec8 x;
vec8 y;
vec8 z;
vec8 w;
};
vec8 dot(const Vec8Lanes& a, const Vec8Lanes& b) {
return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
}The assembly:
dot(Vec8Lanes const&, Vec8Lanes const&):
vmovaps ymm0, ymmword ptr [rdx]
vmovaps ymm1, ymmword ptr [rdx + 32]
vmovaps ymm2, ymmword ptr [rdx + 64]
vmovaps ymm3, ymmword ptr [rdx + 96]
vmulps ymm0, ymm0, ymmword ptr [rcx]
vfmadd231ps ymm0, ymm1, ymmword ptr [rcx + 32]
vfmadd231ps ymm0, ymm2, ymmword ptr [rcx + 64]
vfmadd231ps ymm0, ymm3, ymmword ptr [rcx + 96]
retSame result.
Why doesn’t the compiler do this for us? Because it would have to change the data layout. A compiler will happily vectorize a loop over an array of floats, but turning your array of vec4s into four arrays of floats changes what every other piece of code reading that memory sees, so no compiler will reorder and transpose your data for you. That part is up to us.
A more complex example
Doing many dot products is cool and all, but it doesn’t have much of a practical use on its own. Let’s try applying the same trick to an average frustum culler.
The math for those who are unfamiliar with it: A (camera) frustum consists of six planes with their normals pointing outwards. A plane is stored as (n, d): its unit normal n, and d, its offset along the normal, chosen so that dot(n, p) + d is the signed distance from a point p to the plane. A point is in front of the plane when that distance is positive, and if a box is entirely in front of any one plane, we can cull it.
A box is stored as its center and its half extent: the offset from the center to a corner along each axis, so its corners are center ± halfExtent. A simple way to cull it is to check if all eight corners are in front of a plane.
Sample code, where each plane is a plain vec4 with the normal in xyz and the offset along the normal in w:
struct Box {
vec3 center;
vec3 halfExtent;
};
vec3 abs(vec3 v) {
return __builtin_elementwise_abs(v);
}
bool outsidePlane(vec4 plane, const Box& box) {
for (int c = 0; c < 8; c++) {
vec3 corner = box.center + box.halfExtent * vec3{
c & 1 ? 1.0f : -1.0f,
c & 2 ? 1.0f : -1.0f,
c & 4 ? 1.0f : -1.0f
};
if (dot(plane.xyz, corner) + plane.w < 0.0f) {
return false;
}
}
return true;
}
void cullNaive(uint64_t* results, const vec4* planes, const Box* boxes, uint32_t count) {
for (uint32_t i = 0; i < count; i++) {
bool visible = true;
for (uint32_t p = 0; p < 6; p++) {
if (outsidePlane(planes[p], boxes[i])) {
visible = false;
break;
}
}
const uint32_t shift = i % 64;
results[i / 64] = (results[i / 64] & ~(uint64_t(1) << shift)) | (uint64_t(visible) << shift);
}
}The cost is 13.4 ns per box. We can do better.
A simple optimization
You can achieve the same result by testing only the corner with the lowest signed distance: center - sign(n) * e, one step against the normal on every axis. If even that corner is in front of the plane, the whole box is:
float distance = dot(n, c) - dot(abs(n), e) + d; // outside if distance >= 0
The second term is the distance from the center to that corner, measured along the normal, so one test replaces eight:
void cullNearestCorner(uint64_t* results, const vec4* planes, const Box* boxes, uint32_t count) {
for (uint32_t i = 0; i < count; i++) {
bool visible = true;
for (uint32_t p = 0; p < 6; p++) {
const float distance = dot(planes[p].xyz, boxes[i].center)
- dot(abs(planes[p].xyz), boxes[i].halfExtent) + planes[p].w;
if (distance >= 0.0f) {
visible = false;
break;
}
}
const uint32_t shift = i % 64;
results[i / 64] = (results[i / 64] & ~(uint64_t(1) << shift)) | (uint64_t(visible) << shift);
}
}That’s 3.0 ns per box, so we’re already four times faster than the basic implementation.
The branch predictor
Going from 10,000 boxes to 100,000, the same loop suddenly takes 12.7 ns per box. That’s only 3 MB of data, read front to back, so memory isn’t the problem. To see what is, I timed a 1,000-box scene and cycled through a different number of scenes between runs:
| Scenes | ns per box |
|---|---|
| 1 | 3.0 |
| 16 | 3.0 |
| 18 | 5.2 |
| 20 | 7.3 |
| 32 | 11.3 |
| 200 | 13.3 |
| 200 identical copies | 3.0 |
The last two rows touch the same amount of memory. The difference is that in the last one every scene is a copy, so the early outs come in the same order every time. The CPU predicts each branch before it knows the answer, and Zen 3 can memorize the early outs of about 16,000 boxes. So if your benchmark culls the same small scene in a loop, you’re mostly benchmarking the branch predictor. Here, it made the nearest-corner test look 4.4× faster than it is on new data.
Culling multiple boxes at the same time
Now that our benchmark measures the culler rather than the branch predictor, back to the point of this post: making full use of the available SIMD width. We can just use templates to make the same code work on any type T, so float, double, vecN all work.
template <typename T>
struct Vec3Lanes {
T x, y, z;
};
template <typename T>
T dot(const Vec3Lanes<T>& a, const Vec3Lanes<T>& b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
template <typename T>
Vec3Lanes<T> cross(const Vec3Lanes<T>& a, const Vec3Lanes<T>& b) {
return { a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x };
}And then the boxes can be expanded similarly:
template <typename T>
struct Box3Lanes {
Vec3Lanes<T> center, halfExtent;
T getMinSignedDistance(vec4 plane) const {
const vec3 n = plane.xyz, a = abs(n);
return n.x * center.x - a.x * halfExtent.x
+ n.y * center.y - a.y * halfExtent.y
+ n.z * center.z - a.z * halfExtent.z
+ plane.w;
}
};The boxes are stored in blocks of Box3Lanes, one box per lane, with the last block padded out. The number of lanes W is 8 with AVX2 and 4 with SSE, so the same code fills whatever register size we compile for:
#ifdef __AVX2__
constexpr uint32_t W = 8;
#else
constexpr uint32_t W = 4;
#endif
using floatW = float __attribute__((ext_vector_type(W)));
using intW = int32_t __attribute__((ext_vector_type(W)));A whole block can’t take the same early out, so the branch has to turn into a mask:
template <typename P>
intW visibleMask(const Box3Lanes<floatW>& block, const P* planes) {
intW visible = intW(-1);
for (uint32_t p = 0; p < 6; p++) {
visible &= block.getMinSignedDistance(planes[p]) < 0.0f;
}
return visible;
}
[[gnu::always_inline]] inline void writeLanes(uint64_t* results, intW visible, uint32_t first, uint32_t count) {
for (uint32_t lane = 0; lane < W && first + lane < count; lane++) {
const uint32_t i = first + lane;
const uint64_t held = visible[lane] & 1;
results[i / 64] = (results[i / 64] & ~(uint64_t(1) << (i % 64))) | (held << (i % 64));
}
}
void cullLanes(uint64_t* results, const vec4* planes, const Box3Lanes<floatW>* blocks, uint32_t count) {
for (uint32_t first = 0; first < count; first += W) {
writeLanes(results, visibleMask(blocks[first / W], planes), first, count);
}
}This runs at 1.24 ns per box with AVX2 and 2.87 ns with SSE. On AVX2 that’s 10× faster than the nearest-corner test on fresh data. The SSE number is disappointing, though. Against the nearest-corner test with a trained predictor it’s only 1.3× faster, which isn’t much to show for four lanes. So let’s look at the assembly.
The assembly in question
For one block of eight boxes against six planes, the math compiles to:
- 6 loads;
- 6
vmulpsand 30vfmadd231ps; - 6
vcmpltpsand 5vandpsto combine the masks.
Those 53 instructions are fine. The problem is that the whole loop is about 220.
The frustum is rebuilt in every iteration.
vandps xmm7, xmm3, xmmword ptr [rdx] ; abs(normal) of plane 0
vbroadcastss ymm0, xmm7 ; splat |n.x|
vmovups ymm1, ymmword ptr [rsp + 48] ; reload a permutation...
vpermps ymm5, ymm1, ymm7 ; ...that splats |n.y|Every block gets 30 vbroadcastss, 12 vpermps and 6 vandps, and none of them have anything to do with the block. The compiler does this because six planes need 42 splats and AVX2 only has 16 ymm registers. Its register allocator figures that recomputing a broadcast is cheaper than spilling it and loading it back. For a plain broadcast from memory that’s true, but it redoes the abs and the permutes too.
The results are written one bit at a time.
vpextrd r9d, xmm3, 3 ; take one lane out of the mask
mov r15d, r10d
and r15b, 56 ; the block's first bit within the word
mov ebp, 8
shlx rbp, rbp, r15
andn rbp, rbp, r11 ; clear the lane's bit
or r15b, 3 ; the lane's bit within the word
shlx r11, r9, r15 ; shift the lane to that bit
or r11, rbp
mov qword ptr [rcx + 8*rsi], r11 ; store the word
lea r9, [rax - 4]
cmp r9d, r10d ; was that the last object?
je .LBB12_11That’s thirteen instructions per lane, about 100 per block, which is more than the math itself.
Splatting the planes once
So we do the splatting ourselves, once, before the loop, and keep the results in memory:
template <typename T>
struct PlaneLanes {
Vec3Lanes<T> normal, absNormal;
T constant;
};
// Broadcasts each component of v into its own register
Vec3Lanes<floatW> splat(vec3 v) {
return { floatW(v.x), floatW(v.y), floatW(v.z) };
}
PlaneLanes<floatW> splatPlane(vec4 plane) {
return { splat(plane.xyz), splat(abs(plane.xyz)), floatW(plane.w) };
}
void cullLanesSplatOnce(uint64_t* results, const vec4* frustum, const Box3Lanes<floatW>* blocks, uint32_t count) {
PlaneLanes<floatW> planes[6];
for (uint32_t p = 0; p < 6; p++)
planes[p] = splatPlane(frustum[p]);
for (uint32_t first = 0; first < count; first += W) {
writeLanes(results, visibleMask(blocks[first / W], planes), first, count);
}
}Box3Lanes gets a second overload for the splatted planes:
T getMinSignedDistance(const PlaneLanes<T>& p) const {
return dot(p.normal, center) - dot(p.absNormal, halfExtent) + p.constant;
}Now the compiler reads the planes straight from the stack as memory operands of the FMAs:
vfmadd231ps ymm3, ymm15, ymmword ptr [rsp + 1056]That gets AVX2 down to 0.98 ns per box and SSE to 1.85 ns. SSE benefits more because its blocks only hold four boxes, so it was rebuilding the frustum twice as often.
Writing the block at once
The eight result bits of a block make up exactly one byte, so we can turn the mask into one bit per lane and store all eight at once. With SSE, two blocks of four share a byte, so each one writes its half:
using boolW = bool __attribute__((ext_vector_type(W)));
void writeBlock(uint64_t* results, intW visible, uint32_t first) {
const boolW bits = __builtin_convertvector(visible, boolW); // one bit per lane
uint8_t& byte = reinterpret_cast<uint8_t*>(results)[first / 8];
if constexpr (W == 8) {
byte = __builtin_bit_cast(uint8_t, bits); // store all eight
} else {
const uint32_t shift = first % 8; // two blocks of four share a byte
const uint8_t held = __builtin_bit_cast(uint8_t, bits) & 0xF;
byte = (byte & ~(0xF << shift)) | (held << shift);
}
}
void cullLanesBlockWrite(uint64_t* results, const vec4* planes, const Box3Lanes<floatW>* blocks, uint32_t count) {
uint32_t first = 0;
for (; first + W <= count; first += W) {
writeBlock(results, visibleMask(blocks[first / W], planes), first);
}
if (first < count) { // only the last block, if it isn't full, still goes lane by lane
writeLanes(results, visibleMask(blocks[first / W], planes), first, count);
}
}The partial last block is handled after the loop rather than with a branch inside it. That matters more than it looks: when the mask is created in one basic block and turned into bits in another, Clang no longer knows every lane is either all ones or all zeros, and replaces the single vmovmskps with a chain of shuffles and packs.
The writeback now looks like this:
vmovmskps r10d, ymm2 ; one bit per lane
mov byte ptr [rcx + rax], r10b ; store all eightWith both changes in place, AVX2 takes 0.67 ns per box1 and SSE 1.84 ns. The loop is down to about 95 straight-line instructions per block of eight, or about 3 cycles per box.
Results
All times are in nanoseconds per box, for 10,000 boxes. The numbers in brackets are for 100,000 boxes, which is too many for the branch predictor to memorize.
| Method | AVX2, 8 lanes | SSE4.2, 4 lanes |
|---|---|---|
| The naive method | 13.4 (17.2) | 15.8 (19.1) |
| Nearest corner | 3.04 (12.7) | 3.77 (13.9) |
| One box per lane | 1.24 | 2.87 |
| Splatting the planes once | 0.98 | 1.85 |
| Block writeback, without splatting | 1.04 | 2.69 |
| Splatting once and block writeback | 0.671 | 1.84 |
Splatting once helps both, SSE the most (1.55× on its own, against 1.27× on AVX2), while the block writeback alone does less (1.19× on AVX2, 1.07× on SSE). On AVX2 the two compound, since once the splats are gone, the bit-by-bit writeback is most of what’s left; together they make the one-box-per-lane version 1.85× faster. On SSE the writeback adds nothing once the planes are splatted, and the combination stays at 1.56×. From start to finish, the culler got about 20× faster, or 25× on data the branch predictor hadn’t seen before.
Conclusion
Most of the speedup comes from turning vectors into lanes. Eight dot products cost three instructions, and the early outs become masks. The less obvious part is that the straightforward SIMD version spent more instructions moving data around than on the actual math, and it took two more changes to get rid of that. Precompute whatever the compiler can’t keep in registers, and write the results a whole block at a time. You’ll run into the same two problems in any loop that tests lots of objects against a few shared values.
1 Only when the stack happens to be 32-byte aligned. Windows and Linux only guarantee 16 bytes, and Clang spills the splatted planes’ ymm registers with unaligned stores rather than realigning the stack, so in some processes half of the reloads straddle a cache line, and the same loop takes 0.99 ns. Adding __attribute__((force_align_arg_pointer)) to the function makes Clang realign the stack in its prologue, which gives 0.67 ns every time.