How to find largest value between 4 different Vector128 using SIMD calculations

Question

I am trying to do something with SIMD calculations. I have come quite far in my problem where I then get stuck and wonder how this could be done.

I think the easiest way is to describe this step by step what I have done:

I use Vector128<byte> which then handles 16 bytes at a time

1. I have created a 2 dimensional array(array2D) with 9 columns and 16 rows per column. I have put the numbers in a sequence of: 0 and 2. This means that for example Row: 0 has only 0s. Row: 1 has only 2s etc.

2. Now I Avx.LoadVector128 for each column/dimension which gives: 9 Vector128<byte> which I put in: dimensionLIST

3. Now the task is to count how many of the numbers: 0 and 2 that could be found on EACH ROW. (We have 16 rows). This information is in the end stored in: counts[0]

4. Looking at the result of counts[0] in the MessageBox. Below is shown: MessageBox.Show(counts[0]);

(represents 16 rows)
[0,9,0,9,0,9,0,9,0,9,0,9,0,9,0,9]

9, 2s were found on every other row.

Now the goal is to count how many "9" that were found in:
[0,9,0,9,0,9,0,9,0,9,0,9,0,9,0,9] which is 8.

So somehow we want the the integer 8 as Scalar somehow here?

    public unsafe static void SIMDfunction()
    {
        //Create dummy values
        byte[,] array2D = new byte[9, 16]; byte num = 0;
        for (int i = 0; i < 9; i++)
        {
            for (int i2 = 0; i2 < 16; i2++)
            {
                array2D[i, i2] = num;
                if (num == 0) { num = 2; } else { num = 0; }
            }
        }


        /*----------------------------------------------------------------------------------------*/
        unsafe
        {
            //Below starts SIMD calculations!
            fixed (byte* ptr = array2D)
            {
                //Add all 9 dimensions as Vector128
                List<Vector128<byte>> dimensionLIST = new List<Vector128<byte>>();
                for (int i = 0; i < 9; i++)
                {
                    byte* featuredimension = &*((byte*)(ptr + i * 16)); //This gives the first dimension with start: 0
                    dimensionLIST.Add(Avx.LoadVector128(&featuredimension[0])); //add "featuredimension" as a vector of the 16 next numbers: [0,1,2,3,0,1,2,3,0,1,2,3,0,1,2,3]
                }


                //Now count how many of: 0,1,2,3 are found in total in all "dimensionLIST" together?
                Span<Vector128<byte>> counts = stackalloc Vector128<byte>[1];
                Span<Vector128<UInt64>> sum64 = stackalloc Vector128<UInt64>[1];
                byte nr2 = 2; byte nr3 = 9; 
                for (int i = 0; i < dimensionLIST.Count; i++) //Each column
                {
                    //Compare: dimensionLIST[i] with Vector128 val to find out how many matches of 2 in this loop
                    //[0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2], [2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]
                    var match = Avx.CompareEqual(dimensionLIST[i], Vector128.Create(nr2)); //Create Vector128 for numbers: 2
                    counts[0] = Avx.Subtract(counts[0], match);
                }
                //STEP1: Show result on how many 2s are found == 9 occurences of "2"!
                var result = Avx.CompareEqual(Vector128.Create(nr3), counts[0]); //counts[0]: [0,9,0,9,0,9,0,9,0,9,0,9,0,9,0,9] (In total 9 2s are found on those indexes)


                //result:[0,255,0,255,0,255,0,255,0,255,0,255,0,255,0,255] Puts - 1 where integer == 9
                MessageBox.Show(result.ToString());

                //Now the goal is to count how many "9" that were found in: [0,9,0,9,0,9,0,9,0,9,0,9,0,9,0,9] which is 8.
                //So somehow we want the the integer 8 as Scalar somehow here?
            }
        }
    }

Answer 1

Answer to edited question: counting matches horizontally for one vector

This seems oversimplified, like in practice you won't know that 9 is the value you're looking for. But since you hard-coded that in your source, maybe that's what you want.

