# An Efficient Vectorization Scheme for Stencil Computation

Kun Li  Institute of Computing Technology, Chinese Academy of
SciencesUniversity of Chinese Academy of SciencesBeijingChina
<EMAIL_ADDRESS>,  Liang Yuan  Institute of Computing Technology, Chinese
Academy of SciencesBeijingChina<EMAIL_ADDRESS>,  Yunquan Zhang
Institute of Computing Technology, Chinese Academy of SciencesBeijingChina
<EMAIL_ADDRESS>,  Yue Yue  Institute of Computing Technology, Chinese Academy
of SciencesUniversity of Chinese Academy of SciencesBeijingChina
<EMAIL_ADDRESS>,  Hang Cao  Institute of Computing Technology, Chinese
Academy of SciencesUniversity of Chinese Academy of SciencesBeijingChina
<EMAIL_ADDRESS>and  Pengqi Lu  Institute of Computing Technology, Chinese
Academy of SciencesUniversity of Chinese Academy of SciencesBeijingChina
<EMAIL_ADDRESS>
(2021)

###### Abstract.

Stencil computation is one of the most important kernels in various scientific
and engineering applications. A variety of work has focused on vectorization
and tiling techniques, aiming at exploiting the in-core data parallelism and
data locality respectively. In this paper, the downsides of existing
vectorization schemes are analyzed. Briefly, they either incur data alignment
conflicts or hurt the data locality when integrated with tiling. Then we
propose a novel transpose layout to preserve the data locality for tiling and
reduce the data reorganization overhead for vectorization simultaneously. To
further improve the data reuse at the register level, a time loop unroll-and-
jam strategy is designed to perform multistep stencil computation along the
time dimension. Experimental results on the AVX-2 and AVX-512 CPUs show that
our approach obtains a competitive performance.

Stencil, Vectorization, Data locality, Data alignment conflict

††journalyear: 2021††copyright: acmlicensed

##  1\. Introduction

Stencil is one of the most important kernels widely used across a set of
scientific and engineering applications. It is extensively involved in various
domains from physical simulations to machine learning (Li et al., 2019).
Stencil is also included as one of the seven computational motifs presented in
the Berkeley View (Asanovic et al., 2008; Yuan et al., 2017) and arises as a
principal class of floating-point kernels in high-performance computing.

A stencil contains a pre-defined pattern that updates each point in a
$d$-dimensional spatial grid iteratively along the time dimension. The
stencil’s order (Yuan et al., 2017; Zhao et al., 2019) defines the dependent
relationship in a certain direction. If the order of a symmetric stencil in
one dimension is $r$, the value of one point at time $t$ is a weighted sum of
$(2r+1)$ points at the previous time (Tang et al., 2011). The naive
implementation for a $d$-dimensional stencil contains $d+1$ loops where the
time dimension is traversed in the outmost loop and all grid points are
updated in inner loops. Since stencil is characterized by this regular
computational structure, it is inherently a bandwidth-bound kernel with a low
arithmetic intensity and poor data reuse (Yuan et al., 2019; Krishnamoorthy et
al., 2007).

Performance optimizations of stencils have been exhaustively investigated in
the literature (Dursun et al., 2012; Luo et al., 2015; Dütsch et al., 2014).
Traditional approaches have mainly focused on either vectorization or tiling
schemes, aiming at improving the in-core data parallelism and the data
locality in cache respectively. These two approaches are often regarded as two
orthogonal methods working at different levels. Vectorization seeks to utilize
the SIMD facilities in CPU to perform multiple data processing in parallel,
while tiling tries to increase the reuse of a small set of data fit in cache.
They actually complement each other and can be subtly combined.

Prior work on vectorization of stencil computations primarily falls into two
categories. The first one is based on the associativity of the weighted sums
of neighboring points. Specifically, the execution order of one stencil
computation can be rearranged to exploit common sub-expression or data reuse
at register or cache level (Zhao et al., 2019; Basu et al., 2015; Rawat et
al., [n. d.], 2018). Consequently, the number of load/store operations can be
reduced and the bandwidth usage is alleviated in optimized execution order.
The second one attempts to deal with the data alignment conflicts (Henretty et
al., 2011, 2013), which is the main performance-limiting factor. The data
alignment conflict is a problem caused by vectorization. One milestone
approach to address the data alignment conflict is the DLT method (Dimension-
Lifting Transpose) (Henretty et al., 2011). We will present a deep discussion
on them in the next section.

As one of the crucial transformation techniques to exploit the parallelization
and data locality for stencils, tiling, also known as blocking, has been
widely studied for decades. Since the size of working sets in stencil-based
applications is generally larger than the cache capacity on a processor, the
spatial tiling algorithms are proposed to explore the data reuse by changing
the traversal pattern of grid points in one time step. Generally, a grid point
in cache is utilized to perform stencil computation for all its neighbors
before swapped out cache. Thus, the data transfers between the cache and main
memory could be reduced. However, the improvement of such tiling techniques is
restricted to the size of the neighbor pattern (Yuan et al., 2017;
Krishnamoorthy et al., 2007). Temporal tiling techniques have been developed
to allow more in-cache data reuse across the time dimension.

The aforementioned two approaches of stencil computation optimizations often
have no influence on the implementation of each other. The fundamental reason
is that the vectorization typically applies to the innermost loop. Therefore,
integrating one technique of vectorization with another tiling scheme is often
straightforward. However, the data organization overhead for vectorization may
degrade the data locality. Furthermore, to the best of our knowledge, most of
the prior work only focuses on temporal tiling technique on the cache level.
This only optimizes the data transfer volume between cache and memory and the
high bandwidth demands of CPU-cache communication is still unaddressed or even
worse with vectorization. We will present a deep discussion of these two
problems in the next section.

In this paper, we first design a novel transpose layout to overcome the input
data alignment conflicts of vectorization. The new layout is formed with an
improved in-CPU matrix transpose scheme, which achieves the lower bounds both
on the total number of data organization operations and the whole latency.
Compared with conventional methods, the corresponding computation scheme for
the new layout requires less data organization operations, whose cost can be
further overlapped by arithmetic calculations. To enhance the data reuse on
the register level, we propose an approach to perform multiple time steps for
stencil computations. The in-register data can be reused to perform successive
updates along the time dimension, which has not explored in existing work.
Finally, we integrated the proposed layout with a tiling framework. It only
requires a slight modification of the new vectorization scheme to preserve the
data reuse ability of tiling. The proposed vectorization scheme is evaluated
with AVX-2 and AVX-512 instructions for 1D, 2D, and 3D stencils. The results
show that our approach is obviously competitive with the existing highly-
optimized work(Henretty et al., 2013, 2011; Yuan et al., 2017).

This paper makes the following contributions:

  * •

We propose an efficient transpose layout and corresponding vectorization
scheme for stencil computation. The layout transformation utilizes an improved
matrix transpose of the lowest latency.

  * •

