If you have an input array, and an output array, but you only want to write those elements which pass a certain condition, what would be the most efficient way to do this in AVX2?

I've seen in SSE where it was done like this: (From:https://deplinenoise.files.wordpress.com/2015/03/gdc2015_afredriksson_simd.pdf)

__m128i LeftPack_SSSE3(__m128 mask, __m128 val)
{
 // Move 4 sign bits of mask to 4-bit integer value.
 int mask = _mm_movemask_ps(mask);
 // Select shuffle control data
 __m128i shuf_ctrl = _mm_load_si128(&shufmasks[mask]);
 // Permute to move valid values to front of SIMD register
 __m128i packed = _mm_shuffle_epi8(_mm_castps_si128(val), shuf_ctrl);
 return packed;
}

This seems fine for SSE which is 4 wide, and thus only needs a 16 entry LUT, but for AVX which is 8 wide, the LUT becomes quite large(256 entries, each 32 bytes, or 8k).

I'm surprised that AVX doesn't appear to have an instruction for simplifying this process, such as a masked store with packing.

I think with some bit shuffling to count the # of sign bits set to the left you could generate the necessary permutation table, and then call _mm256_permutevar8x32_ps. But this is also quite a few instructions I think..

Does anyone know of any tricks to do this with AVX2? Or what is the most efficient method?

Here is an illustration of the Left Packing Problem from the above document:

Left.Packing.Problem

Thanks

share|improve this question
1  
You could use VGATHERDPS under the assumption that the src is in memory. Before that you have to create the appropriate indices from the mask. – zx485 Apr 29 '16 at 8:04
1  
It's worse than you think. The AVX2 256-bit VPSHUFB instruction can't move data between the 128-bit vector lanes. You'd need vpermd to do that, which would need a second lookup-table. – EOF Apr 29 '16 at 11:12
1  
@EOF: Thanks for this important addition. That VPSHUFB, (scroll down to 'VEX.256 encoded version') does not operate on a 256-bit vector but instead operates on two separate 128-bit vectors in a YMM is noteworthy. Another major inconsistency in the Intel ISA. – zx485 Apr 29 '16 at 13:03
2  
@zx485: I'll have to disagree with you on the "inconsistency". The separate AVX-lanes are actually fairly consistent, with the few instructions that can cross them being explicitly documented. Also, what other ISA even offers 256-bit vectors at all? Yes, there's a price to pay for compatibility, but AVX2 is a really nice vector instruction set. – EOF Apr 29 '16 at 13:06
2  
@EOF: I'll have to disagree with your preceding elaborations, too, but from my/another point of view. Due to legacy over legacy, the Intel ISA is highly fragmented. IMHO a thorough cleanup would be beneficial. Intel tried that with IA-64, but in a strange way. Some days ago I read a posting of Agner Fog, in which he explains the inconsistencies of the x86/64 architecture proliferated, titled '...a big step forward - but repeating past mistakes!'. – zx485 Apr 29 '16 at 13:23
up vote 14 down vote accepted

AVX2 + BMI2. See my other answer for AVX512. (Update: saved a pdep in 64bit builds.)

We can use AVX2 vpermps (_mm256_permutevar8x32_ps) (or the integer equivalent, vpermd) to do a lane-crossing variable-shuffle.

We can generate masks on the fly, since BMI2 pext (Parallel Bits Extract) provides us with a bitwise version of the operation we need.


For integer vectors with 32-bit or wider elements: Either 1) _mm256_movemask_ps(_mm256_castsi256_ps(compare_mask)).
Or 2) use _mm256_movemask_epi8 and then change the first PDEP constant from 0x0101010101010101 to 0x0F0F0F0F0F0F0F0F to scatter blocks of 4 contiguous bits. Change the multiply by 0xFFU into expanded_mask |= expanded_mask<<4; or expanded_mask *= 0x11; (Not tested). Either way, use the shuffle mask with VPERMD instead of VPERMPS.

For 64-bit integer or double elements, everything still Just Works; The compare-mask just happens to always have pairs of 32-bit elements that are the same, so the resulting shuffle puts both halves of each 64-bit element in the right place. (So you still use VPERMPS or VPERMD, because VPERMPD and VPERMQ are only available with immediate control operands.)


The algorithm:

Start with a constant of packed 3 bit indices, with each position holding its own index. i.e. [ 7 6 5 4 3 2 1 0 ] where each element is 3 bits wide. 0b111'110'101'...'010'001'000.

