10e. SIMD: Unit Strides

Introduction

SIMD improves computational performance by simultaneously executing the same operation on multiple data elements. Technically, this is achieved through the use of specialized vector registers, which can hold several values (e.g., 4 floating-point numbers or 8 integers). They make it possible that operations like multiple additions or multiplications can be completed with a single instruction.

For SIMD to fully exploit this vector-based processing, data should be organized and accessed in memory in specific ways. Two fundamental concepts for understanding this are contiguous memory layout and unit stride access.

Contiguous memory layout means that data elements reside in adjacent memory addresses, therefore not exhibiting gaps. Newly allocated arrays satisfy this condition, enabling the load of entire segments directly into vector registers. In contrast, array views reference the original data structure, and thus don't guarantee contiguity. This can potentially result in highly irregular memory access patterns.

In addition, strides refer to the step size between consecutive memory accesses, with unit strides in particular entailing that elements are accessed sequentially. For example, consider a freshly allocated vector x. Then, accessing its elements through eachindex(x) (or equivalently 1:length(x)) ensures a unit stride, as each access moves sequentially to the next element in memory. This contrast with ranges having a non-unit stride such as 1:2:length(x) or indices with no predictable pattern (e.g., [1, 5, 3, 4] as an index vector).

When it comes to performance, SIMD is most effectively applied when data is stored in a contiguous memory block and accessed in a unit stride pattern. This has significant practical implications for the approach to creating slices, where the options are given by copies or views. Essentially, this choice gives rise to a trade-off, as views reduce memory allocations but at the expense of hindering computational efficiency. This section explores such a decision.

Contiguous Blocks in Memory

In the section discussing decreases in memory allocations, we highlighted the benefits of using views over copies when handling slices. Specifically, views maintain references to the original data, thereby avoiding the cost of additional memory allocation. However, views can lead to irregular memory access patterns if data are too scattered. This is why that section also remarked on that copies could outperform views in some scenarios. We're now in a position to explain in more depth why this occurs.

Creating copies of some data structure involves allocating the information in a new contiguous block of memory. This ensures that all elements are stored sequentially, thus offering two key advantages: a quicker fetching of elements and a more effective use of SIMD instructions.

To illustrate this, consider retrieving books from a library. If every book you need resides on a single shelf, collecting them is straightforward—you move once, grab the entire stack, and proceed. This mirrors contiguous memory access. Conversely, if the books are dispersed across different floors and sections, each retrieval demands additional time and effort, akin to non-contiguous access.

In addition to this, our recollection of books would be even more efficient if we had a cart capable of carrying multiple books at once—SIMD operations act like this cart. Note that even without such a cart, the act of gathering contiguous books remains inherently faster, as the physical effort (or computational cycles) is minimized. In other words, the benefits of SIMD are on top of the faster memory access.

Vector and Boolean Indexing

The following examples highlight the two advantages from contiguous memory access. To isolate their performance effect, we create the slices outside the function to be benchmarked. In this way, the functions we'll present don't entail memory allocations.

For the case of vector indexing:

x       = rand(5_000_000)

indices = sortperm(x)
y       = @view x[indices]

function foo(y)
    output = 0.0

    for a in y
        output += a
    end

    return output
end

  21.481 ms (0 allocations: 0 bytes)

x       = rand(5_000_000)

indices = sortperm(x)
y       = @view x[indices]

function foo(y)
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  20.754 ms (0 allocations: 0 bytes)

x       = rand(5_000_000)

indices = sortperm(x)
y       = x[indices]

function foo(y)    
    output = 0.0

    for a in y
        output += a
    end

    return output
end

  2.281 ms (0 allocations: 0 bytes)

x       = rand(5_000_000)

indices = sortperm(x)
y       = x[indices]

function foo(y)    
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  902.847 μs (0 allocations: 0 bytes)

while for Boolean indexing:

x       = rand(5_000_000)

indices = x .> 0.5
y       = @view x[indices]

function foo(y)
    output = 0.0

    for a in y
        output += a
    end

    return output
end

  2.206 ms (0 allocations: 0 bytes)

x       = rand(5_000_000)

indices = x .> 0.5
y       = @view x[indices]

function foo(y)
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  1.878 ms (0 allocations: 0 bytes)

x       = rand(5_000_000)

indices = x .> 0.5
y       = x[indices]

function foo(y)    
    output = 0.0

    for a in y
        output += a
    end

    return output
end

  1.015 ms (0 allocations: 0 bytes)

x       = rand(5_000_000)

indices = x .> 0.5
y       = x[indices]

function foo(y)    
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  246.526 μs (0 allocations: 0 bytes)

Finally, comparing unit strides with non-unit strides:

x = rand(1_000_000)
y = @view x[1:2:length(x)]

function foo(y)
    output = 0.0

    for a in y
        output += a
    end

    return output
end

  902.479 μs (0 allocations: 0 bytes)

x = rand(1_000_000)
y = @view x[1:2:length(x)]

function foo(y)
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  889.059 μs (0 allocations: 0 bytes)

x = rand(1_000_000)
y = x[1:2:length(x)]

function foo(y)
    output = 0.0

    for a in y
        output += a
    end

    return output
end

  196.497 μs (0 allocations: 0 bytes)

x = rand(1_000_000)
y = x[1:2:length(x)]

function foo(y)
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  38.274 μs (0 allocations: 0 bytes)

Some Remarks

Note that views don't necessarily impede a sequential memory access. If the view consists of sequential elements, then we'd obtain the same performance relative to a copy created outside the function.

x       = rand(1_000_000)

indices = 1:length(x)
y       = @view x[indices]

function foo(y)
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  76.950 μs (0 allocations: 0 bytes)

x       = rand(1_000_000)

indices = 1:length(x)
y       = x[indices]

function foo(y)    
    output = 0.0

    @simd for a in y
        output += a
    end

    return output
end

  76.777 μs (0 allocations: 0 bytes)

A corollary of this example is that views should be employed in cases like this, as they additionally avoid memory allocations.

Another remark is that storing elements contiguously is only a necessary condition for contiguous memory access. To demonstrate this, let's consider a scenario where we access vectors through the ranges.

Specifically, consider a vector x. Moreover, there's another vector y that comprises the same elements as x, but with zeros placed in between. Formally:

x_size = 1_000_000

x = rand(x_size)

y = zeros(eltype(x),x_size * 2)
    temp  = view(y, 2:2:length(y))
    temp .= x

3-element Vector{Float64}:
 0.906299638797481
 0.44349373245960455
 0.7456733811393941

6-element Vector{Float64}:
 0.0
 0.906299638797481
 0.0
 0.44349373245960455
 0.0
 0.7456733811393941

Given this, we'd get the same value whether we add all elements from x or add all elements in y skipping zeros. However, the results below reveal that, despite both x and y being contiguous blocks in memory, the pattern access with unit strides only holds for x. This explains the differences in the performance observed.

function foo(x)
    output = 0.0

    @inbounds @simd for i in 1:length(x)
        output += x[i]
    end

    return output
end

  78.349 μs (0 allocations: 0 bytes)

function foo(y)
    output = 0.0

    @inbounds @simd for i in 2:2:length(y)
        output += y[i]
    end

    return output
end

  188.483 μs (0 allocations: 0 bytes)