We exploit the in-register data reuse by performing multiple time step
computation based upon the new proposed transpose layout.

  * •

An integrated approach is proposed to perform a tiling framework in
conjunction with the vectorization scheme.

  * •

We demonstrate that the proposed approach could achieve superior performance
compared to several highly-optimized stencil benchmarks on multi-core
processors.

The paper is organized as follows. Section 2 presents the relevant background
and elaborates on the addressed problem. Section 3 introduces the proposed
vectorization scheme and the tiling technique. Section 4 provides experimental
results that demonstrate our approach produces a higher performance compared
to the benchmarks. In Section 5, we present the related work and Section 6
concludes the paper.

##  2\. Background

###  2.1. Data Alignment Conflicts of Vectorization

We take the 1D3P stencil as an example to illustrate the fundamental problem
of stencil computations caused by vectorization. Since in most existing work
vectorization is restricted to innermost loops (Datta et al., 2009), the codes
shown in Figure 1 only illustrates the 1D3P stencil of one time step.

In the $i$-th iteration of this scalar code execution, it loads $A[i+1]$ and
$B[i]$ to registers and reuses register data $A[i-1]$ and $A[i]$ referenced by
the previous calculation of $B[i-1]$. Observing the CPU-memory data transfer,
this code is exactly similar to a common array copy code, i.e. the `memcpy`
function (Falke et al., 2013). The computation implementation inside CPU is
straightforward. Loop optimizations like loop unrolling also preserve these
properties.

The vectorization groups a set of data in a vector register and processes them
in parallel. The naive vectorization of the 1D3P stencil code computes
contiguous elements in the output array $B$. Assume the vector register holds
4 elements (vector length $vl=4$), the vectorization code performs the
calculation using vector operations and output $(B[1],B[2],B[3],B[4])$ with
one vector register.

A well-known problem incurred by the vectorization of stencil codes is the
input data alignment conflicts. For example, to compute
$(B[1],B[2],B[3],B[4])$, it requires three vectors: $(A[0],A[1],A[2],A[3])$,
$(A[1],A[2],A[3],A[4])$ and $(A[2],A[3],A[4],A[5])$. The element $A[2]$
appears in all these vector registers but at different positions. We call this
a data alignment conflict. Thus there is no corresponding simple execution as
the scalar code.

Figure 1. Illustration of input data alignment conflicts handling in DLT.

To address the data alignment problem, two common implementations are often
adopted. The first one loads all the needed elements from memory in a vector
form straightforward. Due to the low operational intensity, the stencil
computation is often regarded as a memory-starving application. Compared with
the scalar code, this multiple load vectorization method further increases the
data transfer volume. Moreover, in each iteration of this code, it has at
least two unaligned memory references where the first data address is not at a
32-byte boundary. Since CPU implementations favor aligned data loads and
stores, these unaligned memory references will degrade the performance
considerably.

The second solution is similar to the scalar code in terms of the CPU-memory
data transfer. It loads each input element to vector register only once and
assembles the required vectors via inter-register data permutations
instructions. Compared with the multiple load method, this data permutations
method reduces the memory bandwidth usage and takes the advantages of the rich
set of data-reordering instructions supported by most SIMD architectures.
However, the execution unit for data permutations inside the CPU may become
the bottleneck.

###  2.2. Dimension-Lifting Transpose (DLT)

