How to set bits of a bit vector efficiently in parallel?

Question

Consider a bit vector of N bits in it (N is large) and an array of M numbers (M is moderate, usually much smaller than N), each in range 0..N-1 indicating which bit of the vector must be set to 1. The latter array is not sorted. The bit vector is just an array of integers, specifically __m256i, where 256 bits are packed into each __m256i structure.

How can this work be split efficiently accross multiple threads?

Preferred language is C++ (MSVC++2017 toolset v141), assembly is also great. Preferred CPU is x86_64 (intrinsics are ok). AVX2 is desired, if any benefit from it.

Answer 1

Let's assume you want to divide this work up among T threads. It's a pretty interesting problem since it isn't trivially parallelizable via partitioning and various solutions may apply for different sizes of N and M.

Fully Concurrent Baseline

You could simply divide up the array M into T partitions and have each thread work on its own partition of M with a shared N. The main problem is that since M is not sorted, all threads may access any element of N and hence stomp on each others work. To avoid this, you'd have to use atomic operations such as std::atomic::fetch_or for each modification of the shared N array, or else come up with some locking scheme. Both approaches are likely to kill performance (i.e., using an atomic operation to set a bit is likely to be an order of magnitude slower than the equivalent single-threaded code).

Let's look at ideas that are likely faster.

Private N

One relatively obvious idea to avoid the "shared N" problem which requires atomic operations for all mutations of N is simply to give each T a private copy of N and merge them at the end via or.

Unfortunately, this solution is O(N) + O(M/T) whereas the original single-threaded solution is O(M) and the "atomic" solution above is something like O(M/T)⁴. Since we know that N >> M this is likely to be a poor tradeoff in this case. Still, it's worth noting that the hidden constants in each term are very different: the O(N) term, which comes from the merging step⁰ can use 256-bit wide vpor instructions, meaning a throughput of something close to 200-500 bits/cycle (if cached), while the bit-setting step which is O(M/T) I estimate at closer to 1 bit/cycle. So this approach can certainly be the best one for moderate T even if the size of N is 10 or 100 times the size of M.

Partitions of M

The basic idea here is to partition the indexes in M such that each worker thread can then work on a disjoint part of the N array. If M was sorted, that would be trivial, but it's not, so...

A simple algorithm that will work well if M is smoothly distributed is to first partition that values of M into T buckets, with the buckets having values in the ranges [0, N/T), [N/T, 2N/T], ..., [(T-1)N/T, N). That is, divide N into T disjoint regions and then find the values of M that fall into each of them. You can spread that work across the T threads by assigning each thread an equal size chunk of M, and having them each create the T partitions and then logically merging¹ them at the end so you have the T partitions of M.

The second step is to actually set all the bits: you assign one partition to each thread T which can set the bits in a "single threaded" way, i.e., not worrying about concurrent updates, since each thread is working on a disjoint partition of N².

Both steps O(M) and the second step is identical to the single-threaded case, so the overhead for parallelizing this is the first step. I suspect the first will range from about the same speed as the second to perhaps 2-4 times as slow, depending on implementation and hardware, so you can expect a speedup on a machine with many cores, but with only 2 or 4 it might not be any better.

If the distribution of M is not smooth, such that the partitions created in the first step have very different sizes, it will work poorly because some threads will get a lot more work. A simple strategy is to create say 10 * T partitions, rather than only T and have the threads in the second pass all consume from the same queue of partitions until complete. In this way you spread the work more evenly, unless the array M is very bunched up. In that case you might consider a refinement of the first step which first essentially creates a bucketed histogram of the elements, and then a reduce stage which looks at the combined histogram to create a good partitioning.

Essentially, we are just progressively refining the first stage into a type of parallel sort/partitioning algorithm, for which there is already lots of literature. You might even find that a full (parallel) sort is fastest, since it will greatly help in bit-setting phase, since accesses will be in-order and have the best spatial locality (helping with prefetching and caching, respectively).

---

⁰ ... and also from the "allocate a private array of length N" step, although this is likely to be quite fast.

¹ The conceptually simplest form of merging would be to simply copy each thread's partitions of M such that you have a contiguous partition of all of M, but in practice if the partitions are large you can just leave the partitions where they are and link them together, adding some complexity to the consuming code, but avoiding the compacting step.