Use pext to extract the indices we want into a contiguous sequence at the bottom of an integer register. e.g. if we want indices 0 and 2, our control-mask for pext should be 0b000'...'111'000'111. pext will grab the 010 and 000 index groups that line up with the 1 bits in the selector. The selected groups are packed into the low bits of the output, so the output will be 0b000'...'010'000. (i.e. [ ... 2 0 ])

See the commented code for how to generate the 0b111000111 input for pext from the input vector mask.

Now we're in the same boat as the compressed-LUT: unpack up to 8 packed indices.

By the time you put all the pieces together, there are three total pext/pdeps. I worked backwards from what I wanted, so it's probably easiest to understand it in that direction, too. (i.e. start with the shuffle line, and work backward from there.)

We can simplify the unpacking if we work with indices one per byte instead of in packed 3-bit groups. Since we have 8 indices, this is only possible with 64bit code.

See this and a 32bit-only version on the Godbolt Compiler Explorer. I used #ifdefs so it compiles optimally with -m64 or -m32. gcc wastes some instructions, but clang makes really nice code.

#include <stdint.h>
#include <immintrin.h>

// Uses 64bit pdep / pext to save a step in unpacking.
__m256 compress256(__m256 src, unsigned int mask /* from movmskps */)
{
  uint64_t expanded_mask = _pdep_u64(mask, 0x0101010101010101);  // unpack each bit to a byte
  expanded_mask *= 0xFF;    // mask |= mask<<1 | mask<<2 | ... | mask<<7;
  // ABC... -> AAAAAAAABBBBBBBBCCCCCCCC...: replicate each bit to fill its byte

  const uint64_t identity_indices = 0x0706050403020100;    // the identity shuffle for vpermps, packed to one index per byte
  uint64_t wanted_indices = _pext_u64(identity_indices, expanded_mask);

  __m128i bytevec = _mm_cvtsi64_si128(wanted_indices);
  __m256i shufmask = _mm256_cvtepu8_epi32(bytevec);

  return _mm256_permutevar8x32_ps(src, shufmask);
}

This compiles to code with no loads from memory, only immediate constants. (See the godbolt link for this and the 32bit version).

    # clang 3.7.1 -std=gnu++14 -O3 -march=haswell
    mov     eax, edi                   # just to zero extend: goes away when inlining
    movabs  rcx, 72340172838076673     # The constants are hoisted after inlining into a loop
    pdep    rax, rax, rcx              # ABC       -> 0000000A0000000B....
    imul    rax, rax, 255              # 0000000A0000000B.. -> AAAAAAAABBBBBBBB..
    movabs  rcx, 506097522914230528
    pext    rax, rcx, rax
    vmovq   xmm1, rax
    vpmovzxbd       ymm1, xmm1         # 3c latency since this is lane-crossing
    vpermps ymm0, ymm1, ymm0
    ret

So, according to Agner Fog's numbers, this is 6 uops (not counting the constants, or the zero-extending mov that disappears when inlined). On Intel Haswell, it's 16c latency (1 for vmovq, 3 for each pdep/imul/pext / vpmovzx / vpermps). There's no instruction-level parallelism. In a loop where this isn't part of a loop-carried dependency, though, (like the one I included in the Godbolt link), the bottleneck is hopefully just throughput, keeping multiple iterations of this in flight at once.

This can maybe manage a throughput of one per 3 cycles, bottlenecked on port1 for pdep/pext/imul. Of course, with loads/stores and loop overhead (including the compare, movmsk, and popcnt), total uop throughput can easily be an issue. (e.g. the filter loop in my godbolt link is 14 uops with clang, with -fno-unroll-loops to make it easier to read. It might sustain one iteration per 4c, keeping up with the front-end, if we're lucky, but I think clang failed to account for popcnt's false dependency on its output, so it will bottleneck on 3/5ths of the latency of the compress256 function.)

gcc does the multiply by 0xFF with multiple instructions, using a left shift by 8 and a sub. This takes an extra mov instructions, but the end result is a multiply with a latency of 2. (Haswell handles mov at register-rename stage with zero latency.)


Since all hardware that supports AVX2 also supports BMI2, there's probably no point providing a version for AVX2 without BMI2.

If you need to do this in a very long loop, the LUT is probably worth it if the initial cache-misses are amortized over enough iterations with the lower overhead of just unpacking the LUT entry. You still need to movmskps, so you can popcnt the mask and use it as a LUT index, but you save a pdep/imul/pexp.