One milestone approach to address the data alignment conflict is the DLT
method (Henretty et al., 2011). In DLT the original one-dimensional array of
length $N$ is viewed as a matrix of size $vl$*$(N/vl)$, where $vl$ is the
vector length in vector elements. For example, $vl$=4 for double-precision
floats in a 256-bit vector. It then performs a global transpose. Figure 1
illustrates the DLT method for a one-dimensional array of 28 elements. The DLT
layout overcomes the input data alignment conflicts. For instance, the second
$vl=4$ elements in the transformed layout are formed into one vector (B[1],
B[8], B[15],B[22] and all the three required input vectors: left vector (A[0],
A[7], A[14], A[21]), center vector (A[1],A[8],A[15],A[23]) and right vector
(A[2], A[9], A[15], A[23]) are free of data sharing and stored contiguous in
memory. DLT needs to assemble input vectors for calculating output vectors at
boundary.

DLT has the following disadvantages. First, DLT can be viewed as $vl$
independent stencils if we ignore the boundary processing. Therefore when
incorporated with blocking frameworks, the data reuse decreases $vl$ times.
The reason is that there is no data reuse among the $vl$ independent stencils.
Second, DLT suffers from the overhead of explicit transpose operations
executed before and after the stencil computation. For 1D and 2D stencils in
scientific applications, the number of time loops is often large enough to
amortize the transpose overhead. But for 3D stencils and low-dimensional in
other applications like image processing, the time size is small that makes
the transpose overhead unignorable. Finally, it’s hard to implement the DLT
transpose in-place and it often chooses to use an additional array to store
the transposed data. This increases the space complexity of the code.

##  3\. The Transpose Layout

In this section, we first discuss the drawbacks of existing vectorization
methods. Then we present a new transpose layout and its corresponding stencil
computation scheme. Next, we present several further optimizations on the
transpose layout including the extension to multiple time steps, integration
with a tiling framework and an improved matrix transpose algorithm.

###  3.1. Motivation

From the hardware perspective, the critical approach to boost performance is
to fully utilize the execution units that perform the arithmetic instructions.
Since there is no data dependence in one time-step iteration of stencil
computations, the only bottleneck is data preparation. Equivalently, the key
technique of vectorization is to address the data alignment conflict.

Our starting point is the observation of the disadvantages of existing
methods. The DLT is a promising method that extremely reduces the data
reorganization operations. However, essentially the DLT vectorization format
hurts the locality properties as mentioned above. In particular, the elements
in one vector are distant, thus there is no data reuse among them.

On the contrary, the straightforward multiple load and data reorganization
methods load contiguous element in one vector. They lead to the optimal data
locality when integrated with a temporal tiling scheme.

These two methods seem to be at two extreme ends of a balance between the
number of reorganization operations of data in CPU and the reuse ability of
data in cache. Our scheme seeks to preserve the data locality property and
employs the fundamental idea of DLT to improve the overhead of data
preparation.

Figure 2. Register Transpose Layout for SIMD vector length of 4. Figure 3.
Illustration of stencil computation for transpose layout.

###  3.2. Locally Transpose

To preserve the data locality and reduce the number of data organization
operations, we apply a matrix transpose to a small sub-sequence of contiguous
elements. Specifically, like the dimension-lifting approach in DLT, the one-
dimensional view to the sub-sequence is substituted by a two-dimensional
matrix view. To perform vectorization after a matrix transpose, the column
size of the matrix should be equal to the vector length $vl$. Let the row size
be $m$, the size of the matrix is then $vl*m$. Our locally transpose layout is
equivalent to the DLT when $m=N/vl$ and the original data layout when $m=1$.

After the matrix transpose, it still requires some data reorganizations for
computing the first and last one of the $m$ vectors. For example, if $m=1$,
the original data preparations of the left and right vectors must be done for
computing each output vector. This is the trade-off between data locality and
the number of data preparations explained above.

There are several considerations for deciding the size $m$. First, $m$ should
be large enough to hidden the overhead of the data reorganizations for the
first and last vectors by the actually arithmetic operations of the middle
vectors. Assume the order of a stencil is $r$, then the number of arithmetic
operations of the middle vectors is $(2r+1)*(m-1)+1$. The number of data
operations is $4r$ since the first and last vectors need $2r$ vectors and
assembling each of them requires two reorganization instructions as will be
explained shortly. Thus $m$ should be at least $3$. Notice that this
limitation is irrelevant to the order $r$. Second, to avoid an additional
array that is needed to store the transposed data as in the DLT format, it’s
desirable to complete the matrix transpose in CPU. Thus the $m$ input vectors
and additional auxiliary vectors must be kept in the CPU vector register file.
In this work we always set $m=vl$. The final reason is that transposing a
matrix of size $vl*vl$ is easier to implement on modern CPU products. We will
present a highly efficient algorithm for matrix transpose of size $vl*vl$
later.

Figure 2 illustrates the transpose layout for a one-dimensional stencil with a
vector length of four. The matrix transposes of every sub-sequence of $vl*vl$
length is performed before and after the stencil computation. In the rest of
the paper, we also refer to the $vl$ vectors as a vector set (VS). Note that
in the implementation a vector set is always aligned to a 32-Byte boundary.

The update of one vector set of the 1D3P stencil requires two assembled
vectors. One is the left dependent vector of its first vector and the other is
the right dependent vector of its last vector. Figure 3 describes the data
reorganization of these two vectors. The first vector is $(A,E,I,M)$ and its
left dependent vector is $(Z,D,H,L)$ which is stored in two distant vectors in
the transpose layout, $(*,*,*,Z)$ and $(D,H,L,*)$. These two vectors are
combined by a blend instruction followed by a permute operation to shift the
components to the right circularly.

The stencil computations of the vector set are straightforward as shown in
Figure 3. We then achieve an efficient vectorization scheme by performing
lower-overhead matrix transpose and two data operations per vector set.
Moreover, the proposed vectorization scheme avoids data reloads compared with
the multiple load method and frequent inter-vector permutations compared with
the data reorganization method. The transpose layout could also be applied to
higher-order and multidimensional stencils in the same manner.

###  3.3. Unroll-and-jam the Time Loop

Figure 4.  Illustration of stencil computation for two time steps.

In general, stencil computation is restricted to its input data alignment
conflicts, and all elements are only updated once before the round starts in
the next time step. Although blocking technique (Yuan et al., 2019, 2017;
Bandishti et al., 2012a) can be utilized to decrease the data transfers
between main memory and cache, there is no in-register data reuse between
successive time loops. Therefore the in-CPU flops/byte ratio is limited by the
stencil pattern. To the best of our knowledge, computation for multiple time
steps in registers is not explored in existing work.

We develop an unroll-and-jam strategy of the time loop. It loads one element
at time dimension $t$ and updates it $k$ time steps before store it to memory.
$k$ is called the unrolling factor. The normal execution corresponds to the
case of $k=1$. If $k>1$ the execution is equivalent to unroll the time loop
$k$ times and jam them. Consequently, it improves the in-CPU flops/byte ratio
by $k$ times. A one-dimensional example is illustrated in this subsection and
the strategy is also applicable to multidimensional stencils.

Algorithm 1 Unroll-and-jam the Time Loop

1:function Assemble($\mathbf{v}_{a}$,$\mathbf{v}_{b}$)

2: $\mathbf{v}_{c}$=_mm256_blend_pd($\mathbf{v}_{a}$,$\mathbf{v}_{b}$)

3: $\mathbf{v}_{c}$=_mm256_permute4$\times$64_pd($\mathbf{v}_{c}$)

4: return $\mathbf{v}_{c}$

5:end function

6:function Compute($\mathbf{v}_{left}$, $\mathbf{v}_{1}$, $\mathbf{v}_{2}$,
$\mathbf{v}_{3}$, $\mathbf{v}_{4}$, $\mathbf{v}_{right}$)

7: $\mathbf{v}_{0}\leftarrow$ Assemble($\mathbf{v}_{left}$,$\mathbf{v}_{4}$)

8: $\mathbf{v}_{5}\leftarrow$ Assemble($\mathbf{v}_{1}$,$\mathbf{v}_{right}$)

9: for $i=1\to 4$ do

10: $\mathbf{v}_{i-1}\leftarrow$ Stencil($\mathbf{v}_{i-1}$, $\mathbf{v}_{i}$,
$\mathbf{v}_{i+1}$)

11: end for

12:
$\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3},\mathbf{v}_{4}\leftarrow\mathbf{v}_{0},\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3}$

13:end function

14:function MultipleTimeSteps($\mathbf{VS}_{1:k}$,$\mathbf{vrl}_{0:k-1}$,$k$)

15: for $j=k+1\to N$ do

16: $\mathbf{VS}_{k+1}\leftarrow$ Load the $j$-th Vector Set

17: for $i=k\to 1$ do

18: $\mathbf{vrl}_{i}\leftarrow\mathbf{VS}_{i}[3]$

19: Compute($\mathbf{vrl}_{i-1},\mathbf{VS}_{i}[0:3],\mathbf{VS}_{i+1}[0]$)

20: end for

21: Store $\mathbf{VS}_{1}$

22: for $i=1\to k$ do

23: $\mathbf{VS}_{i}\leftarrow\mathbf{VS}_{i+1}$

24: $\mathbf{vrl}_{i-1}\leftarrow\mathbf{vrl}_{i}$

25: end for

26: end for

27:end function

28:function Main( )

29: while $t<T$ do

30: Booting computation.

31: MultipleTimeSteps($\mathbf{VS}_{1:k}$,$\mathbf{vrl}_{0:k-1}$,$k$)

32: Epilogue computation.

33: $t+=k$

34: end while

35:end function

Overall the algorithm is straightforward. After update one vector set, we keep
the result in registers and process the next neighbor vector set. Then the
current vector set can be forwarded along time dimension one more step using
the new value of the right neighbor.

Algorithm 1 shows the pseudo-code of our multiple time steps updating scheme.
The Compute function receives a set of $vl$ vectors and their dependent
vectors that are assembled by the Assemble function. It computes the elements
in the vector set by one time step. Notice that this is an in-place updating
that the value of last time will be overwritten.

The main function traverses the time loop stepped by the unrolling factor $k$.
For simplicity, we assume $T$ is divisible by $k$. In each iteration of the
while loop, every element is forwarded $k$ steps along the time dimension. The
booting computation prepares the data at head needed by the following
pipelined updating. The top part of Figure 4 illustrates the case of $k=2$
after a booting computation. The vector sets $\mathbf{VS}_{1}$ to
$\mathbf{VS}_{k}$ have been updated $k-1$ to $0$ times, respectively. Due to
the overwriting property of the Compute function, it needs to preserve the
value of the last time of the vector to each vector set’s left, denoted as
$\mathbf{vrl}_{i}$. As the figure shows, $\mathbf{vrl}_{i}$ and
$\mathbf{VS}_{i}[3]$ store the value of the same vector at time $t-1$ and $t$,
respectively.

