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

Question

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++ )
    {
        outputs_.xs[i] = V4_ADD(
            V4_ADD(
                V4_MUL( m00, inputs_.xs[i] ),
                V4_MUL( m01, inputs_.ys[i] )
            ),
            V4_ADD(
                V4_MUL( m02, inputs_.zs[i] ),
                 m03
            )
        );

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

        outputs_.zs[i] = V4_ADD(
            V4_ADD(
                V4_MUL( m20, inputs_.xs[i] ),
                V4_MUL( m21, inputs_.ys[i] )
            ),
            V4_ADD(
                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.

Answer 1

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++ )
{
    outxs[i] = _mm_fmadd_ps(m02, zs[i],
               _mm_fmadd_ps(m01, ys[i],
               _mm_fmadd_ps(m00, xs[i], m03)));

    outys[i] = _mm_fmadd_ps(m12, zs[i],
               _mm_fmadd_ps(m11, ys[i],
               _mm_fmadd_ps(m10, xs[i], m13)));

    outzs[i] = _mm_fmadd_ps(m22, zs[i],
               _mm_fmadd_ps(m21, ys[i],
               _mm_fmadd_ps(m20, xs[i], m23)));
}

The actual code looks pretty good, apart from the fact that it uses short vectors:

$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_

That doesn't really address the SSE2 side of the store directly, but you can do the same thing but with separate add/mul. Proper use of AVX would use the wider vectors. For this code that's a fairly trivial change, just make almost everything wider (except the matrix). If AVX is available, you should use it in this case - it doesn't help for everything, but for this type of code (almost purely vertical SIMD, apart from some broadcasts) it's great. FMA is also a free performance win here (especially on Haswell and Broadwell), well worth using even if it means writing two versions and doing runtime dispatch to also support the Bridges (which are perhaps not old enough yet to completely disregard).

Supporting 32bit machines is annoying. They don't even have enough registers to load all of that matrix, so everything gets bogged down by the extra loads that are suddenly required. I don't see good fixes. Rearranging the multiplication so that the matrix can be held in 4 registers does work, but then a horizontal addition appears which is bad, and the multiplication wastes a lane on padding. I expect it would be worse than the extra loads, on the other hand if we're talking about old hardware like Core2 (which had x64 support, but at the time installing a 64bit OS was a rarity) then extra loads are extra bad since the load throughput used to be only 1/cycle. On the other hand, horizontal addition also used to be much worse than it is now. It just seems like all options are terrible. It depends on your audience of course, but frankly I don't think it's worth expending much energy on.