SIMD Fundamentals. Part III: Implementation & Benchmarks

This is the third part of my series on SIMD Fundamentals, targeted at .NET developers. You may read the other two parts here:

• Part I: From SISD to SIMD
• Part II: AoS, SoA, Gather/Scatter - Oh my!

While the previous two parts where more focused on theory and concepts, we are now going to actually get our hands dirty and write some code to see how different approaches to SIMD processing compare in practice, both from a performance an an ease-of-implementation point of view.

Prelude

In Part II I used F# (pseudo-)code to illustrate the difference between the AoS (array of structures) and SoA (structure of array) approaches. That's why I thought using F# for implementing the benchmarks might be a good idea. So I installed Visual Studio 2015, created a new F# console application project and installed System.Numerics.Vectors in its most recent 4.1.0 incarnation via NuGet. Yet, when I tried to use System.Numerics.Vector<T> IntelliSense wanted to convince me there was no such thing:

Maybe just a problem with the F# language service? I tried to run this little hello world sample

	open System.Numerics
	[<EntryPoint>]
	let main argv =
	let laneWidth = System.Numerics.Vector<float32>.Count
	printfn "%i" laneWidth
	0

view raw fsvectest.fs hosted with ❤ by GitHub

but that didn't work either, because it references the wrong version of System.Numerics.Vectors:

	------ Build started: Project: ConsoleApplication1, Configuration: Release Any CPU ------
	C:\Program Files (x86)\Microsoft SDKs\F#\4.0\Framework\v4.0\fsc.exe -o:obj\Release\ConsoleApplication1.exe --debug:pdbonly --noframework --define:TRACE --doc:bin\Release\ConsoleApplication1.XML --optimize+ --platform:x64 -r:"C:\Program Files (x86)\Reference Assemblies\Microsoft\FSharp\.NETFramework\v4.0\4.4.0.0\FSharp.Core.dll" -r:"C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.6\mscorlib.dll" -r:"C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.6\System.Core.dll" -r:"C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.6\System.dll" -r:"C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.6\System.Numerics.dll" -r:"C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.6\System.Numerics.Vectors.dll" --target:exe --warn:3 --warnaserror:76 --vserrors --validate-type-providers --LCID:1033 --utf8output --fullpaths --flaterrors --subsystemversion:6.00 --highentropyva+ --sqmsessionguid:dbfb1201-d8ea-4f7f-af58-f587b9d5e3cb "C:\Users\Frank\AppData\Local\Temp\.NETFramework,Version=v4.6.AssemblyAttributes.fs" AssemblyInfo.fs Program.fs
	AosVsSoa\ConsoleApplication1\ConsoleApplication1\Program.fs(5,37): error FS0039: The value, constructor, namespace or type 'Vector' is not defined
	Done building project "ConsoleApplication1.fsproj" -- FAILED.

view raw fsbuildfail.txt hosted with ❤ by GitHub

I didn't have luck with manually replacing the "C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.6\System.Numerics.Vectors.dll" reference with the one delivered by the System.Numerics.Vectors NuGet package either.

For now, I thus resorted to using C# as an implementation language instead.

Managed Implementations

For the upcoming benchmarks we will use the same problem we stated in the previous post: Compute the squared L2 norm of a set of 3-vectors. According to the previous post, we want to compare five possible approaches:

• Scalar AoS
• AoS vectorized across each 3-vector
• On-the-fly vectorization (converting AoS to SoA)
• Scalar SoA
• Vectorized SoA

To avoid any memory bottlenecks and keep everything cached, we will use a comparatively small set of vectors. I found that an array of 1024 to 2048 Vec3fs yielded peak performance. Note that this has the nice side effect of being a multiple of the AVX vector lane width (32 byte or 8 single precision floating point values). To get a reliable estimate of the average performance, we repeat the procedure a large number of times: For the benchmarks, all procedures had to compute a total of 32 · 2^30 (approximately 32 billion) dot products.

For further details, you can consult the full benchmark code here. In case you want to run it yourself and change some parameters, make sure that vectorLen is a multiple of laneWidth, as the remaing code assumes.