The MultipleTimeSteps function forwards all the vector sets from right to left
by one time step. Meanwhile, it preserves the old value of their rightmost
vector in $\mathbf{vrl}$. At the end of each iteration, $\mathbf{VS}_{1}$ has
been updated $k$ times and is stored in memory. Then after some data
reassignments, the next loop is ready to execute. Each iteration loads and
stores one vector set of $vl*vl$ elements and performs $k*vl*vl$ stencil
computations. As mentioned above, it increases the in-CPU flops/byte ratio by
$k$ times.

From the algorithm, we see that it needs $k$ vector sets and $k$ additional
vectors to unroll-and-jam the time loop, i.e., total $(vl+1)*k$ registers in
addition to coefficient vector registers. In modern CPUs, the typical number
of available vector registers is $vl*4$, where $vl$ is the capacity of double
precision variables in one register, therefore in this work we always set
$k=2$.

There is another advantage of the algorithm. Conventionally the stencil of
Jacobi style is implemented with two arrays, storing the value at odd and even
time respectively. If we set $k=2$, then the input and output value are all at
the even time. It’s legal to reuse the input data space and make the whole
computation in-place. The space usage is then reduced.

###  3.4. Integrated With Tiling

Vectorization and tiling are two orthogonal methods. They target at different
levels. Vectorization boosts the computation using the data parallelism at the
execution level, while tiling serves to exploit the data reuse at cache
levels. The transpose layout described above identifies a vectorization
technique as the solution to the data alignment conflict for stencils. The
multiple time update further improves the data reuse ability at the CPU vector
register level. In the following, we present the combination of the transpose
layout and a tiling framework.

The tessellation tiling (Yuan et al., 2019) can be viewed as a tessellation in
iteration space by utilizing shaped tiles. Figure 5 (a) and Figure 5 (b)
illustrate the tiling framework for a one-dimensional stencil. The iteration
space is tessellated by triangles and inverted triangles in alternative
stages. Thus, concurrent execution is processed by two stages which are
started in each triangle with a given time range first, followed closely by
the execution of inverted triangles over the same time range concurrently.
Updates in different time steps are distinguished from each other by different
colors, and the state of each element along the time dimension is represented
with a number in Figure 5. For the example in Figure 5 (a), the new state of
each triangle contains (0,1,2,3,4,3,2,1,0) where the center element is updated
four steps and its neighbors are updated fewer steps proportional to the
distance with the center element. To make all elements updated with the same
steps, two half parts from adjacent triangles constitute new inverted
triangles and the elements are updated with the state (4,3,2,1,0,1,2,3,4). As
Figure 5 (c) shows, all elements are updated to four steps by adding the
projection of the triangles with inverted triangles. With the tessellate
tiling strategy, concurrent execution for different tiles is enabled over a
given time range without redundant computation.

Figure 5. Tessellate tiling iteration space for 1D updated with two time steps
on register transpose layout.

The only problem for applying the transpose layout is the calculations at the
two boundaries of each block. The execution of triangles is a ’shrinking’
process, the range of processed elements decreases as the time forwards.
Similarly, an ’expanding’ process occurs in the execution of inverted
triangles. Since the physical neighbor elements are stored apart from each
other in one vector set, the calculations of the vector set that covers a
boundary are too complex to implement. As the basic computing unit in the
transpose layout is a vector set, we convert the vector set at boundary back
to the original format before the computation and employ a simple data
reorganization method to process them. As illustrated in Figure 5 (d), the
shrinking and expanding process could be simplified in this way. When the
boundary slides away, the vector set is transposed again.

Further, the register transpose layout and time loop fusion make it feasible
to achieve multiple time steps computation in registers over the tiles
efficiently without reloading operations.

The tessellate tiling could also be applied for multidimensional stencil
computations. For a $d$-dimensional stencil, tessellation in iteration space
contains $d+1$ stages. Similar to the 1D stencil example in Figure 5, the
spatial space in stage $i$ is tessellated by $tiles_{i}$ ($0\leq i\leq d$).
$tiles_{1}$ is a hypercube (typically a line segment in 1D, square in 2D, cube
in 3D). $tiles_{i+1}$ is built by recombining the sub-tiles split from
adjacent $tiles_{i}$ along some dimensions. Applying the transpose layout to
higher-dimensional stencils is exactly similar to the one-dimensional case
since the layout only affects the unit-stride dimension.

###  3.5. Transpose

Figure 6. Transpose for $double$ type using AVX-2 instructions.

Unlike previous work (Henretty et al., 2011) that performs a global dimension-
lifted transformation, we only need a transpose on-the-fly for each register
set twice throughout the whole process. The lower bound on the memory
operations for completing a matrix transpose of size $vl*vl$ is $vl\log(vl)$,
e.g., 8 data reorganization instructions for $vl=4$.

In modern CPU architectures, these 8 instructions can be launched continuously
in 8 cycles. However, the implementation of existing algorithm adopts lane-
crossing instructions, which increases the overhead by 25%. Figure 6
illustrates our improved version where the long-latency instructions are
hidden by their following single-cycle instructions. In the first stage, pairs
of two vectors with distance $2$, e.g., $(A,B,C,D)$ and $(I,J,K,L)$, exchange
data using the `permute2f128` instruction. In the second stage, the pairs of
two adjacent vectors, e.g., $(A,B,I,J)$ and $(E,F,M,N)$, swap elements by the
`unpackhi` or `unpacklo` instruction. The total cost of the new transpose
scheme is then reduced to 8 cycles. Similarly, the transpose by using AVX-512
instructions contains three stages where the last stage consists of in-lane
instructions.

##  4\. Evaluation

In this section, we evaluate our proposed scheme for 1D, 2D and 3D stencils
with AVX-2 and AVX-512 instructions.

###  4.1. Setup

Table 1. Parameter description for stencils used in experiments Dim | Pts | Problem Size | Blocking Size
---|---|---|---
1D | 3 | 10240000$\times$1000 | 2000$\times$1000
1D | 5 | 10240000$\times$1000 | 2000$\times$500
2D | 5 | 3000$\times$3000$\times$1000 | 200$\times$200$\times$50
2D | 9 | 3000$\times$3000$\times$1000 | 120$\times$128$\times$60
3D | 7 | 128$\times$128$\times$128$\times$1000 | 23$\times$23$\times$10
3D | 27 | 128$\times$128$\times$128$\times$1000 | 23$\times$23$\times$10

