Portable Vectorization with xsimd

In the previous lesson, we opened up the CPU core to look at its wide vector registers. We saw how SIMD (Single Instruction, Multiple Data) allows the processor to crunch 4, 8, or even 32 numbers in a single clock cycle.

However, we also saw that relying on the compiler to automatically vectorize our C++ loops is fragile. A simple if statement, an early return, or the slightest ambiguity about pointer aliasing will scare the optimizer into falling back to slow, sequential scalar instructions.

When we are building high-performance data structures like the Swiss Table later in this chapter, we cannot cross our fingers and hope -O3 does what we want. Our probe loops are full of branches and early exits, which completely paralyze auto-vectorizers. We must take direct, mechanical control of the vector registers ourselves.

In this lesson, we will move from passive auto-vectorization to explicit vectorization. We will integrate the xsimd library into our project, learn how to pack data into hardware batches, and build a vectorized search algorithm capable of scanning 32 elements in a single hardware instruction.

Before we introduce xsimd, let's look at what happens when we try to command vector registers directly using compiler intrinsics.

Every CPU manufacturer provides C-style header files that expose their specific hardware instructions directly to our code. On x86 processors (Intel and AMD), we include <immintrin.h>. On ARM processors (Apple Silicon, mobile devices, AWS Graviton servers), we include <arm_neon.h>.

However, when we try to use these intrinsics, our code becomes incredibly difficult to both read and write:

1#include <immintrin.h>
2#include <cstdint>
3
4uint32_t FindFortyTwo_AVX2(const int8_t* data) {
5  __m256i target = _mm256_set1_epi8(42);
6
7  __m256i chunk = _mm256_loadu_si256(
8    reinterpret_cast<const __m256i*>(data)
9  );
10  
11  __m256i cmp = _mm256_cmpeq_epi8(chunk, target);
12
13  return static_cast<uint32_t>(_mm256_movemask_epi8(cmp));
14}

At a hardware level, this code is fast, but it is incredibly cryptic. The naming conventions read like serial numbers. epi8 stands for "extended packed 8-bit integers," and si256 stands for "256-bit signed integer." Standard arithmetic operators like + or == don't work on __m256i types, forcing us to wrap every basic operation in arcane function calls.

Worse, it locks our code to a specific set of CPUs. If we try to compile this code on an M-series MacBook or an ARM server, the compiler will throw an error. ARM chips do not have the required registers or support AVX2 instructions. Even within the same architecture, an older processor may not support these instructions and will also break.

Enter xsimd

To solve this without losing the performance opportunties of SIMD, the C++ ecosystem relies on abstraction libraries. The two most popular options in modern C++ are xsimd (maintained by QuantStack), and Google's Highway library.

We'll usexsimd here - a zero-overhead, header-only template library. It inspects our compiler flags, detects the widest vector registers our target CPU supports, and exposes a clean, unified C++ type called xsimd::batch. Behind the scenes, this adapts to the most appropriate set of instructions for our target.

When we write batchA == batchB, xsimd overloads the == operator. On an x86 machine compiled with AVX2, that == compiles down to the _mm256_cmpeq_epi8() instruction we used above. If we update our build to target a different family of CPUs, the exact same line of C++ will compile down to the equivalent instructions for that target, requiring no changes to our code.

Remember our zero-cost abstractions philosophy from the very first lesson: at runtime, xsimd completely melts away, leaving behind the same assembly as hand-written intrinsics.

Because xsimd is a header-only library, adding it to our project is straightforward - we just grab the code from its GitHub repository and #include it.

If you're using the CMake build system from our introductory chapter, we can automate this process in the same way we did for Google Benchmark. We can download it using FetchContent and then link it using target_link_libraries:

1# ...
2
3# Fetch xsimd for portable vectorization
4include(FetchContent)
5FetchContent_Declare( 
6  xsimd 
7  GIT_REPOSITORY https://github.com/xtensor-stack/xsimd.git 
8  GIT_TAG 13.0.0 
9) 
10FetchContent_MakeAvailable(xsimd) 
11
12target_link_libraries(dsa_core INTERFACE xsimd) 
13
14# ...

To write SIMD code, it's helpful to stop thinking about variables as single values and start thinking in registers and lanes.

A hardware vector register is simply a fixed-width block of bits sitting inside the CPU core. For example, an AVX2 register is physically 256 bits (32 bytes) wide, while an SSE or ARM NEON register is 128 bits (16 bytes) wide.

Those 256 bits don't care what type of data you put inside them. We can subdivide a 256-bit register into different numbers of lanes depending on the size of our C++ data type.

For example, if we wanted to store 64-bit values like double or int64_t, we'd have 4 lanes ($256/64 = 4$), meaning the register can hold and operate on 4 elements at once. If we use smaller types, we have more lanes and can process more of them in a single clock cycle. For example:

• For 32-bit types like float or int32_t, $256/32 = 8$ lanes, so the register holds 8 elements at once.
• For 16-bit types like int16_t, $256/16 = 16$ lanes, so the register holds 16 elements at once.
• For 8-bit types like int8_t or uint8_t, $256/8 = 32$ lanes, so the register holds 32 elements at once.

Notice the inverse relationship here: the smaller we make our data types, the more items we can process in a single clock cycle. This is yet another reason to ensure our data types are not larger than they need to be, beyond the usual memory and cache efficiency considerations.

For example, if we store some metadata as a 64-bit integer, a 256-bit register only processes 4 items per instruction. If we compress that metadata down to a single 8-bit int8_t byte - just like the 1-byte fingerprints we used in our Cuckoo Filter - that same register now processes 32 items simultaneously.

In xsimd, a hardware vector register is represented by xsimd::batch<T>, where T represents the type of value we want to store within each lane. We can query the resulting lane count at any time using xsimd::batch<T>::size:

1#include <iostream>
2#include <cstdint>
3#include <xsimd/xsimd.hpp>
4
5int main() {
6  std::cout << "Architecture: "
7    << xsimd::default_arch::name() << "\n";
8
9  std::cout << "double lanes:  "
10    << xsimd::batch<double>::size << "\n";
11  std::cout << "float lanes:   "
12    << xsimd::batch<float>::size << "\n";
13  std::cout << "int32_t lanes: "
14    << xsimd::batch<int32_t>::size << "\n";
15  std::cout << "int8_t lanes:  "
16    << xsimd::batch<int8_t>::size;
17}

1Architecture: sse2
2double lanes:  2
3float lanes:   4
4int32_t lanes: 4
5int8_t lanes:  16

The output might be slightly different depending on your hardware. That would be an example of xsimd doing its job. The xsimd::batch type is helping our program automatically adapt to the target CPU.

By default, compilers like GCC, Clang, and MSVC target a lowest-common-denominator CPU from twenty years ago so that your compiled binary can run on ancient hardware. On x86-64, the default baseline is SSE2 (128-bit registers), even if you are compiling on a modern processor with 256-bit AVX2 or 512-bit AVX-512 registers.

To address this, we can pass -march=native on GCC/Clang or /arch:AVX2 on MSVC. This instructs the compiler to query the physical silicon of the host machine and unlock its widest vector instruction sets. The xsimd library then automatically reads these compiler macros and sizes its batches to match your hardware.

If you're using the CMake build system from the first chapter, we already configured this for our dsa_core library when we use Release configurations:

