Prashant Singh

Aravind Sukumaran-Rajam
email: sukumaranrajam.1@osu.edu

Atanas Rountev
email: rountev.1@osu.edu

Fabrice Rastello
email: fabrice.rastello@inria.fr

Louis-Noël Pouchet
email: pouchet@colostate.edu

P Sadayappan
email: sadayappan.1@osu.edu

Register Optimizations for Stencils on GPUs

Keywords:

The recent advent of compute-intensive GPU architecture has allowed application developers to explore high-order 3D stencils for better computational accuracy. A common optimization strategy for such stencils is to expose sufficient data reuse by means such as loop unrolling, with the expectation of register-level reuse. However, the resulting code is often highly constrained by register pressure. While current state-of-the-art register allocators are satisfactory for most applications, they are unable to effectively manage register pressure for such complex high-order stencils, resulting in sub-optimal code with a large number of register spills. In this paper, we develop a statement reordering framework that models stencil computations as a DAG of trees with shared leaves, and adapts an optimal scheduling algorithm for minimizing register usage for expression trees. The effectiveness of the approach is demonstrated through experimental results on a range of stencils extracted from application codes.

Introduction

Stencil computations are an important computational motif in many scientific applications. Typically, a simple stencil computation updates elements of one or more output arrays using elements in the spatial neighborhood from one or more input arrays. The footprint of a stencil is determined by its order, which is the number of input elements from the center read along each dimension. In many scientific applications, the stencil order determines the computational accuracy. For this reason, high-order stencils have been gaining popularity. However, the inherent data reuse within or across statements in such high-order stencils exposes performance challenges that are not addressed by current stencil optimizers.

A significant focus in optimizing stencil computations has been to fuse operations across time steps or across a sequence of stencils in a pipeline [START_REF] Bandishti | Tiling Stencil Computations to Maximize Parallelism[END_REF][START_REF] Grosser | Hybrid Hexagonal/Classical Tiling for GPUs[END_REF][START_REF] Gysi | MODESTO: Data-centric Analytic Optimization of Complex Stencil Programs on Heterogeneous Architectures[END_REF][START_REF] Teja Mullapudi | Poly-Mage: Automatic Optimization for Image Processing Pipelines[END_REF][START_REF] Singh Rawat | Resource Conscious Reuse-Driven Tiling for GPUs[END_REF][START_REF] Wahib | Scalable Kernel Fusion for Memory-bound GPU Applications[END_REF][START_REF] Xue | On Tiling as a Loop Transformation[END_REF]. With high-order stencils, the operational intensity is sufficiently high so that even with just a simple spatial tiling, the computation should theoretically not be memory-bandwidth bound. Consider a GPU with around 300 GBytes/sec global memory bandwidth and a peak double-precision performance of around 1.5 TFLOPS. The required operational intensity to be compute-bound and not memory-bandwidth bound is around 5 FLOPs/byte or 40 FLOPs per double-word. Many high-order stencil computations have much higher arithmetic intensities than 40. For such stencils, achieving a high degree of reuse in cache is very feasible, but high performance is not realized on GPUs. The main hindrance to performance is the high register pressure with such codes, resulting in excessive register spilling and a subsequent loss of performance. As we elaborate in the next section, existing register management techniques in production compilers are not well equipped to address the problem with register pressure for high-order stencils. Addressing this problem in context of GPUs is even more challenging, since most of the widely used GPU compilers like NVCC [START_REF] Nvcc | NVIDIA CUDA Compiler Driver NVCC[END_REF] are closed-source. Even the recent effort by Google (gpucc [START_REF] Wu | gpucc: An Open-source GPGPU Compiler[END_REF]) has only exposes the front-end to the user, and uses the NVCC backend as a black box to perform instruction scheduling and register allocation.

In this paper, we develop an effective pattern-driven global optimization strategy for instruction reordering to address this problem. The key idea behind the instruction reordering approach is to model reuse in high-order stencil computations by using an abstraction of a DAG of trees with shared nodes/leaves, and exploit the fact that optimal scheduling to minimize registers for a single tree with distinct operands at the leaves is well known [START_REF] Sethi | The Generation of Optimal Code for Arithmetic Expressions[END_REF]. We thus devise a statement reordering strategy for a DAG of trees with shared nodes that enables reduction of register pressure to improve performance. The paper makes the following contributions:

• It proposes a framework for multi-statement stencils that reduces register pressure by reordering instructions across statements. • It describes novel heuristics to schedule a DAG of trees that reuse data using a minimal number of registers. • It demonstrates the effectiveness of the proposed framework on a number of register-constrained stencil kernels.

