# Casper: Accelerating Stencil Computation
using Near-cache Processing

Alain Denzler1 Rahul Bera1 Nastaran Hajinazar1,2 Gagandeep Singh1
Geraldo F. Oliveira1 Juan Gómez-Luna1 Onur Mutlu1 1ETH Zürich 2Simon Fraser
University

###### Abstract.

Stencil computation is one of the most used kernels in a wide variety of
scientific applications, ranging from large-scale weather prediction to
solving partial differential equations. Stencil computations are characterized
by three unique properties: (1) low arithmetic intensity, (2) limited temporal
data reuse, and (3) regular and predictable data access pattern. As a result,
stencil computations are typically bandwidth-bound workloads, which only
experience limited benefits from the deep cache hierarchy of modern CPUs.

In this work, we propose Casper, a near-cache accelerator consisting of
specialized stencil compute units connected to the last-level cache (LLC) of a
traditional CPU. Casper is based on two key ideas: (1) avoiding the cost of
moving rarely reused data through the cache hierarchy, and (2) exploiting the
regularity of the data accesses and the inherent parallelism of the stencil
computation to increase the overall performance. With minimal changes in LLC
address decoding logic and data placement, Casper performs stencil
computations at the peak bandwidth of the LLC. We show that, by tightly
coupling lightweight stencil compute units near to LLC, Casper improves
performance of stencil kernels by $1.65\times$ on average, while reducing the
energy consumption by $35$% compared to a commercial high-performance multi-
core processor. Moreover, Casper provides a $37\times$ improvement in
performance-per-area compared to a state-of-the-art GPU.

##  1\. Introduction

A stencil operation (gysi2015modesto, ) defines a computation pattern where
elements in a multidimensional grid are updated based on the values of a fixed
pattern of neighboring points. Computations using stencil operations (called
stencil computations) are a key building block of important HPC applications
(colella2004defining, ) and are used in a wide range of scientific
simulations, including climate modeling (fuhrer2018near, ), seismic imaging
(mcmechan1983migration, ), fluid dynamics, and electromagnetic simulations
(taflove1988review, ). Stencil computations encompass a large percentage of
the overall runtime of such applications (maruyama2014optimizing, ;
fuhrer2013towards, ; olschanowsky2014study, ). For example, stencil
computation represents over 90% and 60% of the overall runtime in
computational fluid dynamics solver and COSMO climate simulation model,
respectively. However, the current compute-centric processing systems, such as
multi-core CPUs and GPUs, fail to fully utilize their on-chip resources (e.g.,
deep cache hierarchy, high throughput floating-point engines) when computing
stencil operations. This results in lower performance and decreased energy
efficiency for stencil computations in current systems.

Figure 1 shows the roofline plot (williams2009roofline, ) of important stencil
kernels on a server-class CPU. The horizontal line demonstrates the peak
floating-point performance of the system. The DRAM and L3 lines show the peak
memory and LLC bandwidth, respectively. We make two observations. First, the
stencil kernels _are not compute bound_ , having low arithmetic intensity
ranging from 0.09 to 0.2 FLOP/Byte. Second, all the kernels are _bounded by
memory resources_ , i.e., located on the left side of the inflection point and
below the memory lines. More importantly, all the kernels are located below
the L3 line and above the DRAM line, which shows that the stencils are bound
by the LLC bandwidth rather than the main memory. Given these observations, we
conclude that the high number of LLC accesses is the main bottleneck for
stencil computations, resulting in such computations to experience only
limited benefits from the deep cache hierarchies in modern CPU architectures.

Prior works (augustin2009optimized, ; datta2008stencil, ;
phillips2010implementing, ) also show that stencil kernels are bottlenecked by
memory due to their low arithmetic intensity, leading to under-utilization of
computational resources in compute-centric platforms (i.e., CPU, GPU). These
prior works demonstrate that stencil kernels can only leverage 21.8%
(augustin2009optimized, ) and 46% (datta2008stencil, ) of the computational
resources of a multi-core CPU and GPU, respectively, even in presence of code
optimizations (e.g., temporal blocking). Therefore, general-purpose processors
cannot deliver high performance and high efficiency (in terms of utilization)
for stencil computations, thus opening up the space for custom stencil-based
accelerators (sano2013multi, ; nero, ).

Figure 1. Roofline plot for a multi-core system running six stencil kernels.

Near-memory processing is a promising paradigm that aims to accelerate memory-
bandwidth-bound workloads with low arithmetic intensity. (seshadri2017ambit, ;
eckert2018neural, ; seshadri2017ambit, ; ahn2015pim, ; li2019pims, ;
Imani2019floatpim, ). The key idea of the NMC paradigm is to move computation
close to (i.e., near-memory) or even into (ı.e., in-memory) the memory devices
(i.e., caches, DRAM, storage) where the data resides, eliminating the need to
move the data to the processor and resulting in increased performance and
reduced energy consumption. Stencil computations are a prime candidate for
acceleration using NMC paradigm. In this work, we explore the opportunity to
improve the performance and energy efficiency of stencil computations in
traditional multi-core CPUs by computing near LLC.

Why near LLC? We exploit LLC as the prime location for computations, as
opposed to offloading the computation to the off-chip main-memory (li2019pims,
) or higher levels of caches (e.g., L2) for three main reasons. First, the
per-thread dataset for stencil kernels in widely-deployed workloads are
normally tiled to fit inside the LLC of traditional workstation-class CPUs
(cosmo_knl, ; fuhrer2018near, ). Hence, placing computation near LLC minimizes
unnecessary data transfers, which eventually provides both higher performance
and higher energy efficiency than computing inside core or main memory.
Second, the on-chip LLC bandwidth is multi-fold higher than a traditional DDR-
based DRAM main-memory bandwidth. Even compared to a high-bandwidth yet costly
HBM interface, the on-chip LLC bandwidth can be readily exploited in any
commodity off-the-shelf processor. Third, though computing in high-level
caches (e.g., L2) can provide higher per-compute bandwidth on paper, moving
data back and forth from LLC to higher-level caches reduces the effective
bandwidth due to cache-protocol-induced additional transactions (e.g., back
invalidations, write backs, etc.). Moreover, pulling data that provides low
reuse (common in stencil computation data) to higher-level caches results in
energy waste.

Thus our goal in this paper is to design a near-LLC accelerator that improves
the performance and energy efficiency of stencil computations by minimizing
the unnecessary data movement.