Our experiments were performed on a machine composed of two Intel Xeon Gold
6140 processors with 2.30 GHz clock speed, which owns 36 physical cores
organized into two sockets. Each core contains a 32KB private L1 data cache, a
1 MB private L2 cache, and a unified 24.75MB L3 cache. AVX-512 instruction set
extension is supported and it’s able to conduct operations for 8 double-
precision floating point data in a SIMD manner, which yields a theoretical
peak performance of 73.6 GFlop/s/core (2649.6 GFlop/s in aggregate).

Since the recent tiling technique proposed by Yuan (Yuan et al., 2017) and the
nested/hybrid tiling technique (denoted as SDSL, which is the name of the
software package.) presented by Henretty (Henretty et al., 2013) outperform
the other stencil research like Pluto (Bondhugula et al., 2008; Bandishti et
al., 2012a) and Pochoir (Tang et al., 2011), we take them as two counterparts
of our work, which are vectorized by data multiple load and DLT methods,
respectively. All programs were compiled using the ICC compiler version
19.0.3, with the ’-O3 -xHost -qopenmp -ipo’ optimization flags.

The detailed parameters for stencils of various orders used in experiments are
described in Table 1, which consists of four star stencils (1D 3-Points, 1D
5-Points, 2D 5-Points, and 3D 7-Points) and two box stencils (2D 9-Points and
3D 27-Points) corresponding to the references (Yuan et al., 2017; Henretty et
al., 2013). The default value of total time steps is 1000 or 200 in the
references. Thus, we fix it as a larger value of 1000 in our experiments.
Other parameters of each stencil are also fine-tuned on the basis of
references work to guarantee that the peak performance for all methods could
be reached exactly. Since the performance is sensitive to the stencil
parameters, significant efforts are required in automatic tuning and this will
be done separately as future work.

###  4.2. Sequential Block-free Results

Figure 7. Absolute performance comparison for tested methods in single-thread
blocking-free experiments. The results are shown separately with different
total time steps.

In this subsection, we present performance results of varied methods across
problem sizes ranging from L1 cache to main memory with a single thread. The
spatial and temporal blocking method is not applied to them for examining the
pure improvements on various storage levels. The multiple loads and data
reorganization methods represent a class of auto-vectorization in modern
compilers (Yuan et al., 2017). DLT is the dimension-lifting transpose strategy
designed by Henretty (Henretty et al., 2011). All the methods are implemented
by hand-written codes optimized with the appropriate strategies such as
alignment and loop unrolling to ensure fairness.

Figure 7 shows the performance comparison of our methods with the other three
methods. The results are illustrated separately in two subfigures on the basis
of the total time steps $T$. It can be seen that our method updating two time
steps outperforms others apparently in both experiments, which demonstrates
the effectiveness of the improvement of the flop/byte ratio. Our method
without time loop unroll-and-jamming also achieves better performance results
than the hand-written DLT in most cases. The performance has a decrease at the
size of 1000 in L1 cache. This can be attributed to the cheaper dimension-
lifting transpose operation in small size for DLT. The multiple loads method
exhibits the worst performance among them due to the overhead caused by
redundant loads.

Table 2. Performance improvements on different storage level in single-thread blocking-free experiments Storage | Data | DLT | Our | Our
---|---|---|---|---
Level | Reorganization | (2 steps)
L1 | 1.28x | 2.06x | 2.16x | 3.13x
L2 | 1.11x | 1.37x | 1.67x | 2.07x
L3 | 1.01x | 0.95x | 2.02x | 2.92x
Memory | 1.00x | 1.01x | 1.97x | 2.96x
Mean | 1.11x | 1.35x | 1.98x | 2.81x

To further investigate the effect of total time steps $T$, we perform a
tenfold increase on the default value to $T=10000$, which is illustrated in
Figure 7 (b). It can be observed that the performance trends of $T=10000$ are
still largely consistent with the results in Figure 7 (a). However, the
performance of our method falls slightly behind the DLT in L1 cache, and this
performance anomaly is primarily due to the diluted dimension-lifting
transpose cost by overly long time steps. Notably, only the performance of DLT
in L1 cache drops gradually as problem size increases for both results in
Figure 7, which is resulted from a costly data layout transformation and
indicates a potential bottleneck for cache-blocking.

###  4.3. Multicore Cache-blocking Experiments

Figure 8. Absolute performance comparison for varied methods in multicore
cache-blocking experiments. The results are shown separately with different
total time steps.

In this subsection, we present the experiment results that exhibit the
benefits of our methods with the temporal blocking and parallelization scheme.
The SDSL employs a split tiling technique (nested tiling in 1D, hybrid tiling
for higher dimensions) to achieve temporal blocking. The tessellate tiling
technique utilized auto-vectorizing supported by the compiler (Yuan et al.,
2017).

The results are shown in Figure 8 (a) and Figure 8 (b) with time steps of
$T=1000$ and $T=10000$ respectively. As can be seen from Figure 8 (a), the
performance drops apparently as the problem size moves from L3 cache to the
memory hierarchy, which is mainly caused by the cost of data transfers. We
also further investigate the influence of the block size on performance. In
the case of L1 blocking, the observed performance is higher than that with L2
blocking overall. Since the smaller stencils could be prefetched into cache
directly, the performance gap between different blocks is further aggravated
when the problem size lies in the memory hierarchy. Surprisingly, our method
with two time steps could still take up a leading position, approximately
3.29x and 3.48x improvements are obtained compared to SDSL with L1 blocking
and L2 blocking respectively. The performance of SDSL is inferior to
tessellation, which is resulted from the blocking technique constrained to its
data layout. Longer time steps of $T=10000$ are evaluated in Figure 8 (b), and
similar performance trends but higher values are observed compared with Figure
8 (a).

Table 3 shows the detailed performance improvements on different storage
levels as before. Our method could obtain better optimization results when the
problem size lies in L3 cache and memory. The speedup ranges from 2.54 to
2.76x with L1 blocking, showing that our method integrated with tiling
provides a significant benefit over others on varied problem size.

###  4.4. Scalability

We also evaluated the scalabilities of our schemes and the counterparts. The
detailed parameters are given in Table 1, where all problem sizes exceed the
L3 cache. Since our tiling framework is the same as the tessellation scheme,
the performance improvements of our method with respect to the tessellation
method are fully derived from the vectorization.

Figure 9 illustrates the results of 1D, 2D and 3D stencils implemented with
AVX-2 and AVX-512 instructions respectively. It can be observed that our
method could obtain the highest performance while the SDSL performs the lowest
performance. In one-dimensional stencils, all these methods achieve nearly
linear scaling on both instruction sets and the proposed time loop fusion
strategy provides a significant improvement. With the increase of the problem
dimension, the scalability for all methods drops as a result of the inherent
complexity for multidimensional stencil computations. Compared to the results
implemented with AVX-2 instructions, the performance of the right half in
Figure 9 shows a slight increase.