You're on the right track with pcmpeqb to find exact matches for that element you're looking for. Then you just need to horizontally sum those matches:

// This is C, sorry I don't know C#
#include <immintrin.h>

// input: a vector like {0,9,0,9,...} and a value like 9 or 0
int count_matches(__m128i v, int val)
{
    __m128i matches = _mm_cmpeq_epi8(v, _mm_set1_epi8(val));  // 0 or -1
    matches = _mm_sub_epi8(_mm_setzero_si128(), matches);     // 0 or +1
    __m128i hsums = _mm_sad_epu8(matches, _mm_setzero_si128());  // psadbw against 0 = hsum of each qword half = count ones
    __m128i high = _mm_unpackhi_epi64(hsums, hsums);      // punpckhqdq to extract high half
    hsums = _mm_add_epi32(hsums, high);                // paddd or paddw would be fine
    unsigned sum = _mm_cvtsi128_si32(hsums);           // movd to extract the low 32-bit scalar

    return sum;
}

(Godbolt - fun fact: clang10 "optimizes" the sub to pand with set1_epi8(1) from memory, even with -march=skylake where vpsubb has the same performance as vpand. Silly compiler.)

i.e. just horizontal sum the number of "true" elements from the pcmpeqb result.

Negating with psubb (or doing pand with set1(1)) before psadbw is more efficient than anything I could come up with for turning sum*255 into sum, if we had just horizontally added the 0 or 255 elements.

-(int8_t)sum or (int8_t)-sum compiles to movsx eax, al / neg eax which is 2 instructions (assuming we need the result as a 32-bit int), worse than vpsubb against a zeroed vector that already exists. Without AVX it might be better though, or if you're bottlenecked on back-end SIMD execution ports instead of the front-end.

sum/255 would obviously be bad, compilers don't have enough info to optimize that, that's why my answer doesn't use it. Another option is (sum + 16) >> 8, which happens to give the right answer for all i*255 values from 0*255 through 16*255. But shifts on Intel CPUs run on port 0 or 6, not any ALU port, so that's probably worse than neg/movsx. neg and movsx can run on any ALU port, so are most flexible at avoiding / not contending with back-end ALU pressure from surrounding code.

vpsubb runs on any of p0, p1, p5 on Intel Skylake and later, but less flexible on earlier CPUs. Without AVX, it might need a movdqa register copy, or redoing an xor-zeroing to make a new zero vector for psadbw.

---

Answer to original / title question, finding which of 4 vectors the max come from in each vertical SIMD element

After you have counts[0..3], if you have spare bits at the top of your counts, left shift by 2 and OR with a 0..3 tag number (to record where it came from) so you can use pmaxub to select the largest count, bringing the tag along with it.

The SIMD MAX operation will act on the whole (counts[i] << 2) | i, so the count part is the most-significant part of the integer value, but the tag part comes with it as part of the integer MAX operation.

The "tag" will act as a tie-breaker, biasing towards the higher i (i.e. toward 3 in your case). If you need equal counts to be treated as vector 0, then XOR with 3 or count down from 3 or something so the tag bits have integer values in the tie-break order you want.

// kind of pseudo-code, I don't know C# so this is more like C with intrinsics
for (int i=0 ; i<4 ; i++){
    counts[i] <<= 2;   // emulate non-existent psllb somehow; 2x paddb or psllw / pand
    counts[i] |= set1(i);  // low 2 bits = tag
}
__m128i m0 = _mm_max_epu8(counts[0], counts[1]);
__m128i m1 = _mm_max_epu8(counts[2], counts[3]);
__m128i max = _mm_max_epu8(m0, m1);
max = _mm_and_si128(max, _mm_set1_epi8(3));  // discard high 6 bits, keep low 2 = tag

---

Or if you don't have room to left shift without losing significant bits, unpack with set1(i) and use _mm_max_epu16 (pmaxuw), with a tree of 2x max -> 1x max separately for the high/low halves. So each integer is (count<<8) | i.