All of the following benchmark code was compiled and run on an Intel Core i5-4570 (Haswell) with 16 GB auf DD3-SDRAM on Windows 10 Pro. The managed implementation was developed using Visual Studio 2015 on .NET 4.6 and targeted x64 (release build, "RyuJIT").

Scalar AoS

This version is probably the easiest to understand and, as long as you don't intend to vectorize, a pretty reasonable one: We have an array of Vector3 structures (here we simply use the one provided by System.Numerics.Vectors) and compute the resulting dot products vector by vector:

	static void Dot3AosScalar(Vector3[] vs, float[] dp) {
	for (var j = 0; j < reps; ++j) {
	for (var i = 0; i < dp.Length; ++i) {
	dp[i] = vs[i].X * vs[i].X + vs[i].Y * vs[i].Y + vs[i].Z * vs[i].Z;
	}
	}
	}

view raw csaosscalar.cs hosted with ❤ by GitHub

The outer loop over j is there to ensure the total number of computed dot products is 32 billion. For the inner loop, the JIT compiler generates the following machine code:

	cmp r9d,r10d
	jae 00007FF98D096F68
	movsxd r11,r9d
	imul r11,r11,3
	vmovss xmm0,dword ptr [rcx+r11*4+10h]
	vmulss xmm0,xmm0,xmm0
	lea r11,[rcx+r11*4+10h]
	vmovss xmm1,dword ptr [r11+4]
	vmulss xmm1,xmm1,xmm1
	vaddss xmm0,xmm0,xmm1
	vmovss xmm1,dword ptr [r11+8]
	vmulss xmm1,xmm1,xmm1
	vaddss xmm0,xmm0,xmm1
	movsxd r11,r9d
	vmovss dword ptr [rdx+r11*4+10h],xmm0
	inc r9d
	cmp r8d,r9d
	jg 00007FF98D096F06

view raw csaosscalar.asm hosted with ❤ by GitHub

As expected, it uses scalar AVX instructions to compute the dot products (VMULSS, VADDSS).

AoS vectorized across each 3-vector

In this version, we still compute a single dot product per iteration, but we compute the squares of the components at once and then determine the (horizontal) sum of those squared components:

	static void Dot3AosVectorDp(Vector3[] vs, float[] dp) {
	for (var j = 0; j < reps; ++j) {
	for (var i = 0; i < dp.Length; ++i) {
	dp[i] = Vector3.Dot(vs[i], vs[i]);
	}
	}
	}

view raw csaosvecdp.cs hosted with ❤ by GitHub

When using Vector3.Dot, the compiler emits code that uses the DPPS (SSE 4.1) or VDPPS (AVX) instructions:

	cmp r9d,r10d
	jae 00007FF98D077A54
	movsxd r11,r9d
	imul r11,r11,3
	lea r11,[rcx+r11*4+10h]
	vmovss xmm1,dword ptr [r11+8]
	vmovsd xmm0,qword ptr [r11]
	vshufps xmm0,xmm0,xmm1,44h
	vmovss xmm2,dword ptr [r11+8]
	vmovsd xmm1,qword ptr [r11]
	vshufps xmm1,xmm1,xmm2,44h
	vdpps xmm0,xmm0,xmm1,0F1h
	movsxd r11,r9d
	vmovss dword ptr [rdx+r11*4+10h],xmm0
	inc r9d
	cmp r8d,r9d
	jg 00007FF98D0779F6

view raw csaosvecdp.asm hosted with ❤ by GitHub

For some reason, it's loading the data twice, first into XMM0 and then into XMM1 (VMOVSS, VMOVSD, VSHUFPS).

On-the-fly vectorization

Because AoS isn't a data layout well suited for vectorization, last time we came up with the idea of reordering the data on the fly. Vector gather instructions would help with that, but System.Numerics.Vector<T> only supports loading consecutive elements for now. The only managed solution I could come with is to first manually gather the required data into temporary arrays and then creating the vector instances from these temporary data structures:

	static void Dot3AosGather(Vector3[] vs, float[] dp) {
	var xtmp = new float[laneWidth];
	var ytmp = new float[laneWidth];
	var ztmp = new float[laneWidth];
	for (var j = 0; j < reps; ++j) {
	for (var i = 0; i < dp.Length; i += laneWidth) {
	for (var k = 0; k < laneWidth; ++k) {
	xtmp[k] = vs[i + k].X;
	ytmp[k] = vs[i + k].Y;
	ztmp[k] = vs[i + k].Z;
	}
	var x = new Vector<float>(xtmp);
	var y = new Vector<float>(ytmp);
	var z = new Vector<float>(ztmp);
	var dpv = x * x + y * y + z * z;
	dpv.CopyTo(dp, i);
	}
	}
	}