Table 3. Performance improvements on different storage level in multicore cache-blocking experiments | Blocking | Tessellation | Our | Our
---|---|---|---|---
| Level | (Two time steps)
L3 Cache | L1 | 1.43x | 2.54x | 2.99x
L2 | 1.21x | 2.58x | 3.01x
Memory | L1 | 1.62x | 2.76x | 3.42x
L2 | 1.39x | 2.92x | 3.58x
Mean | L1 | 1.56x | 2.69x | 3.29x
L2 | 1.32x | 2.79x | 3.48x
Figure 9. Performance comparison for stencils of various orders with different
dimensions in a multicore environment. Table 4. Average performance
improvement for different stencils

| Method | 1D3P | 2d5P | 3D7P | 1D3P | 2d5P | 3D7P | 1D5P | 2d9P | 3D27P | 1D5P | 2d9P | 3D27P
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (AVX-2) | (AVX-2) | (AVX-2) | (AVX-512) | (AVX-512) | (AVX-512) | (AVX-2) | (AVX-2) | (AVX-2) | (AVX-512) | (AVX-512) | (AVX-512)
Speedup | SDSL | 1.00x | 1.00x | 1.00x | - | - | - | 1.00x | 1.00x | 1.00x | - | - | -
over | Tessellation | 1.77x | 1.54x | 1.39x | 1.00x | 1.00x | 1.00x | 1.56x | 1.61x | 1.50x | 1.00x | 1.00x | 1.00x
SDSL | Our | 2.77x | 2.05x | 1.85x | 1.50x | 1.24x | 1.31x | 2.37x | 1.91x | 1.66x | 1.49x | 1.26x | 1.15x
/Tessellation |  Our∗ | 3.52x | 2.26x | 1.97x | 1.86x | 1.34x | 1.38x | 2.92x | 2.12x | 1.76x | 1.98x | 1.39x | 1.24x
Speedup | SDSL | 29.2 | 13.5 | 9.4 | - | - | - | 29.1 | 22.0 | 20.7 | - | - | -
over | Tessellation | 29.8 | 24.9 | 11.5 | 29.0 | 24.1 | 12.4 | 26.2 | 24.1 | 24.2 | 24.5 | 24.0 | 23.6
single | Our | 30.8 | 26.3 | 20.1 | 30.7 | 22.4 | 17.7 | 31.3 | 25.9 | 24.7 | 28.2 | 24.9 | 24.5
core |  Our∗ | 32.1 | 26.7 | 21.2 | 31.4 | 22.3 | 18.2 | 31.8 | 25.7 | 24.9 | 28.7 | 25.6 | 23.8

  * *

Our method updated with two time steps.

The speedups and scalabilities for high-order stencils including 1D5P, 2d9P,
and 3D27P also decrease gradually from 1D to 3D. However, the overall
performance falls behind the corresponding one-order results, which is
resulted from complex data access patterns in high-order stencils. Our method
could also obtain a substantial performance improvement in all experiments.

The average speedup results are given in Table 4. Since the stencil kernels
with AVX-512 instructions are not available in the SDSL(Henretty et al.,
2013), the corresponding position is filled with a hyphen symbol. Thus the
comparison basis is then SDSL or Tessellation for AVX-2 or AVX-512. Taking all
stencils with AVX-2 instructions into account, remarkable performance benefits
are observed from our method updating two time steps, 3.52x and 2.92x
respectively for 1D3P and 1D5P. The performance improvement ranges from 1.66x
to 2.77x with a mean of 2.10x, demonstrating that our vectorization scheme
provides a significant benefit in a large problem size compared to the
referenced work.

The speedup for each method with 36 cores is also given at the bottom of the
table. For scalability, our method obtains a 20.1x speedup while the value of
DLT is only 9.4x for 3D7P, which indicates a sustainable performance for our
method in multidimensional stencils. Additionally, the largest speedup in each
stencil column again corresponds to the performance shown in Figure 9, where
our method outperforms others in most cases.

###  4.5. Discussion

In this subsection, we provide an analysis of the performance on various
configurations in previous experiments to tease out the contributions from
different aspects of our proposed scheme.

Sequential block-free experiments examine a variety of vectorization methods
and demonstrate that our scheme with multiple time steps updating can achieve
an considerable 2.81x improvement on average compared with the multiple loads
method. Subsequently, the performance gains for a larger time steps are still
significant and consistent with the results of the small time steps. Moreover,
the DLT method is more appropriate only on the relatively small size and long
time steps, and this is partly explained by the performance penalty associated
with additional dimension-lifting transpose in memory. Since the problem size
ranges from L1 cache to main memory, clear insights are provided that the
overall performance trends drop consistently with the various memory
hierarchy.

Multicore cache-blocking experiments conduct stencil cases with 36 cores, and
an average 2.69x speedup is obtained by our method on the basis of SDSL. Due
to the reduced data transfers by our time loop unroll-and-jam, our method
updated with two time steps achieve a further 3.29x speedup. We also study the
influence of blocking size, and the results prove that appropriate L1 blocking
or in-cache problem size could contribute to better performance for all
methods. The overall trends are in accord with the sequential block-free
experiments, and our method updated with two time steps outperforms others
obviously.

The scalability experiments demonstrate that our vectorized scheme leveraging
tessellate tiling successfully outperforms the referenced fastest multicore
stencil work to date across a broad variety of configurations. Constrained to
its specific data layout, DLT is slower than other methods. Since
multidimensional or high-order stencils are more compute-intensive, more
dependency data are loaded into cache while they are not fully utilized to
perform their own stencil computation. Thus, the overall performance for each
method falls gradually with the increasing dimensions or orders, and our
method could still obtain a better performance.

##  5\. Related Work

Research on optimizing stencil computation has been intensively studied
(Henretty et al., 2013; Yuan et al., 2019, 2017; Henretty et al., 2011; Kamil
et al., 2005), and it can be broadly classified as optimization methods to
improve the computation performance and enhance the data reuse.

Vectorization by using SIMD instructions is an effective way to improve
computation performance for stencils. Prior work on optimizing the order of
execution instructions could decrease loads/stores operations to relieve the
register pressure, while only the individual element in each vector could be
reused (Zumbusch, 2012). Basu designs a vector code generation scheme to reuse
several vectors in the computation process, and it is constrained to constant-
coefficient and isotropic stencils (Basu et al., 2015). YASK (Yount et al.,
2016) could improve data reuse by using common expression elimination and
unrolling based on their vector-folding methods with fine-grained blocks
(Yount, 2015), which is less feasible for high-order complex stencils (Zhao et
al., 2019). Henretty proposes a new method DLT (Henretty et al., 2011, 2013)
to overcome input data alignment conflicts at the expense of a dimension-
lifting transpose, which makes it infeasible to perfectly utilize the tiling
technique as a result of its spatially separated data elements (Krishnamoorthy
et al., 2007). Essentially DLT can be viewed as the combination of strip-
mining (1-dimensional tiling) and out-loop vectorization (Henretty et al.,
2013). Specifically, the original innermost loop traverses the corresponding
dimension from $1$ to $N$. In DLT the loop is transformed to a depth-2 loop
nest where the size of the outer loop equals the vector length $vl$ and the
inner loop processes each subsequence of length $N/vl$. Note that the strip-
mining was also introduced for vectorization (Allen and Kennedy, 2002).
However, the conventional usage is to make the size of the innermost loop be
the vector length and substitute it by a vector code. In addition, the in-
place matrix transpose involved in our work has also been widely studied and a
kernel of 4$\times$4 matrix transpose consists of two stages basically.
Hormati splits the vector register to some 128-bit lanes (Hormati et al.,
2010), and the lane-crossing instructions for $double$ incur a longer latency,
typically 3 to 4 cycles. Zekri (Zekri, 2014) use in-lane instructions in four
stages only for $float$ type. Springer(Springer et al., 2017) utilize Shuffle
and Permute2f128 instructions for $double$ type in two stages, while it
requires 8 integers as parameters.