² To make it truly disjoint from a threading point of view you want to ensure the partition of N falls on "byte boundaries", and perhaps even cache-line boundaries to avoid false sharing (although the latter is likely not to be a big problem since it only occurs at the edge of each partition, and the order of processing means that you are not likely to get contention).

⁴ In practice, the exact "order" of the baseline concurrent solution using shared N is hard to define because there will be contention so the O(M/T) scaling will break down for large enough T. If we assume N is quite large and T is limited to typical hardware concurrency of at most a dozen cores or so it's probably an OK approximation.

Answer 2

@IraBaxter posted an interesting but flawed idea which can be made to work (at significant cost). I suspect @BeeOnRope's idea of partial-sort / partitioning the M array will perform better (especially for CPUs with large private caches which can keep parts of N hot). I'll summarize the modified version of Ira's idea that I described in comments on his deleted answer. (That answer has some suggestions about how big N has to be before it's worth multi-threading.)

---

Each writer thread gets a chunk of M with no sorting/partitioning.

The idea is that conflicts are very rare because N is large compared to the number of stores that can be in flight at once. Since setting a bit is idempotent, so we can handle conflicts (where two threads want to set different bits in the same byte) by checking the value in memory to make sure it really does have the bit set that we want after a RMW operation like or [N + rdi], al (with no lock prefix).

E.g. thread 1 tried to store 0x1 and stepped on thread 2's store of 0x2. Thread 2 must notice and retry the read-modify-write (probably with lock or to keep it simple and make multiple retries not possible) to end up with 0x3 in the conflict byte.

