High-Performance SIMD in C++

Recently, I have been looking at related content on vectorization, and it seems a bit overwhelming. Taking this opportunity to learn about high-performance SIMD programming.

SIMD stands for Single Instruction, Multiple Data.

In traditional computer architecture, the CPU can only process one data element at a time. However, many tasks involve performing the same operation on a large amount of data, such as addition, multiplication, or logical operations on all elements of an array. SIMD programming provides the CPU with a specialized instruction set that allows it to perform the same operation on multiple data elements simultaneously.

This processing method is particularly suitable for calculations involving vectors, matrices, images, audio, and video data.

Currently, commonly used SIMD instruction sets include SSE, SSE2, AVX128, AVX256, and AVX512.

In this section, we will briefly learn about some operations of AVX512. There are many operations, so we will only introduce a few.

1. Terminology

The first question is that SIMD programming code looks quite different from the usual code, with various underscores and naming conventions that are completely incomprehensible. How do we understand it?

For example:

_mm512_set1_ps

Broadcast single-precision (32-bit) floating-point value a to all elements of dst.

_mm512_set1_epi32

Broadcast 32-bit integer a to all elements of dst.

Thus, I found the following table:

Abbreviation	Full Name	C/C++ Equivalent
ps	Packed Single-Precision	float
ph	Packed Half-Precision	None*
pd	Packed Double-Precision	double
pch	Packed Half-Precision Complex	None*
pi8	Packed 8-Bit Integer	int8_t
pi16	Packed 16-Bit Integer	int16_t
pi32	Packed 32-Bit Integer	int32_t
epi8	Extended Packed 8-Bit Integer	int8_t
epi16	Extended Packed 16-Bit Integer	int16_t
epi32	Extended Packed 32-Bit Integer	int32_t
epi64	Extended Packed 64-Bit Integer	int64_t
epi64x	Extended Packed 64-Bit Integer	int64_t

https://stackoverflow.com/questions/70911872/what-are-the-names-and-meanings-of-the-intrinsic-vector-element-types-like-epi6

For instance:

_mm512_mask_load_ps

_mm512_mask_loadu_ps

Where ‘u’ indicates unordered, meaning that when using _mm512_mask_loadu_ps to load data from memory, no alignment requirements are imposed on the memory address. In contrast, _mm512_mask_load_ps requires a 64-byte alignment.

By comparing and learning this way, one can quickly understand the meaning of each interface.

For relevant APIs, you can check the Intel Intrinsics Guide.

https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html

2. Practical Examples

2.1 Finding Minimum Value

For 512 bits, we can store 16 32-bit floats.

static inline rf_512 load(float* ptr) { return _mm512_loadu_ps(ptr); }
float a[width] = {1.0, 2.0,  3.0,  22.0, 5.0,  17.0, 9.0,  8.0,
                  9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0};
float b[width] = {3.3, 6.2,   5.3,   4.4,   5.5,   6.6,   7.7,  8.8,
                  9.9, 10.10, 21.11, 12.12, 13.13, 14.14, 15.0, 16.16};
rf_512 data = minimum(load(a), load(b));

Thus we can quickly obtain:

1 2 3 4.4 5 6.6 7.7 8 9 10 11 12 13 14 15 16

2.2 Quick Data Shuffling

For the input data:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

We can quickly obtain:

16.16 2 14.14 13.13 5 6 7 8 9 10 11 12 13 14 15 16

The corresponding implementation:

static inline rf_512 permutexvar(ri_512 idx, rf_512 src) {
  return _mm512_permutexvar_ps(idx, src);
}
/*
raw data: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
shuffle data: 16.16 2 14.14 13.13 5 6 7 8 9 10 11 12 13 14 15 16
*/
void print_permutexvar_mask() {
  float a[width] = {1.0, 2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,
                    9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0};
  float b[width] = {1.1, 2.2,   3.3,   4.4,   5.5,   6.6,   7.7,  8.8,
                    9.9, 10.10, 11.11, 12.12, 13.13, 14.14, 15.0, 16.16};
  mask_type mask = make_bit_mask<1, 0, 1, 1, 0>();
  int idx_array[width] = {15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0};
  ri_512 idx = _mm512_loadu_si512((__m512i*)idx_array);
  rf_512 result = permutexvar_mask(load(a), mask, idx, load(b));
  float result_arr[width];
  store(result, result_arr);
  std::cout << "raw data: ";
  for (int i = 0; i < width; i++) {
    std::cout << a[i] << " ";
  }
  std::cout << std::endl;

  std::cout << "shuffle data: ";
  for (int i = 0; i < width; i++) {
    std::cout << result_arr[i] << " ";
  }
  std::cout << std::endl;
}

2.3 Rotation

For an array:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

We can rotate to obtain:

4 5 6 7 8 9 10 11 12 13 14 15 16 1 2 3

14 15 16 1 2 3 4 5 6 7 8 9 10 11 12 13

This means rotating down or rotating up.

Here, we can implement it like this:

rf_512 rotateGeneral(float* arr, int s) {
  int idx_array[width];
  for (int i = 0; i < width; i++) {
    idx_array[i] = (i + s) % width;
  }
  ri_512 idx = _mm512_loadu_si512((__m512i*)idx_array);
  rf_512 result = permutexvar(idx, load(arr));
  return result;
}

And there are other examples; we can find that by using SIMD, we can implement some very interesting algorithms to accelerate the processing of arrays and bulk data.

We will continue to learn SIMD together, let’s keep it up!

Note: The complete example of this section has been updated on the platform, thank you!