More Parallel Patterns

2-D/3-D stencils (heat equation, Jacobi, finite differences, image blur) follow the same pattern with tiles + halos; register tiling along the third dimension for 3-D.

15.1 Histogram (privatized in shared memory)

__global__ void histogram256(const unsigned char* data, int n, unsigned int* hist) {
    __shared__ unsigned int local[256];
    for (int i = threadIdx.x; i < 256; i += blockDim.x) local[i] = 0;
    __syncthreads();
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x)
        atomicAdd(&local[data[i]], 1);                    // fast shared-memory atomics
    __syncthreads();
    for (int i = threadIdx.x; i < 256; i += blockDim.x)
        atomicAdd(&hist[i], local[i]);                    // merge into the global histogram
}

15.2 Stencil (1-D with halo)

2-D/3-D stencils (heat equation, Jacobi, finite differences, image blur) follow the same pattern with tiles + halos; register tiling along the third dimension for 3-D.

#define RADIUS 3
#define BLOCK 256
__global__ void stencil1d(const float* in, float* out, int n) {
    __shared__ float s[BLOCK + 2 * RADIUS];
    int g = blockIdx.x * blockDim.x + threadIdx.x;
    int l = threadIdx.x + RADIUS;
    s[l] = (g < n) ? in[g] : 0.0f;
    if (threadIdx.x < RADIUS) {                            // load halos
        s[l - RADIUS] = (g >= RADIUS) ? in[g - RADIUS] : 0.0f;
        s[l + BLOCK]  = (g + BLOCK < n) ? in[g + BLOCK] : 0.0f;
    }
    __syncthreads();
    if (g < n) {
        float r = 0.0f;
        for (int o = -RADIUS; o <= RADIUS; ++o) r += s[l + o];
        out[g] = r / (2 * RADIUS + 1);
    }
}

15.3 2-D convolution (image filtering) with constant memory

#define K 5
__constant__ float d_filter[K * K];

__global__ void conv2d(const float* img, float* out, int w, int h) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= w || y >= h) return;
    float acc = 0.0f;
    for (int fy = 0; fy < K; ++fy)
        for (int fx = 0; fx < K; ++fx) {
            int ix = min(max(x + fx - K / 2, 0), w - 1);    // clamp at borders
            int iy = min(max(y + fy - K / 2, 0), h - 1);
            acc += img[iy * w + ix] * d_filter[fy * K + fx];
        }
    out[y * w + x] = acc;
}
cudaMemcpyToSymbol(d_filter, h_filter, sizeof(float) * K * K);

15.4 Stream compaction (filter)

Compute flags (keep[i] = pred(in[i])), 2. exclusive scan of flags → output positions, 3. scatter kept elements. Or thrust::copy_if / cub::DeviceSelect::If.

15.5 Sorting

Radix sort — the fastest GPU sort for integers/floats: per-digit histogram + scan + scatter (cub::DeviceRadixSort, thrust::sort).

Bitonic sort / merge sort for small or key-value data; segmented sorts (cub::DeviceSegmentedSort).

15.6 Other patterns

Map, gather/scatter, transpose, SpMV (sparse matrix-vector, CSR/ELL formats), merge, BFS/graph traversal (frontier-based), N-body (all-pairs with shared-memory tiles), Monte Carlo (cuRAND), dynamic programming (wavefront), k-means, matrix factorization.