view raw csaosgather.cs hosted with ❤ by GitHub

That works, in principle, meaning that the compiler can now emit VMULPS and VADDPS instructions to compute 8 dot products at once. Yet, because the JIT compiler doesn't employ VGATHERDPS all this gathering becomes quite cumbersome:

	xor r9d,r9d
	mov r10d,dword ptr [rbx+8]
	movsxd r10,r10d
	cmp r10,8
	setge r10b
	movzx r10d,r10b
	mov r11d,dword ptr [rbp+8]
	movsxd r11,r11d
	cmp r11,8
	setge r11b
	movzx r11d,r11b
	and r10d,r11d
	mov r11d,dword ptr [rax+8]
	movsxd r11,r11d
	cmp r11,8
	setge r11b
	movzx r11d,r11b
	and r10d,r11d
	test r10d,r10d
	je 00007FF98D078059
	mov r10d,dword ptr [rsi+8]
	lea r11d,[r8+r9]
	cmp r11d,r10d
	jae 00007FF98D078144
	movsxd r11,r11d
	imul r14,r11,3
	vmovss xmm0,dword ptr [rsi+r14*4+10h]
	movsxd r14,r9d
	vmovss dword ptr [rbx+r14*4+10h],xmm0
	imul r11,r11,3
	lea r11,[rsi+r11*4+10h]
	vmovss xmm1,dword ptr [r11+4]
	movsxd r14,r9d
	vmovss dword ptr [rbp+r14*4+10h],xmm1
	vmovss xmm2,dword ptr [r11+8]
	movsxd r11,r9d
	vmovss dword ptr [rax+r11*4+10h],xmm2
	inc r9d
	cmp r9d,8
	jl 00007FF98D077FFD
	jmp 00007FF98D0780DF
	mov r10d,dword ptr [rsi+8]
	lea r11d,[r8+r9]
	cmp r11d,r10d
	jae 00007FF98D078144
	movsxd r10,r11d
	imul r11,r10,3
	vmovss xmm0,dword ptr [rsi+r11*4+10h]
	mov r11d,dword ptr [rbx+8]
	cmp r9d,r11d
	jae 00007FF98D078144
	movsxd r11,r9d
	vmovss dword ptr [rbx+r11*4+10h],xmm0
	imul r10,r10,3
	lea r10,[rsi+r10*4+10h]
	vmovss xmm1,dword ptr [r10+4]
	mov r11d,dword ptr [rbp+8]
	cmp r9d,r11d
	jae 00007FF98D078144
	movsxd r11,r9d
	vmovss dword ptr [rbp+r11*4+10h],xmm1
	vmovss xmm2,dword ptr [r10+8]
	mov r10d,dword ptr [rax+8]
	cmp r9d,r10d
	jae 00007FF98D078144
	movsxd r10,r9d
	vmovss dword ptr [rax+r10*4+10h],xmm2
	inc r9d
	cmp r9d,8
	jl 00007FF98D078059
	vmovupd ymm0,ymmword ptr [rbx+10h]
	vmovupd ymm1,ymmword ptr [rbp+10h]
	vmovupd ymm2,ymmword ptr [rax+10h]
	vmulps ymm0,ymm0,ymm0
	vmulps ymm1,ymm1,ymm1
	vaddps ymm0,ymm0,ymm1
	vmulps ymm1,ymm2,ymm2
	vaddps ymm0,ymm0,ymm1
	lea r9d,[r8+7]
	cmp r9d,ecx
	jae 00007FF98D078144
	vmovupd ymmword ptr [rdi+r8*4+10h],ymm0
	add r8d,8
	cmp ecx,r8d
	jg 00007FF98D077FB2

view raw csaosgather.asm hosted with ❤ by GitHub

As you can probably imagine, this code isn't exactly a candidate for the world's fastest dot3 product implementation...

Scalar SoA

