An Extension of OpenACC Directives for
Out-of-Core Stencil Computation with Temporal
Blocking
Nobuhiro Miki, Fumihiko Ino, and Kenichi Hagihara
Graduate School of Information Science and Technology, Osaka University
1-5 Yamadaoka, Suita, Osaka 565-0871, Japan
Email: {n-miki, ino}@ist.osaka-u.ac.jp
Abstract—In this paper, aiming at realizing directive-based
temporal blocking for out-of-core stencil computation, we present
an extension of OpenACC directives and a source-to-source
translator capable of accelerating out-of-core stencil computation
on a graphics processing unit (GPU). Out-of-core stencil computation here deals with large data that cannot be entirely storedin GPU memory. Given an OpenACC-like code, the proposedtranslator generates an OpenACC code such that it decomposeslarge data into smaller chunks, which are then processed in apipelined manner to hide the data transfer overhead needed forexchanging chunks between the GPU memory and CPU memory.
Furthermore, the generated code is optimized with a temporal
blocking technique to minimize the amount of CPU-GPU datatransfer. In experiments, we apply the proposed translator tothree stencil computation codes. The out-of-core performance ona Tesla K40 GPU reaches 73.4 GFLOPS, which is only 13%
lower than the in-core performance. Therefore, we think thatour directive-based approach is useful for facilitating out-of-core
stencil computation on a GPU.

I. I NTRODUCTION
Stencil computation is a class of iterative computation thatfrequently appears in scientific applications of a wide rangeof fields, including computational fluid dynamics [24], computational electromagnetics [1], image processing [5], and soon. This iterative computation typically updates array elements
according to a fixed pattern, called stencil. For example,
finite-difference methods have been widely used for solving
partial differential equations that describe the time evolutionof variables. In general, stencil applications are regarded as
memory-intensive applications [6], [7], [19]. Therefore, many
stencil codes have been implemented on accelerator devices,such as the graphics processing unit (GPU) [16] and Xeon
Phi [8], which typically provide a 5 times higher memory
bandwidth than the CPU.
Accelerator codes are usually developed with a uniqueprogramming language to maximize the performance on accelerator devices. For example, the compute unified device
architecture (CUDA) [16] is widely used for NVIDIA GPUs.
The CUDA code consists of the host code and the device code,
each running on the CPU and the GPU, respectively. The key
for achieving acceleration is to locate the performance bottleneck of the sequential code so that the bottleneck operations
can be offloaded to the accelerator device with device-specific
WACCPD16; Salt Lake City, Utah, USA; November 2016
c978-1-5090-6152-5/16/$31.00 ⃝2016
IEEE

optimization techniques. However, this offloading process istime consuming because code and data structures usually haveto be entirely rewritten to adapt themselves to the device architecture. Moreover, application developers are again enforced torewrite their device-specific code when a new architecture will
be released in the future. From this point of view, scientificcodes should retain performance portability such that the codesachieve high performance on different architectures.
Such tremendous efforts can be minimized by directivebased programming (see Fig. 1), namely an emerging approachfor developing parallel codes on an accelerator device. For
example, OpenACC [20] provides a collection of compiler
directives that are useful for offloading bottleneck workloads
to an accelerator device. In other words, directives represent architecture-specific description, so that the architecturedependent code is clearly separated from the generic code,
which expresses the essentials of computation. Thus, OpenACC significantly lowers the barrier to accelerated computing. However, typical OpenACC-based implementations
assume that the entire data are stored in device memory,which compromises the directive-based approach. Because
the device memory generally has a smaller capacity than thehost memory, OpenACC-based implementations usually fail tosolve a large problem that has been processed by CPU-based
implementations.
To relax this restriction, we developed a directive-based
framework, named pipelined accelerator (PACC) [10], capableof pipelined execution of large-scale stencil computation on aGPU. According to PACC directives, our translator generates
an OpenACC code such that large data are automaticallydecomposed into smaller chunks, which are then processedin a pipelined manner to overlap CPU-GPU data transfer withkernel execution. Thus, PACC directives are useful for solving
a problem of the same size as that processed by CPU-based
implementations. However, data decomposition increased theamount of CPU-GPU data transfer, which limited the performance of stencil applications, particularly those solving time
evolution problems.
In this work, aiming at reducing the amount of data transfer
between the CPU and GPU, we extend PACC directives [10]

# pragma pacc …
for (…) {}PACC codefor (…) {}Sequentialcode# pragma acc …
for (…) {}Translator# pragma acc …
for (…) {for (…){ … }}OpenACC
compiler

OpenACC code

Out-of-core
GPU code

OpenACC
compilerbin

OpenACC code

# pragma omp …
for (…) {}bin

In-coreGPU code
OpenMP
compiler

OpenMP codebin

In-coreCPU code

Fig. 1. Directive-based programming flow.

