# Sorting fixed-length arrays using SIMD

Simply sorting an array is a very simple task, while efficient sorting is very difficult, largely because of the simplicity of the task.

The meaning of this logical paradox is that a complex problem can be solved in many ways, among which you can always try to find an even more effective one. While solving a simple problem is usually possible in only a few ways, it is not so easy to improve, since over decades of work by many people, optimization has achieved a truly impressive result.

The most obvious way to optimize the sorting of an array is to sort its different parts in parallel, but this method faces the problem of high overhead costs for running parallel processes and the disproportionate complexity of solving the problem and the result obtained. As a result, it seems very interesting to use the capabilities of SIMD instructions for sorting, which, even for the “standard” SSE2, can theoretically allow sorting the 4th stream of 32-bit numbers and theoretically achieve a fourfold increase in speed.

Below is a variant of sorting an array of 16, 32-bit unsigned numbers, implemented using the SIMD instructions of the SSE4.1 set. The presented code is essentially a demonstrator of the algorithm and does not imply practical application.

The objections that arrays are rarely of a clearly fixed length, I would like to counter with the remark that large arrays can be divided into sub-arrays of the required length, the algorithm can be easily extended to 32 and 64 members, and arrays of non-multiple lengths can be easily padded” redundant” elements having a maximum/minimum value, which, as a result, during sorting will focus on any of the ends of the array where they can be easily ignored.

The objection that this is a special case only for unsigned 32-bit numbers, I would like to counter with the remark that SIMD instructions allow you to write absolutely similar algorithms for signed integers and single-precision floating-point numbers. The efficiency of working with double precision numbers is no longer so significant, but continues to be preserved.

We will assume that the method takes one single parameter, this is the beginning of the array, which, in accordance with the VECTORCALL calling convention, will be located in the ECX register.

We load all 16 values to be sorted into SIMD registers. This is the first and last access to memory for reading, the next access will be used to write the sorted array. Thus, having received the data once, the algorithm will return an already fully sorted array, minimizing the number of memory accesses.

```
movdqa xmm1, [rcx + xmmword * 0]
movdqa xmm2, [rcx + xmmword * 1]
movdqa xmm4, [rcx + xmmword * 2]
movdqa xmm5, [rcx + xmmword * 3]
```

We sort the array with literally 6 commands, combining the number series into eight ordered pairs.

```
movdqa xmm0, xmm2
pmaxud xmm2, xmm1
pminud xmm1, xmm0
movdqa xmm3, xmm5
pmaxud xmm5, xmm4
pminud xmm4, xmm3
```

The next step is sorting the resulting pairs into tetrads, and if finding the smallest / largest is not difficult, you just need to compare the smallest / largest number from each pair with the same number from the other pair, then finding “intermediate” numbers is not so obvious, but also quite simply. It is necessary to find the smallest of the largest and the largest of the smallest, and then compare and order them. This requires only 9 instructions.

```
movdqa xmm0, xmm1
pminud xmm0, xmm4
movdqa xmm3, xmm2
pmaxud xmm3, xmm5
pminud xmm2, xmm5
pmaxud xmm1, xmm4
movdqa xmm4, xmm1
pminud xmm1, xmm2
pmaxud xmm2, xmm4
```

We transpose the resulting matrix for further sorting. I do not presume to assert that the presented transposition code is the most efficient, perhaps it can be shortened by a couple of cycles, but this is not critical for demonstrating the algorithm.

```
movdqa xmm4, xmm0
punpckhdq xmm4, xmm1
movdqa xmm5, xmm2
punpckldq xmm5, xmm3
punpckldq xmm0, xmm1
movdqa xmm1, xmm0
punpcklqdq xmm0, xmm5
punpckhqdq xmm1, xmm5
punpckhdq xmm2, xmm3
movdqa xmm3, xmm4
punpckhqdq xmm3, xmm2
punpcklqdq xmm4, xmm2
```

Now, after transposition, in each register there are four numbers ordered in ascending order. We will make 4 comparisons and exchanges, during which we will unite 4 tetrads into 2 octaves. It is worth noting that if in the first comparison there is an exchange of 4 pairs of numbers, then in subsequent comparisons the number of exchanges decreases by one each time and in the last comparison there is only one exchange.

```
mov eax, 4
@@:
movdqa xmm2, xmm1
pmaxud xmm1, xmm0
pminud xmm0, xmm2
pshufd xmm1, xmm1, 10010011b
movdqa xmm5, xmm4
pmaxud xmm4, xmm3
pminud xmm3, xmm5
pshufd xmm4, xmm4, 10010011b
dec eax
jnz @b
```

We perform the last series of comparisons, during which we finally arrange the numerical sequence by comparing two octaves.

```
mov eax, 8
@@:
movdqa xmm2, xmm3
movdqa xmm5, xmm4
pmaxud xmm3, xmm0
pmaxud xmm4, xmm1
pminud xmm0, xmm2
pminud xmm1, xmm5
pshufd xmm3, xmm3, 10010011b
pshufd xmm4, xmm4, 10010011b
movss xmm5, xmm3
movss xmm3, xmm4
movss xmm4, xmm5
dec eax
jnz @b
```

Formally, in order for the algorithm to be considered complete, it is necessary to write the received data to memory.

```
movdqa [rcx + xmmword * 0], xmm0
movdqa [rcx + xmmword * 1], xmm1
movdqa [rcx + xmmword * 2], xmm3
movdqa [rcx + xmmword * 3], xmm4
```

I hope you were interested, I will be glad to hear interesting, constructive criticism.

PS The efficiency of the algorithm increases significantly if the target machine has AVX.