Let's move on to a proper SoA data layout. The scalar version is similar to the scalar AoS version, only that we now index the components instead of the array of vectors:

	static void Dot3SoaScalar(float[] xs, float[] ys, float[] zs, float[] dp) {
	for (var j = 0; j < reps; ++j) {
	for (var i = 0; i < dp.Length; ++i) {
	dp[i] = xs[i] * xs[i] + ys[i] * ys[i] + zs[i] * zs[i];
	}
	}
	}

view raw cssoascalar.cs hosted with ❤ by GitHub

The resulting machine code is likewise similar to scalar AoS:

	movsxd rsi,r11d
	vmovss xmm0,dword ptr [rcx+rsi*4+10h]
	vmulss xmm0,xmm0,xmm0
	movsxd rsi,r11d
	vmovss xmm1,dword ptr [rdx+rsi*4+10h]
	vmulss xmm1,xmm1,xmm1
	vaddss xmm0,xmm0,xmm1
	movsxd rsi,r11d
	vmovss xmm1,dword ptr [r8+rsi*4+10h]
	vmulss xmm1,xmm1,xmm1
	vaddss xmm0,xmm0,xmm1
	movsxd rsi,r11d
	vmovss dword ptr [r9+rsi*4+10h],xmm0
	inc r11d
	cmp r10d,r11d
	jg 00007FF98D076BDF

view raw cssoascalar.asm hosted with ❤ by GitHub

Vectorized SoA

The SoA layout makes it easy to use vector arithmetic to compute 8 dot products at once:

	static void Dot3SoaVectorized(float[] xs, float[] ys, float[] zs, float[] dp) {
	for (var j = 0; j < reps; ++j) {
	for (var i = 0; i < dp.Length; i += laneWidth) {
	var x = new Vector<float>(xs, i);
	var y = new Vector<float>(ys, i);
	var z = new Vector<float>(zs, i);
	var d = x * x + y * y + z * z;
	d.CopyTo(dp, i);
	}
	}
	}

view raw cssoavectorized.cs hosted with ❤ by GitHub

This is what the compiler makes of it:

	lea ebp,[r11+7]
	cmp ebp,esi
	jae 00007FF98D06828A
	vmovupd ymm0,ymmword ptr [rcx+r11*4+10h]
	cmp ebp,edi
	jae 00007FF98D06828A
	vmovupd ymm1,ymmword ptr [rdx+r11*4+10h]
	cmp ebp,ebx
	jae 00007FF98D06828A
	vmovupd ymm2,ymmword ptr [r8+r11*4+10h]
	vmulps ymm0,ymm0,ymm0
	vmulps ymm1,ymm1,ymm1
	vaddps ymm0,ymm0,ymm1
	vmulps ymm1,ymm2,ymm2
	vaddps ymm0,ymm0,ymm1
	cmp ebp,r10d
	jae 00007FF98D06828A
	vmovupd ymmword ptr [r9+r11*4+10h],ymm0
	add r11d,8
	cmp r10d,r11d
	jg 00007FF98D068220

view raw cssoavectorized.asm hosted with ❤ by GitHub

Nice, but sprinkled with range (?) checks. I also wonder, why it emits VMOVUPD instead of VMOVUPS instructions.

Unmanaged Implementations

After implementing and running the above variants in C#, I figured it would be useful to have something to compare the results to. Thus, I ported the benchmark code to C++ to see what the Visual C++ optimizer, its auto-vectorizer and SIMD intrinsics can do and how close we can get to the theoretical peak performance of the Haswell CPU. For the "native" implementation I used the Visual C++ 2015 compiler with the following flags:

/GS- /GL /W3 /Gy /Zi /Gm- /Ox /Ob2 /Zc:inline /fp:fast /WX- /Zc:forScope /arch:AVX2 /Gd /Oy /Oi /MD /Ot

Scalar AoS

Again, the code for this version is pretty straightforward:

	void dot3_aos_scalar(const vector<Vec3f>& vs, vector<float>& dp) {
	for (auto j = 0; j < reps; ++j) {
	auto i = vector_len;
	while (i--) {
	dp[i] = vs[i].x * vs[i].x + vs[i].y * vs[i].y + vs[i].z * vs[i].z;
	}
	}
	}