such that both data decomposition and temporal blocking
can be automatically applied to stencil computation on a
CUDA-compatible GPU. Given a sequential C code with
PACC directives, our source-to-source translator generates
an OpenACC code that can efficiently accelerate large-scale
stencil computation with temporal blocking. The generated
code can be rapidly processed on a GPU with several execution
parameters such as the number of data decompositions. We
apply the presented translator to three stencil computation
codes [4] to investigate the impact of execution parameters.
The paper is structured as follows. Section II introduces
related work regarding on acceleration of stencil computation.
Section III presents an overview of temporal blocking and thensummarizes how this optimization technique can be applied to
an OpenACC code. Section IV describes the proposed PACC
directives and translator. Section V shows experimental results.
Finally, Section VI concludes the paper with future work.
II. R ELATED W ORK
Maruyama et al. [12] presented a programming framework,
named Physis, which provided a domain specific language(DSL) based solution to stencil computation. Given a DSLbased code, Physis generates a CUDA code for multi-node
systems, where computing nodes communicate each other with
the Message Passing Interface (MPI) standard [13]. Further,
they extended Physis to automatically apply temporal blocking
to the generated code [9]. This DSL-based approach is similar to our directive-based approach because both approaches
facilitate acceleration of stencil computation. However, thedirective-based approach is more useful for achieving thiscommon goal with less efforts because directive-based approach keeps the original structure of the sequential code; there
is no need to adapt the code to the target DSL.
Endo et al. proposed a run-time library, called hybrid
hierarchical runtime (HHRT) [2], which allows CUDA+MPI
codes to deal with out-of-core data with temporal blocking[3]. The HHRT library virtualized the GPU memory with
automated swapping to the CPU memory. Because the HHRT
library requires a CUDA code as its input, more efforts
are necessary to achieve parallelization on the GPU. Our

directive-based approach avoids this time-consuming process
but execution efficiency can be degraded due to higher-levelof code description. For example, low-level optimization suchas intrinsic functions [16] is not explicitly available from an
OpenACC code.
XcalableACC [15], a hybrid model of the OpenACC and apartitioned global address space (PGAS) language [14], [25],realizes directive-based programming for multi-node accelerator systems. Similar to our approach, their approach is usefulfor implementing a highly-efficient portable code by adding
directives to a sequential code. Although the XcalableACC
code efficiently runs on multi-node systems with minimum
programming efforts, code structure must be modified to use
temporal blocking, which is the key technique to minimize theamount of data transfer between the host and device.
OpenMP [21] provides a collection of directives that are
useful for implementing multithreaded parallel applications onshared memory architectures such as multi-core CPUs. The
latest OpenMP 4.5 supports accelerators with OpenACC-likedirectives. However, the host data cannot be partially mapped
to the device data. Therefore, large problems cannot be solved
due to device memory exhaustion.
Finally, NVIDIA has recently announced their latest architecture, called Pascal (P100) [18], which provides a single,
unified virtual address space for CPU and GPU memory. This
capability is useful to realize out-of-core stencil computation
with a low effort. However, pipelined execution requires
prefetching of chunks to overlap kernel execution with CPUGPU data transfer. It is not clear whether such a software
pipeline can be easily implemented with maintaining theoriginal structure of the sequential code; the P100 architecture
is not available at this time.
III. S TENCIL C OMPUTATION AND T EMPORAL B LOCKING
Temporal blocking is a cache optimization technique forstencil applications that solve time evolution problems. This
technique saves the memory bandwidth by performing reuseof on-cache data. To do so, the computational domain isdecomposed into smaller blocks such that each block evolves
k contiguous time steps at a time before proceeding to the next
block. Hereafter, we call parameter k as the blocking factor.
There are several approaches for realizing temporal blocking.
In this paper, we consider an overlapped tiling approach for
its simplicity. Other temporal blocking approaches such as
wavefront temporal blocking and diamond tiling are useful for
eliminating halos mentioned below (i.e., overheads in terms ofcomputation at block boundaries).
Figure 2 shows an example code of four-point stencil
computation. In this example, a four-point stencil is applied
to X × Y elements stored in a two-dimensional (2-D) array(see Fig. 3); updating an element refers to itself together with
its four neighbors. To realize block-based time evolution, theoriginal t loop in Fig. 2 must be reorganized into double
nested loops; the inner loop evolves k time steps within a
block and the outer loop evolves every k time step for blocks.

1 for (t=0; t<T; t++) { // time evolution2for (x=1; x<X-1; x++) {3for (y=1; y<Y-1; y++)4
q[x][y] = (p[x-1][y] + p[x+1][y] + p[x][y-1] + p[x][yY+1]) * 0.25;

Decomposition

5}6
swap(p, q);
7 }Halo region rk
Block size b
Halo region rk

XChunk
Original data

Fig. 2. Pseudocode of four-point stencil computation.

Fig. 4. A 1-D block decomposition scheme with halo region. Given a stencilof (2r + 1) × (2r + 1) elements, each block requires halos of size rk × Y
to evolve k time steps for all elements within the block.

Y

X

Host
(a)

(b)

buf_p[0](c)Fig. 3. Overview of stencil computation. (a) X × Y elements in thecomputational domain are updated with (b) a four-point stencil stored in a
3 × 3 region. (c) As compared with the original stencil in (b), applying
temporal blocking refers more neighboring elements to compute the target
element (k = 2, in this example).

DeviceStream #0

buf_p[0]Stream #1p

buf_p[1]buf_p[1]Stream #2
Data copy

buf_p[2]buf_p[2]Data transfer

