In all versions of a stencil code, a specific algorithm is applied to every element in a matrix in a sweep. Newly calculated values are piped forward into the next sweep. Cells within the same sweep have no dependencies and are calculated in parallel. Sweeps can also be calculated in parallel if the results from one sweep are continuously piped into the next.

In the 1D representation of a 2D matrix, assuming the length of a row is N and the number of columns is M , the 4-point neighborhood/stencil of any individual point is given by : {x-N,x+N,x-1,x+1}, where all values outside of the boundary of the matrix are assumed to be 0. Assume we have matrix A with values {x₀, x₁, x₂, x₃, …}, where x is in the range 0 to NxM . The matrix initially begins at timestep T₀. The values at T₁ are defined as: T₁(x_i) = 0.25 * (T₀(x_i+1) + T₀(x_i-1) + T₀(x_i+N) + T₀(x_i-N)).

Now assume we have a 2D representation of a 2D matrix B with values {x₀y₀, x₁y₀, x₀y₁, x₁y₁, …} where x in the range of 0 to N, and y is in the range of 0 to M. The matrix initially begins at timestep T₀. The values at T₁ are defined as follows:

T₁(x_iy_j) = 0.25 * (T₀(x_iy_j-1) + T₀(x_iy_j+1) + T₀(x_i-1y_j) + T₀(x_i+1y_j))

The main idea behind this specific stencil pattern is to iteratively update the data within a matrix until the values converge (Jacobi Iteration). Performance is optimized by using caching and feed-forward techniques to reduce to reads from global memory (DDR) as much as possible, ideally to 1. Data is typically retrieved from real world measurements, but for the sake of testing all data in the associated reference design is pseudorandom.

Stencil computations in general are memory-bound applications. The optimizations included in this OpenCL* ^{1 2} reference design seek to leverage the power of an FPGA by both parallelizing as much as possible and using channels/pipes to saturate the bandwidth from on-board DDR to maximize GFLOPS and minimize execution time.

Data is initially read from global memory by the compute unit (CU) named "feeder" and written back to global memory by the "writer"’ CU. Data is read and written in blocks of 16 floats at a time, and the net data bandwidth depends on the frequency of the kernel and supported data read/write rate of the on-board DDR and memory controller.

The following diagram captures the flow of data into and out of the kernel system:

Between the feeder and writer CUs is the main part of the kernel system – a series of chained and replicated calculation CUs. Data flows into the first calculation CU in the same order it was read from DDR, and the calculated data is piped into the next CU in the chain in the same order. Data coming into a CU is cached, boundaries are updated, and the result of each calculation is immediately sent into the next kernel. In this way every single CU, and hence every iteration, can be calculated in parallel once enough data has been sent through the chain.

Within each of these calculation CUs is a local memory system called "cache" that is composed of M20k blocks adjacent to the CUs. The cache size must be large enough to store an entire row of the incoming matrix. The height of the matrix can be as large as device memory allows. Matrix attributes are fed forward through the system before stencil computations begin.

Because of the nature of the computation, only 14 of the 16 floating point values being forwarded can be calculated at a time, so boundary conditions are updated between CUs. This is because each element requires both its left and right neighbors to perform a stencil computation, so the first and last elements in the block cannot be calculated. The entire computation requires a total of 3 blocks of 16 floats to be loaded into private registers, the row being calculated and their top/bottom neighbors. After a calculation is performed one block with 14 new and 2 outdated values is piped to the next CU.

Thirty kernels were chained together in a feed-forward approach in order to perform 30 iterations of the stencil algorithm in parallel. Each individual kernel began execution as soon as it was sent enough information from the previous kernel.

The following heat maps show the convergence of values after running the kernel. The first pair of images represents the unprocessed raw input. The heat maps illustrate 280x280 grids, where the values in the left maps were initially created with a pseudorandom function, and the values in the right maps were created by looping from 1 to 100.

To calculate the execution time for a system requiring more than 30 iterations, divide the desired number of iterations by 30 multiplied by the time it takes 30 iterations to run ($(desirednumberofiterations÷30)\times \left(timefor30iterationtorun\right)$) For example, 300 iterations applied to a 4088x65536 data set should take around 1658 ms to run.

The following heat map shows the result of running the initial values through a kernel with 30 chained calculation CUs. You can see the beginnings of convergence emerge.