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

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

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:

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:

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:

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

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:

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:

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:

The resulting machine code is likewise similar to scalar AoS:

Vectorized SoA

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

This is what the compiler makes of it:

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:

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:

AoS vectorized across each 3-vector

Let's use the DPPS instruction via intrinsics:

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:

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:

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

Scalar SoA

Auto-vectorized SoA

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

Now the compiler even generates FMA instructions:

Vectorized SoA

Of course we can also vectorize manually by using intrinsics:

Vectorized SoA using FMA

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

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.

SIMD Fundamentals. Part II: AoS, SoA, Gather/Scatter - Oh my!

Last time we looked at a very basic example of a data parallel problem: adding two arrays. Unfortunately, data parallelism is not always so easy to extract from existing code bases and often requires considerable effort. Today we will therefore address a slightly more complicated problem and move step by step from a pure scalar to a fully vectorized version.

The problem and its SISD solution

Imagine we wanted to compute the squared ℓ²-norm of a set of m 3-vectors u in R³, what is really just the dot product of each vector x with itself: ‖x‖² = x² + y² + z²

In good OO fashion, you'd probably first define a Vec3f data type, in its simplest form like so (F# syntax)

Our vector field u can then simply be represented by an array of these structs and the code that performs the computation of the dot products is

Load the three components, square them, add up the squares and store the result in dp:

 

For SSE2, that's a total of three MULSS and two ADDSS instructions, so the one dot product we get out is worth 5 FLOPs. Now we would of course like to make use of our shiny SIMD hardware to accelerate the computation. There are a few options of how we can approach this:

Pseudo-SIMD: Bolting SIMD on top of AoS

Our current data layout is an "array of structures" (AoS). As we shall see further below, AoS isn't a particular good choice for vectorization and SIMDifying such code won't yield peak performance. Yet it can be done and depending on the used SIMD instruction set it may even provide some speed-up.

Here's the idea: Instead of loading the three scalar components into separate scalar registers, we store all of them in a single SIMD register. In case of SSE2, n = 4 for single precision floats, so one of the elements in the register is a dummy value that we can ignore.  Now we can multiply the SIMD register with itself to get the squared components (plus the ignored fourth) one. Then we need to horizontally add up the three squares that we are interested in; in pseudo-F#:

Graphically:

Although we were able to replace three MULSS with a single MULPS and therefore express 5 FLOPs through only 3 arithmetic instructions, there are multiple issues with this idea: 

1. We need three scalar loads plus shuffling operations to place the three components in a single SIMD register.
2. We effectively only use 3 of the 4 vector lanes, as we operate on Vec3f, not Vec4f structs. This problem becomes even more severe with wider vector paths.
3. We use vector operations only for multiplication, but not for addition.
4. Horizontal operations like adding the components are comparatively slow, require shuffling the components around and extracting a single scalar (the result).
5. We still only compute one dot product per iteration.

Another example for this kind of pseudo-vectorization can be found in Microsoft's SIMD sample pack (ray tracer sample). Yet, for all the reasons mentioned above, you should avoid this approach whenever possible.

On-the-fly vectorization

If you recall our "Hello, world!" example from Part I, then you may also remember that vectorizing our loop meant that we computed n results per iteration instead of one. And that is what we have to achieve for our current dot-product example as well: We need to compute n dot products in each iteration. How can we accomplish this?

Imagine our data (array of Vec3f structs) to form a 3×m matrix (three rows, m columns). Each column represents a Vec3f, each row the x, y or z components of all Vec3fs. In the previous section, we tried to vectorize along the columns by parallelizing the computation of a single dot product—and largely failed.

The reason is that the data parallelism in this example can really be found along the rows of our imagined matrix, as each dot product can be computed independently and thus in parallel from each other. Furthermore m, the number of Vec3fs, is typically much larger then n and so we don't face problems in utilizing the full SIMD width. The wider the SIMD vector is, the more dot products we can compute per iteration.

As with our example from last time, the vectorized version of the algorithm is actually very similar to the plain scalar one. That's why vectorizing a data-parallel problem isn't really hard once you get the idea. The only difference is that we don't handle individual floats, but chunks of n floats:

1. Load n x values and square them
2. Load n y values and square them
3. Load n z values and square them
4. Add the n x, y and z values
5. Store the n resulting dot products

In pseudo-F# the algorithm is now

Graphically:

That good part about this version is that it uses vector arithmetic instructions throughout, performing n times the work of its scalar counterparts, performing a total of 20 FLOPs each iteration.

The bad part is the one labeled "scalar loads & shuffling" in the picture above: Before we can compute the n results, we have to gather n x, y and z values, but our data is still laid out in memory as an AoS, i.e. ordered [xyzxyzxyzxyzxyz...]. Loading logically successive [x0 x1 x2 x3 ...] values thus requires indexed/indirect load instructions (vector gather). SIMD instruction sets without gather/scatter support, like SSE2, have to load the data using n conventional scalar loads and appropriately place them in the vector registers, hurting performance considerably.

Approaching nirvana: Structure of arrays (SoA)

To avoid this drawback, maybe we should just store the data in a SIMD-friendly way, as a Structure of Arrays (SoA): Instead of modeling the vector field u as a number of Vec3f structs, we consider it to be an object of three float arrays:

This is indeed very similar to switching from a column-major to row-major matrix layout, because it changes our data layout from [xyzxyzxyz...]. to [xxxxx...yyyyy...zzzzz...].

Observe how this implementation differs from the previous one only in that we index the components instead of the vectors:

