Algorithms, data structures, and the line of thinking

A GPU algorithm is not a loop. It is a data-ownership plan: who owns which output, where partial results meet, and how expensive that meeting is. Learn to translate any serial algorithm into independent ownership plus a combine step.

The core reframing

A CPU programmer asks one question: "What loop do I write?" Control is implicit and free — the program counter walks the data and the cache hides latency. On a GPU that question is the wrong one. You have tens of thousands of threads in flight and no cheap global ordering, so you instead ask three questions:

• What work is independent? Independent work is free parallelism. The more arithmetic a thread does before it has to talk to anyone, the better.
• Where do partial results meet? The meeting point is the only place ordering and communication exist. Push it as late and as local as possible.
• How expensive is that meeting? A meeting inside a warp is nearly free (shuffles). Inside a block it costs a barrier. Across blocks it costs an atomic or — more often — a whole new kernel launch.

Reframed this way, "writing the algorithm" becomes "deciding ownership and pricing the combine." The loop is an afterthought the hardware generates for you.

An algorithm pattern map

Almost every parallel kernel is one of seven shapes. Memorize the ownership and the combine column — that pairing is the answer to "how would you parallelize X?"

Notice the pattern: ownership is always "thread / warp / block / grid owns one unit of output," and the combine is always "none / tree / scan / atomic / second launch." Two small vocabularies generate every kernel.

Reduction: many local accumulators, then a tree

Reduction (sum, max, dot product, norm) is the canonical combine. The serial loop for i: acc += x[i] has a single accumulator and a strict dependence chain. The GPU answer inverts that: give every thread its own accumulator, let each chew through many elements with a grid-stride loop (fat independent phase), then fold the thousands of partials with a tree.

The tree has three tiers priced by communication cost:

• Within a warp — __shfl_down_sync. No shared memory, no barrier. Nearly free. This is warp_reduce_sum from lesson 02.
• Within a block — one partial per warp lands in shared memory, then warp 0 reduces those with another warp_reduce_sum. Costs one __syncthreads.
• Across blocks — and here is the catch.

Here is the block-level kernel. It produces exactly one partial per block. warp_reduce_sum is defined in lesson 02 — signature shown for context:

// from lesson 02:
__inline__ __device__ float warp_reduce_sum(float v);   // folds 32 lanes -> lane 0

__global__ void reduce_partials(const float* __restrict__ x, int n,
                                float* __restrict__ partial) {
    // 1. fat independent phase: each thread sums many elements
    float sum = 0.0f;
    for (int i = blockIdx.x * blockDim.x + threadIdx.x;
         i < n;
         i += blockDim.x * gridDim.x) {
        sum += x[i];
    }

    // 2. warp-level combine (free)
    sum = warp_reduce_sum(sum);

    // 3. block-level combine: one slot per warp in shared memory
    __shared__ float warp_sums[32];           // up to 1024 threads = 32 warps
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    if (lane == 0) warp_sums[warp] = sum;
    __syncthreads();

    // warp 0 folds the per-warp partials into one block partial
    if (warp == 0) {
        int n_warps = (blockDim.x + 31) >> 5;
        float v = (lane < n_warps) ? warp_sums[lane] : 0.0f;
        v = warp_reduce_sum(v);
        if (lane == 0) partial[blockIdx.x] = v;   // one partial per block
    }
}

The host then launches the same kernel a second time on the partial array (now small enough for one block) to get the final scalar. Two alternatives avoid the second launch:

• atomicAdd to one global accumulator. Replace the last line with if (...) atomicAdd(out, v);. One atomic per block — cheap because there are few blocks — but the result is nondeterministic in floating point (addition order varies) and you must zero out first.
• Cooperative-groups grid sync. Launch with cudaLaunchCooperativeKernel and call grid.sync() from cooperative_groups. This does give a true grid barrier inside one kernel, but it caps occupancy (all blocks must be co-resident) and needs hardware/driver support. Use it only when relaunch overhead dominates.

Prefix scan (exclusive): "how many valid items came before me?"

An exclusive prefix scan turns an array into the running total excluding the current element. It looks academic until you see what it answers: "for element i, how many qualifying items appeared strictly before me?" That number is a write address. Scan is the engine behind stream compaction, radix sort, run-length encoding, sparse-matrix construction, and token packing in LLM pipelines (packing variable-length sequences into a dense buffer).

index :   0   1   2   3   4   5   6   7
flags :   1   0   1   1   0   0   1   1     (keep this element?)
                |
                v   exclusive scan (sum of all flags before i)
scan  :   0   1   1   2   3   3   3   4
                |
                v   for each kept element, scan[i] is its output slot
out   : [ x0 | x2 | x3 | x6 | x7 ]          total kept = scan[7]+flags[7] = 5

Within a single block, the simplest correct scan is Hillis-Steele: at step d, each element adds the value 2^d positions to its left. After log2(n) steps every element holds its inclusive prefix; shift right by one for exclusive. We double-buffer two shared arrays to avoid a read-after-write race within a step:

// single-block inclusive->exclusive scan, n <= blockDim.x
__global__ void scan_block(const int* __restrict__ in, int* __restrict__ out, int n) {
    extern __shared__ int tmp[];     // size 2*n
    int t = threadIdx.x;
    int pin = 0, pout = 1;           // ping-pong buffer indices

    // load shifted right by one -> makes the result EXCLUSIVE
    tmp[pin * n + t] = (t > 0) ? in[t - 1] : 0;
    __syncthreads();

    for (int offset = 1; offset < n; offset <<= 1) {
        // swap buffers each step so reads and writes never alias
        pout = 1 - pout;
        pin  = 1 - pin;
        int v = tmp[pin * n + t];
        if (t >= offset) v += tmp[pin * n + (t - offset)];
        tmp[pout * n + t] = v;
        __syncthreads();
    }
    out[t] = tmp[pout * n + t];       // exclusive prefix
}

For arrays larger than one block, the recipe is three passes — and it is exactly the "scan, scan the totals, add the offset back" structure:

1. Scan each block locally and write each block's total sum to an array block_totals[blockIdx.x].
2. Scan block_totals (it is small — recurse or use one block). Now block_totals[b] holds the count of everything in blocks 0..b-1.
3. Add the scanned offset back: every element in block b adds block_totals[b]. This stitches the per-block scans into one global scan.

Stream compaction: scan builds a write-address table

Compaction collects the elements that pass a predicate into a dense output. The naive serial version is if (keep(x[i])) out[j++] = x[i]; — but j++ is a serial dependency; every thread would fight over j. The fix: precompute all the j values at once with an exclusive scan of the keep-flags. The scan is the table of output addresses.

// flags[i] = 1 if keep, else 0.  pos = exclusive scan of flags.
__global__ void compact(const float* __restrict__ x,
                         const int*   __restrict__ flags,
                         const int*   __restrict__ pos,   // exclusive scan of flags
                         float* __restrict__ out, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n && flags[i]) {
        out[pos[i]] = x[i];        // pos[i] is a deterministic, collision-free slot
    }
}

Every kept element writes to a unique slot — no atomics, no contention, fully deterministic order. This is the key trick: a scan converts data-dependent branching ("where does this go?") into a precomputed scatter address. The same idea places elements into radix buckets during sort and writes packed tokens into a contiguous KV buffer.

CSR sparse matrix-vector multiply (SpMV)

Most large matrices are mostly zeros, so we store only the nonzeros in Compressed Sparse Row (CSR) format: three arrays where row_ptr marks where each row begins in the flat col_idx / value arrays.

dense A (4x4):              CSR:
  10  0  0  0               row_ptr : [0, 1, 3, 4, 6]   (length rows+1)
   0 20 30  0               col_idx : [0, 1, 2, 3, 0, 3]
   0  0  0 40               value   : [10,20,30,40,50,60]
  50  0  0 60
                            row r's nonzeros: indices row_ptr[r] .. row_ptr[r+1]-1
                            e.g. row 1: cols {1,2} vals {20,30}

SpMV computes y[r] = sum over nonzeros in row r of value * x[col]. Ownership choice is the whole game. The vector (warp-per-row) kernel assigns one warp to one row: the 32 lanes stride across the row's nonzeros (coalesced loads of value and col_idx), then a warp reduction sums their contributions.

__global__ void spmv_csr_vector(const int*   __restrict__ row_ptr,
                                const int*   __restrict__ col_idx,
                                const float* __restrict__ value,
                                const float* __restrict__ x,
                                float* __restrict__ y, int rows) {
    int warp_id = (blockIdx.x * blockDim.x + threadIdx.x) >> 5;
    int lane    = threadIdx.x & 31;
    if (warp_id >= rows) return;

    int start = row_ptr[warp_id];
    int end   = row_ptr[warp_id + 1];

    float sum = 0.0f;
    for (int j = start + lane; j < end; j += 32) {     // lanes stride the row
        sum += value[j] * x[col_idx[j]];
    }
    sum = warp_reduce_sum(sum);                        // lesson 02
    if (lane == 0) y[warp_id] = sum;                   // one writer per row
}

Top-k per row for small k: a data structure in registers

When you need the k largest values per row and k is tiny (say k ≤ 8 — common for attention top-k, beam search, recommendation), a full sort is wasteful. Keep a sorted array of size K in registers and insert each new element with a single insertion step: if it beats the current minimum, shift and drop the loser. That is O(n·k) work in fast registers versus O(n log n) for a full sort touching memory — and for tiny k the constant is unbeatable.

template <int K>
__device__ void topk_insert(float (&best)[K], float v) {
    if (v <= best[K - 1]) return;          // smaller than our worst kept value
    int i = K - 1;
    while (i > 0 && best[i - 1] < v) {       // shift smaller entries right
        best[i] = best[i - 1];
        --i;
    }
    best[i] = v;                           // insert in sorted position
}