To this end, we propose Casper, a novel hardware/software codesign
specifically targeted at stencil computations. We minimize data movement by
placing a set of plSPU near the LLC of a traditional CPU architecture and
provide novel mechanisms to introduce data mapping changes and support
unaligned loads needed for high-performance stencil computations. Computation
is mapped to plSPU so that each SPU operates on the data that is located in
the closest LLC slice. This reduces the overall data access latency and energy
consumption while matching the compute performance to the peak bandwidth of
the LLC.

Placing a stencil accelerator next to the LLC of a CPU introduces two key
challenges. The first challenge is to maximize bandwidth utilization. A
superscalar out-of-order CPU leverages its front-end structures such as
reorder buffers and large instruction windows to exploit instruction-level
parallelism, thus enabling high memory bandwidth utilization. To address this
challenge, Casper introduces the notion of _stream_ to expose the memory level
parallelism that exists in the stencil computation to the SPUs. Each stream
represents a set of consecutive memory accesses with a fixed stride. Knowing
the stride for each stream, Casper can issue the memory accesses of the stream
in parallel. The notion of stream enables Casper to maximize memory bandwidth
utilization without requiring complex structures such as reorder buffers.

The second challenge is to minimize the data movement between different cache
slices. In stencil computations, the neighboring grid points need to be
accessed to compute the stencil operation for each grid point. However,
current systems employ an address mapping scheme that distributes data over
different LLC slices and provides load balance and fairness across CPU cores.
Such a mapping scheme can potentially map neighboring grid points to different
LLC slices, introducing data transfers over the NoC and thereby eliminating
the benefits of near-cache computing. To address this challenge, Casper
includes hardware support to align the data mapping to the grid structure of
the stencil computation and place neighboring grid points into the same slice
of the LLC.

We evaluate Casper using six widely used stencil kernels with up-to
3-dimensional grid domains. Casper outperforms a commercial multi-core CPU on
average by (up-to) $1.65\times$ ($4.16\times$) and reduces the energy
consumption by $35\%$. Compared to a state-of-the-art GPU, Casper improves
performance-per-area by $37\times$ on average.

We make the following contributions:

  * •

We propose Casper, a high-performance near-cache accelerator for stencil
computations. Casper addresses the memory bottleneck in stencil computations
by (1) minimizes the movement of data within the cache hierarchy, and (2)
maximizing the utilization of the LLC bandwidth.

  * •

We provide hardware support to change the mapping of cache lines to LLC slices
to improve spatial data locality for stencil computations.

  * •

We design an _efficient_ mechanism to support unaligned loads directly from
the LLC.

  * •

We evaluate the effectiveness of Casper on six widely used stencil kernels and
demonstrate that Casper improves performance and reduces the energy
consumption by $1.65\times$ and $35$% on average compared to a commercial
multi-core CPU.

##  2\. Background: Stencil Computation

Stencil computations (gysi2015modesto, ) update points in a data grid based on
a fixed pattern of neighboring points. This fixed pattern, called the stencil,
is applied on the complete grid iteratively until either a convergence
criterion or a fixed number of time steps are reached. Stencil computations
are widespread in scientific computing, and are considered one of the 13
dwarfs of scientific computing (colella2004defining, ).

Stencil computations exhibit several common properties. We explain these
properties using the Jacobi stencil as an example (demmel1997applied, ) which
is commonly used to solve discretized partial differential equations. Figure 2
shows the source code and data access pattern of a 2-dimensional Jacobi
stencil. The computation performs arithmetic mean of each point in the grid
and its immediate neighbors in all directions. This implementation uses three
nested loops where the outermost loop iterates over time steps and the two
inner loops sweep over the complete 2D grid.

Figure 2. 2D Jacobi stencil pseudo-code and data access pattern.

From this example, we can identify four common properties of stencil
computations. First, the computation is embarrassingly parallel because the
read and write data sets are disjoint. Second, the computation is regular and
statically analyzable. The data dependencies and the structure of the
computation can be analyzed ahead-of-time and do not depend on dynamic input.
Third, the arithmetic intensity of the stencil is low. Fourth, only a few
types of operations are needed to compute the stencil. For example, in case of
Jacobi stencil, only a floating-point MAC operation is performed for each
input grid point. These properties enable the opportunity for high-performance
stencil acceleration since they provide well-balanced parallelism and
predictable and uniform access patterns. Additionally, being memory-bound
(Figure 1) makes stencils good candidates for near-memory acceleration.

In this work, we focus on stencil computations that apply _one_ stencil to the
complete dataset per time step. Related work on stencil computations, such as
(nero, ; gysi2015modesto, ) focuses on compound stencils, which apply
_several_ stencils using multiple input and intermediate arrays during each
time step. These compound stencils are limited by DRAM bandwidth (nero, ).
While both types of stencils are common in real-world workloads
(gysi2015modesto, ; christen2012patus, ), we focus on the first type which
only applies one stencil per iteration.

##  3\. Casper: Overall Architecture

###  3.1. Overview

The key idea of Casper is to place specialized stencil processing units (SPU)
near the last-level cache (LLC) of a traditional CPU. As a typical LLC is
partitioned into multiple cache slices, Casper places one SPU in each LLC
slice. At a high level, each SPU performs the stencil operation on each grid
point sequentially. The input grid points are accessed in row-major order,
following their layout in memory. For each input grid point, the SPU loads the
data from LLC either from the local or remote LLC slice. Once the data arrives
at the SPU, the SPU’s execution unit performs a single multiply operation with
a predefined constant, whose output is added to the accumulator. The final
accumulated result is stored back into the output grid point only when all
computations for the output grid point have been completed.

To accelerate the accesses to input grid points, Casper abstracts a series of
data accesses using a _stream_. A _stream_ is a sequence of data elements of
the same type located in consecutive physical memory addresses (jouppi_stream,
). A stream is represented by the start address, the data type width, number
of elements, and a position pointer. Initially, the position pointer of a
stream points to the start address. Upon receiving a control signal, the
position pointer is advanced till the final element in the stream. The stream
abstraction facilitates quick iteration over the data instead of loading each
data element using absolute memory addresses. Since a stencil computation
moves through a grid at the same pace maintaining the same relative offsets
for the stencil pattern, SPU uses stream abstraction to accelerate input grid
access. This allows the same stencil code to be reused for all the grid
points, and multiple disjoint data sources to be easily included in the
stencil pattern by referring to different streams.

