Reference Intel algorithm

In this section we show the reference implementation of Intel solution for histogram computation. The implementation presented in the article is the Go. The assembly translation is available in sources.

func HistogramIntelReference(indices []uint32, output []uint32) {
    mem := output

    zeros := vec_set1(0)
    neg1 := vec_set1(0xffffffff)
    for len(indices) > 0 {
        chunk := vec_load(indices)
        conflicts := vec_conflicts(chunk)

        increment := vec_set1(1)

        anyConflict := vec_not_equals(conflicts, zeros)
        if anyConflict != 0 {
            position := vec_lzcnt(conflicts)
            position = vec_sub(vec_set1(31), position)

            active := anyConflict
            for {
                conflicted := vec_permute_masked(increment, position, active)
                conflicted = vec_mask(conflicted, active)

                increment = vec_add_masked(increment, conflicted, active)

                position = vec_permute_masked(position, position, active)

                active = vec_not_equals(position, neg1)
                if active == 0 {
                    break
                }
            }
        }

        oldValues := vec_gather(mem, chunk)
        newValues := vec_add(oldValues, increment)
        vec_scatter(newValues, mem, chunk)

        indices = indices[16:]
    }
}

The main loop of the algorithm works as follows.

1. First we load 16 indices. For simplicity we assume that we already have these bytes expanded into 32-bit values.

   In our example we have a set of only three values: {1, 22, 33}.

   chunk := vec_load(indices)
   

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │   22│   22│    1│   33│    1│   33│   33│    1│    1│   22│   22│   22│   22│    1│    1│    1│ = chunk
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   dword 15                                                                                  dword 0
   

2. Next, we perform the conflict detection on the loaded indices.

   conflicts := vec_conflicts(chunk)
   

   The contents of conflicts is presented in binary.

   conflicts[0]  = 00000000_00000000   // 1
   conflicts[1]  = 00000000_00000001   // 1
   conflicts[2]  = 00000000_00000011   // 1
   conflicts[3]  = 00000000_00000000   // 22
   conflicts[4]  = 00000000_00001000   // 22
   conflicts[5]  = 00000000_00011000   // 22
   conflicts[6]  = 00000000_00111000   // 22
   conflicts[7]  = 00000000_00000111   // 1
   conflicts[8]  = 00000000_10000111   // 1
   conflicts[9]  = 00000000_00000000   // 33
   conflicts[10] = 00000010_00000000   // 33
   conflicts[11] = 00000001_10000111   // 1
   conflicts[12] = 00000110_00000000   // 33
   conflicts[13] = 00001001_10000111   // 1
   conflicts[14] = 00000000_01111000   // 22
   conflicts[15] = 01000000_01111000   // 22
   

3. We initialize the increment vector to all decimal ones. If there are no conflicts, that is a ready to use increment value.

   increment := vec_set1(1)
   

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│ = icrement
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   