template <int K>
__global__ void topk_per_row(const float* __restrict__ a, int rows, int cols,
                             float* __restrict__ out) {   // out: rows * K
    int r = blockIdx.x * blockDim.x + threadIdx.x;
    if (r >= rows) return;

    float best[K];
    #pragma unroll
    for (int i = 0; i < K; ++i) best[i] = -INFINITY;

    const float* row = a + (size_t)r * cols;
    for (int c = 0; c < cols; ++c) topk_insert<K>(best, row[c]);

    #pragma unroll
    for (int i = 0; i < K; ++i) out[(size_t)r * K + i] = best[i];
}

The best[] array, because K is a compile-time template parameter, lives entirely in registers — no shared memory, no global traffic until the final write. If instead a warp cooperates on one row (each lane keeps its own top-k), you do a final merge of the 32 register heaps via shuffles. The line of thinking: the data structure lives in the fastest memory that fits, and only merges at the very end.

Sorting note

You will rarely hand-write a GPU sort, but you must know the two shapes:

• Radix sort — the default GPU sort. It is built directly on prefix scan: for each digit, count occurrences, scan the counts to get bucket offsets, then scatter keys to those offsets (the compaction trick generalized to many buckets). Linear work, no comparisons, stable, bandwidth-bound. This is why scan is the foundational primitive.
• Bitonic sort — a fixed network of compare-exchange steps, O(n log²n) work but completely data-independent (no branching, perfect for warps/blocks). Ideal for sorting small fixed-size groups (e.g. sorting within a block or a warp) where its predictable pattern beats radix's setup cost.

In practice: cub::DeviceRadixSort or thrust::sort. They are tuned per architecture and you should not reimplement them.

Graph BFS: frontier expansion with atomic claims

Breadth-first search over a graph stored in CSR adjacency (same row_ptr / col_idx layout, rows = vertices, entries = neighbors) proceeds in levels. Each level expands the current frontier: for every vertex in it, look at its neighbors, claim the unvisited ones, and append them to the next frontier. Two atomics carry the correctness:

__global__ void bfs_expand(const int* __restrict__ row_ptr,
                           const int* __restrict__ col_idx,
                           const int* __restrict__ frontier, int frontier_n,
                           int* __restrict__ visited,
                           int* __restrict__ next, int* __restrict__ next_n) {
    int t = blockIdx.x * blockDim.x + threadIdx.x;
    if (t >= frontier_n) return;

    int u = frontier[t];
    for (int e = row_ptr[u]; e < row_ptr[u + 1]; ++e) {
        int v = col_idx[e];
        // claim v exactly once: only the thread that flips 0->1 succeeds
        if (atomicCAS(&visited[v], 0, 1) == 0) {
            int slot = atomicAdd(next_n, 1);   // reserve a unique queue slot
            next[slot] = v;
        }
    }
}

The correctness comes entirely from atomicCAS(&visited[v], 0, 1): many threads may reach the same neighbor v in the same level, but the compare-and-swap succeeds for exactly one of them — that one (and only that one) appends v to next. atomicAdd(next_n, 1) hands out collision-free queue slots, the same address-table idea as compaction but computed on the fly.

A GPU hash table

The same atomic-claim idea builds a concurrent hash table. Use open addressing with linear probing: hash the key to a slot; if taken, walk forward until you find an empty slot, claiming it with atomicCAS so concurrent inserters never collide.

__device__ const unsigned EMPTY = 0xffffffffu;

__global__ void hash_insert(unsigned* __restrict__ keys, int capacity,
                            const unsigned* __restrict__ in, int n) {
    int t = blockIdx.x * blockDim.x + threadIdx.x;
    if (t >= n) return;
    unsigned k = in[t];
    unsigned slot = (k * 2654435761u) & (capacity - 1);   // capacity is power of two
    for (int probe = 0; probe < capacity; ++probe) {
        unsigned prev = atomicCAS(&keys[slot], EMPTY, k);
        if (prev == EMPTY || prev == k) return;           // inserted, or already present
        slot = (slot + 1) & (capacity - 1);               // linear probe to next slot
    }
}

When to use libraries

The honest senior answer is "reach for the library first." But you must understand the kernel underneath to choose, size, and debug it.

Interview answer template

When asked "how would you implement / parallelize X on a GPU?", walk this ladder out loud. It signals that you think in ownership and combines, not loops:

1. Start from the data. What is the layout, and is it dense, sparse, or skewed? Where is the reuse? (This decides everything downstream.)
2. Choose ownership. Assign one writer per output unit at the right granularity — thread, warp, block, or grid — matched to the work per output and the access pattern (for coalescing).
3. Name the combine primitive. State it explicitly: none (map), reduction (tree), scan (prefix → scatter), atomic (histogram / queue claim), or a second kernel launch (grid sync).
4. Predict the bottleneck. Usually the combine: global atomics, the launch boundary, load imbalance from skew, or uncoalesced scatter. Say which one and why.
5. Optimize locally or fall back. Collapse the expensive step — privatize atomics into shared memory, push the meeting point into the warp, rebalance ownership for skew — or, if it is a standard primitive, call CUB / cuBLAS / cuSPARSE and explain that beating it is rarely worth it.