Low level optimization of parallel algorithms or SIMD in .NET

Currently, a huge number of tasks require high system performance. Physical limitations do not allow infinitely increasing the number of transistors on the processor chip. The geometrical dimensions of transistors cannot be physically reduced, since when exceeding the possible allowable sizes, phenomena begin to appear that are not noticeable with large sizes of active elements - quantum size effects begin to strongly affect. Transistors do not start to work like transistors.
And Moore’s law has nothing to do with it. This was and remains the law of value, and an increase in the number of transistors on a chip is more likely a consequence of the law. Thus, in order to increase the power of computer systems, one has to look for other ways. This is the use of multiprocessors, multicomputers. This approach is characterized by a large number of processor elements, which leads to independent execution of subtasks on each computing device.

Parallel processing methods:

Source of concurrency	Acceleration	Programmer effort	Popularity
Many cores	2x-128x	Moderate	High
Many cars	1x Infinity	Moderately High	High
Vectorization	2x-8x	Moderate	Low
Graphics Adapters	128x-2048x	High	Low
Coprocessor	40x-80x	Moderately high	Extremely low

There are many ways to improve the efficiency of systems and all are quite different. One of such methods is the use of vector processors, which significantly increase the speed of calculations. Unlike scalar processors, which process one data element per instruction (SISD), vector processors are capable of processing several data elements (SIMD) in one instruction. Most modern processors are scalar. But many of the tasks they solve require a large amount of computation: processing video, sound, graphics, scientific calculations, and much more. To speed up the computation process, processor manufacturers began to integrate additional streaming SIMD extensions into their devices.
Accordingly, with a certain programming approach, it became possible to use vector processing of data in the processor. Existing Extensions: MMX, SSE, and AVX. They allow you to use additional processor capabilities for accelerated processing of large data arrays. At the same time, vectorization allows acceleration without obvious parallelism. Those. it exists from the point of view of data processing, but from the point of view of the programmer it does not require any expenses for the development of special algorithms to prevent the state of a race or synchronization, and the development style does not differ from synchronous. We get acceleration without much effort, almost completely free. And there is no magic in it.

What is SSE?

SSE (Eng. Streaming SIMD Extensions, streaming SIMD-extension of the processor) is a SIMD (Eng. Single Instruction, Multiple Data, One instruction - a lot of data) a set of instructions. The SSE includes eight 128-bit registers and a set of instructions in the processor architecture. SSE technology was first introduced in Pentium III in 1999. Over time, this instruction set has been improved by adding more complex operations. Eight (in x86-64 - sixteen) 128-bit registers were added to the processor: from xmm0 to xmm7.

Initially, these registers could be used only for single precision calculations (i.e., for the float type). But after the release of SSE2, these registers can be used for any primitive data type. Given a standard 32-bit machine this way, we can store and process in parallel:

• 2 double
• 2 long
• 4 float
• 4 int
• 8 short
• 16 char

If you use the AVX technology, then you will be manipulating already 256-bit registers, respectively, more numbers in one instruction. So there are already 512-bit registers.

First, using C ++ as an example (who are not interested, you can skip), we will write a program that will sum two arrays of 8 elements of type float.

C ++ vectorization example

SSE technology in C ++ is implemented by low-level instructions, presented in the form of pseudo-code, which reflect assembly instructions. So, for example, the command __m128 _mm_add_ps (__ m128 a, __m128 b); converted to the assembler instruction ADDPS operand1, operand2 . Accordingly, the command __m128 _mm_add_ss (__ m128 a, __m128 b); will be converted to the ADDSS instruction operand1, operand2 . These two commands do almost the same thing: they add up the elements of the array, but in slightly different ways. _mm_add_ps is completely register-case, so:

• r0: = a0 + b0
• r1: = a1 + b1
• r2: = a2 + b2
• r3: = a3 + b3

Moreover, the entire register __m128 is the set r0-r3. But the _mm_add_ss command adds up only part of the register, so:

• r0: = a0 + b0
• r1: = a1
• r2: = a2
• r3: = a3

Other commands are arranged on the same principle, such as subtraction, division, square root, minimum, maximum and other operations.
To write a program, you can manipulate 128-bit registers like __m128 for float, __m128d for double and __m128i for int, short, char. At the same time, you can not use arrays of type __m128, but use the given pointers of the array float to type __m128 *.
In this case, several working conditions should be taken into account:

