Fast merge of sorted subsets of 4K floating-point numbers in L1/L2

algorithm, assembly, c, c++, performance

Solution

Here's a very naive way to do it. (Please excuse any 4am delirium-induced pseudo-code bugs ;)

//4x sorted subsets
data[4][4] = {
  {3, 4, 5, INF},
  {2, 7, 8, INF},
  {1, 4, 4, INF},
  {5, 8, 9, INF}
}

data_offset[4] = {0, 0, 0, 0}

n = 4*3

for(i=0, i<n, i++):
  sub = 0
  sub = 1 * (data[sub][data_offset[sub]] > data[1][data_offset[1]])
  sub = 2 * (data[sub][data_offset[sub]] > data[2][data_offset[2]])
  sub = 3 * (data[sub][data_offset[sub]] > data[3][data_offset[3]])

  out[i] = data[sub][data_offset[sub]]
  data_offset[sub]++

Edit: With AVX2 and its gather support, we could compare up to 8 subsets at once.

Edit 2: Depending on type casting, it might be possible to shave off 3 extra clock cycles per iteration on a Nehalem (mul: 5, shift+sub: 4)

//Assuming 'sub' is uint32_t
sub = ... << ((data[sub][data_offset[sub]] > data[...][data_offset[...]]) - 1)

Edit 3: It may be possible to exploit out-of-order execution to some degree, especially as K gets larger, by using two or more `max` values:

max1 = 0
max2 = 1
max1 = 2 * (data[max1][data_offset[max1]] > data[2][data_offset[2]])
max2 = 3 * (data[max2][data_offset[max2]] > data[3][data_offset[3]])
...
max1 = 6 * (data[max1][data_offset[max1]] > data[6][data_offset[6]])
max2 = 7 * (data[max2][data_offset[max2]] > data[7][data_offset[7]])

q = data[max1][data_offset[max1]] < data[max2][data_offset[max2]]

sub = max1*q + ((~max2)&1)*q

Edit 4:

Depending on compiler intelligence, we can remove multiplications altogether using the ternary operator:

sub = (data[sub][data_offset[sub]] > data[x][data_offset[x]]) ? x : sub

Edit 5:

In order to avoid costly floating point comparisons, we could simply `reinterpret_cast<uint32_t*>()` the data, as this would result in an integer compare.

Another possibility is to utilize SSE registers as these are not typed, and explicitly use integer comparison instructions.

This works due to the operators `< > ==` yielding the same results when interpreting a float on the binary level.

Edit 6:

If we unroll our loop sufficiently to match the number of values to the number of SSE registers, we could stage the data that is being compared.

At the end of an iteration we would then re-transfer the register which contained the selected maximum/minimum value, and shift it.

Although this requires reworking the indexing slightly, it may prove more efficient than littering the loop with `LEA`'s.

Problem

What is a fast way to merge sorted subsets of an array of up to 4096 32-bit floating point numbers on a modern (SSE2+) x86 processor? Please assume the following: - The size of the entire set is at maximum 4096 items - The size of the subsets is open to discussion, but let us assume between 16-256 initially - All data used through the merge should preferably fit into L1 - The L1 data cache size is 32K. 16K has already been used for the data itself, so you have 16K to play with - All data is already in L1 (with as high degree of confidence as possible) - it has just been operated on by a sort - All data is 16-byte aligned - We want to try to minimize branching (for obvious reasons) Main criteria of feasibility: faster than an in-L1 LSD radix sort. I'd be very interested to see if someone knows of a reasonable way to do this given the above parameters! :)

Original source

Related problems