3. When all words of conflict vector are zeros, it means there are no conflicts. We may skip a conflict resolution loop, and increment counters (step 5). But in our example there are repeated indices, thus we are entering the conflict resolution loop.

   anyConflict := vec_not_equals(conflicts, zeros)
   

   1. The prerequisite to handle the repeated values is computing a permutation vector. The permutation vector contains the indices to the previous occurrence of the given element. The index is -1 when the i-th element is the first occurrence of a value.

      This calculation is possible with LZCNT (Leading Zeros Count) instruction. The previous index equals to 31 - LZCNT(conflicts[i]); if there is no conflict (conflicts[i] = 0), then the result is -1.

      position := vec_lzcnt(conflicts)
      position = vec_sub(vec_set1(31), position)
      

      ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
      │   14│    6│   11│   10│    8│    9│   -1│    7│    2│    5│    4│    3│   -1│    1│    0│   -1│ = position
      └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
      

   2. The next step is to set active elements. We are processing only the repeated items — those having the previous element (position[i] != -1). Initially, the active mask equals to the conflicts mask.

      active := anyConflict
      

      ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
      │   14│    6│   11│   10│    8│    9│   -1│    7│    2│    5│    4│    3│   -1│    1│    0│   -1│ = position
      └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
           1     1     1     1     1     1     0     1     1     1     1     1     0     1     1     0  = active
      
                                                                                   1111_1101_1111_0110  = active
      

   3. We repeat this step until active mask is not zero.

      We use the position vector to pull previous values of indices to the vector conflicted, zeroing non-active elements (lines 1 and 2). Then we update only the active elements of increment (line 3) — in our case it is a single, masked addition. Finally, we pull the previous positions of the position vector itself (line 4), and update active mask accordingly (lines 5 to 8).

      /* 1 */ conflicted := vec_permute_masked(increment, position, active)
      /* 2 */ conflicted = vec_mask(conflicted, active)
      
      /* 3 */ increment = vec_add_masked(increment, conflicted, active)
      
      /* 4 */ position = vec_permute_masked(position, position, active)
      
      /* 5 */ active = vec_not_equals(position, neg1)
      /* 6 */ if active == 0 {
      /* 7 */     break
      /* 8 */ }
      

      Let's see in details how these vectors change across iterations. Because the increment vector is full of ones, we named elements of this vector as i{index}. In the first iteration we have:

         i15   i14   i13   i12   i11   i10    i9    i8    i7    i6    i5    i4    i3    i2    i1    i0
      ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
      │    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│    1│ = increment
      └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
      ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
      │   14│    6│   11│   10│    8│    9│   -1│    7│    2│    5│    4│    3│   -1│    1│    0│   -1│ = position
      └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
      ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
      │    1│    1│    1│    1│    1│    1│    0│    1│    1│    1│    1│    1│    0│    1│    1│    0│ = conflicted
      └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
         i14    i6   i11   i10    i8    i9          i7    i2    i5    i4    i3          i1    i0
      

   After updating the increment vector (line 3) its content is:

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    2│    2│    2│    2│    2│    2│    1│    2│    2│    2│    2│    2│    1│    2│    2│    1│ = increment
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
      i15   i14   i13   i12   i11   i10    i9    i8    i7    i6    i5    i4    i3    i2    i1    i0
     +i14   +i6  +i11  +i10   +i8   +i9         +i7   +i2   +i5   +i4   +i3         +i1   +i0
   

   The permutation of position yields the following vector. Note that now we have twice as much -1 values.

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    6│    5│    8│    9│    7│   -1│   -1│    2│    1│    4│    3│   -1│   -1│    0│   -1│   -1│ = position
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
        1     1     1     1     1     0     0     1     1     1     1     0     0     1     0     0  = active
                                                                                 1111_1001_1110_0100 = active
   

   In the second iteration we have:

      i15   i14   i13   i12   i11   i10    i9    i8    i7    i6    i5    i4    i3    i2    i1    i0
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    2│    2│    2│    2│    2│    2│    1│    2│    2│    2│    2│    2│    1│    2│    2│    1│ = increment
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    6│    5│    8│    9│    7│   -1│   -1│    2│    1│    4│    3│   -1│   -1│    0│   -1│   -1│ = position
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    2│    2│    2│    1│    2│    0│    0│    2│    2│    2│    1│    0│    0│    1│    0│    0│ = conflicted
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
       i6    i5    i8    i9    i7                i2    i1    i4    i3                i0
   

   After updating the increment vector (line 3) its content is:

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    4│    4│    4│    3│    4│    2│    1│    4│    4│    4│    3│    2│    1│    3│    2│    1│ = increment
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
      i15   i14   i13   i12   i11   i10    i9    i8    i7    i6    i5    i4    i3    i2    i1    i0
      +i6   +i5   +i8   +i9   +i7               +i2   +i1   +i4   +i3               +i0
   

   The permutation of position yields the following vector.

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    4│    3│    2│   -1│    1│   -1│   -1│    0│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│ = position
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
        1     1     1     0     1     0     0     1     0     0     0     0     0     0     0     0  = active
                                                                                 1110_1001_0000_0000 = active
   

   In the third iteration we have:

      i15   i14   i13   i12   i11   i10    i9    i8    i7    i6    i5    i4    i3    i2    i1    i0
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    4│    4│    4│    3│    4│    2│    1│    4│    4│    4│    3│    2│    1│    3│    2│    1│ = increment
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    4│    3│    2│   -1│    1│   -1│   -1│    0│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│ = position
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    2│    1│    3│    0│    2│    0│    0│    1│    0│    0│    0│    0│    0│    0│    0│    0│ = conflicted
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
       i4    i3    i2          i1                i0
   

   After updating the increment vector (line 3) its content is:

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    6│    5│    7│    3│    6│    2│    1│    5│    4│    4│    3│    2│    1│    3│    2│    1│ = increment
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
      i15   i14   i13   i12   i11   i10    i9    i8    i7    i6    i5    i4    i3    i2    i1    i0
      +i4   +i3   +i2         +i1               +i0
   

   The permutation of position yields the vector full of -1. This means active mask is zero, and we finished the conflict resolution.

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│   -1│ = position
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   

   Before moving to the final step, let's look closer on increment collated with the input vector chunk.

   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │   22│   22│    1│   33│    1│   33│   33│    1│    1│   22│   22│   22│   22│    1│    1│    1│ = chunk
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
   ┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
   │    6│    5│    7│    3│    6│    2│    1│    5│    4│    4│    3│    2│    1│    3│    2│    1│ = increment
   └─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
        │           │     │
        │           │     └─ the last occurrence of 33: 3 is the count
        │           └─────── the last occurrence of  1: 7 is the count
        └─────────────────── the last occurrence of 22: 6 is the count
   