Thus, temporal blocking exploits temporal locality to takeadvantage of caches, which have a low latency but witha limited capacity compared to the CPU memory. Such a
tradeoff relation can be found between the CPU memory andGPU memory. Consequently, temporal blocking is frequently
used for GPU-based implementations to reuse data blocks thathave been transferred to the GPU memory. Evolving k time
steps within the transferred data blocks reduces the amount of
CPU-GPU data transfer to 1/k.
Notice here that multiple updates within a block require a
halo region around the block (see Fig. 4) because updating
an element refers to its surrounding neighbors. Without thishalo region, blocks cannot be processed independently, andthereby, preventing full parallelization. Suppose that a crossstencil is stored in a (2r+1)×(2r+1) region, i.e., updating anelement refers to itself and its r neighbors in up/down/left/right
directions: r = 1 for a four-point stencil, as depicted inFig. 3(b). In this case, k updates for an element refers to its(2rk + 1) × (2rk + 1) region (Fig. 3(c)). This means that temporal blocking increases the amount of computation because
halo regions of adjacent blocks overlap each other. Thus, there
is a tradeoff relation between the amount of computation andthe degree of available parallelism. Consequently, the blocking
factor k must be optimized to maximize the performance gainof temporal blocking.
Hereafter, we use the term chunk to denote the region that
contains a block and its overlapping halo region (Fig. 4).
A. OpenACC-based implementation
In this section, we present how data decomposition andtemporal blocking can be applied to an OpenACC-based
stencil code. In other words, the OpenACC code presented
here is the output of our translator, which requires a PACC
code as its input.

Kernel execution

Fig. 5. A copy-based scheme for data decomposition. Host and device buffers
are allocated for copying blocks and their halos from the original array.

1) Data decomposition: Data decomposition cannot berealized by simply adding OpenACC directives into the sequential code. Given a large array that exceeds the capacityof the device memory, a different array must be allocated in
the host memory, as depicted in Fig. 5. This is due to theassumption of the OpenACC specification, which allocates thesame variables on both the host memory and device memory.
In this figure, large array p can be successfully processed on
the GPU by copying chunks (i.e., blocks and their halos) todifferent buffers buf_p; memory exhaustion can be avoided
if the total buffer size is smaller than the capacity of device
memory. One drawback of this copy-based scheme is the copy
overhead incurred on the host. However, this overhead can be
hidden by pipelined execution, which overlaps the overhead
with data transfer and kernel execution.
An alternative solution for avoiding device memory exhaustion is a map-based scheme, which maps the allocated device
buffers to the original host array. This direct mapping allows
GPU threads to update values in the original array, so thatthe host buffer (i.e., the copy overhead) can be eliminated.
However, we decided to use the copy-based scheme because
the map-based scheme failed to asynchronously transfer dataon our experimental environment (see Section V); pipelined
execution was not available with this scheme. Another drawback of the map-based scheme is that it may fail to realize multidimensional decomposition, which produces many mapping
points; direct mapping is mainly designed for a contiguous
memory region.
Figure 6 shows a pseudocode of OpenACC-based out-of-

1 Set blocking factor k, number d of decompositions, and number num_queue of queues from environmental variables;
2 Compute block size b from d and data size X;
3 Allocate buf_p[0], buf_p[1], ..., buf_p[num_queue-1] on host memory;
4 #pragma acc create (buf_p[0:num_queue][0:b+2*r*k], ...)
5 Allocate g[0], g[1], ..., g[num_queue-1] on host memory to store completed chunk numbers;
6
7 for (t=0; t<T; t+=k) { // outer loop for time evolution8
for (c=0; c<d; i = (i+1) % num_queue) { // for each chunk910
if (!acc_async_test(i)) // Select i as the ID of an idle queue
11
continue;
1213
if (g[i] != NONE) {14
Copy the g[i]-th chunk from buf_p[i] to p;
15g[i] = NONE;
16
}
17
18
Copy the c-th chunk from p to buf_p[i];19#pragma acc update device (buf_p[i:1][0:b+2*r*k], ...) async (i)2021for (s=0; s<k; s++) { // inner loop for time evolution
22#pragma acc kernels present (buf_p[i:1][0:b+2*r*k], ...) async(i)23{
24
offset = r*(s+1);
25
xsize = b+2*r*(k-1-s);
26
27#pragma acc loop independent28
for (x=offset; x<offset+xsize; x++)29
#pragma acc loop independent30
for (y=1; y<Y-1; y++)31
#pragma acc loop independent32
for (z=1; z<Z-1; z++)33
buf_q[i][x*Y*Z+y*Z+z] += buf_p[i][(x+1)*Y*Z+y*Z+z] + ...;
34}
35
#pragma acc kernels present (buf_p[i:1][0:b+2*r*k], ...) async(i)36{37#pragma acc loop independent38
for (x=offset; x<offset+xsize; x++)39
#pragma acc loop independent40
for (y=1; y<Y-1; y++)41
#pragma acc loop independent42
for (z=1; z<Z-1; z++)43buf_p[i][x*Y*Z+y*Z+z] = buf_q[i][(x+1)*Y*Z+y*Z+z];
44
}
45
}
46
47
#pragma acc update host (buf_p[i:1][0:b+2*r*k], ...) async (i)48g[i] = c;
49
c++;
50
}
51 }
52 for (i=0; i<num_queue; i++) {// Clean up all queues
53#pragma acc wait(i)54
if (g[i] != NONE)
55
Copy the g[i]-th chunk from buf_p[i] to p;
56 }Fig. 6. Pseudocode of out-of-core stencil computation with temporal blocking. This pseudocode applies a cross stencil of (2r + 1) × (2r + 1) elements tothe computational domain of (X − 2r) × Y × Z elements. A 1-D block decomposition scheme is applied to the computational domain, so that the data isdecomposed into d blocks of size (X − 2r)/d × Y × Z.