###  3.2. Stencil Processing Unit

Figure 3 shows the architecture of a SPU. The main building blocks are (1) the
instruction buffer, (2) the load queue, (3) the stream buffer, (4) the
constant buffer, and (5) the execution unit. The complete SPU is pipelined to
maintain single-cycle instruction throughput. The SPU controls the complete
computation, from instruction fetch until retire. Therefore it does not
involve any interaction with the CPU.

Instruction Buffer. The instruction buffer has a capacity to hold $64$
_compressed_ instructions to compute a stencil. Every instruction encodes the
operands and the control signals for one type of stencil operation. The same
sequence of instructions is applied to every stencil output grid point due to
the regular nature of the stencil computation. We introduce the instruction
set in Section 5.

Load Queue. The load queue is responsible for issuing data requests to the
memory subsystem. While the SPU processes instructions in-order, the varying
memory access latency can result in responses that arrive out-of-order. The
load queue acts as a buffer to hold the data until all the previous, longer-
latency memory requests have been completed. This buffer space ensures that
the instruction order is maintained.

Stream and Constant Buffers. The stream buffer holds the current state of the
streams. For every decoded instruction, this state is loaded from the buffer,
and the effective address is calculated to issue the memory requests to LLC.
The constant buffer holds the constant factors needed for the MAC operation.
The stream and constant buffers are initialized by API calls before the SPU
starts the stencil computation.

Execution Unit. As shown in Figure 3, the execution unit of a SPU is comprised
of a $512$-bit vector unit which operates on vectors of $8$ double-precision
floating-point elements. For every input grid point, the SPU loads the data
elements from the neighboring points, multiplies them with a constant factor,
and accumulates the results. This final accumulated result is then stored in
the output grid point. Therefore, each SPU execution unit only computes one
kind of instruction: a double-precision floating-point multiply-accumulate
operation (MAC).

Figure 3. Components of one SPU.

##  4\. Casper: Micro-architectural Support

###  4.1. Support to Unaligned Loads

While the stream abstraction offers efficient memory accesses, the LLC
architecture only supports data accesses aligned to 64-byte boundaries. As the
relative offsets used in streams are not necessarily aligned to 64-byte
cacheline boundary, the loaded data might need to be realigned before it can
be used by MAC compute. For example, in Figure 4 each SPU computes $8$
$64$-bit output grid points (B[i]) for a $3$-point Jacobi-1D stencil. The
access to the center point of the stencil A[i] is correctly aligned to the
$64$-byte boundary (shaded in yellow) such that the data for this point can be
used for computation as soon as it arrives at the SPU. However, to gather
several input grid points at indices $+3$ (A[i+3]) and $-3$ (A[i-3]), the data
coming from the cache needs to be correctly shifted (shaded in orange and
red), and two cache lines need to be combined to assemble the operands for the
computation. As a result, preparing the operand for the MAC unit involves two
loads, one shift, and one combine operation on the data from the two
cachelines. These additional operations due to unaligned data accesses not
only drops the utilization of MAC units (e.g., the MAC utilization in Figure 4
drops from $4$ load/store operations per $3$ MACs to $5$ load/store operations
per $3$ MACs), but also wastes precious cache bandwidth to load the data. To
address the above challenge, we introduce two simple modifications to the LLC
row decoding logic: (1) support to loading data aligned to any 8-byte
boundary, and (2) support to correctly align data for SPU execution. Our
modified LLC row decoding logic reduces SPU’s complexity and area footprint
since by only modifying the row decoding circuitry, we avoid the need for
additional logic to shift and pack the partial cache lines, and/or a large
register file inside SPU.

Figure 4. Unaligned loads occurring during the vectorized execution of a
stencil.

Implementation Challenges. Loading data aligned to 8B boundaries potentially
implies loading data located on two cache lines. This introduces two key
challenges. First, two tags need to be matched to locate the two cache lines.
Second, the correct data from both cache lines must be loaded using only one
supplied address.

To address the first challenge, we enable the cache to match the two tags in
parallel by adding a second read port to the tag array. Since the two cache
lines involved in an unaligned load are always at consecutive addresses, they
are guaranteed to be mapped to different cache sets, and thus, there are no
conflicts during the tag matching process. If one or both cache lines are not
in a readable state, a regular cache miss is generated, and the request is
stalled until the miss is resolved.

To address the second challenge, we make changes into the LLC row decoding
logic to load data either from the requested address or one of the cache lines
located inside the 64B-vicinity of the requested address. More specifically,
we add one 3-to-1 multiplexer for each LLC SRAM row. The inputs to the
multiplexer are set to the row decoder output of the current row and both the
adjacent rows. The multiplexer selects the appropriate row(s) based on the row
selection signal.

Loading Shifted Cache Lines. Figure 5 shows the execution sequence of
unaligned access to the LLC. In this example, we consider the request to cache
line holding elements 8 to 15, shifted to the right by three positions. This
corresponds to the access to input grid points A[i-3] in Figure 4. Let us
first consider the case where the two cache lines involved in the access are
mapped to the _same_ cache way. One cache way consists of 4$\times$32kB data
arrays, each having 2$\times$16kB subarrays. One 16kB subarray holds 64
consecutive bits of a cache line for all the 2048 sets. Thus, each 64-bit
segment of a cache line is stored in a different SRAM subarray.

As shown in Figure 5, all the elements required to build the 64B block of data
are stored in different subarrays of a single cache way (orange and green
elements on the left part of the figure). Therefore, by selecting the correct
data at the subarray level (using shift direction and amount) it is possible
to load all the elements using only one command, while maintaining the same
throughput and latency as a regularly aligned load. To complete the unaligned
access, we need to rotate the data such that the correct element is at the
first position of the 64B block of response data. To accomplish this, we use a
rotate network that performs this operation before the final output. In this
way, the SPU can load data from the LLC aligned to arbitrary 8B boundaries.

If the two cache lines involved in an unaligned load are mapped to _different_
ways, the load sequence is similar. As the data and tag access happen in
parallel inside the cache, the data load is initiated on all cache ways before
the way hit has been confirmed. Thus for any access, all the ways present the
data for the requested set on their output bus. Based on the way hit, the data
from the correct way is selected for the output. We modify the way hit
selection to select each subarray’s output independently, depending on the
shift amount and direction supplied with the request.