view raw cppaosscalar.cpp hosted with ❤ by GitHub

In case you wonder about the inner loop construction: while (i--) turned out to result in slightly faster code than a more traditional for loop.

No surprises regarding the machine code, either:

	lea rax,[rax-0Ch]
	lea rdx,[rdx-4]
	vmovss xmm0,dword ptr [rax-8]
	vmovss xmm2,dword ptr [rax-4]
	vmovss xmm3,dword ptr [rax]
	vmulss xmm1,xmm0,xmm0
	vmulss xmm0,xmm2,xmm2
	vaddss xmm2,xmm1,xmm0
	vmulss xmm1,xmm3,xmm3
	vaddss xmm2,xmm2,xmm1
	vmovss dword ptr [rdx],xmm2
	sub ecx,1
	jne benchmark<<lambda_918319110c26e9fabd8b05fc8f2cd5cd>,<lambda_c092a5680821d9f5b8bc5a7043f59100> >+0B3h (07FF656A82803h)

view raw cppaosscalar.asm hosted with ❤ by GitHub

AoS vectorized across each 3-vector

Let's use the DPPS instruction via intrinsics:

	void dot3_aos_vector_dp(const vector<Vec3f>& vs, vector<float>& dp) {
	// 0000 0000 0111 0001: mul lower three components, store sum in lowest component
	static const auto mask = 0x71;
	for (auto j = 0; j < reps; ++j) {
	const auto pvs = (float*)vs.data();
	auto pdp = (float*)dp.data();
	auto i = vector_len;
	while (i--) {
	// load 16 bytes (xyz|x)
	const auto xyzx = _mm_loadu_ps(pvs + i * 3);
	// compute d = x*x + y*y + z*z + 0*0 -> 000d
	const auto xxyyzz00 = _mm_dp_ps(xyzx, xyzx, mask);
	// store d (lower 4 bytes of dpv)
	_mm_store_ss(pdp + i, xxyyzz00);
	}
	}
	}

view raw cppaosvecdp.cpp hosted with ❤ by GitHub

Notice the little trick of directly loading four consecutive floats instead of scalar loads and shuffling. Strictly speaking, this might go wrong for the last element of the vector, if you try to access unallocated memory... In reality you'd handle that special case separately (or simply allocate a few more bytes). The corresponding machine code is really compact:

	lea rcx,[rcx-0Ch]
	lea rax,[rax-4]
	vmovups xmm0,xmmword ptr [rcx]
	vdpps xmm0,xmm0,xmm0,71h
	vmovss dword ptr [rax],xmm0
	sub edx,1
	jne benchmark<<lambda_d1ac89d5e59a169233af7a419374e043>,<lambda_c092a5680821d9f5b8bc5a7043f59100> >+0C0h (07FF7B83F2100h)

view raw cppaosvecdp.asm hosted with ❤ by GitHub

On-the-fly vectorization

In contrast to the managed version, we can now employ AVX's vector gather instructions to load eight of each 3ed component value into YMM registers:

	void dot3_aos_vector_gather(const vector<Vec3f>& vs, vector<float>& dp) {
	static const auto epi32_one = _mm256_set1_epi32(1);
	static const auto x_offsets = _mm256_setr_epi32(0, 3, 6, 9, 12, 15, 18, 21);
	static const auto y_offsets = _mm256_add_epi32(x_offsets, epi32_one);
	static const auto z_offsets = _mm256_add_epi32(y_offsets, epi32_one);
	for (auto j = 0; j < reps; ++j) {
	const auto pvs = (float*)vs.data();
	auto pdp = (__m256*)dp.data();
	auto i = vector_len / lane_width;
	while(i--) {
	// xyz|xyz|xyz|xyz|xyz|xyz|xyz|xyz -> load 8 * 3 = 24 scattered floats
	const auto base = pvs + i * 3 * lane_width;
	const auto xs = _mm256_i32gather_ps(base, x_offsets, sizeof(float));
	const auto ys = _mm256_i32gather_ps(base, y_offsets, sizeof(float));
	const auto zs = _mm256_i32gather_ps(base, z_offsets, sizeof(float));
	const auto xx = _mm256_mul_ps(xs, xs);
	const auto yyxx = _mm256_fmadd_ps(ys, ys, xx);
	const auto zzyyxx = _mm256_fmadd_ps(zs, zs, yyxx);
	pdp[i] = zzyyxx;
	}
	}
	}