core stencil computation that deploys a 1-D block decomposition scheme and temporal blocking. At line 3, num_queue
buffers, buf_p, are allocated for the original array p, where
num_queue is the number of queues used for pipelining.
Similarly, device buffers with the same variable names areallocated by the create clause at line 4. Note that multiple
buffers are necessary to realize efficient pipelining; a single
buffer causes data dependence between the kernel execution
and host-device data transfer, which avoids overlapped execution. Note also that these buffers are reused to minimize
the allocation overhead during program execution; the buffers
are allocated once at the beginning of an execution. At line 5,
num_queue buffers, g, are allocated to store chunk numbers.

After this allocation, chunks are copied from the originalarray p to the host buffers buf_p, which are then transferredto the device buffers buf_p by using the update clauseat line 19. Using these buffers, a kernel function is invokedat line 22,35 to update elements in the chunks for k time
steps. To do this, a kernels construct is deployed to specify
a code block to be offloaded from the host to the device.
In addition, the present clause at line 22,35 indicates thatchunks to be accessed have already been sent to the device
buffers buf_p, so that additional data transfer can be avoided
at kernel invocation. The updated elements are then transferredback to the host buffers at line 47. Host buffers are copied tothe original array p at line 14 when the same queue is selected1#pragma pacc init
next time.
2
Notice that the code modification mentioned above is nec- 3 #pragma pacc pipeline targetinout(p,q) size([0:X][0:Y])
halo([1:1][1:1]) async
essary to avoid device memory exhaustion, which results in an 4 for (t=0;t<T; t++){#pragma pacc loop dim(2)execution failure. Without this code modification, the original 5for (x=1; x<X-1; x++)code fails to run if array p is too large to fit in the device 67#pragma pacc loop dim(1)for (y=1; y<Y-1; y++)memory. However, as shown in Fig. 6, this modification 8
q[x][y] = (p[x-1][y] + p[x+1][y] + p[x][y-1] + p[xreorganizes the loop structure of the original code (Fig. 2), 9][y+1]) * 0.25;
diminishing the benefits of directives. That is, the modified10#pragma pacc loop dim(2)code degrades the performance portability because it leads1112
for (x=1; x<X-1; x++)#pragma pacc loop dim(1)to inefficient run on other machines equipped with a large-13
14
for (y=1; y<Y-1; y++)capacity memory.
15
p[x][y] = q[x][y];2) Temporal blocking: Efficient temporal blocking can be16 }
implemented by realizing the following points.
Fig. 7. An example of a PACC code.
• Block-based time evolution. The time evolution loop must
be separated into two loops such that one is responsiblefor intra-block and the other is responsible for inter-block,as shown in Fig. 6; the outer t loop at line 7 is responsible be overlapped with the data transfer if chunks are assigned
for processing blocks at every k time steps while the to different queues. A full overlap can be established with
inner s loop at line 21 is responsible for processing a multiple CPU threads, but we leave this issue for a future
block for consecutive k time steps. To facilitate automated work.
parallelization, x, y, and z loops in the code block have
IV. P IPELINED ACCELERATOR (PACC)
a loop construct with an independent clause, which
The proposed PACC directives allow application developers
notifies the OpenACC compiler that iterations can beto specify key information, such as the stencil size and the data
efficiently executed without synchronization.
to be decomposed, which are required to automate applying• Pipelined execution. Asynchronous APIs are deployed
to realize software-based pipelining that overlaps kernel data decomposition and temporal blocking. Moreover, several
execution with host-device data transfer. In Fig. 6, a execution parameters, such as the blocking factor k andchunk is assigned to one of the num_queue asyn- the number d of decompositions, can be flexibly specified
chronous queues. At line 10, a queue is tested whether by environmental variables at the beginning of a program
all associated operations of the queue have completed or execution.
not using an acc_async_test API. That is, chunks
Our translator assumes that the target stencil code satisfies
are transferred at line 19 with an async clause and the following constraints.
its argument i, which specifies the queue ID to be• Constraints on data decomposition. Application developused. As mentioned before, the host and device buffers
ers are prohibited to manually decompose data in thebuf_p[i] are dedicated to the i-th queue for realizing
PACC code. The translator also assumes that the data
concurrent execution. Thus, different queues have no data
are small enough to be stored in the host memory.
dependence between their tasks.
• Constraints on temporal blocking. The number T of total
With respect to the tasks to be queued, each queue istime steps must be fixed at the beginning of a program
responsible for processing the following steps in order.
execution. Accordingly, the target program is prohibitedto use a while loop to terminate its execution according
1) Data copy step. The host copies a chunk from the original
to a user-defined threshold on the error. Furthermore, thearray to the host buffer (line 18).
device evolves data with k time steps at a time, so that the2) Data transfer step (host to device). The host transfers theintermediate values between every k time steps cannot be
copied chunk from the host buffer to the device buffer
accessed from the host.
(line 19).
3) Kernel execution step. The device iteratively executes thekernel to evolve the transferred chunk for k time steps A. PACC directive
The proposed PACC extends OpenACC [20] with three con(lines 21–45).
4) Data transfer step (device to host). The updated chunk structs: the init, pipeline, and loop constructs. Figure 7
is transferred from the device buffer to the host buffer shows an example of a PACC code that implements out-ofcore stencil computation with temporal blocking. Similar to(line 47).
5) Data copy step. The host copies the updated chunk to the OpenACC directives, a PACC directive starts with #pragma
pacc. The extended directives are as follows.
original array (line 14).
• The init construct. This construct allocates host and deNote that the pseudocode code of Fig. 6 cannot overlap thevice buffers for realizing data decomposition (see Section
first step with the last step because it uses a single CPU
III-A). Consequently, the init construct must be placed
thread to copy data. However, the kernel execution step canCode rewrite rulesfor (t=0;t<T;t++)