1# ...
2
3if(CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU")
4  target_compile_options(dsa_core PRIVATE
5    $<$<CONFIG:Release>:-O3 -march=native> 
6  )
7elseif(MSVC)
8  target_compile_options(dsa_core PRIVATE
9    $<$<CONFIG:Release>:/O2 /arch:AVX2>  
10  )
11endif()
12
13# ...

To implement this for code directly written in our dsa_app, we could apply the same configuration to that target too, within the dsa_app/CMakeLists.txt file.

Alternatively, because dsa_app is linking against dsa_core, we could make our dsa_core configuration PUBLIC, meaning it will be inherited by dsa_app as part of that linking process:

1# ...
2
3if(CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU")
4  target_compile_options(dsa_core PRIVATE 
5  target_compile_options(dsa_core PUBLIC 
6    $<$<CONFIG:Release>:-O3 -march=native>
7  )
8elseif(MSVC)
9  target_compile_options(dsa_core PRIVATE 
10  target_compile_options(dsa_core PUBLIC 
11    $<$<CONFIG:Release>:/O2 /arch:AVX2>
12  )
13endif()
14
15# ...

With either approach, once we reconfigure and rebuild our application, we should now see xsimd adapting to take full advantage of the target hardware. If we're running on a modern CPU, our Architecture string will be updated, and our lane counts will be bumped up:

1Architecture: fma3+avx2
2double lanes:  4
3float lanes:   8
4int32_t lanes: 8
5int8_t lanes:  32

Our examples and explanations will generally assume 256-bit registers, but our implementations will be portable - they will adapt to take full advantage of different architectures.

Let's write a simple program that increments every number in a std::vector using SIMD. Our array will be conveniently sized to fit exactly within a single vector register on the target system, but we'll work on more complex examples soon:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <xsimd/xsimd.hpp>
5
6int main() {
7  using Batch = xsimd::batch<int32_t>;
8  constexpr size_t Lanes{Batch::size};
9
10  // Allocate a vector sized to fit one SIMD register
11  std::vector<int32_t> Data(Lanes);
12  for (size_t i{0}; i < Lanes; ++i) {
13    Data[i] = static_cast<int32_t>(i + 1);
14  }
15  
16  std::cout << "Array:      ";
17  for (int32_t val : Data) {
18    std::cout << val << " ";
19  }
20}

Your output may be different depending on your target hardware. On hardware with 256-bit vector registers, we'll have enough lanes for eight 32-bit integers:

1Array:      1 2 3 4 5 6 7 8

Aligned vs. Unaligned Loads

When loading memory into a vector register, we have to be mindful of memory alignment . A 256-bit (32-byte) register prefers to load data from a memory address that is a clean multiple of 32.

However, standard std::vector allocations only guarantee alignment up to the size of the underlying element (e.g., 4-byte alignment for float, or 1-byte alignment for int8_t). If we load a pointer that is not strictly aligned to the vector register's byte width, the CPU will trigger a segmentation fault and crash our program.

Unless we explicitly allocate memory using a custom aligned allocator specifically designed for the vector registers of our target, we should use load_unaligned() and store_unaligned(). On modern CPUs, unaligned vector loads have virtually zero penalty as long as the data doesn't cross a 64-byte cache-line boundary.

The following program loads our array data into a SIMD batch, and then writes the batch back out to the array:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <xsimd/xsimd.hpp>
5
6int main() {
Previous Code Unchanged|Show
19
20  // Pull contiguous RAM into a vector register
21  Batch Register{Batch::load_unaligned(Data.data())};
22
23  // Log the register contents using the batch's
24  // get() method
25  std::cout << "\nRegister A: ";
26  for (size_t i{0}; i < Lanes; ++i) {
27    std::cout << Register.get(i) << " ";
28  }
29
30  // Write the register contents back out to RAM
31  Register.store_unaligned(Data.data());
32
33  std::cout << "\nArray:      ";
34  for (int32_t val : Data) {
35    std::cout << val << " ";
36  }
37}

1Array:      1 2 3 4 5 6 7 8 
2Register A: 1 2 3 4 5 6 7 8
3Array:      1 2 3 4 5 6 7 8

Broadcasting

Now that our data is sitting inside a vector register, we want to increment (or "offset") every element by 1 using the + operator. However, at the silicon level, vector ALU instructions operate strictly between two vector registers of matching width. The processor cannot directly offset an 8-lane vector register using a single scalar integer sitting in a general-purpose register.

To solve this, we must allocate a second vector register and copy our offset (1) into every single lane: [1, 1, 1, 1, 1, 1, 1, 1]. Duplicating a single scalar value across all lanes of a wide register is often called broadcasting.

The xsimd library gives a static helper method to make this easier. The following expression creates a Batch called Offset with 1 in every lane:

1Batch Offset{Batch::broadcast(1)}

Alternatively, we can just pass a single scalar value directly to the Batch constructor, such as Batch Offset(1), which implicitly performs the broadcast.

Either way, once both operands are in vector registers, the + operator multiplies all lanes simultaneously in a single instruction:

1#include <iostream>
2#include <vector>
3#include <cstdint>
4#include <xsimd/xsimd.hpp>
5
6int main() {
Previous Code Unchanged|Show
26
27  // Broadcast 1 into every lane of a second register
28  Batch Offset{Batch::broadcast(1)}; 
29  
30  // Alternatively:
31  // Batch Offset(1);
32
33  std::cout << "\nRegister B: ";
34  for (size_t i{0}; i < Lanes; ++i) {
35    std::cout << Offset.get(i) << " ";
36  }
37
38  // Increment all values in parallel
39  Batch Result{Register + Offset};
40
41  Result.store_unaligned(Data.data());
42
43  std::cout << "\nArray:      ";
44  for (int32_t val : Data) {
45    std::cout << val << " ";
46  }
47}

1Array:      1 2 3 4 5 6 7 8 
2Register A: 1 2 3 4 5 6 7 8
3Register B: 1 1 1 1 1 1 1 1
4Array:      2 3 4 5 6 7 8 9

Reviewing the previous section, to perform vectorized math, every SIMD algorithm follows a three-step cycle:

1. Load / Broadcast: Pull a contiguous chunk of data from RAM into a batch register, or duplicate a single scalar number across all lanes of a batch.
2. Execute: Run parallel ALU instructions across the registers (addition, multiplication, bitwise logic, comparisons).
3. Store / Extract: Write the resulting batch back out to a contiguous array in RAM, or extract a summary result back into a scalar register.

In the previous example, the amount of numbers we needed to offset could perfectly fit in a single vector register, but we're rarely so lucky. In the real world, we need to iterate over a collection containing a dynamic number of elements.

When designing our algorithms, we should always be mindful of the possibility that our collection size is smaller than our vector register. If we have a collection of only 3 elements, and we try to load it into an 8-lane vector register, we would try to read 5 elements past the end of our array, causing an out-of-bounds buffer overflow.

More commonly, our collections will be much larger than what can be stored in a single register. For example, let's imagine we needed to process 27 elements. Each element is 32 bits, and we're building for an architecture with 256 bit vector registers, meaning we can process 8 of these items at once using SIMD. Our algorithm then becomes:

• Use an initial loop that uses SIMD to process 8 items on each iteration. After 3 iterations, this loop will have processed the first 24 items, with 3 items remaining from our original 27
• Follow up with a second loop that processes a single item on each iteration using standard scalar logic. The remaining three items can be processed in 3 iterations of this loop

The following function builds on the previous example, but it implements this two-loop setup to support an input array of any number of elements, represented by the count argument:

1#pragma once
2#include <xsimd/xsimd.hpp>
3#include <cstddef>
4#include <cstdint>
5
6namespace dsa {
7
8inline void Offset(
9  const int32_t* input,
10  int32_t* output,
11  size_t count,
12  int32_t offset
13) {
14  using Batch = xsimd::batch<int32_t>;
15  constexpr size_t Lanes = Batch::size;
16
17  // Broadcast the scalar constant into every lane
18  // of a dedicated vector register
19  Batch offsetBatch(offset); // [offset, offset, offset, ...]
20
21  size_t i{0};
22
23  // Process the array in wide SIMD chunks
24  // (e.g., 8 elements at a time)
25  for (; i + Lanes <= count; i += Lanes) {
26    // Load 8 contiguous elements from RAM into a register
27    Batch data = Batch::load_unaligned(input + i);
28
29    // Offset all 8 lanes in parallel
30    Batch result = (data + offsetBatch);
31
32    // Write all 8 results back to RAM
33    result.store_unaligned(output + i);
34  }
35
36  // Scalar Tail Loop: Clean up any remaining elements if
37  // 'count' is not cleanly divisible by our lane width
38  for (; i < count; ++i) {
39    output[i] = (input[i] + offset);
40  }
41}
42
43}

If Lanes is 8 and our array has 27 elements, the first loop will execute 3 times, processing indices 0..7, 8..15, and 16..23 in three massive gulps. That leaves 3 stragglers at the end of the array (24, 25, and 26) which the scalar tail loop safely cleans up one by one.

Using our new SIMD Offset() function would look something like this:

1#include <iostream>
2#include <vector>
3#include <dsa/SimdAlgorithms.h>
4
5int main() {
6  // 10 elements: processed by the SIMD loop first,
7  // with remaining stragglers handled by the tail loop
8  std::vector<int32_t> Input{
9    1, 2, 3, 4, 5, 6, 7, 8, 9, 10
10  };
11
12  std::vector<int32_t> Output(Input.size());
13
14  dsa::Offset(
15    Input.data(), Output.data(), Input.size(), 1
16  );
17
18  std::cout << "Input:  ";
19  for (int32_t val : Input) {
20    std::cout << val << " ";
21  }
22
23  std::cout << "\nOutput: ";
24  for (int32_t val : Output) {
25    std::cout << val << " ";
26  }
27}

1Input:  1 2 3 4 5 6 7 8 9 10 
2Output: 2 3 4 5 6 7 8 9 10 11

We've now abstracted the SIMD instructions completely - it's all hidden away in our library. A consumer can use this Offset() function without any knowledge of how it works, or even what a vector register is. It's just working behind the scenes, magically making their program run more efficiently.

In this lesson, we moved away from fragile compiler heuristics and took explicit control of the CPU's vector registers:

1. Portable Intrinsics: Writing raw vendor intrinsics locks our code to specific hardware architectures. Using xsimd gives us a clean, zero-overhead C++ API that compiles down to optimal vector instructions for our target architecture.
2. Registers and Lanes: A hardware vector register has a fixed physical bit width that subdivides into parallel lanes. Shrinking our data types from 32-bit integers down to 8-bit bytes quadruples the number of elements we can process per cycle.
3. Unaligned Loads and Stores: Standard heap allocations do not guarantee 32-byte or 64-byte alignment, so we use load_unaligned() and store_unaligned() to safely move data between RAM and vector registers without segmentation faults.
4. Broadcasting and Tail Loops: Scalar constants must be broadcast across every lane of a register before executing vectorized math, and dynamic array sizes require a scalar tail loop to safely clean up trailing elements.

Vectorizing arithmetic is only the beginning. In the next lesson, we will learn how to perform searching and branching inside vector registers, combining SIMD comparison masks with bitwise hardware intrinsics to scan 32 bytes in a single clock cycle.