Designing Scan algorithms and studying circuit diagrams is instructive, but in order to implement Scan for CUDA, we need to map the algorithms onto registers, memory and addressing schemes, and correct synchronization.
The optimal CUDA implementation of Scan depends on the size of the scan being performed; different schemes are best for warp-sized scans, scans that can fit in shared memory, and scans that must spill to global memory. Because blocks cannot reliably exchange data through global memory, scans too large to fit in shared memory must perform multiple kernel invocations3.
Before examining special cases (such as scanning of predicates), we build Scan for CUDA from the bottom up. The warp scan is the primitive at the very bottom; a block scan is built on top of it, and three device-wide strategies on top of that:
Scan-then-fan, a recursive decomposition that serves as a lead-in;
Reduce-then-scan, which trades a global memory pass for a reduction pass (Merrill’s two-pass formulation); and
Single-pass decoupled look-back, the strategy that the CUB and Thrust libraries use today.
The primitive at the bottom of every Scan implementation in this chapter is the warp scan: an inclusive scan across the 32 threads of a warp. Because the threads of a warp execute together and can exchange data without barriers, warp scans are fast, and the block- and device-wide scans that follow all use one as their innermost step.
The warp scan is modeled on the Kogge-Stone circuit (Figure 13-7). Kogge-Stone circuits are work-inefficient – they perform many operations for their small depth – but at the warp level, where the hardware’s execution resources are available regardless whether the developer uses them, that inefficiency is free and the shallow depth wins.
Listing 13-3 gives the warp scan as a __device__ routine that
operates on shared memory, the fastest way for threads to exchange data
with one another. It is packaged not as a free function but as a
policy: a struct, WarpScanShared, with one static
method, inclusive(), that takes a thread’s value and a
pointer into shared scratch and returns the scanned result. The reason
for the struct, rather than a plain function, becomes clear once the
second implementation is in hand.
struct WarpScanShared { // Kogge-Stone in shared, Volta-clean static constexpr const char *name = "WarpScanShared"; template<class T> static inline __device__ T inclusive( T val, volatile T *scr ) { const int lane = threadIdx.x & 31; scr[0] = val; T t = val; #pragma unroll for ( int offset = 1; offset < 32; offset <<= 1 ) { if ( lane >= offset ) t += scr[-offset]; __syncwarp(); scr[0] = t; __syncwarp(); } return t; }};WarpScanShared: warp scan in shared
memory (source on GitHub)Older implementations went to some length to shave instructions off this shared-memory warp scan – for example, interleaving each warp’s data with a block of zeros so the lane-ID conditionals could be removed. Since SM 3.0 (2013), the warp shuffle instruction has made those tricks obsolete: it exchanges registers directly between the threads of a warp, with no shared memory at all. Its “up” and “down” variants implement scan and reverse scan.
Listing 13-4 gives the second policy, WarpScanShuffle.
Each of its five doubling steps calls __shfl_up_sync() to
fetch the partial sum from the lane offset positions below,
and a lane-ID test guards the add so that the bottom lanes, which have
no lower neighbor, keep their own value. The scan touches no shared
memory and lives entirely in registers.
struct WarpScanShuffle { // register-only, __shfl_up_sync static constexpr const char *name = "WarpScanShuffle"; template<class T> static inline __device__ T inclusive( T val, volatile T * /*unused*/ ) { const int lane = threadIdx.x & 31; #pragma unroll for ( int offset = 1; offset < 32; offset <<= 1 ) { T n = __shfl_up_sync( 0xffffffffu, val, offset ); if ( lane >= offset ) val += n; } return val; }};WarpScanShuffle: register-only warp
scan (source on GitHub)Listing 13-5 gives the block-wide scan,
blockScanInclusive(), written once as a function template
parameterized on the warp policy. It scans each warp with
WarpPolicy::inclusive(), records each warp’s total in
shared memory, scans those per-warp totals with a single shuffle pass,
and adds each warp’s exclusive prefix back into its lane. Selecting the
warp scan through a template type parameter – not a runtime flag or a
function pointer – lets it inline with no indirection, and because a
function template cannot be partially specialized, a policy type is the
idiomatic way to choose among implementations of one algorithm. An
exclusive scan needs no code of its own:
warpScanExclusive() returns the inclusive result minus the
thread’s own value.
// Block-wide inclusive scan, composed on the SAME policy -- the warp scan inlines// (no function pointer). scr: blockDim ints of shared, warp-contiguous so the// shared policy's neighbour reads stay within a warp's own run.template<class WarpPolicy, class T>inline __device__ TblockScanInclusive( T val, volatile T *scr ){ const int tid = threadIdx.x, lane = tid & 31, warpid = tid >> 5; const int nwarps = blockDim.x >> 5; __shared__ T warpTotals[32]; T v = WarpPolicy::template inclusive<T>( val, scr + tid ); if ( 31 == lane ) warpTotals[warpid] = v; __syncthreads(); if ( 0 == warpid ) { // exclusive scan of the per-warp totals T w = ( lane < nwarps ) ? warpTotals[lane] : (T)0, inc = w; #pragma unroll for ( int offset = 1; offset < 32; offset <<= 1 ) { T n = __shfl_up_sync( 0xffffffffu, inc, offset ); if ( lane >= offset ) inc += n; } warpTotals[lane] = inc - w; } __syncthreads(); return v + warpTotals[warpid];}blockScanInclusive(): block scan over a
warp policy (source on GitHub)The scan-then-fan approach uses a similar decomposition for global and shared memory. Figure 13-8 shows the approach used to scan a threadblock: a scan is performed on each 32-thread warp, and the reduction of that 32-element subarray is written to shared memory. A single warp then scans the array of partial sums; a single warp is sufficient because CUDA does not support threadblocks with more than 1024 threads. Finally, the base sums are fanned out to each warp’s output elements. Note that Figure 13-8 shows an inclusive scan being performed in step 2), so the first element of its output must be fanned out to the second warp, and so on.
Figure 13-8. Scan-then-fan (shared memory)
The code to implement this algorithm is given in Listing 13-6. The
input array is assumed to have been loaded into shared memory already,
and the parameters sharedPartials and idx
specify the base address and index of the warp to scan, respectively.
(In our first implementation, threadIdx.x is passed as the parameter
idx.) Lines 9-13 implement Step 1 of Figure 13-8; lines
16-21 implement Step 2; and lines 31-45 implement Step 3. The output
value written by this thread is returned to the caller, but used only if
it happens to be the thread block’s reduction.
template<class T>inline __device__ TscanBlock( volatile T *sPartials ){ extern __shared__ T warpPartials[]; const int tid = threadIdx.x; const int lane = tid & 31; const int warpid = tid >> 5; // // Compute this thread's partial sum // T sum = scanWarp<T>( sPartials ); __syncthreads(); // // Write each warp's reduction to shared memory // if ( lane == 31 ) { warpPartials[16+warpid] = sum; } __syncthreads(); // // Have one warp scan reductions // if ( warpid==0 ) { scanWarp<T>( 16+warpPartials+tid ); } __syncthreads(); // // Fan out the exclusive scan element (obtained // by the conditional and the decrement by 1) // to this warp's pending output // if ( warpid > 0 ) { sum += warpPartials[16+warpid-1]; } __syncthreads(); // // Write this thread's scan output // *sPartials = sum; __syncthreads(); // // The return value will only be used by caller if it // contains the spine value (i.e. the reduction // of the array we just scanned). // return sum;}scanBlock(): Block portion of
scan-then-fan for thread blocks (source on GitHub)Figure 13-9 shows how this approach is adapted to global memory: a kernel scans b-element subarrays where b is the block size; the partial sums are written to global memory and another, 1-block kernel invocation scans these partial sums, which are then fanned into the final output in global memory.
Figure 13-9. Scan-then-fan (global memory)
Listing 13-7 gives the CUDA code for the scan kernel of Step 1 in Figure 13-9. It loops over the threadblocks to process, staging the input array into and out of shared memory. The kernel then optionally writes the spine value to global memory at the end. At the bottom level of the recursion, there is no need to record spine values, so the bWriteSpine template parameter enables the kernel to avoid dynamically checking the value of partialsOut.
template<class T, bool bWriteSpine>__global__ voidscanAndWritePartials( T *out, T *gPartials, const T *in, size_t N, size_t numBlocks ){ extern volatile __shared__ T sPartials[]; const int tid = threadIdx.x; volatile T *myShared = sPartials+tid; for ( size_t iBlock = blockIdx.x; iBlock < numBlocks; iBlock += gridDim.x ) { size_t index = iBlock*blockDim.x+tid; *myShared = (index < N) ? in[index] : 0; __syncthreads(); T sum = scanBlock( myShared ); __syncthreads(); if ( index < N ) { out[index] = *myShared; } // // write the spine value to global memory // if ( bWriteSpine && (threadIdx.x==(blockDim.x-1)) ) { gPartials[iBlock] = sum; } }}scanAndWritePartials() (source on GitHub)Listing 13-8 gives the host function that uses Listings 13-6 and 13-7
to implement an inclusive scan on an array in global memory. Note that
the function recurses for scans too large to perform in shared memory;
the first conditional in the function serves both as the base case for
the recursion and to short-circuit scans small enough to perform in
shared memory alone, avoiding any need to allocate global memory. Note
how the amount of shared memory needed by the kernel
(b*sizeof(T)) is specified at kernel invocation time.
For larger scans, the function computes the number of partial sums
needed \(\left\lceil \frac{N}{b}
\right\rceil\), allocates global memory to hold them, and follows
the pattern of Figure 13-9, writing partial sums to the global array for
later use by the scanAndWritePartials() kernel of Listing 13-7.
Each level of recursion reduces the number of elements being processed by a factor of b, so for e.g. b=128 and N= 1048576, two levels of recursion are required: one of size 8192 and a second of size 64.
The partials array is transient scratch – allocated and freed within
a single scanFan call, once per level of the recursion.
scanFan draws it from the stream-ordered allocator of
Section 5.2.3 (cudaMallocAsync/cudaFreeAsync)
rather than cudaMalloc, so each allocate and free is
ordered on the stream instead of synchronizing the device, and the
blocks come from a pool that a repeated scan reuses rather than
re-obtaining from the driver.
template<class T>voidscanFan( T *out, const T *in, size_t N, int b ){ cudaError_t status_cudart; if ( N <= b ) { scanAndWritePartials<T, false><<<1,b,b*sizeof(T)>>>( out, 0, in, N, 1 ); return; } // // device pointer to array of partial sums in global memory // T *gPartials = 0; // // ceil(N/b) // size_t numPartials = (N+b-1)/b; // // one thread block per b-element tile. The kernels are // grid-stride, so any block count within CUDA's limits would // also work; one block per tile keeps the mapping simple. // size_t numBlocks = numPartials; // // The partials array is transient scratch: it lives only for the // duration of this (recursive) call. Allocate it stream-ordered from // the default stream's memory pool so the driver can satisfy the request // by recycling scratch freed at a shallower recursion level, rather than // synchronizing the device on every cudaMalloc/cudaFree. // cuda(MallocAsync( &gPartials, numPartials*sizeof(T), 0 ) ); scanAndWritePartials<T, true><<<numBlocks,b,b*sizeof(T)>>>( out, gPartials, in, N, numPartials ); scanFan<T>( gPartials, gPartials, numPartials, b ); scanAddBaseSums<T><<<numBlocks, b>>>( out, gPartials, N, numPartials ); Error_cudart: if ( gPartials ) cudaFreeAsync( gPartials, 0 );}scanFan() host function. (source on GitHub)Listing 13-9 completes the picture with a very simple kernel to fan-out results from global memory to global memory.
template<class T>__global__ voidscanAddBaseSums( T *out, T *gBaseSums, size_t N, size_t numBlocks ){ const int tid = threadIdx.x; T fan_value = 0; for ( size_t iBlock = blockIdx.x; iBlock < numBlocks; iBlock += gridDim.x ) { size_t index = iBlock*blockDim.x+tid; if ( iBlock > 0 ) { fan_value = gBaseSums[iBlock-1]; } out[index] += fan_value; }}scanAddBaseSums() kernel (source on GitHub)At the highest level of recursion, the scan-then-fan strategy performs 4N global memory operations: the initial scan performs one read and one write, then the fan of Listing 13-9 performs another read and write. We can decrease the number of global memory operations by first computing only reductions on the input array, deferring the scanning to a second pass.
Figure 13-10 shows how this strategy works. As before, an array of \(\left\lceil \frac{N}{b} \right\rceil\) partial sums of the input is computed and scanned to compute an array of base sums; but instead of doing the scan in the first pass, the first pass computes only the reductions. The scan of the final output is then performed, adding the base sum along the way. This performs 3N global memory operations – one read per element in the reduction pass, then a read and a write in the scan pass – as against 4N for scan-then-fan.
Figure 13-10. Reduce-then-scan.
Merrill4 describes a refinement of this
strategy that uses a small, fixed-size number of base sums. The
algorithm is the same as Figure 13-10, except that the array of Step 2)
is a relatively small, fixed size of perhaps a few hundred instead of
\(\left\lceil \frac{N}{b}
\right\rceil\) partial sums. The number of partial sums is the
same as the number of thread blocks to use, both for the reduction pass
and for the Scan pass. Listing 13-10 shows the code to compute these
partial sums, which computes reductions for subarrays of size
elementsPerPartial as opposed to the thread block size.
template<class T, int numThreads>__device__ voidscanReduceSubarray( T *gPartials, const T *in, size_t iBlock, size_t N, int elementsPerPartial ){ extern volatile __shared__ T sPartials[]; const int tid = threadIdx.x; size_t baseIndex = iBlock*elementsPerPartial; T sum = 0; for ( int i = tid; i < elementsPerPartial; i += blockDim.x ) { size_t index = baseIndex+i; if ( index < N ) sum += in[index]; } sPartials[tid] = sum; __syncthreads(); reduceBlock<T,numThreads>( &gPartials[iBlock], sPartials );} /* * Compute the reductions of each subarray of size * elementsPerPartial, and write them to gPartials. */template<class T, int numThreads>__global__ voidscanReduceSubarrays( T *gPartials, const T *in, size_t N, int elementsPerPartial ){ extern volatile __shared__ T sPartials[]; for ( int iBlock = blockIdx.x; iBlock*elementsPerPartial < N; iBlock += gridDim.x ) { scanReduceSubarray<T,numThreads>( gPartials, in, iBlock, N, elementsPerPartial ); }}scanReduceSubarrays() (source on GitHub)Listing 13-11 gives the Scan code, which has been modified to carry over each block’s sum as the Scan of that block is completed.
template<class T>__global__ voidscan2Level_kernel( T *out, const T *gBaseSums, const T *in, size_t N, size_t elementsPerPartial ){ extern volatile __shared__ T sPartials[]; const int tid = threadIdx.x; int sIndex = (threadIdx.x); T base_sum = 0; if ( blockIdx.x && gBaseSums ) { base_sum = gBaseSums[blockIdx.x-1]; } for ( size_t i = 0; i < elementsPerPartial; i += blockDim.x ) { size_t index = blockIdx.x*elementsPerPartial + i + tid; sPartials[sIndex] = (index < N) ? in[index] : 0; __syncthreads(); scanBlock<T>( sPartials+sIndex ); __syncthreads(); if ( index < N ) { out[index] = sPartials[sIndex]+base_sum; } __syncthreads(); // carry forward from this block to the next. base_sum += sPartials[ (blockDim.x-1) ]; __syncthreads(); }}scan2Level_kernel() (source on GitHub)Listing 13-12 gives the host code for Merrill’s two-pass
reduce-then-scan algorithm. Since the number of partials computed is
small and never varies, the host code never has to allocate global
memory in order to perform the scan – instead, we declare a
__device__ array that is allocated at module load time:
__device__ int g_globalPartials[MAX_PARTIALS];
and obtain its address by calling cudaGetSymbolAddress():
status = cudaGetSymbolAddress(
(void **) &globalPartials,
g_globalPartials );
The routine then computes the number of elements per partial and number of threadblocks to use and invokes the three (3) kernels needed to perform the computation.
template<class T>voidscan2Level( T *out, const T *in, size_t N, int b ){ int sBytes = ((b)*sizeof(T)); if ( N <= b ) { return scan2Level_kernel<T><<<1,b,sBytes>>>( out, 0, in, N, N ); } cudaError_t status_cudart; T *gPartials = 0; cuda(GetSymbolAddress( (void **) &gPartials, g_globalPartials )); { // // ceil(N/b) = number of partial sums to compute // size_t numPartials = (N+b-1)/b; if ( numPartials > MAX_PARTIALS ) { numPartials = MAX_PARTIALS; } // // elementsPerPartial has to be a multiple of b // unsigned int elementsPerPartial = (N+numPartials-1)/numPartials; elementsPerPartial = b * ((elementsPerPartial+b-1)/b); numPartials = (N+elementsPerPartial-1)/elementsPerPartial; // // number of CUDA threadblocks to use. The kernels are // blocking agnostic, so we can clamp to any number within // CUDA's limits and the code will work. // const size_t maxBlocks = MAX_PARTIALS; size_t numBlocks = std::min( numPartials, maxBlocks ); scanReduceSubarrays<T>( gPartials, in, N, elementsPerPartial, numBlocks, b ); scan2Level_kernel<T><<<1,b,sBytes>>>( gPartials, 0, gPartials, numPartials, numPartials ); scan2Level_kernel<T><<<numBlocks,b,sBytes>>>( out, gPartials, in, N, elementsPerPartial ); }Error_cudart:;}scan2Level (source on GitHub)Both strategies so far make more than one pass over the input in global memory: scan-then-fan performs 4N operations and reduce-then-scan performs 3N. The extra passes exist because thread blocks cannot share their partial results without returning to the host to launch another kernel. If the blocks could hand their partial sums to one another while the kernel is still running, a single pass – one read and one write per element, 2N operations – would suffice.
Merrill and Garland’s decoupled look-back algorithm does exactly this, and it is the method that CUB and Thrust use to implement device-wide scan. One thread block processes one tile of b elements in a single pass. The blocks are ordered, and each computes the inclusive scan of its own tile in shared memory (using the warp-scan machinery of Section 13.3.1); the last element of that local scan is the tile’s aggregate. To turn its local scan into a global one, a tile needs the exclusive prefix – the sum of all elements in the tiles before it – and computing that prefix is what the look-back does.
The obstacle is the prefix dependency: a tile’s exclusive prefix is its immediate predecessor’s inclusive prefix, so the most literal implementation would have each tile wait for the one before it to finish, serializing the entire grid. The insight that breaks the serialization is that a tile can publish its aggregate – which it computes on its own, with no predecessor – long before it knows its own prefix. A waiting successor can then make progress by summing aggregates rather than blocking on a finished prefix. The two computations are thereby decoupled: prefixes propagate down the grid while aggregates are still being produced.
Each tile advertises its progress through a one-word status descriptor in global memory that moves through three states:
X (invalid): the initial state; the tile has published nothing.
A (aggregate available): the tile has computed its own aggregate but not yet its inclusive prefix.
P (inclusive prefix available): the tile knows the sum of every element up to and including its own tile.
A tile computes its exclusive prefix by looking back over
its predecessors, nearest first (Figure 13-11). The look-back is done
cooperatively by one warp: the 32 lanes read the descriptors of the 32
nearest predecessors at once, a __ballot finds the nearest
lane whose tile has published an inclusive prefix (state
P), and a warp reduction sums the values from the
frontier back to it – the aggregates of the tiles in between, plus that
prefix, which already folds in everything before it. If no lane in the
window holds a P, the warp sums all 32 aggregates and
steps back another 32 predecessors. Because the walk stops at the first
P – and prefixes are continually being filled in behind
it – it is short in practice, usually a single 32-wide step.
Figure 13-11. Single-pass decoupled look-back. Tile 5 has published its aggregate (A) and looks back nearest-first, summing the aggregates of tiles 4 and 3 until it reaches tile 2’s already-computed inclusive prefix (P), which it adds before publishing its own prefix. Tiles ahead (X) have not started.
Two details make the descriptor cheap and the algorithm safe. First, the status flag and its 32-bit value are packed into a single 64-bit word, written and read with one atomic operation. Because a reader that observes the flag necessarily observes the matching value, no separate memory fence is needed between them. This costs one bit of generality: it assumes a 32-bit element type and a commutative, associative operator (here, integer addition). A 64-bit element type needs a two-word descriptor and an explicit fence, which is what CUB implements.
Second, a tile that waits on a predecessor must be certain that
predecessor is resident and running, or the look-back could spin forever
against a block the scheduler has not launched. The kernel therefore
claims tile indices from a global atomic counter on entry rather than
using blockIdx.x. A block that later waits on predecessor
p claimed a larger index than p did, so p
began executing first; the look-back cannot deadlock.
The descriptor’s encoding and its three status codes are given in Listing 13-13.
enum { SCAN_X = 0, SCAN_A = 1, SCAN_P = 2 }; // invalid / aggregate / prefix // 64-bit descriptor word. Deliberately unsigned long long, not uint64_t: CUDA's// 64-bit atomics (atomicOr/atomicExch) are overloaded only for unsigned long// long, and uint64_t is a distinct type (unsigned long on LP64) that won't bind.typedef unsigned long long scanStatus; __device__ __forceinline__ scanStatusscanPackStatus( uint32_t flag, int value ){ return ( (scanStatus) flag << 32 ) | (uint32_t) value;}The kernel itself is given in Listing 13-14: each block claims a tile, inclusive-scans it in shared memory, and its first warp performs the cooperative look-back before the block writes its outputs – all in a single launch.
template<class T>__global__ voidscanDecoupledLookback_kernel( T *out, const T *in, volatile scanStatus *status, // one descriptor per tile (SCAN_X-initialized) uint32_t *tileCounter, // one global counter, 0-initialized size_t N ){ extern __shared__ T s[]; // blockDim.x elements __shared__ uint32_t s_tile; __shared__ T s_exclusive; // // Claim a tile index from a global counter. A block that later waits on // predecessor p is then guaranteed p is already resident and running (p // claimed a smaller index earlier), so the look-back cannot deadlock. // if ( threadIdx.x == 0 ) { s_tile = atomicAdd( tileCounter, 1 ); } __syncthreads(); const uint32_t tile = s_tile; const size_t base = (size_t) tile * blockDim.x; const size_t gidx = base + threadIdx.x; // // Load the tile (zero-filling past the end) and inclusive-scan it in shared // memory (Kogge-Stone). s[blockDim.x-1] then holds the tile's aggregate. // s[threadIdx.x] = ( gidx < N ) ? in[gidx] : (T) 0; __syncthreads(); for ( int off = 1; off < blockDim.x; off <<= 1 ) { T add = ( threadIdx.x >= off ) ? s[threadIdx.x - off] : (T) 0; __syncthreads(); s[threadIdx.x] += add; __syncthreads(); } const T aggregate = s[blockDim.x - 1]; // // Compute this tile's exclusive prefix. Tile 0 has no predecessors; every // other tile advertises its aggregate, then warp 0 looks back cooperatively // over 32 predecessors at a time. // if ( tile == 0 ) { if ( threadIdx.x == 0 ) { atomicExch( (scanStatus *) &status[tile], scanPackStatus( SCAN_P, aggregate ) ); s_exclusive = (T) 0; } } else { // Advertise our aggregate so successors can make progress without us. if ( threadIdx.x == 0 ) atomicExch( (scanStatus *) &status[tile], scanPackStatus( SCAN_A, aggregate ) ); if ( threadIdx.x < 32 ) { const int lane = threadIdx.x; T exclusive = (T) 0; int frontier = (int) tile - 1; // nearest predecessor not yet consumed for ( ;; ) { // Lane L inspects predecessor (frontier - L): the 32 nearest // predecessors at once, lane 0 == nearest. const int pidx = frontier - lane; uint32_t flag; T val; do { // re-read the window until no lane sees X if ( pidx >= 0 ) { scanStatus d = atomicOr( (scanStatus *) &status[pidx], 0ull ); flag = (uint32_t) ( d >> 32 ); val = (T) (int) ( d & 0xffffffffu ); } else { flag = SCAN_P; // out of range: prefix 0, stops the walk val = (T) 0; } } while ( __any_sync( 0xffffffffu, flag == SCAN_X ) ); // Nearest inclusive prefix in the window (lowest lane with P). const uint32_t pmask = __ballot_sync( 0xffffffffu, flag == SCAN_P ); const int firstP = pmask ? ( __ffs( pmask ) - 1 ) : 32; // Sum lanes 0..firstP: aggregates up to it, plus its prefix. T contrib = ( lane <= firstP ) ? val : (T) 0; #pragma unroll for ( int off = 16; off > 0; off >>= 1 ) contrib += __shfl_xor_sync( 0xffffffffu, contrib, off ); exclusive += contrib; if ( firstP < 32 ) // reached a prefix -> done break; frontier -= 32; // whole window was aggregates -> keep walking } if ( lane == 0 ) { s_exclusive = exclusive; atomicExch( (scanStatus *) &status[tile], scanPackStatus( SCAN_P, exclusive + aggregate ) ); } } } __syncthreads(); // // Add the exclusive prefix to the local inclusive scan and write the output. // if ( gidx < N ) { out[gidx] = s_exclusive + s[threadIdx.x]; }}scanDecoupledLookback_kernel() (source on GitHub)Listing 13-15 gives the host function, which allocates the per-tile
status array (zero-initialized to X) and the tile
counter – both stream-ordered from the memory pool, as
scanFan allocates its partials in Section 13.3.2 – then
makes the one kernel launch.
template<class T>voidscanDecoupledLookback( T *out, const T *in, size_t N, int b ){ cudaError_t status_cudart; scanStatus *gStatus = 0; uint32_t *tileCounter = 0; if ( N == 0 ) return; uint32_t numTiles = (uint32_t) ( ( N + b - 1 ) / b ); // // gStatus and tileCounter are transient per-call scratch. Allocate and // zero them stream-ordered on the default stream so repeated calls (e.g. // a timing loop) recycle the same pool memory instead of paying for a // synchronizing cudaMalloc/cudaFree on every invocation. // cuda(MallocAsync( &gStatus, numTiles * sizeof(scanStatus), 0 ) ); cuda(MemsetAsync( gStatus, 0, numTiles * sizeof(scanStatus), 0 ) ); // SCAN_X == 0 cuda(MallocAsync( &tileCounter, sizeof(uint32_t), 0 ) ); cuda(MemsetAsync( tileCounter, 0, sizeof(uint32_t), 0 ) ); scanDecoupledLookback_kernel<T><<<numTiles, b, b * sizeof(T)>>>( out, in, gStatus, tileCounter, N ); Error_cudart: if ( gStatus ) cudaFreeAsync( gStatus, 0 ); if ( tileCounter ) cudaFreeAsync( tileCounter, 0 );}scanDecoupledLookback() host
function. (source on GitHub)Touching global memory only twice per element gives the algorithm the lowest memory traffic of the three, but the kernel above does not yet realize that: with one element per thread and a per-lane serial load, it is bound by low per-thread work rather than by bandwidth. Section 13.3.5 closes that gap.
The decoupled-look-back kernel of the previous section has the least
memory traffic of any scan in this chapter, yet on large inputs it is
slower than the reduce-then-scan of Section 13.3.3. The reason is not
the algorithm but the implementation: one element per thread never gives
the memory system enough in-flight work to reach peak bandwidth, and the
Kogge-Stone block scan spends much of its time in
__syncthreads().
The fix is to increase the amount of work per thread. Listing 13-16
gives an optimized kernel in which each thread sequentially scans a
small contiguous chunk of IPT elements in registers; a
warp-shuffle block scan (Section 13.3.1) then stitches the per-thread
sums together, and the same cooperative tile look-back runs on top.
Loads and stores are coalesced, the register scan has ample
instruction-level parallelism, and the larger tile (blockDim.x ×
IPT elements) means far fewer tiles and far fewer
look-backs.
template<class T, int IPT>__global__ voidscanDecoupledLookback2_kernel( T *out, const T *in, volatile scanStatus *status, uint32_t *tileCounter, size_t N ){ const int B = blockDim.x; extern __shared__ T s[]; // B*IPT elements, logical order __shared__ uint32_t s_tile; __shared__ T s_base; __shared__ T warpsum[32]; if ( threadIdx.x == 0 ) s_tile = atomicAdd( tileCounter, 1 ); __syncthreads(); const uint32_t tile = s_tile; const size_t tileBase = (size_t) tile * B * IPT; // Coalesced striped load; s[k] holds logical element k of the tile. #pragma unroll for ( int i = 0; i < IPT; i++ ) { int k = i * B + threadIdx.x; size_t g = tileBase + k; s[k] = ( g < N ) ? in[g] : (T) 0; } __syncthreads(); // Per-thread sequential inclusive scan of a contiguous chunk (in registers). T chunk[IPT]; const int c0 = threadIdx.x * IPT; T run = (T) 0; #pragma unroll for ( int i = 0; i < IPT; i++ ) { run += s[c0 + i]; chunk[i] = run; } const T threadSum = run; // Block exclusive scan of threadSum via warp shuffles -> offset, aggregate. const int lane = threadIdx.x & 31, wid = threadIdx.x >> 5, numWarps = B >> 5; T x = threadSum; #pragma unroll for ( int off = 1; off < 32; off <<= 1 ) { T y = __shfl_up_sync( 0xffffffffu, x, off ); if ( lane >= off ) x += y; } if ( lane == 31 ) warpsum[wid] = x; __syncthreads(); if ( wid == 0 ) { T w = ( lane < numWarps ) ? warpsum[lane] : (T) 0; #pragma unroll for ( int off = 1; off < 32; off <<= 1 ) { T y = __shfl_up_sync( 0xffffffffu, w, off ); if ( lane >= off ) w += y; } if ( lane < numWarps ) warpsum[lane] = w; } __syncthreads(); const T warpOffset = ( wid == 0 ) ? (T) 0 : warpsum[wid - 1]; const T offset = ( x - threadSum ) + warpOffset; // block-exclusive prefix const T aggregate = warpsum[numWarps - 1]; scanCoopLookback<T>( status, tile, aggregate, s_base ); __syncthreads(); // Add tile prefix + block offset, store coalesced. const T base = s_base + offset; #pragma unroll for ( int i = 0; i < IPT; i++ ) s[c0 + i] = chunk[i] + base; __syncthreads(); #pragma unroll for ( int i = 0; i < IPT; i++ ) { int k = i * B + threadIdx.x; size_t g = tileBase + k; if ( g < N ) out[g] = s[k]; }}scanDecoupledLookback2_kernel():
items-per-thread optimization. (source on GitHub)Table 13-1 shows the effect, measured on an NVIDIA RTX 3060 (peak
memory bandwidth ≈ 360 GB/s) scanning 64M 32-bit integers. Making the
look-back cooperative lifts the teaching kernel from 80 to 125 GB/s;
increasing the work to eight items per thread lifts it to 330 GB/s –
matching CUB’s DeviceScan, which implements the same
algorithm.
| Implementation | GB/s | vs. CUB |
|---|---|---|
| Decoupled look-back, serial look-back | 80 | 0.24× |
| Decoupled look-back, cooperative look-back | 125 | 0.38× |
| Reduce-then-scan (Section 13.3.3) | 177 | 0.53× |
| Decoupled look-back, 8 items/thread | 330 | 1.00× |
CUB
DeviceScan::InclusiveSum |
331 | 1.00× |
Thrust inclusive_scan |
316 | 0.95× |
Table 13-1. Inclusive scan of 64M 32-bit integers on an RTX 3060.
The lesson is the one this chapter opened with: the
decoupled-look-back algorithm is what lets a library reach peak
bandwidth, but realizing that potential takes a carefully engineered
implementation. For production scans, call CUB’s DeviceScan
or Thrust’s inclusive_scan and exclusive_scan,
which implement this algorithm with the full generality – arbitrary
types and operators, 64-bit descriptors – and tuning elided here.
With SM 3.5 hardware, dynamic parallelism can move most of the kernel launches to be “child grids” as opposed to kernel launches initiated by the host.↩︎
Merrill, Duane and Andrew Grimshaw. “Parallel Scan for Stream Architectures.” Technical Report CS2009-14, Department of Computer Science, University of Virginia, December 2009.↩︎