Then you have to re-pack down to just the low byte (tag), probably requiring you to mask away the value byte before _mm_packs_epi16 (packsswb). There's no true inverse to punpcklbw / punpckhbw; pack instructions do signed saturation not truncation.

However, the final mask + pack step is just 2x PAND with set1_epi16(0x00FF) on the inputs to feed one packsswb, not too complicated.

---

You can speed up calculating counts in the first place with Micro Optimization of a 4-bucket histogram of a large array or list - infer counts[3] from set1(total/16) - counts[0..2] (3 SIMD subtractions at the end of the loop, saving a compare/sub every iteration).

Answer 2

I have tried to create a code solution. However the code gives the correct answer: 8 in the end.

Somehow I want to count how many of: 255 are found in:
result:[0,255,0,255,0,255,0,255,0,255,0,255,0,255,0,255]

I then SumAbsoluteDifferences and UnpackHigh.

The thing I wonder about here that I don't know if it is a correct way or OK to do it like this where I divide by: 255? I only did it because it gives the correct answer with a mathematical arithmetic:

int NUM = Convert.ToInt32(Sse2.Add(sum64[0], upper).ToScalar()) / 255;

    public unsafe static void SIMDfunction2()
    {
        //Create dummy values
        byte[,] array2D = new byte[9, 16]; byte num = 0;
        for (int i = 0; i < 9; i++)
        {
            for (int i2 = 0; i2 < 16; i2++)
            {
                array2D[i, i2] = num;
                if (num == 0) { num = 2; } else { num = 0; }
            }
        }


        /*------------------------------------------------------------*/
        unsafe
        {
            //Below starts SIMD calculations!
            fixed (byte* ptr = array2D)
            {
                //Add all 9 dimensions as Vector128
                List<Vector128<byte>> dimensionLIST = new List<Vector128<byte>>();
                for (int i = 0; i < 9; i++)
                {
                    byte* featuredimension = &*((byte*)(ptr + i * 16)); //This gives the first dimension with start: 0
                    dimensionLIST.Add(Avx.LoadVector128(&featuredimension[0])); //add "featuredimension" as a vector of the 16 next numbers: 
                }


                //Now count how many of: 0,2 are found in total in all "dimensionLIST" together?
                Span<Vector128<byte>> matches = stackalloc Vector128<byte>[1];
                Span<Vector128<UInt64>> hsums = stackalloc Vector128<UInt64>[1];
                byte nr2 = 2; byte nr3 = 9;
                for (int i = 0; i < dimensionLIST.Count; i++) //Each column
                {
                    //STEP 0
                    //Compare: dimensionLIST[i] with Vector128 val to find out how many matches of 2 in this loop
                    //[0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2], [2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]
                    var match = Avx.CompareEqual(dimensionLIST[i], Vector128.Create(nr2)); //Create Vector128 for numbers: 2
                    matches[0] = Avx.Subtract(matches[0], match);
                }
                //STEP 1: Show result on how many 2s are found == 9 occurences of "2"!
                var result = Avx.CompareEqual(Vector128.Create(nr3), matches[0]); //matches[0]: [0,9,0,9,0,9,0,9,0,9,0,9,0,9,0,9] (In total 9 2s are found on those indexes)
                //result:[0,255,0,255,0,255,0,255,0,255,0,255,0,255,0,255] Puts 255 where integer == 9

                //STEP 2
                result = Avx.Subtract(Vector128<byte>.Zero, result);

                //STEP 3:
                hsums[0] = Sse2.SumAbsoluteDifferences(result, Vector128<byte>.Zero).AsUInt64();

                //STEP 4:
                Vector128<UInt64> upper = Sse2.UnpackHigh(hsums[0], hsums[0]).AsUInt64(); // punpckhqdq to extract high half

                //STEP 5:
                hsums[0] = Sse2.Add(hsums[0], upper); // paddd or paddw would be fine

                //STEP 6:
                UInt32 SUM = Convert.ToUInt32(hsums[0].ToScalar());

                //This returns: 8 which should be the correct SUM
                MessageBox.Show(SUM.ToString());
            }
        }
    }