You can unpack LUT entries with the same integer sequence I used, but @Froglegs's set1() / vpsrlvd / vpand is probably better when the LUT entry starts in memory and doesn't need to go into integer registers in the first place. (A 32bit broadcast-load doesn't need an ALU uop on Intel CPUs). However, a variable-shift is 3 uops on Haswell (but only 1 on Skylake).

share|improve this answer
1  
I tested it on haswell and it works, nice job! The only issue is that for some reason on MSVC _pdep_u64 and _mm_cvtsi64_si128 are only available if compiling for x64. They get defined out in 32bit builds. – Froglegs Apr 30 '16 at 15:03
1  
Congats on getting this right without having the hardware. I am suprised you have not received more than two (from the OP and me) votes. I added an answer using an instruction LUT. What do you think of this solution? Maybe it's a bad idea. – Z boson May 1 '16 at 9:58
1  
@Christoph: Good catch on the loop increment, thanks. I probably copied it from the avx512 version. Also a good point about sign bits. If you're willing to consider -0.0 as negative, you can even do stuff like XOR floats together and then then use that as the control for a variable blend to do something based on whether their signs are equal or not (e.g. while checking if a point is inside a bounding box), instead of multiplying them or doing two compares. But note that NaNs can have the sign-bit set or not, so this isn't safe wrt to NaNs. – Peter Cordes Jan 30 '17 at 16:01
2  
@Christoph : Correction: On Skylake vpand has latency 1 and throughput 1/3. Note that vpsrlvd is very slow on Haswell: latency 2 and throughput 2. Therefore, on Haswell your solution will be faster. – wim Feb 22 '17 at 15:34
2  
@wim: AMD's new Zen I think still has 128b vector execution units (so 256b ops have half throughput). Doing more in scalar integer will be a win there, if pdep is fast on Zen. (It is supported, but I don't think there are latency numbers yet). I think overall throughput should be more important than latency here, since the loop-carried dependency is only on popcnt and its input. Thanks for the vpmovmskb idea; I'll update my answer with that sometime. (Or feel free to add a paragraph and a godbolt link to the answer yourself; I might not get back to this very soon). – Peter Cordes Feb 22 '17 at 21:27

See my other answer for AVX2+BMI2 with no LUT.

Since you mention a concern about scalability to AVX512: don't worry, there's an AVX512F instruction for exactly this:

VCOMPRESSPS — Store Sparse Packed Single-Precision Floating-Point Values into Dense Memory. (There are also versions for double, and 32 or 64bit integer elements (vpcompressq), but not byte or word (16bit)). It's like BMI2 pdep / pext, but for vectors instead of bits in an integer reg.

The destination can be a vector register or a memory operand, while the source is a vector and a mask register. With a register dest, it can merge or zero the upper bits. With a memory dest, "Only the contiguous vector is written to the destination memory location".

To figure out how far to advance your pointer for the next vector, popcnt the mask.

Let's say you want to filter out everything but values >= 0 from an array:

#include <stdint.h>
#include <immintrin.h>
size_t filter_non_negative(float *__restrict__ dst, const float *__restrict__ src, size_t len) {
    const float *endp = src+len;
    float *dst_start = dst;
    do {
        __m512      sv  = _mm512_loadu_ps(src);
        __mmask16 keep = _mm512_cmp_ps_mask(sv, _mm512_setzero_ps(), _CMP_GE_OQ);  // true for src >= 0.0, false for unordered and src < 0.0
        _mm512_mask_compressstoreu_ps(dst, keep, sv);   // clang is missing this intrinsic, which can't be emulated with a separate store

        src += 16;
        dst += _mm_popcnt_u64(keep);   // popcnt_u64 instead of u32 helps gcc avoid a wasted movsx, but is potentially slower on some CPUs
    } while (src < endp);
    return dst - dst_start;
}

This compiles (with gcc4.9 or later) to (Godbolt Compiler Explorer):

 # Output from gcc6.1, with -O3 -march=haswell -mavx512f.  Same with other gcc versions
    lea     rcx, [rsi+rdx*4]             # endp
    mov     rax, rdi
    vpxord  zmm1, zmm1, zmm1             # vpxor  xmm1, xmm1,xmm1 would save a byte, using VEX instead of EVEX
