Fast affine transformations of many 3D points by one 3×4 matrix
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 &lt; 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.