Inclusive scan takes a binary associative operator ⊕, and an array of length N
[a0, a1, … aN-1]
and returns the array:
[a0, (a0⊕a1), … (a0⊕a1⊕…⊕aN-1)].
Each element of the output depends on the preceding elements in the input.
Exclusive scan is defined similarly, but shifts the output and uses an identity element that has no effect on a value when ⊕ is performed with it (for example, 0 for integer addition, 1 for multiplication, etc.):
[id⊕, a0, a0⊕a1, … a0⊕a1⊕…⊕aN-2].
Inclusive and exclusive scans can be transformed between one another by adding or subtracting the input array element-by-element, as shown in Figure 13-1.
Figure 13-1. Inclusive and Exclusive Scan
Stream compaction is an operation that separates elements in an array according to a criterion; if a predicate (0 or 1) is computed for each element of the input array to determine whether it should be included in the output stream, then an exclusive scan on the predicates computes the indices of the output elements.
A variation of stream compaction, known as stream splitting, writes the compact output separately for each value of the predicate.
Segmented scan is a variation that takes a set of input flags (one per array element) in addition to the array, and performs scans on the subarrays delineated by the flags.
Due to the importance of the Scan primitive, an enormous amount of effort has been put into developing optimized scan implementations for CUDA. A list of references is given at the end of this chapter. Both the CUDPP and Thrust libraries include families of optimized Scan primitives that use templates for the best tradeoff between generality and performance. All that said, however, applications that use Scan as a primitive usually can benefit from custom implementations that take advantage of specific knowledge about the problem.1
A simple implementation in C++ looks like this:
template<class T>TInclusiveScan( T *out, const T *in, size_t N ){ T sum(0); for ( size_t i = 0; i < N; i++ ) { sum += in[i]; out[i] = sum; } return sum;}
For these scalar implementations, the only difference between inclusive and exclusive scan is that the lines
out[i] = sum;
and
sum += in[i];
are swapped2.
template<class T>TExclusiveScan( T *out, const T *in, size_t N ){ T sum(0); for ( size_t i = 0; i < N; i++ ) { out[i] = sum; sum += in[i]; } return sum;}
The serial implementations of Scan are so obvious and trivial that you can be forgiven if you’re wondering what a parallel implementation would look like! The so-called prefix dependency, where each output depends on all of the preceding inputs, may have some wondering if it’s possible. But, upon reflection, you can see that the operations for neighboring pairs (ai⊕ai+1 for 0≤i<N-1) could be computed in parallel; for i=0, ai⊕ai+1 computes a final output of the Scan, and otherwise these pairwise operations compute partial sums that can be used to contribute to the final output, much as we used partial sums in Chapter 12 (Reduction).
Parallel prefix sum has been studied for decades – originally as a hardware problem in the design of fast adders. Blelloch describes a tree-based, two-pass formulation: an upsweep phase that computes the reduction of the array, storing intermediate results along the way, followed by a downsweep that distributes those results back into the leaves to produce the final output. It is the classic work-efficient scan, performing O(N) operations, and early CUDA implementations followed it closely.
As a device-wide GPU algorithm, however, the tree-based scan is a poor fit: a naïve implementation suffers from shared memory bank conflicts, and the addressing schemes that compensate for them incur enough overhead that the costs almost outweigh the benefits. Rather than pursue any one algorithm in the abstract, the next section represents scans as circuits, which lays out the whole design space compactly and makes the tradeoffs visible.
The summed-area table of Section 15.8 is a case in point. A summed-area table is a two-dimensional inclusive scan – an inclusive scan along each row, then down each column – and the fastest implementation is a bespoke decoupled look-back kernel that computes both the first and second moments in one pass over the image.↩︎
As written, the implementation of exclusive scan does
not support an in-place computation. To enable the input and output
arrays to be the same, in[i] must be saved in a temporary
variable.↩︎