Figure 5. Loading unaligned data from the LLC. The two cache lines involved in
the access are mapped to the same cache way. One cache way consists of
4$\times$32kB data arrays, each having 2$\times$16kB subarrays.

###  4.2. LLC Data Mapping

Computing at the LLC’s peak throughput is only possible if all the plSPU in
the system load data from their closest LLC slice, thus avoiding the network-
on-chip (NoC) traffic overhead. However, the existing mapping of memory
addresses to the LLC slices does not guarantee mapping consecutive cache lines
to the same slice to provide better fairness and load balancing
(yarom2015mapping, ). Due to the streaming nature of the stencil computation
and the dependency on a small neighborhood of grid points, scattering data
through the cache slices would increase the number of remote data requests per
SPU, penalizing performance and energy consumption. Thus Casper maps blocks of
consecutive memory addresses to the same LLC slice. To achieve this, we
introduce hardware support for stencil segments (basu2013efficient, ). A
stencil segment is a physically contiguous block of memory that is exposed by
the system to hold stencil data. Casper identifies a stencil segment and
changes the segment mapping to LLC slices such that physically contiguous
blocks of data are mapped to the same slice. At each NoC injection point, the
system checks if the address is part of the stencil segment, and, if valid,
applies a new hash function (i.e., a linear hash that statically maps
contiguous blocks of 128kB to LLC slices in round-robin fashion) to route a
memory request. Otherwise, the existing hash is applied. By deciding which
mapping function to use based on a memory request’s physical address, each
address is mapped to precisely one cache slice, regardless of whether it
contains stencil data or not.

###  4.3. Coherence Support

A SPU loads data from the LLC by directly injecting a load/store request in
the LLC controller’s request queues. Casper does not impact the current cache
coherency mechanism, because requests from the SPU are injected into the same
request queues as conventional requests from the CPU cores and the private
caches. If a write request from an SPU targets a cache line that is not in
writable state in the LLC, the coherency mechanism will trigger necessary
state transitions and invalidations to allow the write to complete.

###  4.4. Concurrent Execution with CPU and Context Switching

To enable concurrent execution with the CPU, we reserve one way of the LLC for
other applications running in the CPU, similar to prior work
(eckert2018neural, ). Additionally, high priority is assigned to the SPU
process to minimize the occurrence of context switches. If a context switch
happens, the state of the SPU remains unchanged, because it is not allowed to
start a new stencil computation before the current computation finishes.

##  5\. Casper: Programming Model Support

###  5.1. Casper ISA

Casper makes use of specialized instructions to compute each point involved in
the stencil operation. Every Casper instruction is $15$-bits wide and
comprises of (1) $4$b constant buffer index, (2) $4$b stream buffer index, (3)
$1$b shift direction, (4) $3$b shift amount, and (5) $3$b control bits. Note
that Casper reuses the instructions for all the grid points, thus reducing
instruction count.

After decoding a new instruction, the SPU uses constant index to load the
double-precision constant value from the constant buffer. The stream index is
used to point the stream buffer to gather memory address of the requested
stream. The shift direction and the shift amount are used to initiate
appropriate load request to the LLC. Three control bits, _clear accumulator_ ,
_enable output_ , and _advance stream_ are used to reset the accumulators in
the execution unit, enable the contents of the accumulator to be stored into
the results array, and to advance the stream pointers, respectively. We also
provide a programming library that allows the user to easily generate Casper
instructions from user-level code.

###  5.2. Casper API

Table 1 shows Casper API functions. Calls to API functions are mapped directly
to ISA instructions, which are broadcasted to all SPUs. These instructions are
integrated into the existing ISA using spare instruction opcodes.

Table 1. Casper Programmer API Function |  Parameters |  Description
---|---|---
initStencilSegment |  int size |  Requests a physically contiguous memory region of the specified size from the system to hold the stencil data
initStencilcode |  addr A, int length |  Takes a pointer to the microcode and the length of the code. After generating the code with helper functions provided by the programming library, the code is then broadcast to the plSPU
initConstant |  double const, int index |  Iinitializes constant values that will be used during the multiplication step of a stencil operation. The function sets the specified constant value at the given index in the constant buffer
initStream |  addr A, int streamID, int accID |  Configures the streams used in the stencil code. The streams are configured per SPU to enable the programmer to tune the data layout to the grid’s structure and minimize the number of off-slice cache accesses
setNElements |  int n, int accID |  To communicates to each SPU how many elements to compute. All the SPUs maintain a counter to keep track of their progress and stop when the desired number of elements has been computed
startAccelerator |  - |  The final function starts with the execution of all the plSPU

##  6\. An Example: Jacobi-2D Stencil

We illustrate how to program Casper using the code in Figure 6 as an example.
The code implements the Jacobi-2D stencil presented in Section 2. We explain
the example using a system consisting of 4 LLC slices and plSPU.

Figure 6. Program code for Jacobi 2D stencil.

First, a stencil segment covering 4MB is allocated (line 4). Then, we define
the start of the arrays A and B such that the same grid point of both arrays
is mapped to the same LLC slice (lines 6-8). The programmer then initializes
the arrays with the stencil data. Furthermore, the constant for the
multiplication (line 12) and the stencil instructions (line 14) are sent to
the plSPU. The streams for all the plSPU are configured inside the loop (lines
22-29). In this example, four streams are configured: three input streams to
load the elements at A[j-1][i], A[j][i] and A[j+1][i] (lines 23-25), and one
output stream to store the result (line 26). As the elements A[j][i-1],
A[j][i], and A[j][i+1] are laid out in consecutive memory addresses, we can
reuse the same stream and leverage the support for unaligned loads by shifting
the access to the left (right) by one element to load the data. Finally, the
number of elements to compute for each SPU (line 28) is configured, and the
computation starts (line 30).

Figure 7 shows the stencil instructions executed on the SPU. As the stencil
loads data from five input points, the instruction sequence consists of five
instructions. All the input points are multiplied by the constant $0.2$. This
means all the instructions encode the same constant factor c0. The
instructions load data from three different streams: the first stream points
to the value A[j-1][i], aligned to a cache line boundary. The values
A[j][i-1], A[j][i], and A[j][i+1] are stored in consecutive memory addresses.
These three accesses all use the same stream, but using a shift amount and
direction included in the request for the elements at indices i+1 and i-1.
This loads the correctly aligned data using the unaligned load mechanism.
Finally, the third stream is configured to load the value at A[j+1][i].