ROSE compiler framework# pragma tpacc …
for (…) {}C / C++
Frontend

PACC codefor (x=1;x<X-1;x++)

Abstract
syntax
tree
(AST)

Code
generatorp = q# pragma acc …
for (…) {for( … ){ … }}for (y=1;y<Y-1; y++)OpenACC codeq[x][y] = p[x][y] + …

(a)Fig. 8. An overview of the proposed translator.

before the pipeline construct and the loop constructmentioned below.
• The pipeline construct. The pipeline constructspecifies the code block to be processed in a pipeline.
In Fig. 7, the pipeline construct is applied to thefor loop at line 4, which is responsible for time evolution. This construct is similar to the data construct of
OpenACC and can have additional clauses such as thetargetin, targetinout, size, and halo clauses.
The targetin and targetinout clauses define readonly and read/write variables, respectively. The size
clause defines an array range specification with startand length for each dimension. For example, Fig. 7
specifies size([0:X][0:Y]) at line 3 to update all
elements in an array of size X × Y . Finally, the haloclause defines the stencil size for each dimension. In
Fig. 7, halo([1:1][1:1]) defines a four-point stencil, which accesses left/right/up/down neighbors to update
an element. Finally, the async clause declares that thecode block must be asynchronously processed to realize
pipelined execution.
• The loop construct. The PACC loop construct is anextension of the OpenACC loop construct. The extended
construct can have an extended clause, dim, whichassociates the loop control variable with the dimensionof the array data. For example, the dim(2) clause atline 5 indicates that the loop control variable x at line 6
corresponds to the second dimension of the array data;
we assume here that the first (last) dimension has theminimum (maximum, respectively) stride between consecutive elements. This association is used for adjusting
the loops for decomposed data. Our translator currently
deploys a 1-D block scheme that decomposes the arraydata in terms of the last dimension. In Fig. 7, the arraydata will be decomposed along the x axis.
B. PACC translator
Our PACC translator was implemented using the ROSE
compiler infrastructure [11], [22], which provides a C/C++
frontend to generate an abstract syntax tree (AST) of the
input code (see Fig. 8). The generated AST is then traversed
to detect PACC directives, i.e., the nodes with the pacc
attribute. During this traverse, detected directives are parsed
to retrieve the key information given by succeeding clauses,
such as stencil size, the variables to be decomposed, andtheir data size. The detected nodes are marked explicitly forInter-block evolutionfor (t=0;t<T;t+=k)for each chunkfor (c=0;c<d;c++)Queueselection

Chunk copyp->buf_p[i]#pragma acc
updatedevicefor (s=0;s<k; s++)#pragma acc
update hostoffset=…
xsize=…

for (x=offset;x<offset+xsize;x++)buf_q[i]=buf_p[i]Chunk copy
buf_p[i]->pIntra-block evolution#pragma acc
kernels …
async

Chunk updatefor(y=1; y<Y-1; y++)buf_q[i][x*Y+y]=…

(b)
Fig. 9. An example of AST transformation. (a) The presented translator
generates an AST from a PACC code. (b) The translator then applies coderewrite rules to obtain an OpenACC code that can deal with large data withdata decomposition and temporal blocking.

modification, so that code rewrite rules are applied to themin the next traversal. Finally, the modified AST is given to acode generator to obtain the pipelined OpenACC code.
The code rewrite rules can be summarized as follows.
1) The rule for the init construct. The AST node that
corresponds to the init construct is replaced with AST
nodes that are responsible for (1) obtaining the blocking
factor k and the number d of data decompositions via
environmental variables, (2) computing the block size bfrom d and the array size X to be decomposed, and (3)allocating host and device buffers.
2) The rule for the pipeline construct. The original for
loop is reorganized into double nested for loops torealize temporal blocking. To achieve this, as shown in
Fig. 9, the code rewrite rule replaces the AST node that
corresponds to the for loop of time evolution with AST
nodes that correspond to (1) the for loop of inter-block
evolution, (2) that of per-chunk operations, and (3) thatof intra-block evolution. Furthermore, several AST nodes
are added as children of the second AST node (i.e., perchunk operations) to select an idle queue, copy chunks on
the host buffers, and exchange them between the host anddevice buffers. Finally, memory accesses to the originalarrays are replaced with those to the buffers by updating
the attributes of AST nodes appropriately.
3) The rule for the loop construct. This rule updates thetransformed double nested loops with new initializationand condition to adapt the loop structure for block-based
evolution. To do this, the rule first locates the for loopthat must be updated due to data decomposition. For

TABLE IS TENCIL COMPUTATION CODES USED FOR EXPERIMENTS . N OTATIONS T , f , l, s, AND G REPRESENT THE NUMBER OF TIME STEPS , THAT OF FLOATING
POINT OPERATIONS PER ELEMENT, THE AMOUNT OF WRITES PER ELEMENT, THE AMOUNT OF READS PER ELEMENT, AND THE ARITHMETIC INTENSITY[23], RESPECTIVELY. T HE ARITHMETIC INTENSITY G IS GIVEN BY G = f /(l + s).
Code
Number of arrays
Array size
Data size (GB) Stencil T (step) f (FLOP) l (B) s (B) G (FLOP/B)
Jacobi2
48, 000 × 48, 000
18.4
4 point 2048
4
16
4
0.20
Himeno14
512 × 512 × 1024
15.0
18 point
256
34
128
4
0.26
CIP
8
22, 000 × 22, 000
15.5
9 point
256
91
120 120.69