• Float data loaded and stored in an __m128 object must have 16-byte alignment
• Some built-in functions require their argument to be of type integer constants, due to the nature of the instruction
• The result of arithmetic operations acting on two NAN arguments is not defined

Such a small digression into theory. However, consider an example program using SSE:

#include "iostream"
#include "xmmintrin.h"	
int main()
{
	const auto N = 8;
	alignas(16) float a[] = { 41982.0,  81.5091, 3.14, 42.666, 54776.45, 342.4556, 6756.2344, 4563.789 };
	alignas(16) float b[] = { 85989.111,  156.5091, 3.14, 42.666, 1006.45, 9999.4546, 0.2344, 7893.789 };
	__m128* a_simd = reinterpret_cast<__m128*>(a);
	__m128* b_simd = reinterpret_cast<__m128*>(b);
        auto size = sizeof(float);
        void *ptr = _aligned_malloc(N * size, 32);
	float* c = reinterpret_cast(ptr);
        for (size_t i = 0; i < N/2; i++, a_simd++, b_simd++, c += 4)
		_mm_store_ps(c, _mm_add_ps(*a_simd, *b_simd));
	c -= N;
	std::cout.precision(10);
	for (size_t i = 0; i < N; i++)
		std::cout << c[i] << std::endl;
	_aligned_free(ptr);
	system("PAUSE");
	return 0;
}