Tiling is one of the most powerful transformation techniques to explore the
data locality of multiple loop nests (Bielecki and Palkowski, 2016; Hartono et
al., 2010; Grosser et al., 2014). Notably work for stencil computations
includes hyper-rectangle tiling (Nguyen et al., 2010; Andonov et al., 2003;
Iooss et al., 2015), time skewed tiling (Jin et al., 2001), diamond tiling
(Bandishti et al., 2012b), cache oblivious Tiling (Frigo and Strumpen, 2005),
split-tiling (Henretty et al., 2013) and tessellating (Yuan et al., 2019).
Wonnacott and Strout present a comparison on the scalability of many existing
tiling schemes (Wonnacott and Strout, 2013). Most of these techniques are
compiler transformation techniques and this paper integrated the new proposed
layout with the tessellation scheme for simplifying the implementation. For
stencil computations, a variety of auto-tuning frameworks (Christen et al.,
[n. d.]; Gysi et al., 2015; Rocha et al., 2017) have been presented by using
varied hyper-rectangular tiles to exploit data reuse alone. However, redundant
computations are involved in these work to resolve the introduced inter-tile
dependencies that hinder the concurrent execution of shaped tiles on different
cores.

##  6\. Conclusion

In this paper, we propose a novel transpose layout to overcome the input data
alignment conflicts efficiently for vectorization. A time loop unroll-and-jam
strategy with in-register multiple time steps processing is designed on the
basis of the proposed transpose layout. Furthermore, we describe how the
proposed vectorization scheme is integrated with a tessellate tiling framework
for enhancing data reuse and concurrency. With the qualitative analysis and
quantitative experiments, we demonstrate that significant performance
improvements are achieved by our vectorization scheme over state-of-the-art
products such as Intel’s ICC and recent work (Henretty et al., 2013; Yuan et
al., 2017). Experimental results provide evidence that our vectorization
scheme incorporated with tessellate tiling could obtain a linear scaling
character and reach a 3.29x improvement over SDSL for one-dimensional
stencils.

##  7\. Acknowledgements

This work is supported by the National Key R&D Program of China under Grant
No.2016YFB0200803.