example, the for loop at line 6 of Fig. 7 is selected
for updating its initialization and condition because thisloop is associated with a dim(2) clause that has thelargest value (2) as the argument of dim clauses. The
new initialization and condition is appropriately given byvariables offset and xsize, as shown in Fig. 6.
V. E XPERIMENTAL R ESULTS
We added PACC directives to three stencil computation
codes and evaluated their performance on an experimental
machine. Our experimental machine had an Intel Xeon E52680 v2 processor, 512 GB main memory, and a Tesla K40
GPU with 12 GB device memory. The Tesla K40 GPU had
two direct memory access (DMA) engines, so that data transferfrom the host to device was overlapped with that of the
opposite direction. We used the PGI Compiler 15.5 [17],CUDA 7.0 [16], and Ubuntu 15.3.
The stencil computation codes used for experiments were
the Jacobi method, Himeno benchmark [4], and constraintinterpolation profile (CIP) method [26] summarized in Table I.
The Jacobi method is an iterative solver for a system of
linear equations. The Himeno benchmark is a linear solver
for 3-D pressure Poisson equations. The CIP method is asolver for hyperbolic partial differential equations. All thethree experimental codes processed at least 15 GB of data,
which could not be stored entirely on the device memory.
A. Performance analysis
We first manually implemented PACC codes by adding
PACC directives appropriately to the sequential codes. The
PACC codes were then given to the proposed translator toautomatically generate out-of-core OpenACC codes. Using thegenerated codes, we measured their effective performance toinvestigate the impact of execution parameters such as theblocking factor k and the block size b.
1) Jacobi method: We solved a problem size of X × Y =48, 000 × 48, 000 elements, which consumed 18.4 GB ofmemory space, with the total time steps of T = 2048.
Figure 10(a) shows the effective performance E given by
E = 4(X − 2)(Y − 2)T /t, where t is the execution time,
which includes the data transfer time between the CPU andGPU.
In Fig. 10(a), the maximum performance of the Jacobi
method was 28.5 GFLOPS, obtained when the block size andthe blocking factor were b = 8000 and k = 32, respectively.
This figure also shows that the performance was determinedby the blocking factor k rather than the block size b.

We then measured the effective performance with varying
the blocking factor k at fixed block size b = 8000 (Fig. 11(a)).
The effective performance did not monotonically increase withthe blocking factor k due to the tradeoff relation mentioned
in Section III. To investigate this behavior, we analyzed thebreakdown of execution time, which we show in Fig. 12(a).
We found that the copy overhead determined the entire performance when k < 32. This overhead was inversely proportional
to k, so that the effective performance increased with kwhen k < 32. In contrast, the effective performance slightly
decreased as we increased k from 32. This performancedegradation was caused by temporal blocking, which increased
kernel execution time due to redundant computation. Thus, thebest tradeoff point was obtained when k = 32, where the datatransfer and copy overheads were fully overlapped with kernel
execution time.
2) Himeno benchmark: We solved the problem size XL(512 × 512 × 1024, 15 GB) with T = 256 time steps. Given
3-D data of X × Y × Z elements, the effective performance
E is given by E = 34(X − 2)(Y − 2)(Z − 2)T /t.
In Fig. 10(b), the effective performance of the Himeno
benchmark monotonically increased with the blocking factor
k. As a result, the maximum performance of 37.5 GFLOPS
was obtained when b = 102 and k = 16, where the buffer
size was maximized; when k > 16, we failed to executethe PACC code due to device memory exhaustion. A larger
memory capacity was required to find the best blocking factor
that might be obtained when k > 16.
Thus, the tradeoff point was not clearly observed for the
Himeno benchmark because we failed to increase the blocking
factor k from 16. In other words, device memory exhaustion occurred though the data were decomposed for saving
memory consumption. This issue can be resolved by realizing
multidimensional decomposition. With our 1-D decomposition
scheme, a chunk consists of (b+2rk)×Y ×Z elements. Therefore, the amount of memory consumption linearly increased
with Y and Z, which restricted the blocking factor k such thatk ≤ 16. Therefore, a multidimensional decomposition scheme
is necessary for our experimental machine to successfully run
the benchmark with k > 16.
As for multidimensional decomposition, data pack andunpack procedures are required to retrieve a small multidimensional array from a large multidimensional array. Therefore,
the PACC translator must be extended such that it can (1)pack several data segments into a host buffer and (2) unpacka device buffer onto the original array. Furthermore, memoryreferences in the kernel must be updated to access the packed

20
15
10
5
0

k=512

k=256

k=128

k=64

k=32

k=16

k=8

k=4

2000

4000
6000
Block size b

40
35
30
25
20
15
10
5
0

8000

k=16
k=8
k=4
k=2
k=1

Effective performance (GFLOPS)

25

0

80

45

Effective performance (GFLOPS)

Effective performance (GFLOPS)

30

60
50
40
30
k=256
k=32
k=4

20
10
0

0 10 20 30 40 50 60 70 80 90 100110
Block size b

(a)

70

0

1000

k=128
k=16

2000 3000
Block size b