• alignas (#) is a standard C ++ portable way of setting custom alignment of variables and user types. It is used in C ++ 11 and supported by Visual Studio 2015. You can use another option - __declspec (align (#)) declarator . Data management tools for alignment during static memory allocation. If alignment with dynamic selection is needed, use void * _aligned_malloc (size_t size, size_t alignment);
• Then convert the pointer to the array a and b to type _m128 * using reinterpret_cast, which allows you to convert any pointer to a pointer of any other type.
  
• After, dynamically allocate the aligned memory using the function _aligned_malloc (N * sizeof (float), 16) already mentioned above ;
• We select the number of necessary bytes based on the number of elements, taking into account the type dimension, and 16 is the alignment value, which should be a power of two. And then we bring the pointer to this piece of memory to another type of pointer, so that we could work with it taking into account the dimension of the float type as an array.
  

Thus, all preparations for SSE work are completed. Further in the loop we summarize the elements of arrays. The approach is based on pointer arithmetic. Since a_simd , b_simd and c are pointers, their increase leads to an offset by sizeof (T) from memory. If we take the dynamic array c for example , then c [0] and * c will show the same value, because c points to the first element of the array. Increment withwill cause the pointer to shift 4 bytes forward and now the pointer will point to 2 elements of the array. Thus, you can move forward and backward through the array by increasing and decreasing the pointer. But at the same time, the dimension of the array should be taken into account, since it is easy to go beyond its boundaries and turn to someone else's memory location. Working with the a_simd and b_simd pointers is similar, only incrementing the pointer will lead to a 128-bit advance, and from the point of view of type float, 4 variables of the array a and b will be skipped. In principle, the pointer a_simd and a , as well as b_simd and b indicate, respectively, one section in memory, with the exception that they are processed differently, taking into account the dimension of the type of pointer:

	for (int i = 0; i < N/2; i++, a_simd++, b_simd++, c += 4)
		_mm_store_ps(c, _mm_add_ps(*a_simd, *b_simd));

Now it’s clear why there are such pointer changes in this loop. At each iteration of the cycle, 4 elements are added and the result is saved at the addresses of the pointer c from the register xmm0 (for this program). Those. as we see this approach does not change the source data, but stores the amount in the register and, if necessary, transfers it to the operand we need. This allows you to improve program performance in cases where you need to reuse operands.
Consider the code that the assembler generates for the _mm_add_ps method :

mov         eax,dword ptr [b_simd]  ;// поместить адрес b_simd в регистр eax(базовая команда пересылки данных, источник не изменяется)
mov         ecx,dword ptr [a_simd]  ;// поместить адрес a_simd в регистр ecx 
movups      xmm0,xmmword ptr [ecx]  ;// поместить 4 переменные с плавающей точкой по адресу ecx в регистр xmm0;            xmm0 = {a[i], a[i+1], a[i+2], a[i+3]}
addps       xmm0,xmmword ptr [eax]  ;// сложить переменные: xmm0 = xmm0 + b_simd
                                    ;// xmm0[0] = xmm[0] + b_simd[0]
                                    ;// xmm0[1] = xmm[1] + b_simd[1]
                                    ;// xmm0[2] = xmm[2] + b_simd[2]
                                    ;// xmm0[3] = xmm[3] + b_simd[3]
movaps      xmmword ptr [ebp-190h],xmm0 ;// поместить значение регистра по адресу в стеке со смещением
movaps      xmm0,xmmword ptr [ebp-190h] ;// поместить в регистр
mov         edx,dword ptr [c] ;// поместить в регистр ecx адрес переменной с
movaps      xmmword ptr [edx],xmm0 ;// поместить значение регистра в регистр ecx или сохранить сумму по адресу памяти, куда указывает (ecx) или с. При этом xmmword представляет собой один и тот же тип, что и _m128 - 128-битовое слово, в   котором 4 переменные с плавающей точкой 

As you can see from the code, one addps instruction processes 4 variables at once, which is implemented and supported by the hardware processor. The system does not take any part in the processing of these variables, which gives a good performance boost without unnecessary external costs.
In this case, I would like to note one feature that in this example the compiler uses the movups instruction , which does not require operands that must be aligned on a 16-byte boundary. From which it follows that one could not align array a . However, array b needs to be aligned, otherwise in the addps operationan error reading the memory will occur, because the register is added with a 128-bit memory location. There may be other instructions in another compiler or environment, therefore it is better in any case for all operands involved in such operations to do border alignment. In any case, to avoid memory problems.
Another reason to do alignment is when we are working with elements of arrays (and not only with them), we actually constantly work with cache lines of 64 bytes in size. SSE and AVX vectors always fall into the same cache line if they are aligned on 16 and 32 bytes, respectively. But if our data is not aligned, then, very likely, we will have to load another "additional" cache line. This process affects the performance quite strongly, and if at the same time we turn to the elements of the array, and therefore to the memory, inconsistently, then everything can be even worse.

SIMD support in .NET

The first mention of JIT support for SIMD technology was announced on the .NET blog in April 2014. Then the developers announced a new preview version of RyuJIT, which provided SIMD functionality. The reason for the addition was the rather high popularity of the support request for C # and SIMD. The initial set of supported types was not large and there were limitations on functionality. Initially, the SSE set was supported, and AVX was promised to be added in the release. Later, updates were released and new types with SIMD support and new methods for working with them were added, which in recent versions represents an extensive and convenient library for hardware data processing.

This approach makes life easier for a developer who does not have to write CPU-dependent code. Instead, the CLR abstracts the hardware by providing a virtual runtime environment that translates its code into machine instructions either at run time (JIT) or during installation (NGEN). Leaving CLR code generation, you can use the same MSIL code on different computers with different processors, without giving up optimizations specific to this particular CPU.
Currently, support for this technology in .NET is represented in the System.Numerics.Vectors namespace and is a library of vector types that can take advantage of SIMD hardware acceleration. Hardware acceleration can lead to a significant increase in performance in mathematical and scientific programming, as well as in graphics programming. It contains the following types:

• Vector - a collection of static convenience methods for working with universal vectors
• Matrix3x2 - represents a 3x2 matrix
• Matrix4x4 - represents a 4x4 matrix
• Plane - represents a three-dimensional plane
• Quaternion - represents a vector used to encode three-dimensional physical rotations
• Vector <(Of <(<'T>)>)> represents a vector of the specified numeric type, which is suitable for low-level optimization of parallel algorithms
• Vector2 - represents a vector with two single precision floating point values
• Vector3 - represents a vector with three single precision floating point values
• Vector4 - represents a vector with four single precision floating point values

The Vector class provides methods for adding, comparing, finding the minimum and maximum, and many other transformations over vectors. At the same time, operations operate using SIMD technology. Other types also support hardware acceleration and contain transformations specific to them. For matrices, this can be multiplication, for vectors, Euclidean distance between points, etc.

C # program example

So, what is needed in order to use this technology? You must first have a RyuJIT compiler and .NET version 4.6. System.Numerics.Vectors via NuGet is not installed if the version is lower. However, even with the library installed, I downgraded the version and everything worked as it should. Then you need to build under x64, for this you need to remove “prefer 32-bit platform” in the project properties and can be built under Any CPU.
Listing:

using System;
using System.Numerics;
    class Program
    {
        static void Main(string[] args)
        {
            const Int32 N = 8;
            Single[] a = { 41982.0F, 81.5091F, 3.14F, 42.666F, 54776.45F, 342.4556F, 6756.2344F, 4563.789F };
            Single[] b = { 85989.111F, 156.5091F, 3.14F, 42.666F, 1006.45F, 9999.4546F, 0.2344F, 7893.789F };
            Single[] c = new Single[N];
            for (int i = 0; i < N; i += Vector.Count) // Count возвращает 16 для char, 4 для float, 2 для      double и т.п.
            {
                var aSimd = new Vector(a, i); // создать экземпляр со смещением i
                var bSimd = new Vector(b, i);
                Vector cSimd = aSimd + bSimd; // или так Vector c_simd = Vector.Add(b_simd, a_simd);
                cSimd.CopyTo(c, i); //копировать в массив со смещением
            }
            for (int i = 0; i < a.Length; i++)
            {
                Console.WriteLine(c[i]);
            }
            Console.ReadKey();
        }
    }

From a general point of view, that C ++ approach that .NET are pretty similar. It is necessary to convert / copy the source data, copy to the final array. However, the approach with C # is much simpler, many things are done for you and you just have to use and enjoy. There is no need to think about data alignment, allocate memory and do it statically or dynamically with certain operators. On the other hand, you have more control over what is happening using pointers, but also more responsibility for what is happening.
And in the loop everything happens the same as in the loop in C ++. And I'm not talking about pointers. The calculation algorithm is the same. At the first iteration, we put the first 4 elements of the source arrays into the aSimd and bSimd structures, then sum and store in the final array. Then, at the next iteration, we enter the next 4 elements using the offset and sum them up. That's how it is done quickly and easily. Consider the code that the compiler generates for this command var cSimd = aSimd + bSimd :

addps       xmm0,xmm1  

The difference from the C ++ version is that both registers are added here, while there was folding of the register with a piece of memory. Placement in registers occurs during initialization of aSimd and bSimd . In general, this approach, when comparing the code of C ++ and .NET compilers, is not very different and gives approximately equal performance. Although the option with pointers will still work faster. I would like to note that SIMD instructions are generated when code optimization is enabled. Those. it will not work to see them in the disassembler in Debug: this is implemented as a function call. However, in Release, where optimization is enabled, you will receive these instructions in an explicit (built-in) form.

In the end

What we have:

• In many cases, vectorization gives a 4-8 × increase in performance
• Sophisticated algorithms will require ingenuity, but without it, nowhere
• System.Numerics.Vectors currently covers only part of the simd instructions. A more serious approach would require C ++
• There are many other ways besides vectorization: the correct use of the cache, multithreading, heterogeneous computing, competent work with memory (so that the garbage collector does not sweat), etc.

During a brief twitter correspondence with Sasha Goldstein (one of the authors of the book “Application Optimization on the .NET Platform”), who examined hardware acceleration in .NET, I asked what SIMD support was implemented in .NET and what it was in Comparison with C ++. To which he replied: “Undoubtedly, you can do more in C ++ than in C #. But you really get cross-processor support in C #. For example, the automatic choice between SSE4 and AVX. " In general, this is good news. With little effort, we can get as much performance from the system as possible, utilizing all possible hardware resources.
For me, this is a very good opportunity to develop productive programs. At least in my thesis on modeling physical processes, I basically achieved efficiency by creating a certain number of threads, as well as using heterogeneous calculations. I use both CUDA and C ++ AMP. The development is carried out on a universal platform for Windows 10, where I am very attracted to WinRT, which allows me to write a program in both C # and C ++ / CX. Basically, on the pros I write the kernel for large calculations (Boost), but in C # I already manipulate the data and develop the interface. Naturally, data transfer via the binary ABI interface for the interaction of the two languages ​​has its price (although not very large), which requires a more reasonable development of the C ++ library.
If you need to manipulate data in C #, I will convert them to .NET types so as not to work with WinRT types, thereby increasing processing performance already in C #. For example, when you need to process several thousand or tens of thousands of elements or processing requirements do not have any special specifications, data can be calculated in C # without using a library (from 3 to 10 million instances of structures are counted in it, sometimes in just one iteration). So a hardware accelerated approach will simplify the task and make it faster.

List of sources when writing an article

• Translation of a .NET blog article about a new SIMD-enabled JIT
• About System.Numerics.Vector on the .NET Blog
• What is alignment and how does it affect your programs?
• Intel Blog on Data Alignment in Loop Vectoring

Well, special thanks to Sasha Goldstein for the help and information provided.