We’ve explored CUDA implementations that operate on one histogram per grid (the output histogram, which we operate on with global atomics) and one histogram per block (which we operate on with shared memory atomics). An even finer granularity – one histogram per thread – is a natural next step to explore, and in fact is the only viable way to parallelize histograms on multicore CPUs, since atomic additions are so expensive compared to incrementing histogram elements that are in each CPU core’s cache hierarchy (each core’s L1 cache has latency of 3-4 clock cycles). A CPU-optimized, multithreaded implementation would spawn M CPU threads for each of the M cores on the system, decompose the problem into M chunks and compute a histogram for each using a fork/join idiom, then combine them into a single histogram at the end. When there is one histogram per thread, the intermediate histograms are known in the literature as “privatized” histograms, presumably because each is private to a thread and not because they are scheduled for divestment by the government.
The CUDA Handbook source code on GitHub includes a multithreaded implementation of the CPU code. On an 8-core AMD Ryzen 7 7700X, the single-threaded version runs at about 885 Mpix/s and the multithreaded version at about 2990 Mpix/s—only about 3.4x faster despite the eight cores. The histogram is memory-bandwidth bound, so a few fast modern cores already saturate the memory system well before the core count would otherwise allow linear scaling. Even at 2990 Mpix/s, this is far slower than the GPU’s per-block implementation, which exceeds 100,000 Mpix/s.
Since privatized, per-thread histograms work well on CPUs, that leads us to wonder if they also might work well in CUDA. The goal is less data-dependent performance; many applications would give up some performance to the per-block or per-grid formulations in exchange for not suffering an order-of-magnitude performance degradation for degenerate inputs.
In CUDA, when allocating one histogram per thread, shared memory is the logical choice since the registers cannot be referenced by index. On the hardware for which this code was first written, shared memory was the binding constraint: with at most 48K per SM, even short histogram elements added up quickly. With 8-bit histogram elements, a 64-thread block uses 256*64=16384 bytes of shared memory—just a bit too much for SM 1.x class hardware (which reserves 256 bytes of the 16K of shared memory for parameter passing), and enough to limit SM 2.0 and later parts to 3 blocks per SM, making for low occupancy. Modern GPUs are far more generous: an Ampere SM can be configured with up to 100K of shared memory, so that same 16K privatized histogram no longer caps occupancy the way it once did, and larger or more numerous per-thread histograms have become practical. (A side benefit of spending shared memory rather than registers is that threads remain free to use many registers without further reducing occupancy.)
Privatization was also, for a time, the best answer on the GPU. Before Maxwell (SM 5.0) gave shared memory a native atomic unit, a shared memory atomic was synthesized from a lock-and-retry sequence whose cost grew with contention, so the per-block formulation degraded on degenerate input while a per-thread formulation stayed level. Section 16.4 reports what the hardware did to that argument: on current parts the per-block implementation is both the faster and the more level of the two. The technique is presented here for the reasoning it contains—packed counters, shared memory layout, and overflow budgeting all apply to any contended accumulator the hardware does not absorb—rather than as the implementation to reach for today.
The 256*NumThreads elements of the privatized histograms may be laid out in shared memory in any number of ways. Figures 16-4 and 5 show two possibilities: a histogram per row (each thread operates on its own row) and a histogram element per row (each thread operates on its own column). These layouts both suffer from poor performance due to using 8-bit memory operands in shared memory.
Figure 16-4. Histogram per row
Figure 16-5. Histogram per column
Figure 16-6 shows a hybrid layout that we settled on: 64 rows that each contain an interleaved set of 4 packed 8-bit counters, one per thread. If we happen to use 64 threads per block, the array of 32-bit integers coincidentally will be square (64x64).
Figure 16-6. Histogram per column (interleaved)
Listing 16-5 shows a device function that increments a privatized histogram element using the layout of Figure 16-6. This operation uses 32-bit shared memory accesses and minimizes bank conflicts for completely uniform data, in which case the threads increment adjacent 32-bit shared memory locations. As a result, performance is less likely to degrade due to contention.
inline __device__ voidincPacked32Element( unsigned char pixval ){ extern __shared__ unsigned int privHist[]; const int blockDimx = 64; unsigned int increment = 1<<8*(pixval&3); int index = pixval>>2; privHist[index*blockDimx+threadIdx.x] += increment;}Once the privatized histograms have been accumulated in shared memory, each histogram element must be reduced across the threads of the block and added to the output. The packed representation makes that reduction cheap: the four 8-bit partial sums in a 32-bit word split into two pairs of 16-bit partial sums with two masks and a shift, as in Listing 16-6. Up to 256 such partial sums can be added together before risking overflow, so a block of fewer than 256 threads is safe.
unsigned int myValue = privHist[i*64+threadIdx.x]; int sum02, sum13; sum02 = myValue & 0xff00ff; myValue >>= 8; sum13 = myValue & 0xff00ff;The full reduction accumulates those pairs down the 64 columns of the
privatized array and fires four atomics per thread into the output
histogram. To avoid bank conflicts in shared memory, the thread ID is
added to the loop index and masked to the range 0..63.
merge64HistogramsToOutput() in
histogram/histogramPerThread64.cuh gives the
implementation.
On its own, the code we’ve discussed will not work correctly if any thread increments the same histogram element more than 256 times. To avoid overflow, there are three basic strategies that can be employed:
Overflow Check: check for overflow after incrementing the histogram element, and fire an atomic into either shared memory or global memory;
Overflow Avoidance: ensure that no thread considers more inputs than any single histogram element can accommodate; or
Periodic Merge: periodically merge the privatized histograms into the output histogram, and zero them out.
We found that the first strategy was too slow – checking for overflow for every input element introduced too much overhead into the inner loop.
The second strategy, overflow avoidance, is easier to implement than
one might expect, because the number of input elements considered by
each thread is proportional to the number of thread blocks in the kernel
launch. If w*h is the number of pixels and
numthreads is the number of threads per block, we can
write:
int numblocks = INTDIVIDE_CEILING( w*h, numthreads*255 );
and no thread will consider more than 255 elements5. If the inner loop is unrolled 4x, as previously described, the number of blocks must be increased accordingly:
int numblocks = INTDIVIDE_CEILING( w*h, numthreads*(255/4) );
The disadvantage of overflow avoidance is that it can be too conservative: every single block must merge its privatized histograms into the final output, and that extra effort is wasted if very few histogram elements ran the risk of overflow.
We found that periodically merging the privatized histograms into the
output, just often enough to prevent any histogram element from
overflowing, was measurably faster than avoiding overflow by launching
more thread blocks. histogram1DPerThread4x64() implements
both policies behind a template parameter, so the two can be measured
against one another.