(b)

k=64
k=8

4000

5000

(c)

Fig. 10. Effective performance with different block size b and blocking factor k. Results for (a) Jacobi method, (b) Himeno benchmark, and (c) CIP method.

1200

24

1000

20

800

16

600

12

400

8

200

4

0

0

8

16

32 64 128 256 512
Blocking factor k(a)Effective Performance
40

400

35

350

70

300

60

250

50

200

40

150

30

100

20

50

10

30
25
20
15
10
5
1

2
4
8
Blocking factor k

(b)

16

0

0

Execution Time48

Effective Performance80

16 32 64 128 256
Blocking factor k

0

Effective performance (GFLOPS)

28

Execution Time

Execution time (s)

32

1400

Execution time (s)

1600

1200
1100
1000
900
800
700
600
500
400
300
200
100
0

Effective performance (GFLOPS)

Effective Performance
36

Execution time (s)Execution Time

Effective performance (GFLOPS)

1800

(c)

Fig. 11. Effective performance with different blocking factor k for fixed block size b. Results for (a) Jacobi method, (b) Himeno benchmark, and (c) CIP
method.

small array accordingly.
Figure 12(b) shows the breakdown of execution time. We
found that the data copy time spent on the host was almost
the same as the kernel execution time when k = 16. Consequently, the blocking factor of k = 16 seems to be the best
configuration for the Himeno benchmark.
3) CIP method: We solved a problem size of 22, 000 ×
22, 000 elements with T = 256 time steps. Given 2-D dataof X × Y elements, the effective performance E is given by
E = 91(X − 2)(Y − 2)T /t.
In Fig. 10, the effective performance reached 73.4 GFLOPS
when b = 2750 and k = 16. As compared with theJacobi method and the Himeno benchmark, the CIP method
maximized its performance with relatively a small blocking
factor, k = 8, due to its high arithmetic intensity [23] (seeTable I). That is, the ratio of the kernel execution time over the
CPU-GPU data transfer time was relatively high, and thereby,
the impact of temporal blocking was rapidly maximized withlow k. Therefore, temporal blocking failed to demonstrate asignificant improvement when k ≥ 8.
Figures 11(c) and 12(c) show the effective performance withdifferent blocking factor k and the breakdown of execution
time, respectively. As we estimated above, the best tradeoff

point was found at k = 8 in Fig. 11(c), where the copy
overhead was close to the kernel execution time, as shownin Fig. 12(c).

B. Comparison with in-core implementation
Finally, we compared our out-of-core implementation withan in-core implementation that processed small data (Table II).
The in-core versions of the Jacobi, Himeno, and CIP methods
solved problem sizes of 24, 000×24, 000 (4.6 GB), 128×128×
256 (0.2 GB), 16, 000 × 16, 000 (8.2 GB), respectively. Eachin-core data was iteratively updated with the same number oftime steps as the out-of-core implementation. Note here thatthe effective performance was derived from the execution timethat included the data transfer time between the CPU and GPU.
The highest in-core performance reached 83.9 GFLOPS,
which was 13% higher than the performance achieved by ourout-of-core implementation. Because data transfer time usuallylimits the performance of GPU applications, we think this 13%
slowdown is acceptable for realizing highly-efficient, out-ofcore stencil computation with a directive-based approach.

1600

Execution time (s)

1400
1200
1000
800
600
400
200
0

8

16

32 64 128 256 512
Blocking factor k

1200
1100
1000
900
800
700
600
500
400
300
200
100
0

Kernel execution
Memcopy on host
Device to host
Host to device

400

Kernel execution
Memcopy on host
Device to host
Host to device

350

Execution time (s)Kernel execution
Memcopy on host
Device to host
Host to device

Execution time (s)

1800

300
250
200
150
100
50

1

(a)

2
4
8
Blocking factor k

16

0

4

8

16 32 64 128 256
Blocking factor k

(b)

(c)

Fig. 12. Breakdown of execution time with different blocking factor k for fixed block size b. Results for (a) Jacobi method, (b) Himeno benchmark, and (c)CIP method.

TABLE IIC OMPARISON OF EFFECTIVE PERFORMANCE OF IN - CORE AND
OUT- OF - CORE IMPLEMENTATIONS .
CodeJacobi
HimenoCIP

In-corep1 (GFLOPS)
32.2
47.5
83.9

Out-of-core
p2 (GFLOPS)
28.5
37.5
73.4

Ratio
p2 /p1 (%)
89
79
87

VI. C ONCLUSION
We presented an extension of OpenACC directives, named
PACC, and its source-to-source translator capable of accelerating out-of-core stencil computation with temporal blocking
on a GPU. Given a PACC code, our translator generates an
OpenACC code such that the code decomposes large data
into smaller chunks, which are then processed in a pipelinedmanner to hide the copy overhead incurred on the CPU.
Furthermore, the generated code is accelerated with temporal
blocking, which reduces the amount of data transfer between
the CPU and GPU.
In experiments, we added PACC directives to three stencil
computation codes, the Jacobi, Himeno, and CIP methods. We
found that the out-of-core performance reached 73.4 GFLOPS
on a Tesla K40 GPU, which was only 13% lower than the incore performance. Thus, PACC directives not only facilitate
out-of-core stencil computation on a GPU but also achieve
high performance on a GPU.
Future work includes an automated framework for finding
the best execution parameters b and k, namely the block size
and the blocking factor.
ACKNOWLEDGMENTS
This study was supported in part by the Japan Society for the
Promotion of Science KAKENHI Grant Numbers 15K12008,
15H01687, 16H02801, and the Japan Science and Technology Agency CREST program, “An Evolutionary Approachto Construction of a Software Development Environment for
Massively-Parallel Computing Systems.” The authors would

