探How 社区

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…

以下正文同步自 Hacker News 最新，版权归原站所有，已转换为易读排版。

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-math lets 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.

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; }

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 ret

Not 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 vmovaps and the memory operand above come from. Marking the function __vectorcall (or compiling with /Gv in clang-cl and MSVC) passes vectors in registers instead, as the System V ABI on Linux already does. /Gv leaves member functions out, so a member function that takes vectors by value needs __vectorcall either 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; }

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] ret

Hardly 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.

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; }

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] ret

Same 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.

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); } }

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.

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); } }

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.

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:

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

正文由探How 从来源站点 RSS 同步，内容未经改写；遇到排版缺失或需要图片、视频时请以原文为准。

Making SIMD Code Scale to Any Vector Width