scan algorithm
**Parallel Scan (Prefix Sum) Algorithms** are the **fundamental parallel computation patterns that compute all prefix reductions of an array — where output[i] = op(input[0], input[1], ..., input[i])** — transforming a seemingly sequential problem into a parallel one with O(N) work in O(log N) steps, serving as the building block for stream compaction, radix sort, histogram computation, and sparse matrix operations across all parallel computing platforms.
**Inclusive vs. Exclusive Scan**
```
Input: [3, 1, 7, 0, 4, 1, 6, 3]
Inclusive scan: [3, 4, 11, 11, 15, 16, 22, 25]
output[i] = sum(input[0..i])
Exclusive scan: [0, 3, 4, 11, 11, 15, 16, 22]
output[i] = sum(input[0..i-1]), output[0] = identity
```
**Naive Parallel Scan (Hillis-Steele)**
```
Step 0: [3, 1, 7, 0, 4, 1, 6, 3] (original)
Step 1: [3, 4, 8, 7, 4, 5, 7, 9] (add offset 1)
Step 2: [3, 4, 11, 11, 12, 12, 11, 14] (add offset 2)
Step 3: [3, 4, 11, 11, 15, 16, 22, 25] (add offset 4)
```
- O(N log N) work, O(log N) steps.
- NOT work-efficient (does more total work than sequential).
- But: All N processors are busy at every step → simple to implement.
**Work-Efficient Scan (Blelloch)**
**Up-Sweep (Reduce) Phase:**
```
[3, 1, 7, 0, 4, 1, 6, 3]
[3, 4, 7, 7, 4, 5, 6, 9] d=1: pairs at distance 1
[3, 4, 7, 11, 4, 5, 6, 14] d=2: pairs at distance 2
[3, 4, 7, 11, 4, 5, 6, 25] d=4: root has total sum
```
**Down-Sweep Phase:**
```
[3, 4, 7, 11, 4, 5, 6, 0] set root to 0
[3, 4, 7, 0, 4, 5, 6, 11] d=4: swap and add
[3, 4, 7, 0, 4, 5, 6, 11] d=2: swap and add at pairs
[0, 3, 4, 11, 11, 15, 16, 22] d=1: final exclusive scan
```
- O(N) work (same as sequential), O(log N) steps.
- Work-efficient → better utilization of parallel resources.
- Two passes (up-sweep + down-sweep) vs. one pass for Hillis-Steele.
**GPU Implementation (Block-Level)**
```cuda
__global__ void scan(float *output, float *input, int n) {
__shared__ float temp[2 * BLOCK_SIZE];
int tid = threadIdx.x;
// Load to shared memory
temp[2*tid] = input[2*tid];
temp[2*tid+1] = input[2*tid+1];
// Up-sweep (reduce)
for (int d = BLOCK_SIZE; d > 0; d >>= 1) {
__syncthreads();
if (tid < d) {
int ai = 2 * d * (tid + 1) - 1;
int bi = ai - d;
temp[ai] += temp[bi];
}
}
// Set root to 0 for exclusive scan
if (tid == 0) temp[2*BLOCK_SIZE - 1] = 0;
// Down-sweep
for (int d = 1; d < 2 * BLOCK_SIZE; d <<= 1) {
__syncthreads();
if (tid < d) {
int ai = 2 * d * (tid + 1) - 1;
int bi = ai - d;
float t = temp[bi];
temp[bi] = temp[ai];
temp[ai] += t;
}
}
__syncthreads();
output[2*tid] = temp[2*tid];
output[2*tid+1] = temp[2*tid+1];
}
```
**Multi-Block Scan (Large Arrays)**
1. Each block scans its portion → produces block prefix sums.
2. Scan the block totals (small array) → produces block offsets.
3. Each block adds its offset to all its elements.
- Three kernel launches for arbitrary-size arrays.
- CUB library: cub::DeviceScan::ExclusiveSum() handles all this.
**Applications of Prefix Sum**
| Application | How Scan Is Used |
|------------|------------------|
| Stream compaction | Scan flag array → compute output indices |
| Radix sort | Scan digit histograms → compute scatter addresses |
| Sparse matrix (SpMV) | Scan row pointers → determine output ranges |
| Histogram | Segmented scan over sorted keys |
| Memory allocation | Scan sizes → compute offsets for variable-size outputs |
| Run-length encoding | Scan to compute output positions |
Parallel scan is **the algorithmic foundation that makes irregularity tractable in parallel computing** — by converting the inherently sequential prefix computation into a parallel-friendly tree reduction, scan enables efficient parallelization of workloads where output positions depend on input data, making it the second most important parallel primitive after reduction and the key enabler of parallel sorting, compaction, and sparse data processing on GPUs.