view raw cppaosgather.cpp hosted with ❤ by GitHub

_mm256_fmadd_ps results in FMA3 (fused multiply add) instructions, combining multiplication and addition/accumulation in one instruction:

	lea rax,[rax-60h]
	lea r8,[r8-20h]
	vpcmpeqb ymm2,ymm2,ymm2
	vmovups ymm5,ymm0
	vgatherdps ymm5,dword ptr [rax+ymm6*4],ymm2
	vpcmpeqb ymm2,ymm2,ymm2
	vmovups ymm4,ymm0
	vgatherdps ymm4,dword ptr [rax+ymm7*4],ymm2
	vpcmpeqb ymm2,ymm2,ymm2
	vmovups ymm1,ymm0
	vgatherdps ymm1,dword ptr [rax+ymm8*4],ymm2
	vmulps ymm2,ymm5,ymm5
	vfmadd231ps ymm2,ymm4,ymm4
	vmovups ymm0,ymm2
	vfmadd231ps ymm0,ymm1,ymm1
	vmovups ymmword ptr [r8],ymm0
	sub edx,1
	jne dot3_aos_vector_gather+190h (07FF624C81470h)

view raw cppaosgather.asm hosted with ❤ by GitHub

Scalar SoA

	void dot3_soa_scalar(const vector<float>& xs, const vector<float>& ys, const vector<float>& zs, vector<float>& dp) {
	for (auto j = 0; j < reps; ++j) {
	auto i = vector_len;
	while(i--) {
	dp[i] = xs[i] * xs[i] + ys[i] * ys[i] + zs[i] * zs[i];
	}
	}
	}

view raw cppsoascalar.cpp hosted with ❤ by GitHub

	lea rax,[rax-4]
	vmovss xmm0,dword ptr [rax]
	vmovss xmm2,dword ptr [rdx+rax]
	vmovss xmm3,dword ptr [r8+rax]
	vmulss xmm1,xmm0,xmm0
	vmulss xmm0,xmm2,xmm2
	vaddss xmm2,xmm1,xmm0
	vmulss xmm1,xmm3,xmm3
	vaddss xmm2,xmm2,xmm1
	vmovss dword ptr [r9+rax],xmm2
	sub ecx,1
	jne dot3_soa_scalar+50h (07FF6F2F31330h)

view raw cppsoascalar.asm hosted with ❤ by GitHub

Auto-vectorized SoA

In order for the auto-vectorizer to kick in, we need to use a for-loop for the inner iteration:

	void dot3_soa_autovec(const vector<float>& xs, const vector<float>& ys, const vector<float>& zs, vector<float>& dp) {
	for (auto j = 0; j < reps; ++j) {
	for (auto i = 0; i < vector_len; ++i) {
	dp[i] = xs[i] * xs[i] + ys[i] * ys[i] + zs[i] * zs[i];
	}
	}
	}

view raw cppsoaautovec.cpp hosted with ❤ by GitHub

Now the compiler even generates FMA instructions:

	lea rax,[rcx+rbx]
	vmovups ymm2,ymmword ptr [r10+rcx]
	vmovups ymm3,ymmword ptr [rcx]
	vmovups ymm0,ymmword ptr [r11+rcx]
	lea rcx,[rcx+20h]
	vmulps ymm1,ymm2,ymm2
	vfmadd231ps ymm1,ymm3,ymm3
	vfmadd231ps ymm1,ymm0,ymm0
	vmovups ymmword ptr [rax+r8],ymm1
	sub r9,1
	jne dot3_soa_autovec+190h (07FF7B8661510h)

view raw cppsoaautovec.asm hosted with ❤ by GitHub

Vectorized SoA

Of course we can also vectorize manually by using intrinsics:

	void dot3_soa_vectorized(const vector<float>& xs, const vector<float>& ys, const vector<float>& zs, vector<float>& dp) {
	for (auto j = 0; j < reps; ++j) {
	const auto px = (__m256*)xs.data();
	const auto py = (__m256*)ys.data();
	const auto pz = (__m256*)zs.data();
	auto pd = (__m256*)dp.data();
	auto i = vector_len / lane_width;
	while (i--) {
	pd[i] = _mm256_add_ps(_mm256_add_ps(_mm256_mul_ps(px[i], px[i]), _mm256_mul_ps(py[i], py[i])), _mm256_mul_ps(pz[i], pz[i]));
	}
	}
	}