like to thank the anonymous reviewers for helpful comments
to improve their paper.
R EFERENCES[1] S. Adams, J. Payne, and R. Boppana. Finite difference time domain (FTDT) simulations using graphics processors. In Proc. High
Performance Computing Modernization Program Users Group Conf.
(HPCMP-UGC’07), pages 334–338, June 2007.
[2] T. Endo and G. Jin. Software technologies coping with memory
hierarchy of GPGPU clusters for stencil computations. In Proc. 16th
IEEE Int. Conf. Cluster Computing (CLUSTER’14), pages 132–139,
Sept. 2014.
[3] T. Endo, Y. Takasaki, and S. Matsuoka. Realizing extremely large-scale
stencil applications on GPU supercomputers. In IEEE Trans. Paralleland Distributed Systems, pages 625–632, Dec. 2015.
[4] R. Himeno. Himeno benchmark, 2015. http://accc.riken.jp/en/supercom/
himenobmt/.
[5] K. Ikeda, F. Ino, and K. Hagihara. Efficient acceleration of mutual
information computation for nonrigid registration using CUDA. IEEE
J. Biomedical and Health Informatics, 18(3):956–968, May 2014.
[6] T. Ikuzawa, F. Ino, and K. Hagihara. Reducing memory usage by thelifting-based discrete wavelet transform with a unified buffer on a gpu.
J. Parallel and Distributed Computing, 93/94:44–55, July 2016.
[7] F. Ino, Y. Munekawa, and K. Hagihara. Sequence homology searchusing fine grained cycle sharing of idle GPUs. IEEE Trans. Paralleland Distributed Systems, 23(4):751–759, Apr. 2012.
[8] Intel Corporation. Intel Xeon Phi product family.
[9] G. Jin, M. Wahib, N. Maruyama, T. Endo, and S. Matsuoka. Locality
optimizations for stencil computations: Algorithms and implementations. In 1st Workshop Programming Abstractions for Data Locality(PADAL’14), Apr. 2014.
[10] T. Kato, F. Ino, and K. Hagihara. PACC: An extension of OpenACC
for pipelined processing of large data on a GPU. In Poster 27th Int’l
Conf. High Performance Computing, Networking, Storage and Analysis(SC’14), Nov. 2014.
[11] C. Liao, D. J. Quinlan, R. Vuduc, and T. Panas. Effective sourceto-source outlining to support whole program empirical optimization.
In Proc. 23rd Int’l Workshop Languages and Compilers for Parallel
Computing (LCPC’10), pages 308–322, Oct. 2010.
[12] N. Maruyama, T. Nomura, K. Sato, and S. Matsuoka. Physis: An
implicitly parallel programming model for stencil computations onlarge-scale GPU-accelerated supercomputers. In Proc. Int’l Conf. High
Performance Computing, Networking, Storage and Analysis (SC’11),
Nov. 2011. 12 pages.
[13] Message Passing Interface Forum. MPI: A message-passing interface
standard. Int’l J. Supercomputer Applications and High Performance
Computing, 8(3/4):159–416, 1994.

[14] H. Murai and M. Sato. An efficient implementation of stencil communication for the XcalableMP PGAS parallel programming language.
In Proc. 7th International Conference on PGAS Programming Models(PGAS’13), pages 142–156, Oct. 2013.
[15] M. Nakao, H. Murai, T. Shimosaka, A. Tabuchi, T. Hanawa, Y. Kodama,
T. Boku, and M. Sato. XcalableACC: Extension of XcalableMPPGAS language using OpenACC for accelerator clusters. In Proc. 2nd
Workshop Accelerator Programming using Directives (WACCPD’14),
pages 27–36, Nov. 2014.
[16] NVIDIA Corporation. CUDA C Programming Guide Version 7.0, March
2015.
[17] NVIDIA Corporation. PGI compiler, 2015.
[18] NVIDIA Corporation. NVIDIA Tesla P100, June 2016.
[19] D. Okada, F. Ino, and K. Hagihara. Accelerating the Smith-Waterman
algorithm with an interpair pruning method for all-pairs comparison ofbase sequences. BMC Bioinformatics, 16(321), Oct. 2015. 15 pages.

[20] OpenACC-Standard.org. The OpenACC application programming interface, version 2.5, Oct. 2015.
[21] OpenMP Architecture Review Board. OpenMP application programming interface, version 4.5, Nov. 2015.
[22] rosecompiler.org.
ROSE compiler infrastructure, 2015.
http://
rosecompiler.org/.
[23] S. W. Williams, A. Waterman, and D. A. Patterson. Roofline: An
insightful visual performance model for floating-point programs andmulticore architectures. Communications of the ACM, 52(4):65–76, Apr.
2009.
[24] E. Wu, Y. Liu, and X. Liu. An improved study of real-time fluid simulation on GPU. Computer Animation and Virtual Worlds, 15(3/4):139–146,
July 2004.
[25] xcalablemp.org. XcalableMP. http://www.xcalablemp.org/.
[26] T. Yabe, F. Xiao, and T. Utsumi. The constrained interpolation profile
method for multiphase analysis. J. Computational Physics, 169(2):556–
593, May 2001.