# Fast affine transformations of many 3D points by one 3×4 matrix

I wrote a function to batch-transform 3D vectors by a single 3x4 matrix using SSE2:

struct alignas(16) Matrix3x4f
{
union
{
Vector4 r[3];   // (r[0].w, r[1].w, r[2].w) contains the translation
float   m[3][4];
};
// (0, 0, 0, 1) is the forth row of the full 4x4 matrix
};

struct Positions_SoA
{
Vector4 *   xs;
Vector4 *   ys;
Vector4 *   zs;
};

void transform_points_SSE2(
const Matrix3x4f& transform_,
const Positions_SoA& inputs_,
const int number_of_packets_,
Positions_SoA &outputs_
)
{
//FIXME: woudn't this cause too much register pressure and memory loads?
// For AVX/AVX2 we'll have to splat 8/16 values.
const Vector4 m00 = SPLAT_X( transform_.r[0] );
const Vector4 m01 = SPLAT_Y( transform_.r[0] );
const Vector4 m02 = SPLAT_Z( transform_.r[0] );
const Vector4 m03 = SPLAT_W( transform_.r[0] ); // translation X

const Vector4 m10 = SPLAT_X( transform_.r[1] );
const Vector4 m11 = SPLAT_Y( transform_.r[1] );
const Vector4 m12 = SPLAT_Z( transform_.r[1] );
const Vector4 m13 = SPLAT_W( transform_.r[1] ); // translation Y

const Vector4 m20 = SPLAT_X( transform_.r[2] );
const Vector4 m21 = SPLAT_Y( transform_.r[2] );
const Vector4 m22 = SPLAT_Z( transform_.r[2] );
const Vector4 m23 = SPLAT_W( transform_.r[2] ); // translation Z

for( int i = 0; i < number_of_packets_; i++ )
{
V4_MUL( m00, inputs_.xs[i] ),
V4_MUL( m01, inputs_.ys[i] )
),
V4_MUL( m02, inputs_.zs[i] ),
m03
)
);

V4_MUL( m10, inputs_.xs[i] ),
V4_MUL( m11, inputs_.ys[i] )
),
V4_MUL( m12, inputs_.zs[i] ),
m13
)
);

V4_MUL( m20, inputs_.xs[i] ),
V4_MUL( m21, inputs_.ys[i] )
),
V4_MUL( m22, inputs_.zs[i] ),
m23
)
);
}
}


The complete, runnable code is on rextester.

1. When I run the tests and compare the numbers with the reference (scalar) version, the results of the SSE version appear slightly wrong in half of all the cases, although the printed numbers look the same. Is that due to rounding errors? (FPU internally uses 80-bit precision, and SSE units compute with 32-bit precision.) Should I ignore those small discrepancies? Is my 'optimized' function even correct?

2. Why is the speedup due to SSE so small? I expected 3-3.5 speed increase, but got less than 2. Should I use AVX?

3. I plan to transform up to 64 points, is it still worth to use SSE/AVX? Will performance be dominated by the cost of loading into registers? I'm worried about the 12 shuffle/broadcast/replicate instructions.

4. Finally, I'll appreciate comments on my programming style, my mistakes in using SSE, a note on best practices and performance pitfalls.

• Is this meant for a 32bit target or 64bit? Commented Nov 26, 2017 at 17:03
• @harold The code is meant to run on 64-bit targets for best performance, but I'd also like to achieve good performance on older, 32-bit machines. Commented Nov 26, 2017 at 18:46

Small differences due to precision are expected and can usually be ignored.

12 shuffles like that are a bit much, though not necessarily avoidable, depending on whether there is AVX support. With AVX, it is better to literally broadcast from memory, rather than emulate broadcasting with a load and shuffles. Even though this means there will be more loads, the shuffles are a bigger problem: loads (even broadcasting loads) typically have a throughput of 2/cycle while the shuffles typically have a throughput of 1/cycle. This means that _mm_set1_ps(float) is actually the better option, as long as you compile with AVX support. That wouldn't use the wider width yet though.

Also annoying is that MSVC apparently likes to reload those pointers inputs_.xs etc continuously. There is really no spare load-throughput for that, despite the presence of 9 FP additions (which are also bad though, since they're 1/cycle pre-Skylake), there are 21 loads there, limiting the throughput of the loop to once every 10.5 cycles (which is more than 9, if barely). Fortunately it can very easily be convinced to stop doing that, simply by copying them into local variables before the loop.

Next, I don't know what MSVC is up to there with its refusal to contract addition and multiplication into FMA, for processors that support it it should definitely be used. Using the corresponding intrinsic does the trick of course, but that makes it harder to compile for pre-FMA targets (Ivy Bridge and older). Anyway if I write the high level code like this:

Vector4 *   xs = inputs_.xs;
Vector4 *   ys = inputs_.ys;
Vector4 *   zs = inputs_.zs;
Vector4 *   outxs = outputs_.xs;
Vector4 *   outys = outputs_.ys;
Vector4 *   outzs = outputs_.zs;

for( int i = 0; i < number_of_packets_; i++ )
{


$LL4@transform_: lea rcx, QWORD PTR [rdx+rax] lea rax, QWORD PTR [rax+16] vmovups xmm0, xmm7 vfmadd231ps xmm0, xmm4, XMMWORD PTR [rax-16] vfmadd231ps xmm0, xmm5, XMMWORD PTR [r10+rax-16] vfmadd231ps xmm0, xmm6, XMMWORD PTR [r11+rax-16] vmovups XMMWORD PTR [rcx+rbx], xmm0 vmovups xmm0, xmm3 vfmadd231ps xmm0, xmm8, XMMWORD PTR [rax-16] vfmadd231ps xmm0, xmm9, XMMWORD PTR [r10+rax-16] vfmadd231ps xmm0, xmm10, XMMWORD PTR [r11+rax-16] vmovups XMMWORD PTR [rcx], xmm0 vmovups xmm0, xmm1 vfmadd231ps xmm0, xmm11, XMMWORD PTR [rax-16] vfmadd231ps xmm0, xmm12, XMMWORD PTR [r10+rax-16] vfmadd231ps xmm0, xmm13, XMMWORD PTR [r11+rax-16] vmovups XMMWORD PTR [r9+rax-16], xmm0 sub r8, 1 jne SHORT$LL4@transform_