view raw cppsoavectorized.cpp hosted with ❤ by GitHub

	lea rax,[rax-20h]
	vmovups ymm0,ymmword ptr [rbx+rax]
	vmovups ymm2,ymmword ptr [rax]
	vmovups ymm3,ymmword ptr [rdi+rax]
	vmulps ymm1,ymm0,ymm0
	vmulps ymm0,ymm2,ymm2
	vaddps ymm2,ymm1,ymm0
	vmulps ymm1,ymm3,ymm3
	vaddps ymm2,ymm2,ymm1
	vmovups ymmword ptr [rsi+rax],ymm2
	sub r11d,1
	jne dot3_soa_vectorized+52h (07FF7183817D2h)

view raw cppsoavectorized.asm hosted with ❤ by GitHub

Vectorized SoA using FMA

This is the same as above, but it additionaly makes use of FMA:

	void dot3_soa_vectorized_fma(const vector<float>& xs, const vector<float>& ys, const vector<float>& zs, vector<float>& dp) {
	for (auto j = 0; j < reps; ++j) {
	const auto px = (__m256*)xs.data();
	const auto py = (__m256*)ys.data();
	const auto pz = (__m256*)zs.data();
	auto pd = (__m256*)dp.data();
	auto i = vector_len / lane_width;
	while (i--) {
	pd[i] = _mm256_fmadd_ps(pz[i], pz[i], _mm256_fmadd_ps(py[i], py[i], _mm256_mul_ps(px[i], px[i])));
	}
	}
	}

view raw cppsoafma.cpp hosted with ❤ by GitHub

	lea rax,[rax-20h]
	vmovups ymm0,ymmword ptr [rbx+rax]
	vmovups ymm2,ymmword ptr [rax]
	vmovups ymm3,ymmword ptr [rdi+rax]
	vmulps ymm0,ymm0,ymm0
	vfmadd231ps ymm0,ymm2,ymm2
	vfmadd231ps ymm0,ymm3,ymm3
	vmovups ymmword ptr [rsi+rax],ymm0
	sub r11d,1
	jne dot3_soa_vectorized+52h (07FF7DF9913D2h)

view raw cppsoafma.asm hosted with ❤ by GitHub

Results

The following figure displays the performance in GFLOP/s of the different versions. The dashed line is at 51.2 GFLOP/s, the theoretical peak performance of a single Haswell core (single precision):

First of all, both, all the AoS variants and the scalar SoA version, don't even come close to the vectorized SoA versions. Second, any attempts at accelerating the original AoS version failed (C#) or only provide insignificant performance gains (C++). Even vector gather can't save the day and in fact further impairs performance. In any event, the gains don't justify the more complicated code.

If you really need the performance SIMD can provide, you have to switch to a SoA layout: While Visual C++'s auto-vectorizer may relieve you of writing SIMD intrinsics directly, it still requires SIMD-friendly—that is: SoA—code. As long as it works, it provides the most accessible way of writing high-performance code. The second-best way, from a usability stand point, is probably C# and System.Numerics.Vectors, which enables (explicit) SIMD programming via a comparably easy-to-use interface.

Yet, the plot above also shows that non of the managed solutions is really able to keep up with any of the vectorized C++ versions. One reason for that is the inferior code generation of the JIT compiler compared to the C++ optimizer. Others are more intrinsic to the managed programming model (null-pointer checks, range checks). But also System.Numerics.Vectors is far from being complete: For instance, there is no support for FMA or scatter/gather operations. A "Vector8" type could help treating a float[] as AVX-sized chunks.

Conclusions

Want speed? Use a vectorized SoA approach. Want more speed? Use C++. That's what I learned from this little investigation. That and—until we have an AI driven optimizing compiler, that's smart enough to understand, what we are trying to achieve and then automatically transforms our code to the most efficient form—that writing high-performance code will remain both an art and hard work for the time being.