The CUDA Handbook ▸ Chapter 16. Histograms and Radix Sort

16.1 One Global Histogram, and Contention

For CUDA, the tradeoffs are more complex, since the hardware supports more threads and provides mechanisms for enforcing mutual exclusion in both global and shared memory. Listing 16-2 gives a CUDA C implementation that looks similar to the serial CPU implementation of Listing 16-1; it uses global memory atomics to operate directly on the output histogram².

__global__ voidhistogram1DPerGrid(    unsigned int *pHist,    const unsigned char *base, size_t N ){    for ( size_t i = blockIdx.x*blockDim.x+threadIdx.x;                 i < N;                 i += blockDim.x*gridDim.x ) {        atomicAdd( &pHist[ base[i] ], 1 );    }}

Global memory atomics were introduced in SM 1.1, but they were so slow as to be almost unusable. For randomly distributed data, this kernel on flagship chips for Tesla (GeForce GTX 280), Fermi (Tesla M2050) and Kepler (GeForce GRID K520) yields the performance results summarized in Table 16-2.

Community medians; full data and raw logs: /benchmarks/histogram/.

Table 16-2. Global Atomics Performance, SM 1.0 through Ampere

You read that correctly: across the SM 1.x-3.x era, SM 2.0 increased performance by 26x over SM 1.3, and SM 3.0 increased performance again by another 7x³. Apparently NVIDIA thought developers wanted fast atomics! (and NVIDIA would be correct.)

The Ampere entry breaks the trend, though: on a mid-range GeForce RTX 3060, this per-grid kernel manages only about 3500 Mpix/s, which is actually much slower than the Kepler part⁴. The shortfall is not an artifact of feeding the GPU too little work; the throughput holds essentially constant from one million to 256 million pixels, so the kernel is genuinely limited by contention for the 256 shared histogram elements rather than by the GPU’s ability to keep its cores busy. Heavily contended global atomics are a pathological case that faster hardware cannot rescue—which is exactly why the rest of this chapter develops per-block and per-thread strategies that avoid the contention in the first place.

One downside of a histogram implementation that uses global atomics is data-dependent behavior: if all the input values are the same, the kernel of Listing 16-2 will perform N atomic adds on the same 32-bit memory location. In this case, the hardware facilities that ensure mutual exclusion for atomics suffer from contention. Our test program can measure the effects of contention by reducing the number of possible random values (specified by the --random parameter to the test program).

With threads contending for the hardware facilities that enable atomicity in global memory, performance degrades as the number of available memory locations goes down, as shown in Table 16-3.

Table 16-3. Performance v. # of Values (Mpix/s; Ampere column is a GeForce RTX 3060)

It’s interesting to note that for Fermi and Kepler, the performance doesn’t decrease monotonically: on Fermi, images with 32 values run faster than images with 64 values, while on Kepler, images with 16 and 8 values exhibit anomalous performance. It is also interesting to note the slowdown in the degenerate case (only a single value in the input image): Tesla was 14.5x slower, Fermi was 22.5x slower, and Kepler was 11x slower. Ampere degrades the least of all—only about 2x from 256 values down to 1—and, like Fermi and Kepler, it does so non-monotonically. It is great that the newer architectures don’t degrade quite as much as the earlier ones, but such data-dependent performance is difficult to avoid in the face of the realities of hardware implementation.

16.1.1 Loop Unrolling

The CUDA kernel of Listing 16-2 uses 8-bit memory operations, which gives much lower performance than 32- or 64-bit memory operations. The time-honored (and, in CUDA, very effective) strategy of loop unrolling can be employed to good effect. The resulting kernel is shown in Listing 16-3. This kernel is more than 20% faster than the kernel of Listing 16-2. Note: N must be a multiple of 4.

This combination of loop unrolling and larger memory operands is used throughout the rest of this chapter; it is even more effective when the operations being unrolled expose more instructions for the compiler to schedule.

__global__ voidhistogram1DNaiveAtomic(    unsigned int *pHist,    const unsigned char *base, size_t N ){    for ( size_t i = blockIdx.x*blockDim.x+threadIdx.x;                 i < N/4;                 i += blockDim.x*gridDim.x ) {        unsigned int value = ((unsigned int *) base)[i];        atomicAdd( &pHist[ value & 0xff ], 1 ); value >>= 8;        atomicAdd( &pHist[ value & 0xff ], 1 ); value >>= 8;        atomicAdd( &pHist[ value & 0xff ], 1 ); value >>= 8;        atomicAdd( &pHist[ value ] , 1 );    }}

The performance of Listings 16-2 and 3 is the same for degenerate input data, as the effects of contention become dominant.