.L2:
    vmovups zmm0, ZMMWORD PTR [rsi]
    add     rsi, 64
    vcmpps  k1, zmm0, zmm1, 29           # AVX512 compares have mask regs as a destination
    kmovw   edx, k1                      # There are some insns to add/or/and mask regs, but not popcnt
    movzx   edx, dx                      # gcc is dumb and doesn't know that kmovw already zero-extends to fill the destination.
    vcompressps     ZMMWORD PTR [rax]{k1}, zmm0
    popcnt  rdx, rdx
    ## movsx   rdx, edx         # with _popcnt_u32, gcc is dumb.  No casting can get gcc to do anything but sign-extend.  You'd expect (unsigned) would mov to zero-extend, but no.
    lea     rax, [rax+rdx*4]             # dst += ...
    cmp     rcx, rsi
    ja      .L2

    sub     rax, rdi
    sar     rax, 2                       # address math -> element count
    ret
share|improve this answer
    
that makes it much easier, avx512 looks very nice! – Froglegs Apr 30 '16 at 1:08
    
You can get rid of the sign-extension, by using _mm_popcnt_u64(). – EOF Apr 30 '16 at 11:11
    
@EOF: Good point. Hopefully AMD's AVX2 CPUs will run 64bit popcnt at the same speed at 32bit popcnt (unlike Steamroller where it has one per 4c throughput instead of one per 2c). Agner's tables don't report popcnt_u64 being slower on other CPUs, just Bulldozer-family. – Peter Cordes Apr 30 '16 at 13:23
    
I dunno why, but msvc also doesn't provide popcnt_u64 unless compiling in x64(disabled in 32 bit builds) – Froglegs Apr 30 '16 at 15:12
    
@Froglegs: of course it doesn't. It's a 64bit integer instruction! Stop wasting your time with 32bit builds, IMO. – Peter Cordes Apr 30 '16 at 15:24

I came up with this method, which uses a compressed LUT, which is 768(+1 padding) bytes, instead of 8k. It requires a broadcast of a single scalar value, which is then shifted by a different amount in each lane, then masked to the lower 3 bits, which provides a 0-7 LUT.

Here is the intrinsics version, along with code to build LUT.

//Generate Move mask via: _mm256_movemask_ps(_mm256_castsi256_ps(mask)); etc
__m256i MoveMaskToIndices(int moveMask) {
    u8 *adr = g_pack_left_table_u8x3 + moveMask * 3;
    __m256i indices = _mm256_set1_epi32(*reinterpret_cast<u32*>(adr));//lower 24 bits has our LUT

    __m256i m = _mm256_sllv_epi32(indices, _mm256_setr_epi32(29, 26, 23, 20, 17, 14, 11, 8));

    //now shift it right to get 3 bits at bottom
    __m256i shufmask = _mm256_srli_epi32(m, 29);
    return shufmask;
}

u32 get_nth_bits(int a) {
    u32 out = 0;
    int c = 0;
    for (int i = 0; i < 8; ++i) {
        auto set = (a >> i) & 1;
        if (set) {
            out |= (i << (c * 3));
            c++;
        }
    }
    return out;
}
u8 g_pack_left_table_u8x3[256 * 3 + 1];

void BuildPackMask() {
    for (int i = 0; i < 256; ++i) {
        *reinterpret_cast<u32*>(&g_pack_left_table_u8x3[i * 3]) = get_nth_bits(i);
    }
}

Here is the assembly generated by VS2015:

lea eax, DWORD PTR [rcx+rcx*2]
movsxd  rcx, eax
lea rax, OFFSET FLAT:?g_pack_left_table_u8x3@@3PAEA ; g_pack_left_table_u8x3
vpbroadcastd ymm0, DWORD PTR [rcx+rax]
vpsllvd ymm0, ymm0, YMMWORD PTR __ymm@000000080000000b0000000e0000001100000014000000170000001a0000001d
vpsrld  ymm0, ymm0, 29
share|improve this answer
    
It looks like you're inventing your own function names for a lot of stuff. e.g. is shift_right_zero_extend_v supposed to be _mm_srlv_epi32? I don't really like Intel's intrinsic names (too much typing, their asm mnemonics are much nicer), but it's clearer. Using auto doesn't even make it clear that the result is a vector, esp. since m & 7 is taking advantage of GNU C vector syntax where a scalar can be implicitly broadcast when used with a vector. – Peter Cordes Apr 30 '16 at 1:21
    