And yet it saves us a lot of loads and shuffling, resulting in pure, 100% vectorized SIMD code:

As this version really only replaces the scalar instructions with vector equivalents but executes 20 FLOPs in each iteration, we should now indeed get about 4x the performance of the scalar version and a much more favourable arithmetic-to-branch/load/store ratio.

In fact, once the hard part of the work is done (switching from AoS to an SoA data layout), many compilers can even automatically vectorize loops operating on that kind of data.

Conclusions

Vectorizing OO code needs some getting used to, the SoA view of the world may be a bit alien to non-APL programmers. If you can foresee the need for vectorizing your code at some point in the future, it may be wise to use an SoA layout from the get-go instead of having to rewrite half of your code later on. Experience certainly helps in identifying data-parallel problems; but as a rule of thumb good candidates for extracting data parallelism are those problems where thinking in terms of operations that apply to a whole array/field of items comes naturally.

In part III we are going to discuss concrete implementations of the above concepts using F#, RyuJIT and System.Numerics.Vectors and compare the performance of the different versions—once they sorted out this issue.

SIMD Fundamentals. Part I: From SISD to SIMD

For a long time, .NET developers didn't have to care about SIMD programming, as it was simply not available to them (except maybe for Mono's Mono.Simd). Only recently Microsoft introduced a new JIT compiler RyuJIT for .NET that, in conjunction with a special library System.Numerics.Vectors, offers access to the SIMD hardware of modern CPUs. The goal of this three-part series of articles is therefore to introduce the fundamental ideas of SIMD to .NET developers who want to make use of all the power today's CPUs provide.

Parallelism in modern CPUs

In recent years CPUs gained performance mainly by increasing parallelism on all levels of execution. The advent of x86-CPUs with multiple cores established true thread level parallelism in the PC world. Solving the "multi-core problem", the question of how to distribute workloads between different threads in order to use all that parallel hardware goodness—ideally (semi-)automatically, suddenly became and continues to be one of the predominant goals of current hard- and software related research.

Years before the whole multi-core issue started, CPUs already utilized another kind of parallelism to increase performance: Instruction level parallelism (ILP) was exploited via pipelining, super-scalar and out-of-order execution and other advanced techniques. The nice thing about this kind of parallelism is that your average Joe Developer doesn't has to think about it, it's all handled by smart compilers and smart processors.

But then in 1997 Intel introduced the P55C and with it a new 64-bit-wide SIMD instruction set called MMX ("multimedia extension"; after all, "multimedia" was "the cloud" of the 90s). MMX made available a level of parallelism new to most PC developers: data level parallelism (DLP). Contrary to ILP however, DLP requires specifically designed code to be of any good. This was true for MMX and it remains to be true for its successors like Intel's 256-bit-wide AVX2. Just as with multi threading, programmers need to understand how to use those capabilities properly in order to exploit the tremendous amounts of floating point horse power.

From SISD to SIMD

Whenever I learn new concepts, I like to first think of examples and generalize afterwards (inductive reasoning). In my experience that's true for most people, so let us therefore start with the "Hello, world!" of array programming: Suppose you have two floating point arrays, say xs and ys of length m and, for some reason, you want to add those arrays component-wise and store the result in an array zs. The conventional "scalar" or SISD (Single Instruction Single Data, see Flynn's taxonomy) way of doing this looks like this (C# syntax)

for (var i = 0; i < m; ++i) {
    zs[i] = xs[i] + ys[i];
}

Fetch two floats, add them up and store the result in zs[i]. Easy:

For this specific example, adding SIMD (Single Instruction Multiple Data) is almost trivial: For simplicity, let us further assume that m is evenly divisible by n, the vector lane width. Now instead of performing one addition per iteration, we add up xs and ys in chunks of n items, in pseudo-C# with an imagined array range expression syntax

for (var i = 0; i < m / n; i += n) {
    zs[i:(i+n-1)] = xs[i:(i+n-1)] + ys[i:(i+n-1)];
}

We fetch n items from xs and n items from ys using vector load instructions add them up using a single vector add and store the n results in zs. Compared to the scalar version, we need to decode fewer instructions and perform n times the work in each iteration. Think of each SIMD register as a window, n values wide into an arbitrarily long stream (array) of scalar values. Each SIMD instructions modifies n scalar values of that stream at once (that's where the speed-up comes from) and moves on to the next chunk of n values:

So SIMD really is just a form of array processing, where you think in terms of applying one operation (single instruction) to a lot of different elements (multiple data), loop-unrolling on steroids. And that's already really all you need to know about SIMD, basically.

Yet, this example is perhaps the most obvious data parallel case one could think of. It gets a whole lot more interesting once you add OOP and the data layout this paradigm propagates to the mix. We will look into the issue of vectorizing typical OOP code in the next part of this series.

NativeInterop 2.4.0

NativeInterop v2.4.0 is ready and awaiting your download from NuGet. This version brings a revised version of Buffer.Copy tuned for .NET 4.6/RyuJIT: By using a 16 byte block copy, the JITter can generate movq (x86) or movdqu (x64) instructions for further improved performance. Enjoy!

Now I just have to figure out, where the larger performance drop for large data sizes comes from compared to the other methods. memcpy somehow reaches approx. 20 GB/s beyond the L3 cache.

Btw., I have a further post on SIMD programming with .NET almost ready, but I can't publish it yet due to problems with the current version of System.Numerics.Vectors in combination with .NET 4.6/VS 2015 RC. Hope that'll get fixed soon!