Figure 7. Instruction sequence for the Jacobi-2D stencil.

The final 3 bits of each instruction hold control information for the SPU. For
example, the first bit clears the accumulator before the accumulation step.
This bit must be set by the first instruction of a grid point (line 2). The
next bit enables the output. It is set in the last instruction of the sequence
(line 10) and generates a store request. The final bit signals the SPU to
advance the stream. It has to be set in the last instruction consuming data
from each stream (lines 2, 8, 10).

##  7\. Methodology

###  7.1. Experimental Setup

We simulate the performance of Casper using the gem5 simulator in syscall
emulation mode (binkert2011gem5, ). Table 2 describes our system configuration
in details. We use the ARM ISA and the Ruby memory model. The system is based
on a generic modern server-class CPU, consisting of 16 out-of-order cores and
three levels of cache, having a sliced LLC with 2MB per slice. Our baseline
configuration uses the same system configuration without the plSPU and the LLC
changes we propose. Also, we include stride prefetchers at all levels of the
cache hierarchy. We evaluate the performance and energy benefit of Casper
against the baseline CPU architecture, an NVIDIA Titan V GPU(titanV, ). We use
an energy model based on CACTI 7.0 (balasubramonian2017cacti, ) and energy
models proposed by prior works (tsai2018adaptive, ; tsai2017jenga, ). We use
the area model presented in (shao2014aladdin, ) scaled down to 22nm to
estimate the area of the SPU. For the GPU performance/area comparisons, we use
the complete die size of 815mm² of the Titan V (nvidiaTitanDie, ) because
typical GPU-accelerated systems also need a host CPU to perform the
computation. As a result, the total area for the end-user is either the
complete GPU or the Casper hardware modifications, added to the area of the
existing host CPU.

Table 2. Simulation Parameters

Component | Configuration
---|---
Casper | 16 plSPU, 1 SIMD unit/SPU (512bits wide)
| 10 entry load queue, 0.016 nJ/instruction
Core | 16 out-of-order cores, 2 GHz, 8-wide issue,
| 72 entry load queue, 64 entry store queue,
| 1 SIMD unit/core (512bits wide)
| 224 entry ROB, 0.08nJ per instruction
L1I/D Cache | 32kB, private, 8-way, 16 MSHRs,
| 4 cycle round-trip latency, 2 load ports, 1 store port
| 15/33 pJ per hit/miss (tsai2018adaptive, )
L2 Cache | 256kB, private, 8-way, 16 MSHRs,
| 12 cycle round-trip latency, 1 load port, 1 store port
| 46/93 pJ per hit/miss (tsai2018adaptive, )
L3 Cache | 32MB, shared, 16-way, 16 slices, 32 MSHRs/slice
| 36 cycle round-trip latency, 1 load/store port per slice
| 945/1904 pJ per hit/miss (tsai2018adaptive, )
Coherence Protocol | MESI
Replacement Policy | LRU replacement
Hardware Prefetchers | Stride prefetchers at all levels of the cache
On-Chip Network | mesh, XY-routing, 64B/cycle per direction
Main Memory | 16GB, DDR4, 4 channels, 160nJ per read/write (tsai2017jenga, )

###  7.2. Benchmarks

We evaluate Casper on six stencil benchmarks, including up to 3-dimensional
stencils with varying data reuse characteristics, ranging from 3 input points
(Jacobi 1D) up to 33 points for the 33-point 3D stencil. We use Jacobi 1D,
-2D, 7-point 3D (heat diffusion) from Polybench (pouchet2012polybench, ), a
5$\times$5 Gaussian blur filter (canny1986computational, ), the 7-point 1D
kernel from (holewinski2012high, ), and a 33-point 3D kernel to represent
higher-order scientific simulations(datta2008stencil, ; datta2009auto, ). All
the benchmarks use Jacobi-style stencils with disjoint read and write data
sets. The benchmarks are elementary stencils that can be conjugated to form
more complex stencils occurring in real-world applications. In addition to the
standard data set sizes (which fit inside the LLC of the CPU), we also
evaluate Casper on data sets that 1) exceed the size of the LLC (named DRAM),
and 2) fit within the private L2 caches of the CPU (L2). Table 3 lists the
domain sizes. We report the runtime for one timestep (T = 1 in Figure 2),
measured at the third timestep of the benchmarks. Iteration runtimes are
stable, starting from the third timestep in our benchmarks. Because of $T$ is
usually large enough to amortize the cost of moving the data from memory, we
do not consider this in the evaluations.

Table 3. Domain size used for evaluations Level |  1D |  2D |  3D
---|---|---|---
L2 |  131,072 |  $512\times 256$ |  $64\times 64\times 32$
L3 |  1,048,576 |  $1024\times 1024$ |  $128\times 128\times 64$
DRAM |  4,194,304 |  $2048\times 2048$ |  $256\times 256\times 64$

##  8\. Results

###  8.1. Performance

Figure 8 shows the speedup of Casper compared to a 16-core baseline system.
For the datasets that fit within the LLC, which represent typical data set
sizes for stencil computations, we observe an average speedup of 1.65$\times$.
The 1- and 2-dimensional stencils achieve speedups between 1.66$\times$ and
3.0$\times$. However, the 3-dimensional stencils cannot achieve the same
performance and even experience a slowdown in the 33-point stencil case. The
reason for this performance is twofold. First, the 3-dimensional stencils need
to load a significant part of their input data from remote LLC slices, which
introduces longer access latencies and lowers the throughput of the plSPU.
Second, the 33-point 3D stencil has good L1 cache behavior in the baseline,
achieving a hit-rate of 95%. Thus, this stencil is well-suited for execution
on a traditional CPU with regards to performance. We conclude that the
performance benefits of Casper are larger on lower-dimensional stencils that
load most of their input data from the local LLC slice.

The average performance improvement for the smaller data set that fits within
the L2 caches of the CPU is 1.89$\times$. Even though the data is stored
closer to the core and does not need to travel through the entire cache
hierarchy, the speedups are similar to the larger data sets that fit within
the LLC. This is the case because the access latency from the CPU to the L2
cache is similar to the latency between SPU and the closest LLC slice (12 vs 8
cycles load-to-use). The 3-dimensional benchmarks even reach better
performance improvements than on larger datasets. This is the case because the
baseline performance is limited by cache lines ping-ponging between private
caches via the LLC. Casper avoids this problem by computing directly on the
shared LLC.