Anyway, nice job with using a smaller LUT, I think, but the code might as well be pseudo-code. (Or maybe it's intended that way?) Also, you don't show the popcnt of the mask, and this functions doesn't return the movmskps result. So I guess the caller will also have to movemask, and then it's up to the compiler to CSE that away? Poor design IMO: have this function take an unsigned int mask. – Peter Cordes Apr 30 '16 at 1:25
    
It is real C++ code that compiles actually-- I don't use the intrinsics manually, I have my own wrappers. This was just a test function, for usage I'd have to take the movemask out of the body of the function, you are correct. The m & 7 bit broadcasts 7 to a register and calls _mm256_and_si256. – Froglegs Apr 30 '16 at 2:32
1  
My point was that writing it the boring / annoying way with Intel's really long function names will make it a better answer, since it makes it clearer exactly what steps are taken. I think your LUT has shuffle masks packed into 3 bytes. And you decompress with pmovzx or something, then vpsrlv, then mask away high garbage in each element? Or are broadcasting one 32b element and then using a variable shift to extract eight 3b elements? I think the latter. Feel free to copy/paste my text description of what you do. – Peter Cordes Apr 30 '16 at 2:38
    
Ya, perhaps I should post it with raw intrinsics then, I'll convert it over and post it again. I can post the table gen code also – Froglegs Apr 30 '16 at 2:39

In case anyone is interested here is a solution for SSE2 which uses an instruction LUT instead of a data LUT aka a jump table. With AVX this would need 256 cases though.

Each time you call LeftPack_SSE2 below it uses essentially three instructions: jmp, shufps, jmp. Five of the sixteen cases don't need to modify the vector.

static inline __m128 LeftPack_SSE2(__m128 val, int mask)  {
  switch(mask) {
  case  0:
  case  1: return val;
  case  2: return _mm_shuffle_ps(val,val,0x01);
  case  3: return val;
  case  4: return _mm_shuffle_ps(val,val,0x02);
  case  5: return _mm_shuffle_ps(val,val,0x08);
  case  6: return _mm_shuffle_ps(val,val,0x09);
  case  7: return val;
  case  8: return _mm_shuffle_ps(val,val,0x03);
  case  9: return _mm_shuffle_ps(val,val,0x0c);
  case 10: return _mm_shuffle_ps(val,val,0x0d);
  case 11: return _mm_shuffle_ps(val,val,0x34);
  case 12: return _mm_shuffle_ps(val,val,0x0e);
  case 13: return _mm_shuffle_ps(val,val,0x38);
  case 14: return _mm_shuffle_ps(val,val,0x39);
  case 15: return val;
  }
}

__m128 foo(__m128 val, __m128 maskv) {
  int mask = _mm_movemask_ps(maskv);
  return LeftPack_SSE2(val, mask);
}
share|improve this answer
1  
If you're going to branch on the mask, you might as well hard-code the popcnt in each case. Return it in an int * parameter or something. (popcnt came after pshufb, so if you have to fall back to an SSE2 version, you don't have hardware popcnt either.) If SSSE3 pshufb is available, a (data) LUT of shuffle masks may be better if the data is unpredictable. – Peter Cordes May 1 '16 at 15:56
    
Since the pshufb masks have a known relationship inside each group of 4B, they could be compressed from [ D+3 D+2 D+1 D | C+3 ... ] down to just 4B [ D C B A ], and unpacked with punpcklbw same,same / punpcklwd same,same / paddb x, [ 3 2 1 0 | 3 2 1 0 | ... ]. That's 3 shuffles and an add instead of just one pshufb, though. Or unpack the mask with a pshufb, so it's 2 shuffles and a paddb. Anyway, that makes the LUT only 16 * 4B = 64B = one cache line, at the cost of needing two other 16B constants in registers, or as memory operands. – Peter Cordes May 1 '16 at 16:06
    
@PeterCordes, good point about popcnt. I thought of that but the OP did not request it even though it seems like something the OP would need. When I look at the jump table in GCC I can't make out yet how GCC created it. I would expect it jump in the order I listed the codes 0,1,2,3...15 (or maybe reverse order) but it seems random. It's 0xe, 0x34, 0x1, 0x38, 0x2, 0x8,0x9, 0x39, 0x3, 0xc, oxd. If you look at the codes above you can see that that's not at all in order. Do you know why this is? – Z boson May 1 '16 at 18:17
1  
Maybe it started to order it for a decision-tree of branches before deciding on a jump-table strategy. It amuses me that when making PIC code, it decided on a table of 4B displacements that it loads with movsx. If it's going to movsx anyway, might as well use 1B displacements for a smaller table. It also doesn't know that the input will always be 0..15, so it checks for outside that range and returns zero :/ – Peter Cordes May 1 '16 at 19:31
1  
re: hex: you mean like this Godbolt feature-request? Having gcc do it internally would probably be ideal, maybe submitting a patch to gcc would be better than having godbolt post-process the output. Esp. because it would be useful outside of godbolt.org! – Peter Cordes May 2 '16 at 14:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.