(a) Stencil with lexicographical sweeps for (i=2; i<N-2; i++) for (j=2; j<N-2; j++) { out [i][j] = 0; for (ii=2; ii>=-2; ii--) for (jj=-2; jj<=2; jj++) out Register Allocation and Instruction Scheduling A compiler has several optimization passes, register allocation and instruction scheduling being two of them. Passes before register allocation manipulate an intermediate representation with an unbounded number of temporary variables. The goal of register allocation is to assign those temporaries to physical storage locations, favoring the few but fast registers to the slower but larger memory.

For a fixed schedule, a common approach to perform register allocation is to build an interference graph of the program, which captures the intersection of the live-ranges of temporaries at any program point. Register assignment is then reduced to coloring the interference graph, where each color represents a distinct register [START_REF] Briggs | Improvements to Graph Coloring Register Allocation[END_REF][START_REF] Chaitin | Register Allocation & Spilling via Graph Coloring[END_REF]. Interfering nodes in the interference graph are assigned different colors due to their adjacency. The number of registers needed by the coloring algorithm is lower-bounded by the maximum number of intersecting live-ranges at any program point (MAXLIVE). If MAXLIVE is more than the number of physical registers, spilling of registers and the consequent load/store operations from/to memory are unavoidable.

Register pressure can sometimes be alleviated by reordering the schedule of dependent instructions to reduce the MAXLIVE. Reordering independent instructions is often used to enhance the amount of instruction-level parallelism (ILP), for hiding memory access latency. Thus, there is a complex interplay between instruction scheduling and register allocation, affecting instruction-level parallelism and register pressure, and the associated optimization problem is highly combinatorial. Production compilers generally use heuristics for increasing ILP, with a best-effort greedy control on register pressure. For typical application codes, the negative effect on register pressure is not very significant. However, for high-order stencil codes with a large number of operations and a lot of potential register-level reuse, the impact can be very high, as illustrated by an example below.

Illustrative Example Consider an unrolled version of the double-precision 2D Jacobi stencil computation (Figure 1a) from [START_REF] Stock | A Framework for Enhancing Data Reuse via Associative Reordering[END_REF]. NVCC interleaves the contribution from each input point to different output points to increase instruction level parallelism (ILP). The interleaving performed to increase ILP also has the serendipitous effect of reducing the live range of the register data, and a consequent reduction in register pressure. Nvprof [START_REF] Nvprof | NVIDIA Profiler[END_REF] profiling data on a Tesla K40c device shows that under maximum occupancy, this version performs 3.73E+06 spill transactions, achieving 467 GFLOPS.

Figure 1b shows the same stencil computation after changing the order of accumulation. Exactly the same contributions are made to each result array element, but the order of the contributions has been reversed. With this access pattern for the code in Figure 1b, NVCC fails to perform the same interleaving despite allowing reassociation via appropriate compilation flags. In fact, the register pressure is now exacerbated by the consecutive scheduling of independent operations to increase ILP. For this version, 1.58E+08 spill transactions were measured, with performance dropping to 51 GFLOPS.

This example illustrates a problem with register allocation when the computation has a specific reuse pattern, characteristic of high-order stencil computations. The problem stems from the fact that for most compliers the register allocation and instruction scheduling algorithms that operate at a basic-block level have a peephole view of the computation -they make scheduling/allocation decisions without a global perspective, and thus sometimes work antagonistically. Meanwhile, stencil computations typically have a very regular access pattern. With a better understanding of the pattern, and a global perspective on the computation, it is feasible to devise an instruction reordering strategy to alleviate register pressure.

Solution Approach

In this paper, we circumvent the complexity of the general optimization problem of instruction reordering and register allocation by devising a pattern-specific optimization strategy. Stencil computations involve accumulation of contributions from array data elements in a small neighborhood around each element. The additive contributions to a data element may be viewed as an expression tree. Thus, for multi-statement stencils, we have a DAG of expression trees. Due to the fact that an element may contribute to several result elements, the trees within the DAG can have many shared leaves.

Given a single tree without any shared leaves, it is well known [START_REF] Sethi | The Generation of Optimal Code for Arithmetic Expressions[END_REF] how to schedule its operations in order to minimize the number of registers needed. We use this as the basis for developing heuristics to schedule the operations from the DAG of trees with shared leaves. In contrast to the problem of reordering an arbitrary sequence of instructions to minimize register pressure, a structured approach of adapting the optimal schedule for isolated trees to the case of DAG of trees with shared leaves results in an efficient and effective algorithm that we develop in the next two sections.

Scheduling DAG of Expression Trees

Stencil computations are often succinctly represented using a domain-specific language (DSL). Listing 1 shows a 7-point Jacobi stencil expressed in an illustrative DSL, similar in spirit to stencil computation DSLs such as SDSL [START_REF] Henretty | A Stencil Compiler for Shortvector SIMD Architectures[END_REF] and Forma [START_REF] Ravishankar | Forma: A DSL for Image Processing Applications to Target GPUs and Multi-core CPUs[END_REF]. The core computation is shown in lines 2-4. As with similar DSLs, the user can specify unroll factors for loop iterators (line 9). Loop unrolling, or thread coarsening on GPUs, is often used to exploit register-level reuse in the code. The computation is automatically unrolled as a preprocessing step, before the code is generated and optimized.

It is important to note that using a DSL is not a prerequisite for using the scheduling techniques proposed in this work. As described shortly, our approach works on a DAG of expression trees. This DAG can be automatically extracted either from the DSL representation or from C/Fortran code.

A stencil statement can be defined by the stencil shape (as in lines 2-4) and the input/output data (as in line 8). Each such stencil statement can be represented by a labeled expression tree. For example, the tree corresponding to the computation in Listing 1 has array element out[k, j, i] as its root, scalars a, b, c and accesses to elements of array in as its leaves, and arithmetic operators * and + as inner nodes.

An expression tree for a stencil computation has three types of nodes: (1) nodes n ∈ N mem representing accesses to memory locations, (2) nodes n ∈ N op representing binary/unary arithmetic operators, and (3) leaf nodes representing constants. All leaf nodes in N mem correspond to reads of array elements (e.g., in[k + 1, j, i]) or scalars. The root of the expression tree is also in N mem and corresponds to a write to an array element (e.g., out[k, j, i]) or a scalar. We associate a unique label with each read/written memory location, and assign to each node in N mem the corresponding label. The remaining tree nodes are in N op . Figure 2b shows the expression tree for an illustrative expression.

In a preprocessing step, we introduce k-ary nodes for associative operators. For example, for the tree in Figure 2b, the chain of + nodes is replaced with a single "accumulation" + node. Figure 2c shows the resulting expression tree; the numbers on the nodes will be described shortly. The semantics of an accumulation node is as expected: the value is initialized as appropriate (e.g., 0 for +, 1 for * ) and the contributions of the children are accumulated in arbitrary order.

We often consider a sequence of stencil computations-for example, in image processing pipelines [START_REF] Ravishankar | Forma: A DSL for Image Processing Applications to Target GPUs and Multi-core CPUs[END_REF]. Each computation in the sequence will be represented by a separate expression tree. Similarly, unrolling will result in distinct Listing 1. The input representation in the DSL expression trees for each unrolled instance. For example, after unrolling along dimensions k and j in Listing 1, there will be a sequence of four expression trees. In some cases the output of a tree is used as an input to a later tree in the sequence. In such a case, there is a flow dependence: the root of the producer tree has the same label as some leaf node of the consumer tree (without an in-between tree that writes to that label). In the input to our analysis, this flow is represented by a dependence edge from the root node to the leaf node. Thus, the entire computation is represented as a DAG of expression trees. Throughout the paper, we make the following two assumptions: (1) the assembly instructions generated for the DAG of trees after register allocation are of the form r 1 ← r 2 op r 3 , where r 1 and r 2 can be the same; (2) each operand/result requires exactly one register. This condition is only enforced to simplify the presentation of the next two sections, and can be very easily relaxed [START_REF] Appel | Generalization of the Sethi-Ullman Algorithm for Register Allocation[END_REF]. Our objective is to schedule the computations in the DAG so that the register pressure is reduced.

1 function j3d7pt ( out , in , a , b , c ) { 2 out [ k ][ j ][ i ] = a *( in [ k +1][ j ][ i ]) + b *( in [ k ][j -1][ i] + 3 in [ k ][ j ][ i -1] + in [ k ][ j ][ i ] + in [ k ][ j ][ i +1] + 4 in [ k ][ j +1][ i ]) + c *( in [k -1][ j ][ i ]); 5 } 6 parameter L , M , N ; 7 iterator k , j , i ; 8 double in [ L ][ M ][ N ] , out [ L ][ M ][ N ] , a , b , c ; 9 unroll k =2 , j =2;

Sethi-Ullman Scheduling

We will use "data sharing" to refer to cases where the same memory location is accessed at multiple places. There are two types of data sharing: (1) within a tree: several nodes from N mem have the same label; and (2) across trees: in a DAG of trees, nodes from distinct trees have the same label.

A classic result, due to Sethi and Ullman [START_REF] Sethi | The Generation of Optimal Code for Arithmetic Expressions[END_REF], applies to a single expression tree without data sharing (i.e., each n ∈ N mem has a unique label), and with binary/unary operators. They present a scheduling algorithm that minimizes the number of registers needed to evaluate such an expression tree under a spill-free model. 1 Each tree node n is assigned an Ershov number [START_REF] Aho | Compilers: Principles, Techniques, and Tools[END_REF]; we will refer to them as "Sethi-Ullman numbers" and denote them by su(n). They are defined as

su(n) =                1 n is a leaf su(n 1 ) n has one child n 1 max (su(n 1 ), su(n 2 )) su(n 1 ) su(n 2 ) 1 + su(n 1 ) su(n 1 ) = su(n 2 ) (1)
The last two cases apply to a binary op node n with children n 1 and n 2 . Intuitively, su(n) is the smallest possible number of registers used for the evaluation of the subtree rooted at n. The first two cases are self-explanatory. For a binary op node n, if one child n ′ has a higher register requirement (case 3), this "big" child should be evaluated first. The result of n ′ will be stored in a register, which will be alive while the second ("small") child is being evaluated. The remaining su(n ′ ) -1 registers used for n ′ are available (and enough) to evaluate the small child. Finally, the register of n ′ can be used to store the result for n, meaning that su(n) is equal to su(n ′ ). If the order of evaluation were reversed, the result of the small child would have to be kept in a register while n ′ is being evaluated, which would lead to sub-optimal su(n) = 1 + su(n ′ ). In the last case in equation ( 1), both children have the same register needs; thus, their relative order of evaluation is irrelevant and one extra register is needed for n. Of course, under the definitions in equation [START_REF] Aho | Compilers: Principles, Techniques, and Tools[END_REF], su(n) is the same as MAXLIVE for the tree rooted at n.

It is straightforward to generalize Sethi-Ullman numbering to trees containing accumulation nodes (as in Figure 2c). Each such accumulation node n has children n i for 1 ≤ i ≤ k. Let mx = max i {su(n i )}. If there is a single child n j with su(n j ) = mx, this child is scheduled for evaluation first, and therefore su(n) = mx. If two or more children n j have su(n j ) = mx, one of them is scheduled first; however, in this case su(n) = 1 + mx. In both cases, the order of evaluation of the remaining children is irrelevant. Figure 2c shows the Sethi-Ullman numbers for the sample expression tree.

Note that the schedules produced by this approach perform atomic evaluation of subexpressions: one of the children is evaluated completely before the other ones are considered. For a tree without data sharing, this restriction does not affect the optimality of the result. In the presence of data sharing, atomic evaluation may not be optimal.

Since stencils read values from a limited spatial neighborhood, data sharing often manifests in the DAG of expression trees. For example, in Listing 1, in[k][j][i] will be an input to all four expression trees corresponding to the unrolled stencil statements. One can also find other nodes in Listing 1 that will be shared across multiple expression trees. For such DAGs, the Sethi-Ullman algorithm cannot be directly applied to obtain an optimal schedule. In Section 3.2, we present an approach to compute an optimal atomic schedule for a DAG of expression trees with data sharing. In cases when finding an optimal evaluation can be prohibitively expensive, Section 3.3 presents heuristics to trade off optimality in favor of pruning the exploration space. Finally, restricting the evaluation to be atomic can generate sub-optimal schedules. Section 4 presents a remedial slice-and-interleave algorithm that performs interleaving on the output schedule generated by the approach presented in Section 3.2.

Scheduling a Tree with Data Sharing

Figure 3a shows an expression tree with data sharing. For illustration, nodes with the same label are connected in the figure . Recall that we assume a spill-free model, therefore a shared label loaded once into a register will remain live for all its uses. With data sharing, there is a possibility that (1) a label is already live before we begin the recursive evaluation of a subtree that has its subsequent use, and/or (2) a label must remain live even after the evaluation of the subtree in which it is used. The optimal schedule of a subtree is affected by the labels that are live before and after the evaluation of the subtree. Therefore, we need to add live-in/out states as parameters to the computation of the optimal schedule of a subtree. In this section, we present an approach to optimally schedule a tree with data sharing, under the model of atomic evaluation of children; we defer the interleaving of computation across subtrees to Section 4. For a node n, let uses(n) be the set of labels used in the subtree rooted at n. Figure 3a shows uses(n) for each internal node n. The live-in set for a node n, denoted by in(n), contains all labels that are live before the subtree rooted at n is evaluated. The live-out set is

out (n) = (in(n) ∪ uses(n)) \ kill(n) (2)
where kill(n) is the set of labels that have their last uses in the subtree rooted at n. Note that kill(n) is context-dependent, i.e., the set will vary depending on the order in which the node is evaluated. The kill sets can be computed on the fly by maintaining the number of occurrences of each label l in the current schedule, and comparing it with the total number of occurrences of l in the entire DAG. We now show how to compute a modified Sethi-Ullman number, su ′ for each node n, when provided with an "evaluation context" in terms of live-in and live-out labels. Consider a node n with some in and out state. Just before the evaluation of n begins, |in(n)| registers are live. Similarly, just after the evaluation of n finishes, |out (n)| registers will be live. During the evaluation of n, additional registers may become live, while some of the other live registers may be released. Now su ′ (n, in, out) represents the maximum number of registers that were simultaneously live at any point during the evaluation of n. We also define su ′ (π , in, out), where π is a sequence of the children nodes of n. This value will represent the maximum number of registers that were simultaneously live at any point during the evaluation of n with its children ordered in the sequence described by π .

For simplicity we will use su ′ (n) instead of su ′ (n, in, out), but the definitions will use the live-in/out sets in(n) and out (n).

For a leaf node n ∈ N mem with |in(n

)| = α, su ′ (n) =      α + 1 label(n) in(n) α label(n) ∈ in(n)
The first case implies that a new register must be reserved for the label of n if it is not already live before the evaluation of n. The second case is self-explanatory.

To compute su ′ for a k-ary (binary or accumulation) node n with children n 1 . . . n k , we need to explore all k! evaluation orders of the children. Let π be any permutation of the children of n representing their evaluation order. Then

su ′ (n) = min π su ′ (π ).
For the purpose of explanation, suppose the permutation π = ⟨n 1 , n 2 ⟩ is one particular evaluation order for a binary node n with children n 1 , n 2 . To compute su ′ (π ), first we determine the live-in and live-out sets for nodes in π as follows: in(n 1 ) = in(n), in(n 2 ) = out (n 1 ), and out (n 2 ) = out (n); here out (n 1 ) is as defined in equation 2. This provides the required context to compute su ′ (n 1 ) and su ′ (n 2 ). Let mx = max i {su ′ (n i )}, so that mx equals the maximum number of simultaneously live registers at any time during the evaluation of π . Then,

su ′ (π ) =      1 + mx n i ∈ N mem & label(n i ) ∈ out (n) mx otherwise In case 2, if n 1 ∈ N op , or if n 1 ∈ N mem but label(n 1 ) out (n),
then the result of the computation, identified by the label of the node n, can reuse the register of n 1 (similarly for n 2 ). However, in case 1, where both n 1 and n 2 are leaf nodes in N mem and both must be live after evaluating n, we need an additional register to hold the result.

For an accumulation node with k children, consider permutation π = ⟨n 1 , n 2 . . . , n k ⟩. Suppose we have computed all su ′ (n i ) and let mx = max i {su ′ (n i )}. Then,

su ′ (π ) =                      mx su ′ (n 1 ) = mx & n 1 ∈ N mem & label(n 1 ) out (n 1 ) & su ′ (n j ) mx, 2 ≤ j ≤ k mx su ′ (n 1 ) = mx & n 1 ∈ N op & su ′ (n j ) mx, 2 ≤ j ≤ k 1 + mx otherwise
Just like the generalization of su(n) for accumulation nodes in Section 3.1, su ′ = mx when the following two conditions hold: (1) n 1 requires the maximum number of simultaneously live registers, and the rest of the nodes in π can be completely evaluated using the registers released by n 1 , and (2) the register holding n 1 can be reused by n, i.e., either n 1 ∈ N op (case 2), or n 1 is a leaf node that is not live beyond this point (case 1). In all other scenarios, we need mx + 1 registers (case 3).

The computation of su ′ (n) for a tree without an evaluation context is shown in Figure 3b, and with an evaluation context is shown in Figure 4a. For the same tree, Figure 4b shows the permutation with minimum su ′ . In all three figures, the children of a node are ordered left-to-right, which defines the corresponding permutation.

In some cases, exhaustively exploring all permutations of the children may be unnecessary. In the tree of Figure 4a, there are two subtree operands of the accumulation node that share no data (the first and the third subtree operands in the left-to-right order). Therefore, even though the scheduling within those two subtrees may be influenced by the evaluation context, they do not influence each other's evaluation. Let passthrough denote the labels that are live both before and after the evaluation of node n:

passthrough(n) = in(n) ∩ out (n).
Then, for a k-ary node n, any two of its children n i and n j do not influence each other's evaluation if

(uses(n i ) ∩ uses(n j )) \ passthrough(n) = ∅ (3)
In such a scenario, for the node n, we can create maximal clusters of its children that share data, but the only data shared among clusters is the passthrough labels of n. For example, if children t 1 and t 2 share label l 1 , and children t 2 and t 3 share label l 2 , where l 1 , l 2 are non-passthrough labels of the parent node, then {t 1 , t 2 , t 3 , t 4 } must belong to the same cluster. We extend the intuition behind Sethi-Ullman scheduling algorithm to establish that different clusters cannot influence each other's evaluation. Then, for each cluster c i , we can independently compute su ′ (c i ) with in(c i ) = in(n). We only need to explore all permutations within the nonsingleton clusters. We propose the following theorems to establish an evaluation order for the clusters.

Theorem 3.1. For k clusters c i 1 ≤ i ≤ k such that |in(c i )| ≤ |out (c i )|,
the one with larger su ′ (c i ) -|out (c i )| will be prioritized for evaluation over others in the optimal schedule. In the special case where all the clusters have the same su ′ (c i ) -|out (c i )|, they can be evaluated in any order without affecting MAXLIVE.

This result is a direct consequence of the Sethi-Ullman algorithm. The cluster with larger su ′ (c i ) -|out (c i )| will release more registers, which can be reused by the next cluster. The special case too is a direct consequence of the Sethi-Ullman algorithm, where two sibling nodes with the same su can be evaluated in any order (case 4 of equation 1). We prove the result by contradiction. Suppose that c 2 is evaluated before c 1 in the optimal schedule. Since the schedule is optimal, su ′ (c 2 ) ≥ su ′ (c 1 ). Now we change this optimal schedule by moving the evaluation of c 1 before c 2 . Evaluating c 1 earlier will release |in(c 1 )| -|out (c 1 )| (i.e., ≥ 1) registers, which can then be used in the evaluation of c 2 . Based on the previous equations, the su ′ (c 2 ) will either decrease or remain the same, depending on whether the number of registers released by c 1 is greater than, or equal to 1. This modified schedule therefore either has the same, or has lower su ′ than the optimal schedule, making it an optimal schedule. Again, we prove the result by contradiction. Suppose that su ′ (c 1 ) < su ′ (c 2 ), and c 2 is scheduled before c 1 in the optimal schedule. We change this schedule by moving the evaluation of c 1 before that of c 2 . From Theorem 3.2, su ′ (c 2 ) after this change will either remain the same, or decrease. Thus, su ′ for the new schedule will either be the same, or reduce if su ′ (c 2 ) was the maximum, making it an optimal schedule.

Based on these theorems, Algorithm 1 summarizes the evaluation an optimal schedule for a tree with data sharing.

Heuristics for Tractability

For a non-singleton cluster c, the algorithm presented in Section 3.2 can become prohibitively expensive if |c | is large. For example, when |c | changes from 7 to 8, the permutations explored increase from 5040 to 40320. We now present some heuristics that trade off optimality for tractability, and a caching technique to further speed up the algorithm.

Pruning Heuristics We begin by establishing, for any node n, the bounds on su ′ (n). When n is evaluated with non-empty contexts in and out, the bounds are:

su ′ (n, ∅, ∅) ≤ su ′ (n, in, out ) ≤ su ′ (n, ∅, ∅) + |in ∪ out |
We prove the lower bound by contradiction; the proof for upper bound is omitted due to space constraints.

su ′ (n, ∅, ∅) ≤ su ′ (n, in, out ): Assume to the contrary that su ′ (n, in, out ) < su ′ (n, ∅, ∅). We will modify the schedule S corresponding to su ′ (n, in, out ) as follows: prepend a stage to S that loads the labels ∈ in(n) into |in| registers, and make in(n) = ∅. Append a state to S that stores all the labels ∈ out (n) from the respective registers into memory, and make out (n) = ∅. This modified schedule corresponds to su ′ (n, ∅, ∅), and hence, su ′ (n, ∅, ∅) = su ′ (n, in, out ).

With the bounds established, instead of exploring all permutations, we can sacrifice optimality and stop further exploration when we are close to the optimal schedule. We use a tunable parameter d, and stop trying the permutations any further when su ′ (n, in, out ) -su ′ (n, ∅, ∅) ≤ d. For the experimental evaluation in Section 5, we set d to 1.

For a cluster c with |c | > 8, we also apply a partitioning heuristic, which recursively partitions the subtrees in c into sub-partitions where each sub-partition can be of a maximum size p, with p < 8. The partitioning is based on either of the two criteria:

on "label affinity": the subtrees that share the maximum labels are greedily assigned to the same subpartition as long as the size of the sub-partition is less than p. Such partitioning is based on the notion that evaluating subtrees with maximum uses together will potentially reduce passthrough labels, and MAXLIVE. on "release potential": the subtrees that have the last uses of some labels are placed in a sub-partition, and that sub-partition is eagerly evaluated. This partitioning is based on the notion that the released registers can be reused by the next partition. Once the sub-partitions are created, we only exhaustively explore all permutations of subtrees within a sub-partition. If the number of sub-partitions created is less than 8, then we also try all the permutations of the sub-partitions themselves. For example, if |c | = 8, and the partitioning heuristic creates two sub-partitions p 1 and p 2 of size 4 each, then our exploration space will be {p 1 , p 2 } and {p 2 , p 1 }, while exploring all 4! permutations of subtrees within p 1 and p 2 each -a total of 2 × 4! × 4! permutations instead of 8! permutations.

We also let the user externally specify a threshold that upper-bounds the total number of permutations for a tree.

Memoization For a node n, a lot of permutations of its children will differ in only a few positions. In such cases, we end up recomputing su ′ for a child multiple times, even when the live-in/out context for the child remains unchanged.

These recomputations can be avoided by a simple memoization, where for a node n, we map su ′ (n) as a function of a minimal context. The minimal context strips away labels that are not in uses (n), but are in passthrouдh(n). The su ′ (n) with minimal context can be suitably adjusted to get su ′ (n) with a different context that has some passthrough labels added to the minimal live-in/out. For example, suppose that su ′ (n) is 3 when the minimal in(n) = {a,b} and the minimal out (n) = ∅. Then su ′ (n) when evaluating it with in(n) = {a, b, c}, out (n) = {c} and c uses (n) will be 2+1=3, and the optimal schedule will remain unchanged. Memoization greatly reduces the total evaluation time, thereby enabling the exploration of a large number of permutations.

Scheduling a DAG of expression trees

For each stencil statement that is mapped to an expression tree, Section 3.2 described a way to schedule it. This section ties everything together for a multi-statement stencil by describing how to schedule a DAG of expression trees. For optimal scheduling, one needs to explore all topological orders for the trees in the DAG, and then evaluate all the trees independently for each topological order. This may be practical if the size of the DAG is small. Otherwise, we must sacrifice optimality for tractability, and fix the evaluation order of the trees in the DAG before the trees are individually evaluated.

Algorithm 2: Schedule-DAG (D, R)

Input : D: DAG of expression trees, R: Per-thread register limit Output : An optimal schedule S for the D We use the greedy heuristic described by Rawat et al. [START_REF] Singh Rawat | Resource Conscious Reuse-Driven Tiling for GPUs[END_REF] to fix the evaluation order of trees in the DAG. At each step, the heuristic tries to fix the evaluation order of two nodes in the DAG. We begin by computing, for each pair of transitive dependence-free trees p i in the DAG, a metric M i that encodes: (a) the number of labels shared between them, and (b) the number of common input arrays read by them. Among the computed M i , we choose the one that has the highest non-zero value, and fix the evaluation order of its tree pair to be contiguous to enhance reuse proximity in the final schedule. The DAG is updated by fusing the nodes corresponding to the two trees into a "macro node". Post fusion, we update the dependence edges to and from the macro node, and recompute the metrics for the next step. The algorithm terminates when no more nodes can be fused.

Once the algorithm terminates, we perform a topological sort of the final DAG, and expand the DAG macro nodes to their tree sequences. For these ordered trees, we can generate code versions with different degree of splits. One extreme would be a version where all the trees are in a single kernel (max-fuse), and another extreme would be a version where each tree is a distinct kernel (max-split) [START_REF] Bondhugula | A Practical Automatic Polyhedral Parallelizer and Locality Optimizer[END_REF][START_REF] Wahib | Automated GPU Kernel Transformations in Large-Scale Production Stencil Applications[END_REF]. For computeintensive stencils with many-to-many reuse, a single kernel can have extremely high register pressure, sometimes causing spills despite allowing for the maximum permitted registers per thread. For such cases, performing kernel fission instead of generating a single kernel for the entire computation might improve performance. The split kernels will incur additional data transfers from global memory, but the register pressure per kernel will be much lower, giving the user an opportunity to further enhance register-level For each split version created, the tree sequence in it is evaluated using Algorithm 1. The returned schedule is the one that gives maximum performance. Algorithm 2 outlines the entire process.

Interleaving Expressions

At this point, we have a schedule for the entire DAG of trees, but with atomic evaluation enforced. However, interleaving within/across trees can be instrumental in reducing MAXLIVE. For example, in the unrolled stencil of Listing 1, there is no reuse within a stencil statement, but plenty of reuse across stencil statements. We will see later in Section 5 that relaxing the constraint of atomic evaluation, and performing interleaving is imperative for performance in such stencils. A compiler optimization that performs some interleaving is common subexpression elimination (CSE). However, we require a more general interleaving that works at the granularity of common labels instead of common subexpressions. For example, Figure 5a shows an expression tree where su ′ (S ) is the largest, and the operands of the accumulation node are evaluated in order from left to right in the final schedule. Also, {c[i], b} uses (S ). The fact that {c[i],b} ∈ passthrouдh(S ) adds to su ′ (S ). By slicing the expression (e[i] * b)/c[i] and placing it after the expression b * c[i] as shown in Figure 5b,c[i] and b will no longer be in in(S ). Instead of those two labels, a temporary label holding the value of the sliced expression will be added to in(S ), and hence su ′ (S ) will reduce by 1. Note that this is not CSE, but a more general optimization aimed to reduce MAXLIVE. This slice-and-interleave optimization slices a target expression, and interleaves it with a source expression, so that su ′ at a program point reduces. It subsumes CSE if the source and target expressions are the same.

We perform the slice-and-interleave optimization at two levels: (a) within an expression tree, where the source and target expressions belong to the same tree; and (b) across the expression trees in the DAG, where source and target Algorithm 3: slice-and-interleave (T , in, out) Input :T : an input tree with schedule S and contexts in and out Output : S : The schedule after applying slice-and-interleave 1 L ilv ← ∅; 2 min_exprs ← sequence of minimal expressions extracted from S whose operands are leaf nodes;

3 for each expression s in min_exprs do 4

L ilv ← all the labels seen in the schedule till s ∪ labels with single occurrence in T; expressions belong to different trees. For a chosen source expression e s rooted at node n, we compute a set of labels, L ilv , which is a union of all the labels that were observed in the schedule till now, with the labels that have a single occurrence in the DAG. We now try to find a set of target expressions operating on just the labels from L ilv . To find the target expressions, we start with minimal expressions, i.e., the simplest expressions whose operands are leaf nodes ∈ N mem . Once we find a minimal expression e m that operates only on the labels ∈ L ilv , we find the root node r of e m , and grow e m to the expression rooted at parent (r ). We continue to grow the expression till we have a maximal expression that only operates on the labels ∈ L ilv . For each target expressions thus discovered, we check if slicing and placing it between the source expression e s and subtree t s immediately following e s in the schedule decreases |in(t s )|. If it does, then slice-and-interleave is performed.

Illustrative Example Let b * c[i] be the source expression in the tree of figure 5a. One of the explored target expressions will be e[i] * b, since it only uses nodes ∈ L ilv . Now we try to grow the target expression by changing the root from * to /, making (e[i] * b)/c[i] the new target expression. All the labels used in the grown target expression also belong to L ilv . A further attempt to grow will change the target expression to

((e[i] * b)/c[i]) * f [i]. However, f [i] L ilv .
Therefore, we backtrack and finalize (e[i] * b)/c[i] as a target expression, since it is the maximal expression with all the labels in L ilv . Placing the target expression after the source expression decreases in(S ) by 1. Therefore, we perform the slice-and-interleave optimization. Algorithm 3 outlines the slice-and-interleave algorithm that tries out different source expressions, and continuously finds the target expressions within the tree to interleave in order to reduce the live ranges. The slice-and-interleave across the trees in a DAG is similar.

Our framework parses stencil statements written in a subset of C: the array access indices in the stencil statements must be an affine function of the surrounding loop iterators, program parameters, and constants; loop iterators and parameters must be immutable in the stencil statements. The framework supports auto-unrolling along different dimensions to expose spatial reuse across stencil statements.

To ensure a tight coupling, several prior efforts on guiding register allocation or instruction scheduling were implemented as a compiler pass in research/prototype compilers [START_REF] Berson | Integrated Instruction Scheduling and Register Allocation Techniques[END_REF][START_REF] Domagala | Register Allocation and Promotion Through Combined Instruction Scheduling and Loop Unrolling[END_REF][START_REF] Ramaswamy Govindarajan | Minimum Register Instruction Sequence Problem: Revisiting Optimal Code Generation for DAGs[END_REF][START_REF] Poletto | Linear Scan Register Allocation[END_REF][START_REF] Rong | Tree Register Allocation[END_REF], or open-source production compilers [START_REF] Ryan | A Global Progressive Register Allocator[END_REF][START_REF] Sarkar | Extended Linear Scan: An Alternate Foundation for Global Register Allocation[END_REF]. However, like some other recent efforts [START_REF] Basu | Compiler-Directed Transformation for Higher-Order Stencils[END_REF][START_REF] Jin | Libra: An Automated Code Generation and Tuning Framework for Registerlimited Stencils on GPUs[END_REF][START_REF] Stock | A Framework for Enhancing Data Reuse via Associative Reordering[END_REF], we implement our reordering optimization at source level for the following reasons: (1) it allows external optimizations for closed-source compilers like NVCC; (2) it allows us to perform transformations like exposing FMAs using operator distributivity, and performing kernel fusion/fission, which can be performed more effectively and efficiently at source level; and (3) it is input-dependent, not machine-or compilerdependent -with an implementation coupled to compiler passes, it would have to be re-implemented across compilers with different intermediate representation. Our framework massages the input to a form that is more amenable to further optimizations by any GPU compiler, and we use appropriate compilation flags whenever possible to ensure that our reordering optimization is not undone by the compiler passes.

We evaluate our framework for the benchmarks listed in Table 1 on a Tesla K40c GPU (peak double-precision performance 1.43 TFLOPS, peak bandwidth 288 GB/s) with NVCC-8.0 [START_REF] Nvcc | NVIDIA CUDA Compiler Driver NVCC[END_REF] and LLVM-5.0.0 compiler (previously gpucc [START_REF] Wu | gpucc: An Open-source GPGPU Compiler[END_REF]). The first five benchmarks are stencils typically used in iterative processes such as solving partial differential equations [START_REF] Hpgmg | High-Performance Geometric Multigrid[END_REF]. The remaining three are representative of complex stencil operations extracted from applications. hypterm is a routine from the ExpCNS Compressible Navier-Stokes mini-application from DoE [START_REF] Exact | ExaCT: Center for Exascale Simulation of Combustion in Turbulence: Proxy App Software[END_REF]; the last two stencils are from the Geodynamics Seismic Wave SW4 application code [START_REF]Seismic Wave Modelling (SW4) -Computational Infrastructure for Geodynamics[END_REF]. For each benchmark, the original version is as written by application developers without any loop unrolling; the unrolled version has the loops unrolled explicitly; and the reordered version is the output from our code generator. On an Intel i7-4770K processor, the code generator generated each reordered version under 4 seconds. When the net unrolling factor is limited to 4, the size of each reordered version is under 600 lines of code. The read-only arrays are annotated with the restrict keyword in all the versions to allow efficient loads via the texture pipeline. All the stencils are double-precision, compiled with NVCC flags '-use_fast_math Xptxas "-dlcm=ca"', and LLVM flags '-O3 -ffast-math -ffp-contract=fast'. Since none of the versions use shared memory, using 'dlcm=ca' for NVCC enhances performance by caching the global memory accesses at L1. However, we 1. Benchmark characteristics notice no discernible performance difference when compiling without the '-use_fast_math' flag in NVCC. We explore different instruction schedulers implemented in LLVM (default, list-hybrid, and list-burr) for the original and unrolled code, and report the best performance. To minimize instruction reordering for our reordered code, we use LLVM's default instruction scheduler, and do not use the -ffast-math option during compilation. We test all the versions against a standard C implementation of the benchmarks for correctness: the difference in each computed output value with that of the C implementation must be less than a tolerance of 1E-5.

Loop Unrolling For the experiments, the iterative kernels were unrolled along a single dimension to expose spatial reuse. Loop unrolling offers the compiler an opportunity to exploit ILP, but scheduling independent instructions contiguously may increase register pressure. Consider an unrolled version of 2d25pt, compiled with 32 registers. From table 1, it is clear that the unrolled code has a high degree of reuse. Listing 2 shows the SASS (Shader ASSembler) snippet generated using NVCC for the unrolled version of 2d25pt after register allocation. The instructions not relevant to the discussion are omitted in Listing 2 (and 3), leading to non-contiguous line numbers. The lines highlighted in red show the instructions involving the same memory location -line 1 loads a value from global memory into register R4, and spills it in line 2 without using R4 in any of the intermediate instructions. Such wasteful spills are a characteristic of register-constrained codes. The same value is reloaded from local memory into R4 in line 4, and R4 is subsequently used in lines 5 and 8. The uses of R4 are placed far apart in SASS, adding to the register pressure. Interspersed with these instructions are the load (line 3) and subsequent uses of register R12. The interleaving increases ILP, but the uses of R12 are also placed very far apart. A better schedule can perhaps achieve the same ILP with less register pressure and fewer spills.

Listing 3 shows the SASS snippet for the reordered code generated by our code generator. Using operator distributivity, the multiplication of the coefficient to the additive contributions is converted by our preprocessing pass into fused multiply-adds. Notice that all the uses of register R20 (highlighted in red) are tightly coupled. The same holds for Listing 2. SASS snippet for the unrolled code For the 3d125pt stencil, table 2 shows some profiling metrics gathered by Nvprof with NVCC. The texture throughput for the original code indicates that the stencil performance is limited by the texture cache bandwidth. Loop unrolling halves the accesses to texture cache and the executed load instructions, but results in a significant drop in IPC due to lowered occupancy. To better expose reordering opportunities after unrolling, the preprocessing pass of the reordering framework exploits operator distributivity and converts all the contributions in an individual statement to FMA operations. Therefore, instead of 130 FLOPs per stencil point, the reordered version performs 250 FLOPs. As measured by Nvprof, we incur a 2× increase in floating point operations, but achieve significant reuse in registers at a higher occupancy, which consequently improves the IPC and execution time.

Register Pressure Sensitivity In GPUs, the number of registers per thread can be varied at compile time by trading off the occupancy. Many auto-tuning efforts have recently been proposed to that end [START_REF] Hayes | Orion: A Framework for GPU Occupancy Tuning[END_REF][START_REF] Li | Critical points based register-concurrency autotuning for GPUs[END_REF]. Table 4 shows the performance, and the local memory transactions reported by Nvprof, with varying register pressure. Due to space constraints, we only present the data for NVCC compiler. We make the following observations: (a) our optimization strategy reduces the register pressure for all the thread configurations; (b) increasing registers per thread for codes exhibiting very high spills results in better performance, e.g., 8× for rhs4th3fort; and (c) for low spills, better performance can be achieved by either increasing occupancy (e.g., reordered code for 3d125pt and hypterm), or maximizing registers per thread (e.g., all the codes for rhs4th3fort).

Finding a right balance between register pressure and occupancy is non-trivial, and an active research field [START_REF] Hayes | Orion: A Framework for GPU Occupancy Tuning[END_REF][START_REF] Li | Critical points based register-concurrency autotuning for GPUs[END_REF][START_REF] Unkule | Automatic Restructuring of GPU Kernels for Exploiting Inter-thread Data Locality[END_REF][START_REF] Xie | Enabling Coordinated Register Allocation and Thread-level Parallelism Optimization for GPUs[END_REF]. We perform a simple auto-tuning by varying the tile sizes by powers of 2, and varying registers per thread [START_REF] Li | Critical points based register-concurrency autotuning for GPUs[END_REF]. The best performance in GFLOPS for the auto-tuned code with NVCC and LLVM compilers is shown in figure 6. Unlike the case with 32 and 64 registers per thread, the unrolled code outperforms the original code for all benchmarks, highlighting the importance of loop unrolling and register-level reuse. Our reordering optimization improves the performance by (a) producing a code version that uses fewer registers, and hence can achieve higher occupancy; and (b) helping expose and schedule independent FMAs together for simple accumulation stencils, thereby hiding latency.

Kernel fission

From table 1, we select the last three multistatement, compute-intensive stencils, for which we anticipate high volume of spills in the max-fuse form, and expect kernel fission to be beneficial. For these three stencils, we generate versions with varying degree of splits (Section 3.4). Some splits require promoting the storage from scalars to global arrays, while others require recomputations due to dependence edges in the DAG. Table 3 shows the Nvprof metrics with NVCC for two reordered codes: a version with maximum fusion (max-fuse), and a version with split kernels. Note that in each case, even though the DRAM reads increase going from max-fuse to split version, the IPC also increases. This is because the register pressure per kernel is much lower in the split version, and hence we can unroll the computation to further exploit register-level reuse. This increase in register-level reuse is reflected in the reduced L2 read transactions. We observe nearly 10% performance improvement for the split version over max-fuse version for all three stencils. Prior works have noted the importance of kernel fusion for bandwidth-bound stencils [START_REF] Gysi | MODESTO: Data-centric Analytic Optimization of Complex Stencil Programs on Heterogeneous Architectures[END_REF][START_REF] Singh Rawat | Resource Conscious Reuse-Driven Tiling for GPUs[END_REF][START_REF] Wahib | Scalable Kernel Fusion for Memory-bound GPU Applications[END_REF], and trivial kernel fission to aid fusion by reducing shared . Performance on Tesla K40c with benchmarks tuned for tile size and register limit memory usage [START_REF] Wahib | Automated GPU Kernel Transformations in Large-Scale Production Stencil Applications[END_REF]. Such a kernel fission has very limited applicability in stencils with significant many-to-many reuse across statements. However, our motivation for kernel fission is orthogonal to prior efforts -we use kernel fission as a means to reduce the register usage of the max-fuse kernel, and then improve the register reuse for split kernels by ample unrolling and instruction reordering.

With our reordering optimizations applied to the benchmarks, we achieve speedups in the range of 1.22×-2.34× for NVCC, and 1.15×-2.08× for LLVM. We finally discuss the effect of optimizations discussed in Section 3.3. The benchmark derivative has a large number of independent trees; each tree is an accumulation. The framework takes 3.42 seconds to generate code for it, and the memoization function is invoked 1.42E+05 times. Without memoization, the code generation time increases to 5.71 seconds. Our framework is well suited to enhance the performance of "optimize once, execute multiple times" stencils found in production codes where the compilation/optimization time is amortized over the stencil execution.

Related Work

Register allocation has been extensively studied: from the seminal work of Chaitin [START_REF] Chaitin | Register Allocation & Spilling via Graph Coloring[END_REF] on using graph coloring, to the more recent SSA based schemes that exploit the possible decoupling of the allocation phase to the assignment phase [START_REF] Colombet | Graph-coloring and treescan register allocation using repairing[END_REF][START_REF] Mössenböck | Linear Scan Register Allocation in the Context of SSA Form and Register Constraints[END_REF]. Many extensions/improvements have been applied [START_REF] Lueh | Fusion-based Register Allocation[END_REF][START_REF] Magno | Register Allocation by Puzzle Solving[END_REF][START_REF] Michael | A Generalized Algorithm for Graph-coloring Register Allocation[END_REF], but almost all the existing works are restricted to the register allocation for fixed schedule. However, it is folk knowledge that improvements to register allocation alone does not bring much performance gain. The interplay of register allocation and scheduling becomes quite important for architectures with ILP, since there is a tradeoff between increasing ILP and exposing locality. Hence, prior work on hyper/super-block list or modulo pre-pass scheduling [START_REF] Codina | A unified modulo scheduling and register allocation technique for clustered processors[END_REF][START_REF] Motwani | Combining Register Allocation and Instruction Scheduling[END_REF][START_REF] Wang | Software Pipelining with Register Allocation and Spilling[END_REF] were extended to account for register pressure. Other works proposed a reverse scheme where register allocation was made sensitive to not change the minimum initiation interval for the scheduler to expose sufficient ILP [START_REF] Touati | Early Periodic Register Allocation on ILP Processors[END_REF]. Most of the current mainstream open-source compilers [START_REF] Lattner | LLVM: A Compilation Framework for Lifelong Program Analysis & Transformation[END_REF][START_REF] Stallman | Using The GNU Compiler Collection: A GNU Manual For GCC Version 4[END_REF] Bench.

Reg 4. Spill metrics and performance in GFLOPS with NVCC on K40c device for different register configurations have adopted the first approach: when the register pressure is too high, the pre-pass list scheduling heuristics prioritize scheduling instructions that reduce the register pressure. However, such algorithms lack a global view, focusing only on the local register pressure at the current scheduled point. The associated optimization problem is NP-hard, and it is known that heuristics implemented in production compilers perform quite poorly on long straight line code [START_REF] Kral | FFT Compiler Techniques[END_REF], such as loop-body of highly unrolled loops. This observation motivated developers of auto-tuned libraries [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF][START_REF] Goto | Anatomy of Highperformance Matrix Multiplication[END_REF] to consider specific properties of the computational DAG to generate codelets that expose good locality for register reuse at source level: for certain domain-specific applications like FFT, a scheduling that minimizes the spill volume is well understood [START_REF] Jia-Wei | I/O Complexity: The Red-blue Pebble Game[END_REF]. The problem addressed in this paper belongs to a similar category, where one has to optimize register reuse for long straight-line code arising from domain-specific kernels. The main difference is that for the computational DAGs considered by our framework, optimal or nearly-optimal scheduling is unknown. Our proposed heuristic is a solution to address the optimization of such DAGs. Prior work on code generation for expression trees [START_REF] Aho | Optimal Code Generation for Expression Trees[END_REF][START_REF] Aho | Code Generation for Expressions with Common Subexpressions[END_REF][START_REF] Sethi | The Generation of Optimal Code for Arithmetic Expressions[END_REF] were discussed in details in Section 3. We now discuss some related work on combined register allocation and instruction scheduling, and register optimizations for high-order stencils.

Integrated Register Allocation and Instruction Scheduling Motwani et al. [START_REF] Motwani | Combining Register Allocation and Instruction Scheduling[END_REF] show that integrated register allocation and instruction scheduling is NP-hard. They propose a combined heuristic that provides relative weights for controlling register pressure and instruction parallelism. For instructions where the register of an operand can be used for the result, Govindrajan et al. [START_REF] Ramaswamy Govindarajan | Minimum Register Instruction Sequence Problem: Revisiting Optimal Code Generation for DAGs[END_REF] try to generate an instruction sequence from data dependence graph that is optimal in register usage. Berson et al. [START_REF] Berson | Integrated Instruction Scheduling and Register Allocation Techniques[END_REF] use register reuse DAGs to identify instructions whose parallel scheduling will require more resources than available, and optimize them to reduce their resource demands. Pinter [START_REF] Shlomit | Register Allocation with Instruction Scheduling[END_REF] describes an algorithm that colors a parallel interference graph to obtain a register allocation that does not introduce false dependences, and therefore exploits maximal parallelism. Norris et al. [START_REF] Norris | A scheduler-sensitive global register allocator[END_REF] propose an algorithm that constructs an interference graph with all feasible schedules, and then removes interferences for schedules that are to be least likely followed.

While these efforts consider integrated register management and instruction scheduling, the goals are very different and the contexts quite dissimilar. Prior work in this category has focused on maximizing parallelism without significantly increasing the MAXLIVE. In our context, the main reason for reordering instructions is to effectively exploit the significant potential for many reuses of values held in registers, while reducing the MAXLIVE.

Register Optimizations for High-Order Stencils Stock et al. [START_REF] Stock | A Framework for Enhancing Data Reuse via Associative Reordering[END_REF] identify a performance issue with register reuse for iterative stencils by noting that even though the computation becomes less memory-bound with increase in the stencil order, their register pressure worsens. They use a generalized version of [START_REF] De | Introducing the Semi-stencil Algorithm[END_REF] to interleave the additive contributions in an unrolled computation. Their approach uses retiming to homogenize stencil accesses and reduce register pressure, whereas our approach is based on a very different computational abstraction, and is more broadly applicable.

Basu et al. [START_REF] Basu | Compiler-Directed Transformation for Higher-Order Stencils[END_REF] propose a partial sum optimization implemented within the CHiLL compiler [START_REF] Hall | Loop Transformation Recipes for Code Generation and Auto-tuning[END_REF]. The partial sums are computed over planes for 3D stencils, and redundant computation is eliminated by performing array common subexpression elimination (CSE) [START_REF] Deitz | Eliminating Redundancies in Sum-of-product Array Computations[END_REF]. However, this optimization is only applicable to stencils with constant and symmetrical coefficients. While their work does not claim to reduce register pressure, it may do so as a consequence of array CSE. Based on the concepts of rematerialization [START_REF] Briggs | Rematerialization[END_REF], Jin et al. [START_REF] Jin | Libra: An Automated Code Generation and Tuning Framework for Registerlimited Stencils on GPUs[END_REF] propose a code generation framework that trades off recomputations for reduction in register pressure. It uses dynamic programming to iteratively determine the minimum amount of recomputations required to reduce the register consumption by one. This approach is limited to the stencils described in [START_REF] Jin | Libra: An Automated Code Generation and Tuning Framework for Registerlimited Stencils on GPUs[END_REF], where the recomputation of an expression does not increase the live ranges of the values involved in it.

In summary, existing approaches targeting register optimization for high-order stencils do not generalize well. Unlike these approaches, our framework optimizes both iterative and more general multi-statement stencils.

Conclusion

Despite a rich literature on register allocation and instruction scheduling, current production compilers are not sufficiently effective in reducing register pressure for compute-intensive high-order stencil codes. For such codes, register spills are a major performance limiter. Unfortunately, the compiler fails to perform an instruction reordering that can relieve register pressure, and the reordering it does perform to increase ILP often increases register pressure.

This paper presents a register optimization framework for multi-statement, high-order stencils, which views such computations as a collection of trees with significant data reuse across nodes, and systematically attempts to reduce register pressure by decreasing the simultaneously live ranges. Just as pattern-specific optimization techniques have demonstrably been more beneficial over traditional compiler optimizations for stencil computations, we show through this work that a specialized register management framework can be highly beneficial for stencil computations.

Figure 1 .

 1 Figure 1. Comparing same stencil computation with different sweeping order 2 Background and Motivation

  10 j3d7pt ( out , in , a , b , c ); 11 return out ; out = a + (b * c[i]) + d[i] + ((e[i] * f) / 2.3); (a) Illustrative stencil statement (b) Expression tree (c) Expression tree with accumulations

Figure 2 .

 2 Figure 2. Expression tree example

( a )Figure 3 .

 a3 Figure 3. Scheduling a tree with data sharing

( a )

 a Tree with context at root (b) Optimal schedule

Figure 4 .

 4 Figure 4. Example: computing su ′ (π )

Theorem 3 . 2 .

 32 For two clusters c 1 and c 2 such that |in(c 1 )| > |out (c 1 )| and |in(c 2 )| ≤ |out (c 2 )|, c 1 must be evaluated before c 2 in the optimal schedule.

Theorem 3 . 3 .

 33 For two clusters c 1 and c 2 such that |in(c 1 )| > |out (c 1 )| and |in(c 2 )| > |out (c 2 )|, the one with smaller su ′ must be prioritized for evaluation in the optimal schedule.

Algorithm 1 : 8 else 9 compute 10 π

 18910 Schedule-Tree (n, in, out) Input : A tree rooted at n with live-in/out contexts in and out Output : An optimal schedule S for the tree1 sched_cost ← ∅, S ← ∅; 2 C ← cr eat e_max imal _clust er s (n); (Sec. 3.2) 3 for each cluster c in C do 4 if |c | is 1 then 5 in(c ) ← in(n);6 out (c ) ← computed using equation 2; 7 sched_cost[c] ← su ′ (c ); in and out for each tree in c; (Sec. 3.2) ← all permutations of the trees in c; 11 sched_cost[c] ← su ′ (π ); 12 P ← sequence clusters using sched_cost and Thms. 3.1, 3.2, 3.3; 13 for each subtree t s in the sequence described by P do 14 append the schedule for t s in S ; 15 return S ;

Figure 5 .

 5 Figure 5. Example: interleaving to reduce MAXLIVE

  1 D ′ ← D; fusion_feasible ← t r ue; tree_order ← ∅; S ← ∅; 2 while fusion_feasible do 3 for each pair of transitive dep-free nodes t i ,t j in D ′ do 4 M ∪ = comput e_met r ic (D ′ , t i , t j ); (sec.3.4)

	5	sor t _descendinд (M );
	6	(t p , t q , fusion_feasible) ← f ind_f usion_candidat e (M );
	7	fuse t p and t q ;
	8	updat e_dependence_edдes (D ′ , t p , t q );
	9 for each node d in D ′ do
	10	append the tree sequence of d in tree_order;
	11 split_versions ← cr eat e_spl it _ver sions (tree_order);
	12 for each split in split_versions do
	13	S ′ ← ∅;
	14	for each kernel k in split do
	15	for each tree t in k do

16 compute in and out for t ; 17 append output of Schedule-Tree(t, in, out ) to S ′ ; 18 execute S ′ after compiling it with register limit R; 19 S ← S ′ if S ′ is a faster schedule than S , or if S is ∅; 20 return S ;

Table 2 .

 2 Metrics for 3d125pt for tuned configurations

	Metrics	rhs4th3fort	hypterm	derivative
		maxfuse	split-3	maxfuse	split-3	maxfuse	split-2
	Inst. Exec.	8.52E+9	8.25E+8	7.48E+8	7.71E+8	8.79E+8	8.96E+8
	IPC	1.07	1.11	0.97	1.06	1.02	1.14
	DRAM reads	9.07E+7	1.65E+8	1.57E+8	1.77E+8	1.34E+8	2.47E+8
	ldst exec.	1.55E+8	1.08E+8	1.27E+8	1.46E+8	1.45E+8	1.30E+8
	FLOPs	1.73E+10 1.81E+10	9.66E+9	9.36E+9	1.28E+10 1.34E+10
	tex txn.	1.11E+9	8.24E+8	9.72E+8	1.06E+9	1.14E+9	1.01E+9
	l2 read txn.	4.64E+8	3.79E+8	6.52E+8	5.90E+8	4.97E+8	4.51E+8
	GFLOPS	237.16	274.52	140.71	155.02	168.27	182.83

Table 3 .

 3 Metrics for reordered max-fuse vs. split versions

	registers R22,R30, and the remaining instructions. Indepen-
	dent FMAs are scheduled together without increasing the
	MAXLIVE. This reduces register pressure without compro-
	mising ILP. Therefore, even though the unrolled version
	performs fewer FLOPs than the reordered version, we incur
	less spill LDL/STL instructions per thread (101 for unrolled
	vs. 7 for reordered).

In a spill-free model of the computation, a data element is loaded only once into a register for all its uses/defs.

Acknowledgments

We thank the anonymous reviewers for their feedback and suggestions that helped improve the paper. This work was supported in part by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration, and by the U.S. National Science Foundation (NSF) through awards 1440749 and 1513120.

A Artifact appendix

Submission and reviewing guidelines and methodology: http://cTuning.org/ae/submission.html

A.1 Abstract

The artifact comprises an automated framework to perform reordering optimization on stencil computations; the algorithmic details of the framework are described in the PPoPP ' Hat Enterprise Linux Server release 6.7 operating system. • Hardware: We recommend a linux platform, and a GPU device with compute capability >= 3.5. • Output: GFLOPS for all the input benchmarks.

• Publicly available?: Yes.

A.3 Description

A.3.1 How delivered

The framework is open-source, and is available for download from the git repositoryhttps://github.com/pssrawat/ppopp-artifact. The downloaded package comprises the source code, the benchmarks, and the evaluation instructions and scripts. All the files in the repository are licensed to The Ohio State University.

A.3.2 Hardware dependencies

The framework has been tested on Ubuntu 16.04 and Red Hat Enterprise Linux Server release 6.7. The generated code can be executed on any NVidia device with compute capability >=3.5. For reproducing the results reported in the paper, we suggest using Kepler K40c device.

A.3.3 Software dependencies

flex version >= 2.6.0 (2.

A.4 Installation

First clone the artifact source to a local machine:

$. git clone https://github.com/pssrawat/ppopp-artifact Then compile the source code to create the executable test: $. cd ppopp-artifact $. make all Now set some environmental variables that are necessary for benchmarking. The first variable is to indicate the path to CUDA installation, and the second variable is to identify the compute capability of the GPU device:

$. export CUDAHOME=path-to-CUDA-installation $. export CAPABILITY=capability-of-GPU-card The compute capability of the K40c device is 35; for a Pascal device, it is 60.

To run the benchmarks: $. cd examples $. ./run-benchmarks.sh

The computed GFLOPS for all the benchmarks will be redirected to the output file output.txt in the examples directory.

A.5 Evaluation and expected result

The performance for each stencil benchmark in GFLOPS will be in output.txt after the evaluation script successfully finishes.

A.6 Experiment customization

The unrolling factors can be changed in the input .cu file, along with the launch paramaters, and the reordered versions regenerated by using the reorder.sh script provided with each benchmark. The documentation in docs/ folder provides additional details about adding new benchmarks, and optimizing them with the framework.