Our results show that for large data sets that exceed the size of the LLC,
Casper improves performance by 1.4$\times$, on average. The highest speedups
are achieved by the 2-dimensional stencils, with Blur 2D reaching
4.16$\times$. We explain this by the fact that the baseline has a very low LLC
hit-rate (only 2%), and thus the number of memory accesses is 4$\times$ higher
than in Casper. The main reason for this low hit-rate is that the prefetchers
are interfering with demand accesses, evicting cache lines out of the LLC
before they are used. The remaining stencils perform similar to the baseline
or even experience slowdowns. The reason for this is the fact that memory
bandwidth is the main bottleneck. The 33-point 3D stencil again experiences a
slowdown because it is well-suited for computation with smaller private caches
for each core. We conclude that even though Casper cannot alleviate the memory
bandwidth bottleneck for large data sets, it does not lead to significant
slowdowns for most stencils.

Figure 8. Speedup compared to the baseline multi-core system.

###  8.2. Energy Consumption

Figure 9 shows the normalized energy consumption of Casper compared to the
baseline CPU. For the data sets that fit into the LLC, Casper reduces energy
consumption by 55%, on average. Casper reduces energy consumption by 49% even
for the 33-point 3D stencil, whose performance is slower in Casper than in the
baseline. The reduction in energy consumption is higher for simpler stencils
(Jacobi 1D/2D, 7-point 3D), reaching up to as 65% for 7-point 3D. The reason
for this is that more complex stencils perform more LLC accesses and since the
plSPU are situated close to the LLC slice, accessing the LLC is not as energy-
efficient as accessing the smaller L1 cache. This results in lower energy
savings for the more complex stencils because of higher L1 reuse in baseline
CPUs.

For the smaller and larger data sets, Casper reduces energy consumption by to
26% and 23%, on average respectively. We make two observations: First, Casper
increases the energy consumption of the 1-dimensional benchmarks (Jacobi 1D
and 7-point 1D) when compared to the baseline, for both small and large data
sets. For the large data sets, this is the case because the CPU cores can be
idle for most of the runtime, waiting for memory. For the smaller data set,
the baseline’s energy consumption is very low because there are very few LLC
accesses. Since Casper loads the data from the shared LLC, it increases energy
consumption in such cases. Second, for all the other benchmarks, Casper
reduces energy consumption significantly. Casper is more energy-efficient even
for the benchmarks that the CPU baseline outperforms since our SPU design is
more energy-efficient than the CPU baseline. Thus, we conclude that Casper
achieves significant reductions in energy consumption compared to an out-of-
order CPU.

Figure 9. Normalized energy consumption compared to the 16-core baseline.

###  8.3. Comparison with GPU

Figure 10 shows Casper’s performance and performance-per-area normalized to an
NVIDIA Titan V GPU. GPU outperforms Casper by $3.71\times$, $2.89\times$, and
$36.64\times$ on average across all stencils that fit inside L2, LLC, and DRAM
respectively. However, in all stencil kernels Casper provides higher
performance-per-area, up to $190\times$ compared to GPU. We observe that the
L2- and LLC-sized data sets achieve performance-per-area improvements of
47$\times$ and 60$\times$, respectively. For these data set sizes, Casper has
the advantage of its tight integration into the large LLC. At the same time,
the data does not fit into the GPU caches. For the large DRAM-sized data sets,
however, the average improvement in performance/area is only 4.78$\times$. In
this case, GPU improves relative performance because of its higher memory
bandwidth.

Figure 10. Performance/area compared to an NVIDIA Titan V GPU for three
different domain sizes.

###  8.4. Effect of Individual Optimizations

The SPUs in Casper take advantage of two key optimizations: 1) a custom _data
mapping_ in LLC, which creates stencil segments to maximize the accesses of an
SPU to the own LLC slice, and 2) the _near-cache_ (near-LLC) location of the
SPUs, which minimizes data access latency and leverages the peak bandwidth of
the LLC. In this section, we evaluate the contribution of each optimization to
the overall performance. The baseline for this analysis is a system where the
SPUs are located next to the private L1 caches of CPU cores. The baseline LLC
data mapping places consecutive cache lines in consecutive LLC slices (similar
to (yarom2015mapping, )). First, we apply the data mapping optimization and
compare the performance to the baseline. Second, we apply both the data
mapping optimization and the near-cache optimization.

We make two observations from the results in Figure 11. First, computing near-
cache (green portion of the bars) is the major contributor to the speedup.
Second, the custom data mapping (blue portion of the bars) produces up to 30%
of the speedup (Jacobi 1D with LLC data set), but its effect is negligible (or
even negative) in several cases (1D and 2D benchmarks with L2/DRAM data sets,
7-point 3D with LLC data set). In these cases, the custom data mapping results
in a number of accesses to remote LLC slices which is similar to the baseline
data mapping.

Figure 11. Contribution of custom data mapping and near-cache SPU location to
the speedup of Casper over a baseline system with SPUs near L1 for
L2-/LLC-/DRAM-sized data sets.

###  8.5. Hardware Cost

Stencil Processing Unit. The area of one SPU scaled to 22nm technology is
0.146 mm2. The most significant contributors to this area are the execution
unit and the SRAM array used to buffer complete memory requests.

Unaligned Loads. Our hardware mechanism to support unaligned loads consumes an
additional 0.14 mm2 compared to the baseline 2MB SRAM cache slice. This
amounts to a 5% increase in area per LLC slice. Almost the complete area
overhead is attributed to the second read port of the tag array, which
consumes 0.12 mm2 of space. The remaining hardware overhead of one 3:1
multiplexer per SRAM row, the barrel shifter to rotate the final output, and
the split multiplexers for way selection are minimal compared to the tag array
overhead.

Address to LLC slice mapping. Identifying the stencil segment needs two
registers to store the start and the length of the segment. The comparison
needs one adder and one comparator. The new mapping is a simple bit-select
from the physical address, and thus requires minimal additions. This hardware
is introduced at every NoC injection point.

In summary, the hardware additions proposed in this paper require an
additional 4.65 mm2 of die area for a system using 16 SPUs. This amounts to a
0.77% area increase compared to the Marvell ThunderX 2 CPU (thunderx2, ), a
server-class ARM CPU implemented in 16nm hosting 32MB of LLC.