5. The last step is histogram update. First, we gather the counters for the bytes we loaded in the step 1 (line 1). Then, we update all counters with the increment vector calculated earlier (line 2). Finally, we scatter the updated values back (line 3).

   The real counts of repeated values are stored at position corresponding the last occurrence of given value. The vector newValues contains garbage for earlier occurrences of a repeated value. It is not a problem, as the scatter instruction guarantees that the value written to the memory will be the correct, last one.

   /* 1 */ oldValues := vec_gather(mem, chunk)
   /* 2 */ newValues := vec_add(oldValues, increment)
   /* 3 */ vec_scatter(newValues, mem, chunk)
   

   Note: due to high latency of the gather instructions, in the assembly code the gather instruction is issued before the conflict resolution loop. Here we put it closer just to make the algorithm more readable.

Downsides of Intel approach

• The core algorithm has to be repeated twice -- in the conflict resolution loop and in the main loop. The histogram algorithm requires just a simple addition, thus that is not that visible. But when we need more operations, it becomes more significant.
• Likewise, when the algorithm needs more than one variable, then all variables have to be permuted in the same way we permute the increment vector.
• Unpredictable number of iterations.
• Underuse of available resources: in each iteration we process fewer and fewer items. The following table shows how many active elements were processed in the conflict resolution loop for various inputs.

Alternative algorithm

The solution we propose completely rids of the conflict resolution loop. The key observation is that in a single iteration of the main loop we have up to 16 conflicts. Each input item might have assigned the unique repetition id (from 0 to 15) within the subset of the same values.

By extending the algorithm state by the factor of 16, we can update without conflicts all values at once. While the original algorithm updates counters at indices chunk[i], we update a sub-counters at chunk[i] * 16 + subbucket, where subbucket is a number from 0 to 15.

After the main loop completes, we finalize computations by merging these sixteen sub-states into one.

Let's look how this approach would assign subbucket numbers in our example input vector.

┌─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┬─────┐
│   22│   22│    1│   33│    1│   33│   33│    1│    1│   22│   22│   22│   22│    1│    1│    1│ = chunk
└─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┴─────┘
     │     │     |     │     │     │     │     │     │     │     │     │     │     │     │     │
     │     │     6     │     5     │     │     4     3     │     │     │     │     2     1     0
     5     4           │           │     │                 3     2     1     0
                       2           1     0

The following Go code shows how the histogram problem can be solved with this new approach.

func HistogramV2Reference(indices []uint32, output []uint32) {
    mem := output

    for len(indices) > 0 {
        chunk := vec_load(indices)
        conflicts := vec_conflicts(chunk)

        subbuckets := vec_lzcnt(conflicts)
        subbuckets = vec_sub(vec_set1(32), subbuckets)

        offsets := vec_sll(chunk, 4)
        offsets = vec_add(offsets, subbuckets) // 16 * index + subbucket

        oldValues := vec_gather(mem, offsets)
        newValues := vec_add(oldValues, vec_set1(1))
        vec_scatter(newValues, mem, offsets)

        indices = indices[16:]
    }

    // finalize
    for i := 0; i < bins; i++ {
        pos := i * 16
        sum := uint32(0)
        for j := 0; j < 16; j++ {
            sum += output[pos+j]
        }

        output[i] = sum
    }
}