## References

  * (1)
  * Allen and Kennedy (2002) Randy Allen and Ken Kennedy. 2002.  _Optimizing compilers for modern architectures: a dependence-based approach_.  Taylor & Francis US.
  * Andonov et al. (2003) R. Andonov, S. Balev, S. Rajopadhye, and N. Yanev. 2003\.  Optimal semi-oblique tiling.  _IEEE Transactions on Parallel and Distributed Systems_ 14, 9 (2003), 944–960.
  * Asanovic et al. (2008) Krste Asanovic, Ras Bodik, James Demmel, et al. 2008\.  The parallel computing laboratory at UC Berkeley: A research agenda based on the Berkeley view.  _University of California, Berkeley, Tech. Rep_ (2008).
  * Bandishti et al. (2012a) Vinayaka Bandishti, Irshad Pananilath, and Uday Bondhugula. 2012a.  Tiling stencil computations to maximize parallelism. In _SC’12_. IEEE, 1–11.
  * Bandishti et al. (2012b) V. Bandishti, I. Pananilath, and U. Bondhugula. 2012b.  Tiling stencil computations to maximize parallelism _(SC ’12)_. 1–11.
  * Basu et al. (2015) P. Basu, M. Hall, S. Williams, B. V. Straalen, L. Oliker, and P. Colella. 2015\.  Compiler-Directed Transformation for Higher-Order Stencils. In _2015 IEEE International Parallel and Distributed Processing Symposium_. 313–323.
  * Bielecki and Palkowski (2016) Wlodzimierz Bielecki and Marek Palkowski. 2016.  Tiling arbitrarily nested loops by means of the transitive closure of dependence graphs.  _AMCS : International Journal of Applied Mathematics and Computer Science_ 26, 4 (2016), 919–939.  https://doi.org/10.1515/amcs-2016-0065
  * Bondhugula et al. (2008) Uday Bondhugula, Albert Hartono, Jagannathan Ramanujam, and Ponnuswamy Sadayappan. 2008.  A practical automatic polyhedral parallelizer and locality optimizer. In _PLDI 08_.
  * Christen et al. ([n. d.]) Matthias Christen, Olaf Schenk, and Helmar Burkhart. [n. d.].  Patus: A code generation and autotuning framework for parallel iterative stencil computations on modern microarchitectures. In _IPDPS 2011_. IEEE.
  * Datta et al. (2009) Kaushik Datta, Shoaib Kamil, Samuel Williams, Leonid Oliker, John Shalf, and Katherine Yelick. 2009\.  Optimization and performance modeling of stencil computations on modern microprocessors.  _SIAM review_ 51, 1 (2009), 129–159.
  * Dursun et al. (2012) Hikmet Dursun, Manaschai Kunaseth, Ken-ichi Nomura, Jacqueline Chame, RobertF. Lucas, Chun Chen, Mary Hall, RajivK. Kalia, Aiichiro Nakano, and Priya Vashishta. 2012\.  Hierarchical parallelization and optimization of high-order stencil computations on multicore clusters.  _The Journal of Supercomputing_ 62, 2 (2012), 946–966.  https://doi.org/10.1007/s11227-012-0764-z
  * Dütsch et al. (2014) Fabian Dütsch, Karim Djelassi, Michael Haidl, and Sergei Gorlatch. 2014.  HLSF: A High-Level; C++-Based Framework for Stencil Computations on Accelerators _(WOSC ’14)_. 41–4.
  * Falke et al. (2013) Stephan Falke, Florian Merz, and Carsten Sinz. 2013\.  Extending the Theory of Arrays: memset, memcpy, and Beyond. In _Working Conference on Verified Software: Theories, Tools, and Experiments_. Springer, 108–128.
  * Frigo and Strumpen (2005) Matteo Frigo and Volker Strumpen. 2005.  Cache oblivious stencil computations _(ICS ’05)_. 361–366.
  * Grosser et al. (2014) Tobias Grosser, Sven Verdoolaege, Albert Cohen, and P. Sadayappan. 2014.  The Relation Between Diamond Tiling and Hexagonal Tiling.  _Parallel Processing Letters_ 24, 03 (2014).
  * Gysi et al. (2015) Tobias Gysi, Tobias Grosser, and Torsten Hoefler. 2015\.  Modesto: Data-centric analytic optimization of complex stencil programs on heterogeneous architectures. In _ICS 2015_. 177–186.
  * Hartono et al. (2010) A. Hartono, M. M. Baskaran, J. Ramanujam, and P. Sadayappan. 2010.  DynTile: Parametric tiled loop generation for parallel execution on multicore processors _(IPDPS ’10)_. 1–12.
  * Henretty et al. (2011) Tom Henretty, Kevin Stock, Louis-Noël Pouchet, Franz Franchetti, J Ramanujam, and P Sadayappan. 2011.  Data layout transformation for stencil computations on short-vector simd architectures. In _International Conference on Compiler Construction_. Springer, 225–245.
  * Henretty et al. (2013) Tom Henretty, Richard Veras, Franz Franchetti, Louis-Noël Pouchet, J. Ramanujam, and P. Sadayappan. 2013.  A Stencil Compiler for Short-vector SIMD Architectures _(ICS ’13)_. 13–24.
  * Hormati et al. (2010) Amir H Hormati, Yoonseo Choi, Mark Woh, Manjunath Kudlur, Rodric Rabbah, Trevor Mudge, and Scott Mahlke. 2010.  MacroSS: macro-SIMDization of streaming applications.  _ACM SIGARCH computer architecture news_ 38, 1 (2010), 285–296.
  * Iooss et al. (2015) Guillaume Iooss, Sanjay Rajopadhye, Christophe Alias, and Yun Zou. 2015.  _Mono-parametric Tiling is a Polyhedral Transformation_.  Research Report.
  * Jin et al. (2001) Guohua Jin, John Mellor-Crummey, and Robert Fowler. 2001\.  Increasing Temporal Locality with Skewing and Recursive Blocking _(SC ’01)_. 43–43.
  * Kamil et al. (2005) Shoaib Kamil, Parry Husbands, Leonid Oliker, John Shalf, and Katherine Yelick. 2005.  Impact of Modern Memory Subsystems on Cache Optimizations for Stencil Computations _(MSP ’05)_. 36–43.
  * Krishnamoorthy et al. (2007) Sriram Krishnamoorthy, Muthu Baskaran, Uday Bondhugula, J. Ramanujam, Atanas Rountev, and P Sadayappan. 2007.  Effective Automatic Parallelization of Stencil Computations _(PLDI ’07)_. 235–244.
  * Li et al. (2019) Kun Li, Honghui Shang, Yunquan Zhang, et al. 2019\.  OpenKMC: a KMC design for hundred-billion-atom simulation using millions of cores on Sunway Taihulight. In _SC’19_.
  * Luo et al. (2015) Yulong Luo, Guangming Tan, Zeyao Mo, and Ninghui Sun. 2015\.  FAST: A Fast Stencil Autotuning Framework Based On An Optimal-solution Space Model _(ICS ’15)_. 187–196.
  * Nguyen et al. (2010) A. Nguyen, N. Satish, J. Chhugani, C. Kim, and P. Dubey. 2010.  3.5-D Blocking Optimization for Stencil Computations on Modern CPUs and GPUs _(SC ’10)_. 1–13.
  * Rawat et al. ([n. d.]) Prashant Singh Rawat, Fabrice Rastello, Aravind Sukumaran-Rajam, Louis-Noël Pouchet, Atanas Rountev, and P Sadayappan. [n. d.].  Register optimizations for stencils on GPUs. In _PPOPP’2018_. 168–182.
  * Rawat et al. (2018) P. S. Rawat, A. Sukumaran-Rajam, A. Rountev, F. Rastello, L. Pouchet, and P. Sadayappan. 2018.  Associative Instruction Reordering to Alleviate Register Pressure. In _SC’18_. 590–602.
  * Rocha et al. (2017) Rodrigo C. O. Rocha, Alyson D. Pereira, Luiz Ramos, and Luís F. W. Góes. 2017.  TOAST: Automatic tiling for iterative stencil computations on GPUs.  _Concurrency and Computation: Practice and Experience_ (2017).
  * Springer et al. (2017) Paul Springer, Jeff R Hammond, and Paolo Bientinesi. 2017\.  TTC: A high-performance compiler for tensor transpositions.  _ACM Transactions on Mathematical Software (TOMS)_ 44, 2 (2017), 1–21.
  * Tang et al. (2011) Yuan Tang, Rezaul Alam Chowdhury, Bradley C. Kuszmaul, Chi-Keung Luk, and Charles E. Leiserson. 2011\.  The Pochoir Stencil Compiler. In _SPAA’11_. ACM, New York, NY, USA, 12.  https://doi.org/10.1145/1989493.1989508
  * Wonnacott and Strout (2013) David G Wonnacott and Michelle Mills Strout. 2013.  On the scalability of loop tiling techniques.  _IMPACT 2013_ (2013).
  * Yount (2015) Charles Yount. 2015\.  Vector Folding: improving stencil performance via multi-dimensional SIMD-vector representation. In _2015 IEEE 17th International Conference on High Performance Computing and Communications, 2015 IEEE 7th International Symposium on Cyberspace Safety and Security, and 2015 IEEE 12th International Conference on Embedded Software and Systems_. IEEE, 865–870.
  * Yount et al. (2016) Charles Yount, Josh Tobin, Alexander Breuer, and Alejandro Duran. 2016.  YASK—Yet another stencil kernel: A framework for HPC stencil code-generation and tuning. In _WOLFHPC 2016_. IEEE, 30–39.
  * Yuan et al. (2019) Liang Yuan, Shan Huang, Yunquan Zhang, and Hang Cao. 2019\.  Tessellating Star Stencils. In _Proceedings of the 48th International Conference on Parallel Processing_. 1–10.
  * Yuan et al. (2017) Liang Yuan, Yunquan Zhang, Peng Guo, and Shan Huang. 2017\.  Tessellating Stencils. In _SC’17_. ACM, New York, NY, USA, Article Article 49, 13 pages.  https://doi.org/10.1145/3126908.3126920
  * Zekri (2014) Ahmed Sherif Zekri. 2014\.  Enhancing the matrix transpose operation using intel avx instruction set extension.  _International Journal of Computer Science & Information Technology_ 6, 3 (2014), 67.
  * Zhao et al. (2019) Tuowen Zhao, Protonu Basu, Samuel Williams, Mary Hall, and Hans Johansen. 2019.  Exploiting reuse and vectorization in blocked stencil computations on CPUs and GPUs. In _SC’19_. 1–44.
  * Zumbusch (2012) Gerhard Zumbusch. 2012\.  _Vectorized Higher Order Finite Difference Kernels_.  343–357.  https://doi.org/10.1007/978-3-642-36803-5_25