##  9\. Related Work

To the best of our knowledge, we present the first work that tightly
integrates specialized compute units into the last-level cache of a CPU to
perform stencil computations. In this section, we succinctly compare prior
proposals against Casper.

A wide body of research has focused on studying and analyzing stencil
computations (li2019pims, ; gysi2015modesto, ; ragan2013halide, ; nguyen20103,
; stengel2015quantifying, ; fuhrer2018near, ; armejach2018stencil, ). PIMS
(li2019pims, ) is the most closely-related work to Casper. PIMS exploits high-
bandwidth provided by 3D-stacking (e.g., HMC, HBM) to accelerate stencils.
Casper, being a near-LLC accelerator, can be integrated with any commodity
processor without the need of costly interfacing using through-silicon vias.
Szustak et al. (szustak2013using, ) accelerate the MPDATA stencil kernel on
multi-core CPU while Thaler et al. (cosmo_knl, ) port weather stencil kernles
to a many-core system. GPUs (anjum2019efficient, ; fuhrer2018near, ) have been
shown to increase performance due to the high degree of parallelism present in
the computation. Wahib et al. (wahib2014scalable, ) develop an analytical
performance model for choosing an optimal GPU-based execution strategy for
various scientific stencil kernles. Gysi et al. (gysi2015modesto, ) provide
guidelines for optimizing stencil kernels for CPU–GPU systems.

More recently, the use of FPGAs to accelerate stencils has been proposed
(van2019coherently, ; singh2019narmada, ; chi2018soda, ; nero, ). Augmenting
general-purpose cores with specialized FPGA accelerators is a promising
approach to enhance overall system performance. However, data still needs to
be moved to these off-chip external accelerators. Moreover, taking full
advantage of FPGAs for accelerating a workload is not a trivial task as it
requires sufficient FPGA programming skills to map the workload and optimize
the design for the FPGA microarchitecture. In contrast, processing close to
memory tightly integrates compute units close to the memory. This integration
avoids unnecessary data-movement.

##  10\. Conclusion

We present Casper, a novel near-cache acceleration system for stencil
computations. Our throughput-optimized stencil processing units and their
placement near the LLC enable a high performance and energy-efficient
execution for stencil computations. The design performs stencil computations
close to where data resides and orchestrates data accesses to minimize data
movements. With the addition of hardware support to map consecutive blocks of
data to the same LLC slice and by exploiting the memory level parallelism, we
reach our goal to match the system’s compute potential to the bandwidth
offered by the LLC with an area overhead of less than $1\%$. Our evaluations
show that Casper improves performance by 65% on typical data sets used for
stencil computations. Compared to a state-of-the-art GPU, performance per area
improves on average by 37$\times$. Across all our benchmarks, Casper reduces
energy consumption by 35% compared to the CPU baseline.