We need an mfence instruction before the read-back. Otherwise store-forwarding will give us the value we we just wrote before other threads see our store. In other words, a thread can observe its own stores earlier than they appear in the global order. x86 does have a Total Order for stores, but not for loads. Thus, we need mfence to prevent StoreLoad reordering. (Intel's "Loads Are not Reordered with Older Stores to the Same Location" guarantee is not as useful as it sounds: store/reload isn't a memory barrier; they're just talking about out-of-order execution preserving program-order semantics.)

mfence is expensive, but the trick that makes this better than just using lock or [N+rdi], al is that we can batch operations. e.g. do 32 or instructions and then 32 read-back. It's a tradeoff between mfence overhead per operation vs. increased chance of false-sharing (reading back cache lines that had already been invalidated by another CPU claiming them).

Instead of an actual mfence instruction, we can do the last or of a group as a lock or. This is better for throughput on both AMD and Intel. For example, according to Agner Fog's tables, mfence has one per 33c throughput on Haswell/Skylake, where lock add (same performance as or) has 18c or 19c throughput. Or for Ryzen, ~70c (mfence) vs. ~17c (lock add).

If we keep the amount of operations per fence very low, the array index (m[i]/8) + mask (1<<(m[i] & 7)) can be kept in registers for all the operations. This probably isn't worth it; fences are too expensive to do as often as every 6 or operations. Using the bts and bt bit-string instructions would mean we could keep more indices in registers (because no shift-result is needed), but probably not worth it because they're slow.

Using vector registers to hold indices might be a good idea, to avoid having to reload them from memory after the barrier. We want the load addresses to be ready as soon as the read-back load uops can execute (because they're waiting for the last store before the barrier to commit to L1D and become globally visible).

Using single-byte read-modify-write makes actual conflicts as unlikely as possible. Each write of a byte only does a non-atomic RMW on 7 neighbouring bytes. Performance still suffers from false-sharing when two threads modify bytes in the same 64B cache-line, but at least we avoid having to actually redo as many or operations. 32-bit element size would make some things more efficient (like using xor eax,eax / bts eax, reg to generate 1<<(m[i] & 31) with only 2 uops, or 1 for BMI2 shlx eax, r10d, reg (where r10d=1).)

Avoid the bit-string instructions like bts [N], eax: it has worse throughput than doing the indexing and mask calculation for or [N + rax], dl. This is the perfect use-case for it (except that we don't care about the old value of the bit in memory, we just want to set it), but still its CISC baggage is too much.

In C, a function might look something like

/// UGLY HACKS AHEAD, for testing only.

//    #include <immintrin.h>
#include <stddef.h>
#include <stdint.h>
void set_bits( volatile uint8_t * restrict N, const unsigned *restrict M, size_t len)
{
    const int batchsize = 32;

    // FIXME: loop bounds should be len-batchsize or something.
    for (int i = 0 ; i < len ; i+=batchsize ) {
        for (int j = 0 ; j<batchsize-1 ; j++ ) {
           unsigned idx = M[i+j];
           unsigned mask = 1U << (idx&7);
           idx >>= 3;
           N[idx] |= mask;
        }

        // do the last operation of the batch with a lock prefix as a memory barrier.
        // seq_cst RMW is probably a full barrier on non-x86 architectures, too.
        unsigned idx = M[i+batchsize-1];
        unsigned mask = 1U << (idx&7);
        idx >>= 3;
        __atomic_fetch_or(&N[idx], mask, __ATOMIC_SEQ_CST);
        // _mm_mfence();

        // TODO: cache `M[]` in vector registers
        for (int j = 0 ; j<batchsize ; j++ ) {
           unsigned idx = M[i+j];
           unsigned mask = 1U << (idx&7);
           idx >>= 3;
           if (! (N[idx] & mask)) {
               __atomic_fetch_or(&N[idx], mask, __ATOMIC_RELAXED);
           }
        }
    }
}

This compiles to approximately what we want with gcc and clang. The asm (Godbolt) could be more efficient in several ways, but might be interesting to try this. This is not safe: I just hacked this together in C to get the asm I wanted for this stand-alone function, without inlining into a caller or anything. __atomic_fetch_or is not a proper compiler barrier for non-atomic variables the way asm("":::"memory") is. (At least the C11 stdatomic version isn't.) I should probably have used the legacy __sync_fetch_and_or, which is a full barrier for all memory operations.

It uses GNU C atomic builtins to do atomic RMW operations where desired on variables that aren't atomic_uint8_t. Running this function from multiple threads at once would be C11 UB, but we only need it to work on x86. I used volatile to get the asynchronous-modification-allowed part of atomic without forcing N[idx] |= mask; to be atomic. The idea is to make sure that the read-back checks don't optimize away.

I use __atomic_fetch_or as a memory barrier because I know it will be on x86. With seq_cst, it probably will be on other ISAs, too, but this is all a big hack.

Answer 3

There are a couple of operations involved in sets (A,B = set, X = element in a set):

Set operation           Instruction
---------------------------------------------
Intersection of A,B     A and B
Union of A,B            A or B
Difference of A,B       A xor B
A is subset of B        A and B = B     
A is superset of B      A and B = A       
A <> B                  A xor B <> 0
A = B                   A xor B = 0
X in A                  BT [A],X
Add X to A              BTS [A],X
Subtract X from A       BTC [A],X

Given the fact that you can use the boolean operators to replace set operations you can use VPXOR, VPAND etc.
To set, reset or test individual bits you simply use

mov eax,BitPosition
BT [rcx],rax

You can set if a set is (equal to) empty (or something else) using the following code

vpxor      ymm0,ymm0,ymm0       //ymm0 = 0
//replace the previous instruction with something else if you don't want
//to compare to zero.
vpcmpeqqq  ymm1,ymm0,[mem]      //compare mem qwords to 0 per qword
vpslldq    ymm2,ymm1,8          //line up qw0 and 1 + qw2 + 3
vpand      ymm2,ymm1,ymm2       //combine qw0/1 and qw2/3
vpsrldq    ymm1,ymm2,16         //line up qw0/1 and qw2/3
vpand      ymm1,ymm1,ymm2       //combine qw0123, all in the lower 64 bits.
//if the set is empty, all bits in ymm1 will be 1.
//if its not, all bits in ymm1 will be 0.     

(I'm sure this code can be improved using the blend/gather etc instructions) From here you can just extend to bigger sets or other operations.

Note that bt, btc, bts with a memory operand is not limited to 64 bits.
The following will work just fine.

mov eax,1023
bts [rcx],rax   //set 1024st element (first element is 0).