parallel stencil computation

**Parallel Stencil Computation** is the **high-performance computing pattern for algorithms where each output element is computed from a fixed neighborhood (stencil) of input elements** — the core computational pattern in finite difference methods for PDEs, image processing kernels, cellular automata, and lattice physics simulations. Stencil computations are memory-bandwidth bound and require carefully designed tiling and communication strategies to achieve high efficiency on multi-core CPUs, GPUs, and distributed clusters. **What Is a Stencil?** - A stencil defines which neighboring grid points contribute to the computation of each output point. - **1D 3-point stencil**: `out[i] = a*in[i-1] + b*in[i] + c*in[i+1]`. - **2D 5-point stencil (von Neumann)**: `out[i,j] = in[i-1,j] + in[i+1,j] + in[i,j-1] + in[i,j+1] − 4*in[i,j]`. - **3D 7-point stencil**: Laplacian in 3D → heat equation, Poisson's equation. - **High-order stencil**: 25-point or 49-point stencil → more neighbors → higher accuracy, more flops, more memory traffic. **Arithmetic Intensity of Stencil Computations** - 3D 7-point stencil: 7 loads + 1 store per output point, 7 FMAs → AI ≈ 7/8 ≈ 0.88 FLOP/byte. - Modern GPUs: Ridge point ~30 FLOP/byte → stencil is extremely memory bandwidth bound. - CPU cache blocking: Reuse data from cache → AI increases → approaches cache bandwidth limit. **GPU Stencil Implementation** ```cuda __global__ void stencil_3d_7pt(float* out, const float* in, int N) { int i = blockIdx.x*blockDim.x + threadIdx.x; int j = blockIdx.y*blockDim.y + threadIdx.y; int k = blockIdx.z*blockDim.z + threadIdx.z; if (i>0 && i0 && j0 && k