## References

  * (1) J. Ahn, S. Yoo, O. Mutlu, and K. Choi, “PIM-enabled instructions: a low-overhead, locality-aware processing-in-memory architecture,” in _ISCA_ , 2015.
  * (2) O. Anjum, S. Garcia de Gonzalo, M. Hidayetoglu, and W.-M. Hwu, “An Efficient GPU Implementation Technique for Higher-Order 3D Stencils,” in _HPCC_ , 2019\.
  * (3) A. Armejach, H. Caminal, J. M. Cebrian, R. González-Alberquilla, C. Adeniyi-Jones, M. Valero, M. Casas, and M. Moretó, “Stencil codes on a vector length agnostic architecture,” in _PACT_ , 2018.
  * (4) W. Augustin, V. Heuveline, and J.-P. Weiss, “Optimized stencil computation using in-place calculation on modern multicore systems,” in _ECPP_ , 2009\.
  * (5) R. Balasubramonian, A. B. Kahng, N. Muralimanohar, A. Shafiee, and V. Srinivas, “CACTI 7: New tools for interconnect exploration in innovative off-chip memories,” _TACO_ , 2017.
  * (6) A. Basu, J. Gandhi, J. Chang, M. D. Hill, and M. M. Swift, “Efficient virtual memory for big memory servers,” in _ISCA_ , 2013.
  * (7) N. Binkert, B. Beckmann, G. Black, S. K. Reinhardt, A. Saidi, A. Basu, J. Hestness, D. R. Hower, T. Krishna, S. Sardashti _et al._ , “The gem5 simulator,” _ACM SIGARCH Computer Architecture News_ , 2011.
  * (8) J. Canny, “A computational approach to edge detection,” _IEEE Transactions on pattern analysis and machine intelligence_ , 1986.
  * (9) Y. Chi, J. Cong, P. Wei, and P. Zhou, “SODA: Stencil with Optimized Dataflow Architecture,” in _ICCAD_ , 2018.
  * (10) M. Christen, O. Schenk, and Y. Cui, “Patus for convenient high-performance stencils: evaluation in earthquake simulations,” in _SC_ , 2012.
  * (11) P. Colella, “Defining software requirements for scientific computing,” 2004.
  * (12) K. Datta, M. Murphy, V. Volkov, S. Williams, J. Carter, L. Oliker, D. Patterson, J. Shalf, and K. Yelick, “Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures,” in _SC_ , 2008\.
  * (13) K. Datta and K. A. Yelick, _Auto-tuning stencil codes for cache-based multicore platforms_. University of California, Berkeley, 2009.
  * (14) J. W. Demmel, _Applied numerical linear algebra_. SIAM, 1997.
  * (15) C. Eckert, X. Wang, J. Wang, A. Subramaniyan, R. Iyer, D. Sylvester, D. Blaauw, and R. Das, “Neural cache: Bit-serial in-cache acceleration of deep neural networks,” in _ISCA_ , 2018.
  * (16) O. Fuhrer, T. Chadha, T. Hoefler, G. Kwasniewski, X. Lapillonne, D. Leutwyler, D. Lüthi, C. Osuna, C. Schär, T. C. Schulthess _et al._ , “Near-global climate simulation at 1 km resolution: establishing a performance baseline on 4888 GPUs with COSMO 5.0,” _Geoscientific Model Development_ , 2018.
  * (17) O. Fuhrer, C. Osuna, X. Lapillonne, T. Gysi, M. Bianco, and T. Schulthess, “Towards gpu-accelerated operational weather forecasting,” in _The GPU Technology Conference_ , 2013.
  * (18) T. Gysi, T. Grosser, and T. Hoefler, “Modesto: Data-centric analytic optimization of complex stencil programs on heterogeneous architectures,” in _SC_ , 2015.
  * (19) J. Holewinski, L.-N. Pouchet, and P. Sadayappan, “High-performance code generation for stencil computations on GPU architectures,” in _SC_ , 2012\.
  * (20) M. Imani, S. Gupta, Y. Kim, and T. Rosing, “FloatPIM: In-memory Acceleration of Deep Neural Network Training with High Precision,” in _ISCA_ , 2019.
  * (21) N. P. Jouppi, “Improving Direct-Mapped Cache Performance by the Addition of a Small Fully-Associative Cache and Prefetch Buffers,” in _ISCA_ , 1990.
  * (22) J. Li, X. Wang, A. Tumeo, B. Williams, J. D. Leidel, and Y. Chen, “PIMS: a lightweight processing-in-memory accelerator for stencil computations,” in _ISMS_ , 2019.
  * (23) N. Maruyama and T. Aoki, “Optimizing stencil computations for nvidia kepler gpus,” in _Proceedings of the 1st international workshop on high-performance stencil computations, Vienna_ , 2014, pp. 89–95.
  * (24) Marvell, “ThunderX2 CPU.” [Online]. Available: https://en.wikichip.org/wiki/cavium/thunderx2
  * (25) G. A. McMechan, “Migration by extrapolation of time-dependent boundary values,” _Geophysical Prospecting_ , vol. 31, no. 3, pp. 413–420, 1983.
  * (26) A. Nguyen, N. Satish, J. Chhugani, C. Kim, and P. Dubey, “3.5-D blocking optimization for stencil computations on modern CPUs and GPUs,” in _SC_ , 2010.
  * (27) Nvidia, “NVIDIA TITAN V.” [Online]. Available: https://www.nvidia.com/en-us/titan/titan-v/
  * (28) Nvidia, “The nvidia titan v preview,” Dec 2017. [Online]. Available: https://www.anandtech.com/show/12170/nvidia-titan-v-preview-titanomachy
  * (29) C. Olschanowsky, M. M. Strout, S. Guzik, J. Loffeld, and J. Hittinger, “A study on balancing parallelism, data locality, and recomputation in existing PDE solvers,” in _SC_ , 2014.
  * (30) E. H. Phillips and M. Fatica, “Implementing the Himeno benchmark with CUDA on GPU clusters,” in _IPDPS_ , 2010.
  * (31) L.-N. Pouchet, “Polybench: The polyhedral benchmark suite,” _URL: http://www. cs. ucla. edu/pouchet/software/polybench_ , 2012.
  * (32) J. Ragan-Kelley, C. Barnes, A. Adams, S. Paris, F. Durand, and S. Amarasinghe, “Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines,” in _PLDI_ , 2013.
  * (33) K. Sano, Y. Hatsuda, and S. Yamamoto, “Multi-FPGA accelerator for scalable stencil computation with constant memory bandwidth,” _TPDS_ , 2013.
  * (34) V. Seshadri, D. Lee, T. Mullins, H. Hassan, A. Boroumand, J. Kim, M. A. Kozuch, O. Mutlu, P. B. Gibbons, and T. C. Mowry, “Ambit: In-memory accelerator for bulk bitwise operations using commodity DRAM technology,” in _MICRO_ , 2017\.
  * (35) Y. S. Shao, B. Reagen, G.-Y. Wei, and D. Brooks, “Aladdin: A pre-rtl, power-performance accelerator simulator enabling large design space exploration of customized architectures,” in _ISCA_ , 2014.
  * (36) G. Singh, D. Diamantopoulos, C. Hagleitner, J. Gómez-Luna, S. Stuijk, O. Mutlu, and H. Corporaal, “NERO: A Near High-Bandwidth Memory Stencil Accelerator for Weather Prediction Modeling,” in _FPL_ , 2020.
  * (37) G. Singh, D. Diamantopoulos, C. Hagleitner, S. Stuijk, and H. Corporaal, “NARMADA: Near-Memory Horizontal Diffusion Accelerator for Scalable Stencil Computations,” in _FPL_ , 2019.
  * (38) H. Stengel, J. Treibig, G. Hager, and G. Wellein, “Quantifying performance bottlenecks of stencil computations using the execution-cache-memory model,” in _ICS_ , 2015.
  * (39) L. Szustak, K. Rojek, and P. Gepner, “Using Intel Xeon Phi Coprocessor to Accelerate Computations in MPDATA Algorithm,” in _PPAM_ , 2013.
  * (40) A. Taflove, “Review of the formulation and applications of the finite-difference time-domain method for numerical modeling of electromagnetic wave interactions with arbitrary structures,” _Wave Motion_ , 1988.
  * (41) F. Thaler, S. Moosbrugger, C. Osuna, M. Bianco, H. Vogt, A. Afanasyev, L. Mosimann, O. Fuhrer, T. C. Schulthess, and T. Hoefler, “Porting the COSMO Weather Model to Manycore CPUs,” in _PASC_ , 2019.
  * (42) P.-A. Tsai, N. Beckmann, and D. Sanchez, “Jenga: Software-defined cache hierarchies,” in _ISCA_ , 2017.
  * (43) P.-A. Tsai, C. Chen, and D. Sanchez, “Adaptive scheduling for systems with asymmetric memory hierarchies,” in _MICRO_ , 2018.
  * (44) J. van Lunteren, R. Luijten, D. Diamantopoulos, F. Auernhammer, C. Hagleitner, L. Chelini, S. Corda, and G. Singh, “Coherently Attached Programmable Near-Memory Acceleration Platform and its application to Stencil Processing,” in _DATE_ , 2019.
  * (45) M. Wahib and N. Maruyama, “Scalable Kernel Fusion for Memory-Bound GPU Applications,” in _SC_ , 2014.
  * (46) S. Williams, A. Waterman, and D. Patterson, “Roofline: an insightful visual performance model for multicore architectures,” _Communications of the ACM_ , vol. 52, no. 4, pp. 65–76, 2009.
  * (47) Y. Yarom, Q. Ge, F. Liu, R. B. Lee, and G. Heiser, “Mapping the intel last-level cache.” _IACR Cryptology ePrint Archive_ , vol. 2015, p. 905, 2015\.