1201.3496v1

## Optimizing the Performance of Streaming Numerical Kernels on the IBM Blue Gene/P PowerPC 450 Processor Tareq M. Malas 1 , Aron J. Ahmadia 1 , Jed Brown 2 , John A. Gunnels 3 , David E. Keyes 1 1 King Abdullah University of Science and Technology Thuwal, Saudi Arabia 2 Argonne National Laboratory 9700 South Cass Avenue, Argonne, IL 60439, USA IBM T.J. Watson Research Center 1101 Kitchawan Avenue, Yorktown Heights, NY 10598, USA ## Abstract Several emerging petascale architectures use energy-efficient processors with vectorized computational units and in-order thread processing. On these architectures the sustained performance of streaming numerical kernels, ubiquitous in the solution of partial differential equations, represents a challenge despite the regularity of memory access. Sophisticated optimization techniques are required to fully utilize the Central Processing Unit (CPU). We propose a new method for constructing streaming numerical kernels using a highlevel assembly synthesis and optimization framework. We describe an implementation of this method in Python targeting the IBM Blue Gene/P supercomputer's PowerPC 450 core. This paper details the high-level design, construction, simulation, verification, and analysis of these kernels utilizing a subset of the CPU's instruction set. Wedemonstrate the effectiveness of our approach by implementing several three-dimensional stencil kernels over a variety of cached memory scenarios and analyzing the mechanically scheduled variants, including a 27-point stencil achieving a 1.7x speedup over the best previously published results. keywords : High Performance Computing, Performance Optimization, Code Generation, SIMD, Blue Gene/P ## 1 Introduction ## 1.1 Motivation As computational science soars past the petascale to exascale, a large number of applications continue to achieve disappointingly small fractions of the sustained performance capability of emerging architectures. In many cases this shortcoming in performance stems from issues at the single processor or thread level. The 'many-core' revolution brings simpler, slower, more power-efficient cores with vectorized floating point units in large numbers to a single processor. Performance improvements on such systems should be multiplicative with improvements derived from processor scaling. The characteristics of Blue Gene/P that motivate this work will persist in the processor cores of exascale systems, as one of the fundamental challenges of the exascale relative to petascale is electrical power [13, 25]. An optimistic benchmark (goal) for a petascale system is the continuous consumption of about a MegaWatt of electrical power, which represents the average continuous power consumption of roughly one thousand people in an OECD country (1.4 kW per person). Power reductions relative to delivered flop/s of factors of one to two orders of magnitude are expected en route to the exascale, which means more threads instead of faster-executing threads, power growing roughly as the cube of the clock frequency. It also means much less memory and memory bandwidth per thread because the movement of data over copper interconnects consumes much more power than the operations on data in registers. Mathematical formulations of problems and algorithms to implement them will be rewritten to increase arithmetic intensity in order to avoid data movement. Implementations in hardware will have to be made without some of today's popular powerintensive performance optimizations. At each stage of such a design, tradeoffs are made that complicate performance optimization for the compiler and programmer in exchange for improved power and die efficiency from the hardware. For example, out-of-order execution logic allows a microprocessor to reorder instruction execution on the fly to avoid pipeline and data hazards. When out-of-order logic is removed to save silicon and power, the responsibility for avoiding these hazards is returned to the compiler and programmer. In the same vein, a wide vector processing unit can significantly augment the floating point performance capabilities of a processing core at the expense of the efficient single-element mappings between input and output in scalar algorithms. Such vector units also incur greater bandwidth demands for a given level of performance, as measured by percentage of theoretical peak. Graphics processing units provide a set of compute semantics similar to a traditional Single Instruction Multiple Data (SIMD) vector processing unit with Single Instruction Multiple Thread (SIMT), but still execute in-order and require vectorized instruction interlacing to achieve optimal performance. We observe a broad trend to improve efficiency of performance with wider vector units and inorder execution units in the architectures of the IBM Blue Gene/P PowerPC 450 [38], the Cell Broadband Engine Architecture [10], Intel's MIC architecture [27], and NVIDIA's Tesla Graphics Processing Unit (GPU) [29]. In general terms, we may expect less memory per thread, less memory bandwidth per thread, and more threads per fixed data set size, creating an emphasis on strong scaling within a shared-memory unit. We also foresee larger grain sizes of SIMDization and high penalization of reads and writes from main memory as we move towards exascale. ## 1.2 Background We define streaming numerical kernels as small, cycle-intensive regions of a program where, for a given n bytes of data accessed in a sequential fashion, O ( n ) computational operations are required. Streaming numerical kernels are generally considered to be memory-bandwidth bound on most common architectures due to their low arithmetic intensity. The actual performance picture is substantially more complicated on high performance computing microprocessors, with constraints on computational performance stemming from such disparate sources as software limitations in the expressiveness of standard C and Fortran when targeting SIMD processors and a host of hardware bottlenecks and constraints, from the number of available floating point registers, to the available instructions for streaming memory into SIMD registers, to the latency and throughput of buses between the multiple levels of the memory hierarchy. In this paper we focus upon stencil operators, a subset of streaming kernels that define computations performed over a local neighborhood of points in a spatial multi-dimensional grid. Stencil operators are commonly found in partial differential equation solver codes in the role of finite-difference discretizations of continuous differential operators. Perhaps the most well-known of these is the 7-point stencil, which usually arises as a finite difference discretization of the Laplace kernel on structured grids. Although adaptive numerical methods and discretization schemes have diminished this operator's relevance for many problems as a full-fledged numerical solver, it is still a cycle-intensive subcomponent in several important scientific applications such as Krylov iterations of Poisson terms in pressure corrections, gravitation, electrostatics, and wave propagation on uniform grids, as well as block on adaptive mesh refinement methods [4]. We also target the 7-point stencil's 'boxier' cousin, the 27-point stencil. The 27-point stencil arises when cross-terms are needed such as in the NASA Advanced Supercomputing (NAS) parallel Multi Grid (MG) benchmark, which solves a Poisson problem using a V-cycle multigrid method with the stencil operator [3]. Finally, we examine the 3-point stencil, the one-dimensional analogue to the 7-point stencil and an important sub-kernel in our analysis. We concentrate on the Blue Gene/P architecture for a number of reasons. The forwardlooking design of the Blue Gene series, with its power-saving and ultra-scaling properties exemplifies some of the characteristics that will be common in exascale systems. Indeed, successor BlueGene/Q now tops the GREEN500 [16] list. Blue Gene/P has SIMD registers, a slow clock rate (850 MHz), an in-order, narrow superscalar execution path, and is constructed to be highly reliable and power-efficient, holding 15 of the top 25 slots of the GREEN500 list as recently as November 2009. Blue Gene continues to generate cutting-edge science, as evidenced by the continued presence of Blue Gene systems as winners and finalists in the Gordon Bell competition [1, 34, 18]. In the work presented here our simulations and performance enhancements focus on using code synthesis and scheduling to increase arithmetic intensity with unroll-and-jam, an optimization technique that creates tiling on multiply nested loops through a two-step procedure as in [5, 6]. Unroll-and-jam combines two well-known techniques, loop unrolling on the outer loops to create multiple inner loops, then loop fusion, or 'jamming,' to combine the inner loops into a single loop. Figure 1 shows an example of unroll-and-jam applied to a copy operation between two three-dimensional arrays, where each of the two outer most loops is unrolled once. This technique can work well for three-dimensional local operators because it promotes register reuse and can increase effective arithmetic intensity, though it requires careful register management and instruction scheduling to work effectively. ``` ``` - (a) Regular loop - (b) Unrolled-and-jammed loop ``` ``` Figure 1: Unroll-and-jam example for a three-dimensional array copy operation ## 2 Related Work Several emerging frameworks are facilitating the development of efficient high performance code without having to go down to the assembly level, at least not directly. These frameworks are largely motivated by the difficulties involved in utilizing vectorized Floating Point Unit (FPU) and other advanced features in the processor. CorePy [30], a Python implementation similar to our approach, provides a code synthesis package with an API to develop high performance applications by utilizing the low-level features of the processor that are usually hidden by the programming languages. Intel has introduced a new dynamic compilation framework, Array Building Blocks (ArBB) [31], which represents/enables a high level approach to automatically using the SIMD units on Intel processors. In addition, several techniques are developed in the literature to address the alignment problems in utilizing the SIMD capabilities of modern processors. In [15], an algorithm is introduced to reorganize the data in the registers to satisfy the alignment constraints of the processor. A compilation technique for data layout transformation is proposed in [20] that reorganizes the data statically in memory for minimal alignment conflicts. As stencil operators in particular have been identified as an important component of many scientific computing applications, a good deal of effort has been spent attempting to improve their performance through optimization techniques. Christen et al. optimized a 7-point stencil in three-dimensional grids on the CELL BE processor and a GPU system in [9]. On CELL BE, they reduced bandwidth requirements through spacial and temporal blocking. To improve the computation performance, they utilized the SIMD unit with shuffling intrinsics to handle alignment issues. They utilized optimization techniques including preloading, loop unrolling, and instruction interleaving. Rivera and Tseng [35] utilized tiling in space to improve spatial locality and performance of stencil computations. A domain-specific technique is time skewing [28] [26]. Unfortunately, time skewing is not generalizable because no other computations between time steps are allowed to occur. Recently, Nguyen et al. [32] optimized the performance of a 7-point stencil and a Lattice Boltzmann stencil on CPUs and GPUs over three-dimensional grids by performing a combination of spatial and temporal blocking, which they dubbed 3.5D blocking, to decrease the memory bandwidth requirements of their memory bound problems. Wellein performed temporal blocking on the thread level to improve the performance of stencil computations on multicore architectures with shared caches in [40]. Kamil conducted a study on the impact of modern memory subsystems on three-dimensional stencils [24] and proposed two cache optimization strategies for stencil computations in [23], namely cache-oblivious and cache-aware techniques. Several other recent studies in performance improvements in stencil operators, including multilevel optimization techniques [14, 33]. There are several approaches to obtaining performance improvements in stencil computations that are more automated, requiring less manual intervention and resulting in more generality: Christen et al. introduced a stencil code generation and auto-tuning framework [8], PATUS, targeting CPUs and GPUs. Williams presented an auto-tuning frameworks in [42]. They performed their optimization work on a Lattice Boltzmann application over several architectures. Kamil built an auto-tuning framework for code generation in [22]. Their framework accepts the stencil's kernel in Fortran and then converts it to a tuned version in Fortran, C, or Compute Unified Device Architecture (CUDA). In [37], a software synthesis approach was proposed to generate optimized stencil code. Machine Learning strategies were proposed in [17] to tune the parameters of the optimization techniques for the 7- and the 27-point stencil. Recently, Tang introduced Pochoir stencil compiler in [39]. Their framework aims to simplify programming efficient stencil codes, utilizing parallel cache oblivious algorithms. We call special attention to a comprehensive work by Datta [12], who constructed an auto-tuning framework to optimize the 7- and the 27-point stencils and the Gauss-Seidel Red-Black Helmholtz kernel. Datta's work was performed on diverse multicore architectures modern at the time. He achieved impressive performance by employing a search over a variety of algorithmic and architecture-targeting techniques including common subexpression elimination and Non-Uniform Memory Access (NUMA)-aware allocations. It is interesting to note that aside from register blocking and common subexpression elimination, Datta's techniques were ineffective in improving the performance of stencil operators on the Blue Gene/P platform. This was attributed in part to the difficulties in achieving good performance for the in-order-execution architecture of the PowerPC 450 processor. ## 3 Implementation Considerations ## 3.1 Stencil Operators A three-dimensional stencil is a linear operator on R MNP , the space of scalar fields on a Cartesian grid of dimension M × N × P . Apart from some remarks in Section 6, we assume throughout this paper that the operator does not vary with the location in the grid, as is typical for problems with regular mesh spacing and space-invariant physical properties such as constant diffusion. The input, A , and output, R , are conventionally stored in a onedimensional array using lexicographic ordering. We choose a C-style ordering convention so that an entry a i,j,k of A has flattened index ( iN + j ) P + k , with zero-based indexing. Further, we assume that the arrays A and R are aligned to 16-byte memory boundaries. Formally, the 3-point stencil operator defines a linear mapping from a weighted sum of three consecutive elements of A to one element in R : <!-- formula-not-decoded --> We further assume certain symmetries in the stencils that are typical of self-adjoint problems, e.g. w 0 = w 2 . The effect of this assumption is that fewer registers are required for storing the w coefficients of the operator, allowing us to unroll the problem further. This assumption is not universally applicable to numerical schemes such as upwinding, where adaptive weights are used to capture the direction of flow. The 7-point stencil (Figure 2) defines the result at r i,j,k as a linear combination of the input a i,j,k and its six three-dimensional neighbors with Manhattan distance one. The 27point stencil uses a linear combination of the set of 26 neighbors with Chebyshev distance one. The boundary values r i,j,k with i ∈ { 0 , M -1 } , j ∈ { 0 , N -1 } , or k ∈ { 0 , P -1 } are not written as is standard for Dirichlet boundary conditions. Other boundary conditions would apply a different one-sided stencil at these location, a computationally inexpensive consideration that we do not regard for our experiments. Figure 2: 7-point stencil operator <details> <summary>Image 1 Details</summary>  ### Visual Description ## 3D Coordinate System Diagram: Spatial Orientation and Vector Representation ### Overview The image depicts a 3D Cartesian coordinate system with labeled axes (i, j, k) and discrete points positioned along the axes. A yellow plane is oriented perpendicular to the k-axis, intersecting the origin. Points are color-coded and labeled with numerical identifiers (0, 1, 2, 3), suggesting a focus on vector positioning and spatial relationships. --- ### Components/Axes - **Axes**: - **i-axis**: Vertical axis (labeled "i"), oriented upward. - **j-axis**: Horizontal axis (labeled "j"), oriented rightward. - **k-axis**: Horizontal axis (labeled "k"), oriented leftward (orthogonal to i and j). - **Origin**: Labeled "0" at the intersection of all three axes. - **Yellow Plane**: A shaded plane perpendicular to the k-axis, spanning the i-j plane (analogous to the xy-plane in standard 3D systems). - **Points**: - **Point 0**: White dot at the origin (0,0,0). - **Point 1**: Gray dot at (1,0,0) on the i-axis. - **Point 2**: Red dot at (0,1,0) on the j-axis. - **Point 3**: Blue dot at (0,0,1) on the k-axis. --- ### Detailed Analysis - **Axis Labels**: All axes are explicitly labeled with lowercase letters (i, j, k), deviating from the conventional (x, y, z) notation. This suggests a domain-specific or pedagogical context. - **Point Placement**: - Point 1 (gray) is positioned 1 unit along the i-axis. - Point 2 (red) is positioned 1 unit along the j-axis. - Point 3 (blue) is positioned 1 unit along the k-axis. - All points are unit vectors aligned with their respective axes. - **Plane Orientation**: The yellow plane is anchored at the origin and extends infinitely along the i and j axes, confirming its role as the i-j plane. --- ### Key Observations 1. **Color-Coding**: Points 1 and 2 share the same red/gray color scheme, while Point 3 is distinctly blue, potentially indicating a categorical distinction (e.g., magnitude, type, or group). 2. **Unit Vectors**: All labeled points (1, 2, 3) represent unit vectors, as they are equidistant (1 unit) from the origin along their respective axes. 3. **Plane Significance**: The yellow plane emphasizes the i-j plane, which is critical for understanding 2D projections in 3D space. --- ### Interpretation This diagram illustrates the foundational structure of a 3D coordinate system, emphasizing: - **Orthogonality**: The i, j, and k axes are mutually perpendicular, forming a right-handed system. - **Vector Representation**: Points 1, 2, and 3 correspond to the standard basis vectors **i**, **j**, and **k**, respectively. These vectors form the basis for expressing any 3D vector as a linear combination. - **Pedagogical Purpose**: The use of non-standard axis labels (i, j, k) and color-coded points suggests this diagram may be used to teach vector algebra, spatial reasoning, or coordinate transformations. The absence of numerical scales beyond unit increments implies a focus on conceptual understanding rather than quantitative analysis. The yellow plane reinforces the importance of 2D subspaces within 3D geometry. </details> The 27-point stencil (Figure 3) can be seen as the summation over nine independent 3-point stencil operators into a single result. We assume symmetry along but not between the three dimensions, leading to 8 unique weight coefficients. The symmetric 7-point stencil operator has 4 unique weight coefficients. ## 3.2 PowerPC 450 Designed for delivering power-efficient floating point computations, the nodes in a Blue Gene/P system are four-way Symmetric Processors (SMPs) comprised of PowerPC 450 processing cores [21]. The processing cores possess a modest superscalar architecture capable of issuing a SIMD floating point instruction in parallel with various integer and load/store instructions. Each core has an independent register file containing 32 4-byte general purpose registers and 32 16-byte SIMD floating point registers which are operated on by a pair of fused floating point units. A multiplexing unit on each end of this chained floating point pipeline enables a rich combination of parallel, copy, and cross semantics in the SIMD floating point operation set. This flexibility in the multiplexing unit can enhance computational efficiency by replacing the need for independent copy and swap operations on the floating point registers with single operations. To provide backward compatibility as well as some additional functionality, these 16-byte SIMD floating point registers are divided into independently addressable, 8-byte primary and secondary registers wherein non-SIMD floating point instructions operate transparently on the primary half of each SIMD register. An individual PowerPC 450 core has its own 64KB L1 cache, divided evenly into a 32KB instruction cache and a 32KB data cache. The L1 data cache uses a round-robin (First In First Out (FIFO)) replacement policy in 16 sets, each with 64-way set-associativity. Each L1 cache line is 32 bytes in size. Every core also has its own private prefetch unit, designated as the L2 cache, 'between' the L1 and the L3. In the default configuration, each PowerPC 450 core can support up to 5 'deep fetching' streams or up to 7 shallower streams. These values stem from the fact that the L2 prefetch cache has 15 128-byte entries. If the system is fetching two lines ahead (settings are configured on job startup), each stream occupies three positions, one current and two 'future,' while a shallower prefetch lowers the occupancy to two per stream. The final level of cache is the 8MB L3, shared among the four cores. The L3 features a Least Recently Used (LRU) replacement policy, with 8-way set associativity and a 128-byte line size. Figure 3: 27-point stencil operator <details> <summary>Image 2 Details</summary>  ### Visual Description ## 3D Cube Diagram: Labeled Nodes and Colored Data Points ### Overview The image depicts a 3D cube structure divided into three sections along the **i-axis** (x-axis), each containing a 2x2 grid of points. The cube is oriented with axes labeled **i** (vertical), **j** (horizontal), and **k** (depth). Vertices are labeled with numbers 0–7, and colored points (red, blue, green, purple, brown, black) are distributed across the cube. A central white point is marked at the intersection of the three axes. ### Components/Axes - **Axes**: - **i-axis (vertical)**: Labeled with an upward arrow. - **j-axis (horizontal)**: Labeled with a rightward arrow. - **k-axis (depth)**: Labeled with a forward arrow. - **Vertices**: - Labeled with numbers **0–7** at the corners of the cube. - Example: Vertex (0,0,0) = 0, (2,0,0) = 1, (2,2,0) = 2, (0,2,0) = 3, (0,0,2) = 4, (2,0,2) = 5, (2,2,2) = 6, (0,2,2) = 7. - **Colored Points**: - **Red**: (0,0,0), (2,0,0), (2,2,0), (0,2,0) - **Blue**: (0,0,2), (2,0,2), (2,2,2), (0,2,2) - **Green**: (0,0,4), (2,0,4), (2,2,4), (0,2,4) - **Purple**: (0,0,6), (2,0,6), (2,2,6), (0,2,6) - **Brown**: (0,0,8), (2,0,8), (2,2,8), (0,2,8) - **Black**: (0,0,10), (2,0,10), (2,2,10), (0,2,10) - **White**: Central point at (2,1,2). ### Detailed Analysis - **Axis Labels**: - The **i-axis** (vertical) spans from 0 to 4, with the cube divided into three sections (i=0, 2, 4). - The **j-axis** (horizontal) spans 0–2, with grid lines dividing the cube into 2x2 sections. - The **k-axis** (depth) spans 0–10, with points at k=0, 2, 4, 6, 8, 10. - **Vertex Labels**: - Numbers 0–7 are assigned to the 8 vertices of the cube, likely representing unique identifiers or positions. - The central white point (2,1,2) is not labeled with a number but is distinct in color. - **Color Distribution**: - Colors are distributed along the **k-axis**, with each color corresponding to a specific k-value: - Red (k=0), Blue (k=2), Green (k=4), Purple (k=6), Brown (k=8), Black (k=10). - This suggests a categorical or sequential relationship between color and the k-axis. ### Key Observations 1. **Structural Division**: The cube is divided into three equal sections along the i-axis, each containing a 2x2 grid of points. </details> On this architecture the desired scenario in highly performant numerical codes is the dispatch of a SIMD floating point instruction every cycle (in particular, a SIMD Floating point Multiply-Add (FMA)), with any load or store involving one of the floating point registers issued in parallel, as inputs are streamed in and results are streamed out. Floating point instructions can be retired one per cycle, yielding a peak computational throughput of one (SIMD) FMA per cycle, leading to a theoretical peak of 3.4 GFlops/s per core. Blue Gene/P's floating point load instructions, whether they be SIMD or non-SIMD, can be retired every other cycle, leading to an effective read bandwidth to the L1 of 8 bytes a cycle for aligned 16-byte SIMD loads (non-aligned loads result in a significant performance penalty) and 4 bytes a cycle otherwise. As a consequence of the instruction costs, no kernel can achieve peak floating point performance if it requires a ratio of load to SIMD floating point instructions greater than 0.5. It is important to ensure packed 'quad-word' SIMD loads occur on 16-byte aligned memory boundaries on the PowerPC 450 to avoid performance penalties that ensue from the hardware interrupt that results from misaligned loads or stores. An important consideration for achieving high throughput performance on modern floating point units is pipeline latency, the number of cycles that must transpire between accesses to an operand being written or loaded into in order to avoid pipeline hazards (and their consequent stalls). Floating point computations on the PowerPC 450 have a latency of 5 cycles, whereas double-precision loads from the L1 require at least 4 cycles and those from the L2 require approximately 15 cycles [38]. Latency measurements when fulfilling a load request from the L3 or DDR memory banks are less precise: in our performance modeling we assume an additional 50 cycle average latency penalty for all loads outside the L1 that hit in the L3. In the event of an L1 cache miss, up to three concurrent requests for memory beyond the L1 can execute (an L1 cache miss while 3 requests are 'in-flight' will cause a stall until one of the requests to the L1 has been fulfilled). Without assistance from the L2 cache, this leads to a return of 96 bytes (three 32-byte lines) every 56 cycles (50 cycles of memory latency + 6 cycles of instruction latency), for an effective bandwidth of approximately 1.7 bytes/cycle. This architectural characteristic is important in our work, as L3-confined kernels with a limited number of streams can effectively utilize the L2 prefetch cache and realize as much as 4.5 bytes/cycle bandwidth per core, while those not so constrained will pay the indicated bandwidth penalty. The PowerPC 450 is an in-order unit with regards to floating point instruction execution. An important consequence is that a poorly implemented instruction stream featuring many non-interleaved load/store or floating point operations will suffer from frequent structural hazard stalls with utilization of only one of the units. Conversely, this in-order nature makes the result of efforts to schedule and bundle instructions easier to understand and extend. ## 3.3 Instruction Scheduling Optimization We wish to minimize the required number of cycles to execute a given code block composed of PowerPC 450 assembly instructions. This requires scheduling (reordering) the instructions of the code block to avoid the structural and data hazards described in the previous section. Although we use greedy heuristics in the current implementation of the code synthesis framework to schedule instructions, we can formulate the scheduling problem as an Integer Linear Programming (ILP) optimization problem. We base our formulation on [7], which considers optimizations combining register allocation and instruction scheduling of architectures with multi-issue pipelines. To account for multi-cycle instructions, we include parts of the formulation in [41]. We consider two minor extensions to these approaches. First, we consider the two separate sets of registers of the PowerPC 450, the General Purpose Registers (GPR), and the Floating Point Registers (FPR). Second, we account for instructions that use the Load Store Unit (LSU) occupying the pipeline for a varying number of cycles. We begin by considering a code block composed of N instructions initially ordered as I = { I 1 , I 2 , I 3 , ..., I N } . A common approach to represent the data dependencies of these instructions is to use a Directed Acyclic Graph (DAG). Given a DAG G ( V, E ), the nodes of the graph ( V ) represent the instructions and the directed edges ( E ) represent the dependencies between them. Figure 4 shows an example of a DAG representing a sequence of instructions with the edges representing data dependencies between instructions. Each read-after-write data dependency of an instruction I j on an instruction I i is represented by a weighted edge e ij . This weighted edge corresponds to the number of cycles needed by an instruction I i to produce the results required by an instruction I j . The weights associated with write-after-read and write-after-write data dependencies are set to one. Figure 4: Example of DAG representation of the instructions <details> <summary>Image 3 Details</summary>  ### Visual Description ## Diagram: Instruction Sequence and Data Flow Graph (DAG) Representation ### Overview The image contains two components: 1. **(a) Instruction Sequence**: A table defining 8 instructions (I₁–I₈) with labels and operations. 2. **(b) DAG Representation**: A directed acyclic graph (DAG) with 8 nodes (1–8) connected by edges labeled with instruction numbers (1–8), illustrating data flow dependencies. --- ### Components/Axes #### (a) Instruction Sequence Table - **Structure**: 8 rows (I₁–I₈), each with: - **Label**: Instruction identifier (e.g., `ld`, `sub`, `mul`, `st`). - **Operations**: Comma-separated operands (e.g., `fa,rx,ry`). - **Key Details**: - `ld`: Load operations (e.g., `I₁: ld fa,rx,ry`). - `sub`: Subtract (`I₃: sub fd,fa,fb`). - `mul`: Multiply (`I₅: mul fe,fb,fd`). - `st`: Store (`I₇: st fc,fd,fc`). #### (b) DAG Representation - **Nodes**: 8 nodes labeled 1–8. - **Edges**: Directed arrows with instruction numbers (1–8) indicating transitions. - Example: Node 1 → Node 2 (edge labeled `1`), Node 1 → Node 3 (edge labeled `4`). - **Legend**: Edge numbers correspond to instruction numbers in the table (e.g., `1` = `I₁`, `4` = `I₄`). --- ### Detailed Analysis #### (a) Instruction Sequence Table | Instruction | Label | Operations | |-------------|-------|---------------------| | I₁ | ld | fa, rx, ry | | I₂ | ld | fb, rx, ry | | I₃ | sub | fd, fa, fb | | I₄ | ld | fb, rt, ry | | I₅ | mul | fe, fb, fd | | I₆ | st | fe, rz, ry | | I₇ | mul | fc, fd, fc | | I₈ | st | fc, rv, ry | #### (b) DAG Representation - **Node Connections**: - **Node 1**: Edges to 2 (1), 3 (4), 4 (4), 7 (5). - **Node 2**: Edges to 3 (2), 5 (3), 6 (6). - **Node 3**: Edges to 4 (1), 5 (7). - **Node 4**: Edges to 5 (8). - **Node 5**: Edges to 6 (5). - **Node 6**: Edges to 7 (4). - **Node 7**: Edges to 8 (7). - **Node 8**: No outgoing edges. --- ### Key Observations 1. **Instruction Usage**: - `I₁` (load) and `I₄` (load) are the most frequently used (4 edges each). - `I₅` (mul) and `I₇` (mul) are used for intermediate computations. - `I₈` (st) is the final store operation. 2. **Data Flow**: - Nodes 1–5 represent computation stages, while nodes 6–8 handle storage. - Node 5 (final computation) feeds into Node 6 (storage). 3. **Dependencies**: - Node 1 initiates the flow, branching to multiple nodes via different instructions. - Node 8 is a terminal node with no outgoing edges. --- ### Interpretation - **Purpose**: The DAG visualizes the execution path of instructions, showing how data dependencies are resolved. - **Critical Path**: The sequence `I₁ → I₄ → I₅ → I₇ → I₈` represents the primary computation and storage flow. - **Anomalies**: - Node 1 has multiple outgoing edges (4), indicating parallelism or branching logic. - Node 3 and Node 4 both feed into Node 5, suggesting convergence of intermediate results. - **Significance**: The diagram highlights the interplay between load/store operations (`ld`, `st`) and arithmetic operations (`sub`, `mul`), emphasizing the importance of instruction scheduling in data flow optimization. --- **Note**: No numerical values or trends are present in the DAG; the focus is on structural relationships and instruction usage. </details> The PowerPC 450 processor has one LSU and one FPU. This allows the processor to execute a maximum of one floating point operation per cycle and one load/store operation every two cycles (each operation using the LSU consumes at least 2 cycles). The instructions use either the LSU ( I LSU ∈ LSU ) or the other computation resources of the processor including the FPU ( I FPU ∈ FPU ). A critical path, C , in the DAG is a path in the graph on which the sum of the edge weights attains the maximum. We can compute the lower bound ( L bound ) to run the instructions I as follows: <!-- formula-not-decoded --> We can compute an upper bound U bound ( I ) by simulating the number of cycles to execute the scheduled instructions. If U bound ( I ) = L bound ( I ) we can generate an optimal schedule. We define the Boolean optimization variables x j i to take the value 1 to represent that the instruction I i is scheduled at a given cycle c j and 0 otherwise. These variables are represented in the array shown in Figure 5, where M represents the total number of the cycles. Figure 5: A Boolean array of size N × M representing the scheduling variables | | c 1 | c 2 | c 3 | . | c M | |-----|-------|-------|-------|-----|-------| | I 1 | x 1 1 | x 2 1 | x 3 1 | . | x M 1 | | I 2 | x 1 2 | x 2 2 | x 3 2 | . | x M 2 | | I 3 | x 1 3 | x 2 3 | x 3 3 | . | x M 3 | | . | . | . | . | . | . | | I N | x 1 N | x 2 N | x 3 N | . | x M N | The first constraint in our optimization is to force each instruction to be scheduled only once in the code block. Formally: <!-- formula-not-decoded --> The PowerPC 450 processor can execute a maximum of one floating point operation every cycle and one load operation every two cycles (store instructions are more complicated, but for this derivation we assume two cycles as well). This imposes the following constraints: <!-- formula-not-decoded --> <!-- formula-not-decoded --> Finally, to maintain the correctness of the code's results we enforce the dependency constraints: <!-- formula-not-decoded --> The instruction latency to write to a GPR is 1 cycle ( e ij = 1). The latency for writing to a FPR is 5 cycles for I i ∈ I FPU and at least 4 cycles for I i ∈ I LSU . Load instructions have higher latency when the loaded data is present in the L3 cache or the RAM. This can be considered in future formulations to maximize the number of cycles between loading the data of an instruction I i and using it by maximizing e ij in the objective. We also wish to constrain the maximum number of allocated registers in a given code block. All the registers are assumed to be allocated and released within the code block. An instruction I i scheduled at a cycle c j allocates a register r by writing on it. The register r is released (deallocated) at a cycle c k by the last instruction reading from it. The life span of the register r is defined to be from cycle c j to cycle c k inclusive. We define two sets of Boolean register optimization variables g j i and f j i to represent the usage of the GPR and the FPR files, respectively. Each of these variables belongs to an array of the same size as the array in Figure 5. The value of g j i is set to 1 during the life span of a register in the GPR modified by the instruction I i scheduled at the cycle c j and last read by an instruction scheduled at the cycle c k , that is g z i = 1 when j ≤ z ≤ k and zero otherwise. The same applies to the variables f j i to represent the FPR register allocation. To compute the values of g j i and f j i , we define the temporary variables ˆ g p i and ˆ f p i . Let K i be the number of instructions reading the from the register allocated by the instruction I i . These temporary variables are computed as follows: <!-- formula-not-decoded --> <!-- formula-not-decoded --> Our optimization variables g j i and f j i will equal to 1 only when ˆ g j i > 0 and ˆ f j i > 0, respectively. This can be formulated as follow: <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> Now we can constrain the maximum number of used registers FPR max and GPR max by the following: <!-- formula-not-decoded --> <!-- formula-not-decoded --> Our optimization objective is to minimize the required cycles to execute the code ( C run ). The complete formulation of the optimization problem is: <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> <!-- formula-not-decoded --> The general form of this optimization problem is known to be NP-complete [19], although many efficient implementations exist for solving integer problems exist; no known algorithms can guarantee global solutions in polynomial time. In our implementation, we use a greedy algorithm that yields an approximate solution in practical time as we describe in 4.4. ## 4 Implementation ## 4.1 C and Fortran We implement the three streaming kernels in C and Fortran and utilize published results as a benchmark to gain a better understanding of the relative performance of the general class of streaming numerical kernels on the IBM PowerPC 450. Our first challenge is in expressing a SIMDized mapping of the stencil operator in C or Fortran. Neither language natively supports the concept of a unit of SIMD-packed doubles, though Fortran's complex type comes very close. Complicating matters further, the odd cardinality of the stencil operators necessitates careful tactics to efficiently SIMDize the code. Regardless of the unrolling strategy employed, one of every two loads is unaligned or requires clever register manipulation. It is our experience that standard C and Fortran implementations of our stencil operators will compile into exclusively scalar instructions that cannot attain greater than half of peak floating point performance. For example, with a na¨ ıve count of 53 flops per 27-point stencil, peak performance is 62 Mstencil/s and we observe 31.5 Mstencil/s on a single core of the PowerPC 450. Register pressure and pipeline stalls further reduce the true performance to 21.5 Mstencil/s. We note that Datta [11] was able to improve this to 25 Mstencil/s using manual unrolling, common subexpression elimination, and other optimizations within C. The XL compilers support intrinsics which permit the programmer to specify which instructions are used, but the correspondence between intrinsics and instructions is not exact enough for our purposes. The use of intrinsics can aid programmer productivity, as this method relies upon the compiler for such tasks as register allocation and instruction scheduling. However, our methods require precise control of these aspects of the produced assembly code, making this path unsuitable for our needs. Our code generation framework is shown in Figure 6. High-level Python code produces the instruction sequence of each code block. First, registers are allocated by assigning them variable names. Next, a list of instruction objects is created that utilize the variables to address registers. Enabling the instruction scheduler allows the simulator to execute the instructions out of order; otherwise they will be executed in order. The instruction simulator uses virtual GPR, FPR, and memory to simulate the pipeline execution and to simulate the expected results of the given instruction sequence. A log is produced by the instruction simulator to provide the simulation details, showing the cycles at which the instructions are scheduled and any encountered data and structural hazards. Also, the log can contain the contents of the GPR, the FPR, and the memory, at any cycle, to debug the code for results correctness. The C code generator takes the simulated instructions and generates their equivalent inline assembly code in a provided C code template using the inline assembly extension to C standard provided by GCC. For added clarity, we associate each generated line with a comment showing the mapping between the used registers numbers and their corresponding variables names in the Python code. Figure 6: The components of the code generation framework in this work <details> <summary>Image 4 Details</summary>  ### Visual Description ## Flowchart: System Architecture for PowerPC 450 Simulation and Code Generation ### Overview The diagram illustrates a workflow for simulating PowerPC 450 instructions using Python code and generating documented C code. The process involves register allocation, instruction scheduling, and memory management within a simulator, followed by C code generation. ### Components/Axes 1. **Python Code Block**: - **Register allocation**: Light blue box with dark blue border. - **Instructions (list of objects)**: Light blue box with dark blue border. 2. **PowerPC 450 Simulator Block**: - **GPR/FPR/Memory**: Light blue box with dark blue border. - **Instruction Scheduler and simulator**: Light blue box with dark blue border. - **C code generator**: Light blue box with dark blue border. 3. **Outputs**: - **Simulation log and debugging information**: Light gray box with dark gray border. - **Documented C code**: Light gray box with dark gray border. 4. **Inputs**: - **C code template**: Light gray box with dark gray border. ### Flow Direction - Arrows indicate the following sequence: 1. Python code → PowerPC 450 simulator (GPR/FPR/Memory → Instruction Scheduler). 2. Simulator output → C code generator. 3. C code generator output → Documented C code and Simulation log. ### Key Observations - The system bridges Python and C ecosystems via PowerPC simulation. - The simulator acts as an intermediary, processing Python-defined instructions and register allocations. - The C code generator produces both executable code and documentation, suggesting automated documentation generation. ### Interpretation This architecture demonstrates a toolchain for: 1. **Simulation**: Translating Python-based PowerPC instructions into executable cycles via GPR/FPR/Memory management. 2. **Code Generation**: Converting simulated logic into optimized C code with embedded documentation. 3. **Debugging**: Providing simulation logs to trace execution paths and identify bottlenecks. The bidirectional arrow between "Instruction Scheduler and simulator" and "C code generator" implies iterative refinement, where simulation results inform code optimization. The separation of "Simulation log and debugging information" suggests post-execution analysis capabilities, critical for low-level hardware emulation workflows. </details> ## 4.2 Kernel Design The challenges in writing fast kernels in C and Fortran motivate us to program at the assembly level, a (perhaps surprisingly) productive task when using our code synthesis framework: an experienced user was able to design, implement, and test several efficient kernels for a new stencil operator in one day using the framework. Much of our work is based on the design of two small and na¨ ıvely scheduled, 3-point kernels, dubbed mutate-mutate and load-copy, that we introduce in 4.2. The kernels distinguish themselves from each other by their relative balance between load/store and floating point cycles consumed per stencil. We find that efficiently utilizing SIMD units in stencil computations is a challenging task. To fill the SIMD registers, we pack two consecutive data elements from the fastest moving dimension, k , allowing us to compute two stencils simultaneously, as in [2]. Computing in this manner is semantically equivalent to an unrolling by 2 in the k direction. As a connventional notation, since i and j are static for any given stream, we denote the two values occupying the SIMD register for a given array by their k indices, e.g., SIMD register a 34 contains A i,j, 3 in its primary half and A i,j, 4 in its secondary half. Many stencil operators map a subset of adjacent A elements with odd cardinality to each R element, as is illustrated in the left half of Figure 7, which depicts the SIMD register contents and computations mid-stream of a 3-point kernel. The odd cardinality prevents a straightforward mapping from SIMD input registers to SIMD output registers. We note that aligned SIMD loads from A easily allow for the initialization of registers containing a 23 and a 45 . Similarly, the results in r 34 can be safely stored using a SIMD store. The register containing a 34 , unaligned elements common to the aligned registers containing adjacent data, requires a shuffle within the SIMD registers through the use of additional load or floating point move instructions. We introduce two kernels, which we designate as mutate-mutate (mm) and load-copy (lc), as two different approaches to form the packed data into the unaligned SIMD registers while streaming through A in memory. A 'mutate' operation is defined as the replacment of one operand of a SIMD register by a value loaded from memory. A 'load' operation is defined as the replacement of both operands in a SIMD register by a SIMD load from memory. A 'copy' operation is the replacement of one operand of a SIMD register by an operand from another SIMD register. Mutates and loads utilize the LSU, copies utilize the FPU. The mutate-mutate kernel replaces the older half of the SIMD register by with the next element of the stream. In our example in Figure 7 we start with a 23 loaded in the SIMD register, then after the first computation we load a 4 into the primary part of the SIMD register so that the full contents are a 43 . The load-copy kernel instead combines the unaligned values in a SIMD register from two consecutive aligned quad-words loaded in two SIMD registers by copying, the primary element from the second register to the primary element of the first. In our example a 23 and a 45 are loaded in two SIMD registers. Then, after the needed computations involving a 23 have been dispatched, a floating point move instruction replaces a 2 with a 4 to form a 43 . The two kernels use an identical set of floating point instructions, visually depicted in the right half of Figure 7, to accumulate the computations into the result registers. The two needed weight coefficients are packed into one SIMD register. The first floating point Figure 7: SIMD stencil computations on one-dimensional streams <details> <summary>Image 5 Details</summary>  ### Visual Description ## Diagram: Contribution Computation Process ### Overview The diagram illustrates a three-step computational process for calculating contributions from specific elements (a23, a43, a45) using weighted operations. Each step involves selecting numerical ranges, applying transformations (W0/W2, W1), and combining results through arithmetic operations (*=, +=). The process is visualized through horizontal rows (A_ij, R_ij), highlighted boxes, and directional arrows. ### Components/Axes - **Rows**: - **A_ij**: Top row with indices 0–6, labeled with "k" (horizontal axis). - **R_ij**: Bottom row with indices 0–6, also labeled with "k". - **Highlighted Boxes**: - **(a)**: Boxes around indices 2–3 in A_ij and 3–4 in R_ij. - **(b)**: Boxes around indices 3–4 in A_ij and 3–4 in R_ij. - **(c)**: Boxes around indices 4–5 in A_ij and 3–4 in R_ij. - **Legend (Right Side)**: - **Symbols**: - `*`: Represents multiplication (W0/W2). - `+=`: Represents addition (W1). - **Labels**: - **W0/W2**: Weighted operation for multiplication. - **W1**: Weighted operation for addition. - **A**: Intermediate result (A_ij). - **R**: Final result (R_ij). - **Operations**: - `*=`: Multiplication assignment. - `+=`: Addition assignment. ### Detailed Analysis #### Step (a): Compute contributions from a23 - **Highlighted Ranges**: - A_ij: Indices 2–3 (values 2, 3). - R_ij: Indices 3–4 (values 3, 4). - **Arrows**: - Blue arrows from A_ij[2–3] point to W0/W2 (multiplication). - Red arrows from R_ij[3–4] point to W1 (addition). - **Operations**: - A_ij[2–3] undergoes `*=` with W0/W2. - R_ij[3–4] undergoes `+=` with W1. #### Step (b): Compute contributions from a43 - **Highlighted Ranges**: - A_ij: Indices 3–4 (values 3, 4). - R_ij: Indices 3–4 (values 3, 4). - **Arrows**: - Blue arrows from A_ij[3–4] point to W0/W2. - Red arrows from R_ij[3–4] point to W1. - **Operations**: - A_ij[3–4] undergoes `*=` with W0/W2. - R_ij[3–4] undergoes `+=` with W1. #### Step (c): Compute contributions from a45 - **Highlighted Ranges**: - A_ij: Indices 4–5 (values 4, 5). - R_ij: Indices 3–4 (values 3, 4). - **Arrows**: - Blue arrows from A_ij[4–5] point to W0/W2. - Red arrows from R_ij[3–4] point to W1. - **Operations**: - A_ij[4–5] undergoes `*=` with W0/W2. - R_ij[3–4] undergoes `+=` with W1. ### Key Observations 1. **Sequential Processing**: Each step processes distinct ranges of A_ij and R_ij, with overlapping contributions in R_ij (indices 3–4 appear in all steps). 2. **Weighted Operations**: - Multiplication (`*=` with W0/W2) applies to A_ij ranges. - Addition (`+=` with W1) applies to R_ij ranges. 3. **Cumulative Effect**: R_ij[3–4] is modified in all three steps, suggesting it accumulates results across iterations. ### Interpretation The diagram represents a **weighted aggregation algorithm** where contributions from specific element pairs (a23, a43, a45) are computed by: 1. Selecting relevant indices in A_ij and R_ij. 2. Applying element-wise multiplication (W0/W2) to A_ij ranges. 3. Accumulating results via addition (W1) to R_ij ranges. 4. The repeated modification of R_ij[3–4] implies it acts as a **central accumulator** for contributions from overlapping element pairs. This process could model scenarios like: - **Graph-based computations** (e.g., edge weight aggregation in networks). - **Signal processing** (e.g., feature weighting in time-series data). - **Machine learning** (e.g., attention mechanisms in transformers). The use of W0/W2 and W1 suggests **conditional weighting** based on element properties, while the directional arrows emphasize **data flow dependencies** between steps. </details> operation is a cross copy-primary multiply instruction, multiplying two copies of the first weight coefficient by the two values in a 23 , then placing the results in the SIMD register containing r 34 (Figure 7a). Then, the value of a 23 in the SIMD register is modified to become a 43 , replacing one data element in the SIMD register either through mutate or copy. The second floating point operation is a cross complex multiply-add instruction, performing a cross operation to deal with the reversed values (Figure 7b). Finally, the value a 45 , which has either been preloaded by load-copy or is created by a second mutation in mutate-mutate, is used to perform the last computation (Figure 7c). We list the resource requirements of the load-copy and the mutate-mutate kernels in Table 1. The two kernels are at complementary ends of a spectrum. The mutate-mutate kernel increases pressure exclusively on the load pipeline while the load-copy kernel incurs extra cycles on the floating point unit. The two strategies can be used in concert, using mutate-mutate when the floating point unit is the bottleneck and the load-copy when it is not. Table 1: Resource usage per stencil of mutate-mutate and load-copy | Kernel | Operations | Operations | Cycles | Cycles | Registers | Registers | |---------------|--------------|--------------|----------|----------|-------------|-------------| | | ld-st | FPU | ld-st | FPU | Input | Output | | mutate-mutate | 2-1 | 3 | 6 | 3 | 1 | 1 | | load-copy | 1-1 | 4 | 4 | 4 | 2 | 1 | ## 4.3 Unroll-and-Jam The 3-point kernels are relatively easy to specify in assembly, but the floating point and load/store instruction latencies will cause pipeline stalls if they are not unrolled. Further unrolling in the k -direction beyond two is a possible solution that we do not explore in this paper. Although this solution would reduce the number of concurrent memory streams, it would also reduce data reuse for the other stencils studied in this paper. Unrolling and jamming once in transverse directions provides independent arithmetic operations to hide instruction latency, but interleaving the instructions by hand produces large kernels that are difficult to understand and modify. To simplify the design process, we constructed a synthetic code generator and simulator with reordering capability to interleave the jammed FPU and load/store instructions to minimize pipeline stalls. The synthetic code generator also gives us the flexibility to implement many general stencil operators, including the 7-point and 27-point stencils using the 3-point stencil as a building block. Unroll-and-jam serves a second purpose for the 7-point and 27-point stencil operators due to the overlapped data usage among adjacent stencils. Unroll-and-jam eliminates redundant loads of common data elements among the jammed stencils, reducing pressure on the memory subsystem by increasing the effective arithmetic intensity of the kernel. This can be quantified by comparing the number of input streams, which we refer to as the 'frame size,' with the number of output streams. For example, an unjammed 27-point stencil requires a frame size of 9 input streams for a single output stream. If we unroll once in i and j , we generate a 2 × 2 jam with 4 output streams using a frame size of 16, improving the effective arithmetic intensity by a factor of 9 4 . We used mutate-mutate and load-copy kernels to construct 3-, 7-, and 27-point stencil kernels over several different unrolling configurations. Table 2 lists the register allocation requirements for these configurations and provides per-cycle computational requirements. In both the mutate-mutate and load-copy kernels, the 27-point stencil is theoretically FPUbound because of the high reuse of loaded input data elements across streams. As can be seen in Table 1, the mutate-mutate kernel allows more unrolling for the 27-point stencil than the load-copy kernel because it uses fewer registers per stencil. The number of allocated registers for input data streams at the mutate-mutate kernel is equal to the number of the input data streams, while the load-copy kernel requires twice the number of registers for its input data streams. ## 4.4 PowerPC 450 Simulator The high-level code synthesis technique generates as many as hundreds of assembly instructions with aggressive unrolling. For optimal performance, we greedily schedule the assembly Table 2: Computational requirements | 3-lc-2x4 | 3-lc-2x3 | 3-lc-2x2 | 3-lc-2x1 | 3-lc-1x1 | 7-lc-2x3 | 7-mm-2x3 | 27-mm-2x3 | 27-mm-2x2 | 27-mm-1x3 | 27-mm-1x2 | 27-mm-1x1 | | Kernel | |------------|------------|------------|------------|------------|------------|------------|-------------|-------------|-------------|-------------|-------------|---------------------|---------------------| | 8 | 6 | 4 | 2 | 1 | 16 | 16 | 20 | 16 | 15 | 12 | 9 | Frame | Configurations | | 16 | 12 | 8 | 4 | 2 | 12 | 12 | 12 | 8 | 6 | 4 | 2 | Stencils/ Iteration | | | 16 | 12 | 8 | 4 | 2 | 22 | 16 | 20 | 16 | 15 | 12 | 9 | | Input | | 8 | 6 | 4 | 2 | 1 | 6 | 6 | 6 | 4 | 3 | 2 | 1 | Result | Registers | | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 4 | 4 | 4 | 4 | 4 | | Weight | | 8-8 | 6-6 | 4-4 | 2-2 | 1-1 | 16-6 | 22-6 | 40-6 | 32-4 | 30-3 | 24-2 | 18-1 | ld-st | Count | | 32 | 24 | 16 | 8 | 4 | 48 | 42 | 162 | 108 | 81 | 54 | 27 | FPU | | | 16-16 | 12-12 | 8-8 | 4-4 | 2-2 | 32-12 | 44-12 | 80-12 | 64-8 | 60-6 | 48-4 | 36-2 | ld-st | Instructions Cycles | | 32 | 24 | 16 | 8 | 4 | 48 | 42 | 162 | 108 | 81 | 54 | 27 | FPU | | | 100 | 100 | 100 | 100 | 100 | 91.7 | 100 | 56.8 | 66.7 | 81.5 | 96.3 | 100 | ld-st | Utilization | | 100 | 100 | 100 | 100 | 100 | 100 | 75 | 100 | 100 | 100 | 100 | 71.1 | FPU | % | | 16 | 16 | 16 | 16 | 16 | 29.3 | 29.3 | 34.7 | 40 | 48 | 56 | 80 | stencil | Bandwidth Bytes/ | instructions for in-order execution using the constraints outlined in 15. First, we produce a list of non-redundant instructions using a Python generator. The simulator reflects an understanding of the constraints by modeling the instruction set, including semantics such as read and write dependencies, instruction latency, and which execution unit it occupies. It functions as if it were a PowerPC 450 with an infinite-lookahead, greedy, out-of-order execution unit. On each cycle, it attempts to start an instruction on both the load/store and floating point execution units while observing instruction dependencies. If this is not possible, it provides diagnostics about the size of the stall and what dependencies prevented another instruction from being scheduled. The simulator both modifies internal registers that can be inspected for verification purposes and produces a log of the reordered instruction schedule. The log is then rendered as inline assembly code which can be compiled using the XL or GNU C compilers. ## 5 Performance Code synthesis allows us to easily generate and verify the performance of 3-, 7-, and 27-point stencil operators against our predictive models over a range of unrolling-and-jamming and inner kernel options. We use individual MPI processes mapped to the four PowerPC 450 cores to provide 4-way parallelization and ascertain performance characteristics out to the shared L3 and DDR memory banks. The generated assembly code, comprising the innermost loop, incurs a constant computational overhead from the prologue, epilogue, and registers saving/restoring. The relative significance of this overhead is reduced when larger computations are performed at the innermost loop. This motivated us to decompose the problem's domain among the four cores along the outermost dimension for the 3- and 7-point stencils, resulting in better performance. However, the 27-point stencil has another important property, dominating the innermost loop overhead cost. Its computations inhibit relatively high input data points sharing among neighbor stencils. Splitting the innermost dimension allows the shared input data points to be reused by consecutive middle dimension iterations, where the processor will likely have them in the L1 cache, resulting in faster computations. Conversely, if the computation is performed using a large innermost dimension, the shared input data points will no longer be in the L1 cache at the beginning of the next iteration of the middle loop. All three studies were conducted over a range of cubic problems from size 14 3 to 362 3 . The problem sizes were chosen such that all variations of the kernels can be evaluated naively without extra code for cleanup. We computed each sample in a nested loop to reduce noise from startup costs and timer resolution, then select the highest performing average from the inner loop samples. There was almost no noticeable jitter in sample times after the first several measurements. All performance results are given per-core, though results were computed on an entire node with a shared L3 cache. Thus, full node performance can be obtained by multiplying by four. During the course of our experiments on the stencils, we noticed performance problems for many of the stencil variants when loading data from the L3 cache. The large number of concurrent hardware streams in the unroll-and-jam approach overwhelms the L2 streaming unit, degrading performance. This effect can be amplified in the default optimistic prefetch mode for the L2, causing wasted memory traffic from the L3. We made use of a boot option that disables optimistic prefetch from the L2 and compare against the default mode where applicable in our results. We distinguish the two modes by using solid lines to indicate performance results obtained in the default mode and dashed lines to indicate results where the optimistic prefetch in L2 has been disabled. ## 5.1 3-Point Stencil Computations We begin our experiments with the 3-point stencil (Figure 8), the computational building block for the other experiments in our work. For a more accurate streaming bandwidth peak estimation, we considers 3.7 bytes/s read, from the DRAM, and 5.3 bytes/s write bandwidth at the stencil computations. Also, we compute the L3 bandwidth peak using 4.7 bytes/s read and 5.3 bytes/s write bandwidth at the stencil computations. The 3-point stencil has the lowest arithmetic intensity of the three stencils studied, and unlike its 7-point and 27-point cousins, does not see an increase in effective arithmetic intensity when unrolland-jam is employed. It is clear from Section 4.2 that the load-copy kernel is more efficient in bandwidth-bound situations, so we use it as the basis for our unroll-and-jam experiments. We see the strongest performance in the three problems that fit partially in the L1 cache (the peak of 224 Mstencil/s is observed at 26 3 ), with a drastic drop off as the problem inputs transition to the L3. The most robust kernel is the 2 × 1 jam, which reads and writes to two streams simultaneously, and can therefore engage the L2 prefetch unit most effectively. The larger unrolls (2 × 2, 2 × 3, and 2 × 4), enjoy greater performance in and near the L1, but then suffer drastic performance penalties as they exit the L1 and yet another performance dip near 250 3 . Disabling optimistic L2 prefetch does not seem to have any large effect on the 2 × 1 kernel, though it unreliably helps or hinders the other kernels. The 3-point kernel seems to be an ideal target on the PowerPC 450 for standard unrolling in the fastest moving dimension, k , a technique we did not attempt due to its limited application to the larger problems we studied. Unroll-and-jam at sufficient sizes to properly cover pipeline hazards overwhelms the L2 streaming unit due to the large number of simultaneous streams to memory. Unrolling in k would cover these pipeline hazards without increasing the number of streams. ## 5.2 7-Point Stencil Computations Our next experiment focuses on the performance of the 7-point stencil operator (Figure 9). We compare the mutate-mutate and load-copy kernels using the same unroll configurations. We note that the mutate-mutate kernel can support a slightly more aggressive unroll-andjam on this problem with a compressed usage of general purpose registers that was only implemented for the 27-point stencil. Once again we notice strong performance within the L1, then a dropoff as the loads start coming from the L3 instead of the L1. The performance drop near 256 3 is caused when the 2 × 3 kernel's frame size of (2 + 2)(3 + 2) -4 = 16 multiplied by the length of the domain exceeds the size of L1. For smaller sizes, neighbors in the j direction can reside in L1 between consecutive passes so that only part of the input frame needs to be supplied by streams from memory. With up to 16 input streams and 2 · 3 = 6 output streams, there is no hope of Figure 8: 3-point stencil performance with load-copy kernel at various unroll-and-jams (ixj). L3 bandwidth peak is 265 Mstencil/s <details> <summary>Image 6 Details</summary>  ### Visual Description ## Line Graph: Streaming Bandwidth Peak Performance ### Overview The graph illustrates the performance of different streaming bandwidth (BW) configurations over varying input lengths, measured in Mstencil/s. Key elements include a legend in the top-right corner, a dashed "BW peak" reference line at 240 Mstencil/s, and six data series representing different BW multipliers (2x4, 2x3, 2x2, 2x1, and 1x1). ### Components/Axes - **X-axis (Input length)**: Ranges from 0 to 400 in increments of 50. - **Y-axis (Mstencil/s)**: Ranges from 60 to 240 in increments of 20. - **Legend**: Positioned in the top-right corner, with color-coded labels: - `2x4`: Cyan (solid line) - `2x3`: Magenta (solid line) - `2x2`: Green (solid line) - `2x1`: Red (solid line) - `1x1`: Black (solid line) - `BW peak`: Dashed black line at 240 Mstencil/s. ### Detailed Analysis 1. **2x4 (Cyan)**: - Starts at ~185 Mstencil/s at input length 0. - Peaks sharply at ~220 Mstencil/s near input length 20. - Drops rapidly to ~100 Mstencil/s by input length 50. - Stabilizes with minor fluctuations (~90–110 Mstencil/s) afterward. 2. **2x3 (Magenta)**: - Begins at ~190 Mstencil/s. - Peaks at ~210 Mstencil/s near input length 20. - Declines to ~110 Mstencil/s by input length 50. - Stabilizes with oscillations (~95–115 Mstencil/s). 3. **2x2 (Green)**: - Initial value ~175 Mstencil/s. - Peaks at ~200 Mstencil/s near input length 20. - Drops to ~120 Mstencil/s by input length 50. - Stabilizes with fluctuations (~100–125 Mstencil/s). 4. **2x1 (Red)**: - Starts at ~160 Mstencil/s. - Peaks at ~180 Mstencil/s near input length 20. - Declines to ~120 Mstencil/s by input length 50. - Stabilizes with minor variations (~110–125 Mstencil/s). 5. **1x1 (Black)**: - Flat line at ~60 Mstencil/s across all input lengths. 6. **BW Peak (Dashed Black)**: - Horizontal reference line at 240 Mstencil/s. - No data series reaches this value. ### Key Observations - **Initial Peak**: All BW multipliers exhibit a sharp performance spike near input length 20, followed by a rapid decline. - **Stabilization**: After input length 50, performance stabilizes with minor oscillations for all multipliers except 1x1. - **1x1 Baseline**: The 1x1 line remains constant, suggesting it represents a control or baseline condition. - **BW Peak Unattained**: No configuration reaches the theoretical BW peak of 240 Mstencil/s. ### Interpretation The graph demonstrates that higher BW multipliers (e.g., 2x4) achieve greater initial performance gains but experience rapid degradation as input length increases. The stabilization phase suggests diminishing returns or system saturation beyond input length 50. The unattained BW peak implies hardware or algorithmic limitations preventing optimal performance. The 1x1 line’s flatness indicates a consistent baseline, possibly representing a non-scaled or reference configuration. The sharp initial drop may reflect resource contention or inefficiencies at higher BW settings. </details> effectively using the L2 prefetch unit. The load-copy kernel shows better performance than the mutate-mutate kernel, as it is clear here that load/store cycles are more constrained than floating point cycles. We also notice that performance of the load-copy kernel improves with the L2 optimistic prefetch disabled slightly within the L3, and drastically when going to the DDR banks. This is likely due to the kernel's improved performance, and therefore increased sensitivity with regards to latency from the memory subsystem. It is likely that the 7-point stencil could attain better results by incorporating cache tiling strategies, though we note that without any attempts at cache tiling the performance of this result is commensurate with previously reported results for the PowerPC 450 that focused on cache tiling for performance tuning. ## 5.3 27-Point Stencil Computations The 27-point stencil should be amenable to using a large number of jammed unrolls due to the high level of reuse between neighboring stencils. Indeed, we see nearly perfect scaling in Figure 10 as we increase the number of jams from 1 to 6 using the mutate-mutate kernel. Although there is a gradual drop off from the peak of 54 Mstencil/s (85% of arithmetic peak) Figure 9: 7-point stencil performance comparison between mutate-mutate (mm) and loadcopy (lc) kernels at a fixed unroll-and-jam (i=2 and j=3). Experiments with L2 optimistic prefetch ends with 'o' <details> <summary>Image 7 Details</summary>  ### Visual Description ## Line Chart: Mstencil/s Performance Across Input Lengths ### Overview The chart displays the performance of multiple configurations (measured in Mstencil/s) as input length increases from 0 to 400. Key trends include sharp initial drops followed by stabilization or fluctuation, with distinct patterns for each configuration. ### Components/Axes - **X-axis (Input length)**: Ranges from 0 to 400 in increments of 50. - **Y-axis (Mstencil/s)**: Ranges from 40 to 160 in increments of 20. - **Legend**: Located in the top-right corner, with the following entries: - Solid black: `mm_2x3o` - Dotted black: `mm_2x3` - Solid green: `lc_2x3o` - Dotted green: `lc_2x3` - Dashed black: `Streaming` - Dash-dot green: `BW peak` - Dotted black: `L3 BW peak` ### Detailed Analysis 1. **`mm_2x3o` (Solid black)**: - Starts at ~140 Mstencil/s at input length 0. - Drops sharply to ~60 by input length 50. - Fluctuates between 60–70 for input lengths 100–350. 2. **`lc_2x3o` (Solid green)**: - Starts at ~150 Mstencil/s at input length 0. - Drops to ~80 by input length 50. - Fluctuates between 60–80 for input lengths 100–350. 3. **`lc_2x3` (Dotted green)**: - Starts at ~80 Mstencil/s at input length 0. - Drops to ~60 by input length 50. - Fluctuates between 60–80 for input lengths 100–350. 4. **`Streaming` (Dashed black)**: - Starts at ~120 Mstencil/s at input length 0. - Drops to ~50 by input length 50. - Fluctuates between 50–60 for input lengths 100–350. 5. **`BW peak` (Dash-dot green)**: - Starts at ~140 Mstencil/s at input length 0. - Drops to ~60 by input length 50. - Fluctuates between 60–70 for input lengths 100–350. 6. **`L3 BW peak` (Dotted black)**: - Starts at ~140 Mstencil/s at input length 0. - Drops to ~60 by input length 50. - Fluctuates between 60–70 for input lengths 100–350. ### Key Observations - **Initial Drop**: All configurations experience a sharp performance drop (~50% reduction) between input lengths 0–50. - **Stabilization**: After input length 50, performance stabilizes or fluctuates within narrow ranges. - **Highest Initial Performance**: `lc_2x3o` (solid green) achieves the highest starting value (~150 Mstencil/s). - **Lowest Sustained Performance**: `Streaming` (dashed black) maintains the lowest values (~50–60 Mstencil/s) after input length 50. - **Similar Patterns**: `BW peak` and `L3 BW peak` lines exhibit nearly identical trends, suggesting shared characteristics. ### Interpretation The chart likely represents computational or data processing performance metrics across different hardware/software configurations. The initial drop may reflect setup overhead or resource allocation delays. Configurations with higher initial performance (`lc_2x3o`, `mm_2x3o`) degrade less sharply than others, indicating better scalability or optimization. The `Streaming` configuration’s consistently low values suggest it prioritizes throughput over peak performance. The `BW peak` and `L3 BW peak` lines imply bandwidth thresholds or limits influencing performance. The stabilization phase after input length 50 highlights a transition from initialization to steady-state operation. </details> as the problem sizes increase to the point that there is little reuse from the L1 cache, the kernel consistently sustains an average of 45 Mstencil/s (70% of arithmetic peak), even when the problem sizes greatly exceed the L3 cache. Despite the L2-overwhelming frame size ((2 + 2)(3 + 2) = 20 input streams and 2 · 3 = 6 output streams), the jammed stencil achieves good performance with no blocking largely due to the high level of reuse of input data afforded by the unrolls in i and j . ## 5.4 Model Validation As we utilized a simulator which incorporates a model of the architecture's performance characteristics to produce our kernels, we sought to validate our performance model by comparing the implicit predictions of our generative system to the empirical results seen in Table 3. Since performance within a core is often considerably easier to predict than when one must go beyond the core for memory access, we divide our comparisons into those on-core (L1) and those that go off the core, to L3 or main memory (streaming). Figure 10: 27-point stencil performance with mutate-mutate kernel at various unroll-andjams (ixj). Best L3 bandwidth peak is 118 Mstencil/s and best streaming bandwidth peak is 97.5 Mstencil/s, both corresponds to 2x3 unroll-and-jam <details> <summary>Image 8 Details</summary>  ### Visual Description ## Line Graph: Mstencil/s vs Input Length ### Overview The image is a line graph comparing the performance of different computational configurations (e.g., 2x3, 2x2, 1x3, etc.) in terms of "Mstencil/s" (measured on the y-axis) against "Input length" (x-axis). A dashed horizontal line labeled "Compute peak" at ~60 Mstencil/s serves as a reference. The graph includes five data series with distinct colors and labels, showing varying performance trends across input lengths. --- ### Components/Axes - **X-axis (Input length)**: Ranges from 0 to 400, with increments of 50. Labeled "Input length." - **Y-axis (Mstencil/s)**: Ranges from 0 to 70, with increments of 10. Labeled "Mstencil/s." - **Legend**: Located in the bottom-right corner, with the following mappings: - **Black dashed line**: "Compute peak" (~60 Mstencil/s). - **Black solid line**: "2x3" (highest performance). - **Green solid line**: "2x2" (moderate performance). - **Purple solid line**: "1x3" (lower performance). - **Cyan solid line**: "1x2" (lower performance). - **Red solid line**: "1x1" (lowest performance). --- ### Detailed Analysis 1. **Black solid line (2x3)**: - Starts at ~42 Mstencil/s at input length 0. - Peaks at ~52 Mstencil/s around input length 50. - Fluctuates between ~40–50 Mstencil/s for input lengths 100–400. - **Trend**: Highest performance, with a sharp initial rise and gradual stabilization. 2. **Green solid line (2x2)**: - Starts at ~34 Mstencil/s at input length 0. - Peaks at ~42 Mstencil/s around input length 50. - Stabilizes between ~38–40 Mstencil/s for input lengths 100–400. - **Trend**: Moderate performance, with a peak followed by steady values. 3. **Purple solid line (1x3)**: - Starts at ~28 Mstencil/s at input length 0. - Peaks at ~34 Mstencil/s around input length 50. - Fluctuates between ~28–30 Mstencil/s for input lengths 100–400. - **Trend**: Lower performance, with a peak and gradual decline. 4. **Cyan solid line (1x2)**: - Starts at ~19 Mstencil/s at input length 0. - Peaks at ~23 Mstencil/s around input length 50. - Stabilizes between ~22–23 Mstencil/s for input lengths 100–400. - **Trend**: Lower performance, with a peak and steady values. 5. **Red solid line (1x1)**: - Remains flat at ~10 Mstencil/s across all input lengths. - **Trend**: Constant, lowest performance. 6. **Dashed line (Compute peak)**: - Horizontal line at ~60 Mstencil/s, above all data series. - **Trend**: Theoretical maximum, not achieved by any configuration. --- ### Key Observations - The **2x3 configuration** consistently outperforms others, with the highest Mstencil/s values. - The **1x1 configuration** is the baseline, showing no improvement with input length. - All configurations except the 1x1 show a peak around input length 50, followed by stabilization. - The **Compute peak** (~60 Mstencil/s) is unattainable by any configuration, suggesting a theoretical limit. --- ### Interpretation The data demonstrates that larger computational configurations (e.g., 2x3) achieve higher performance, but with variability. The 1x1 configuration serves as a baseline, while the Compute peak represents an idealized target. The fluctuations in performance (e.g., 2x3's peaks and troughs) may indicate sensitivity to input length or resource allocation. The absence of any configuration reaching the Compute peak suggests potential bottlenecks or inefficiencies in the system. This graph highlights the trade-off between configuration complexity and performance, with diminishing returns as input length increases. </details> Our modeling of the 27-point stencils can be seen to be highly accurate in Table 3. Inside the L1 cache the disparity between predicted and actual performance is consistently less than 1%. Shifting our attention to the streaming predictions, our accuracy can be seen to be considerably degraded. This is not surprising, given that our simulator was largely targeted to model the L1 domain. However, the relative error is less than 15% in all cases; this appears to be sufficient for producing highly efficient code. This shortcoming appears to stem directly from the level of detail with which we model the shared L3 cache and main memory subsystem and we are working to address this in our simulator. The match between predicted and witnessed performance for the 7-point stencil shows the same pattern. When modeling performance inside the L1 our relative error is less than 10%, but when extending our prediction to the components of the system shared by all four cores, our error is as great as 17.5%. Our greatest error in this instance is an underprediction that is probably attributable to a fortuitous alignment of the continuous vectors in the k-direction, as staggering these carefully often result in bandwidth benefits on the order of 10-15%. Finally, we assess our model for the 3-point stencil. Again we see good agreement between the model and observed performance within the L1, though prediction accuracy degrades for problem sizes requiring streaming. From some further experimentation, we are reasonably certain that the chief reason for our lack of accuracy in predicting performance outside of the L1 stems from the bandwidth that must be shared between the multiple write streams and our failure to account for this in our model. It is most apparent in the 3-point stencil predictions as the ratio of write streams to either read streams or floating point operations is highest in this case. ## 6 Concluding Remarks ## 6.1 Conclusion The main contribution of this work is effective register and instruction scheduling for constant coefficient linear operators on power-efficient processors. The loads of the input vector elements and stores of the output vector elements are minimized and the fraction of multiplyadds among all cycles is maximized. This is achieved by using two novel 3-point streamcomputation sub-kernels designed for the PowerPC 450's instruction set, mutate-mutate and load-copy. Both kernels were possible without data layout reordering because of the extensively multiplexed SIMD-like floating point units implemented in the PowerPC 450 core. Recommendations for the research agenda for computational software libraries in the exascale domain include the the fusion of library routine implementations as well as the creation of frameworks that enable the optimal instantiation of a given routine when supplied with architectural information [13]. We feel that the work presented here is a contribution to that end. Further, the nature of our simulator allows us to optimize our code as measured by other metrics such as bandwidth or energy consumption, given a simple model of the cost. We also demonstrate performance comparable to advanced cache tiling approaches for the 7-point stencil, despite the fact that we make no effort to optimize for cache reuse. ## 6.2 Future Work While the three problems considered (3-point stencils in one dimension and 7-point and 27-point stencils in three dimensions, with constant coefficients and symmetry within each spatial dimension, but not across them) are heavily used in applications, there are numerous generalizations. The suitability of our approach can be characterized by the arithmetic intensity associated with each generalization. We elaborate on two that tend to increase the arithmetic intensity, higher-order stencils and chained iterative passes over the vectors, and two that tend to decrease arithmetic intensity, irregular stencils and spatial varying coefficients. Higher-order stencils expand the number of adjacent input vector elements that enter into a single output vector element, in successive steps of semi-width one in each of the spatial dimensions i , j , and k , with an additional weight coefficient corresponding to each additional increment of semi-width in each dimension. This is a modest generalization. Higher-order Table 3: Predictions vs. observations for in-L1 and streaming performance, in Megastencils per second | 3-lc-2x4 | 3-lc-2x3 | 3-lc-2x2 | 3-lc-2x1 | 3-lc-1x1 | 7-lc-2x3 | 7-mm-2x3 | 27-mm-2x3 | 27-mm-2x2 | 27-mm-1x3 | 27-mm-1x2 | 27-mm-1x1 | Kernel | |------------|------------|------------|------------|------------|------------|------------|-------------|-------------|-------------|-------------|-------------|---------------------| | 425 | 425 | 425 | 425 | 425 | 212.5 | 182.14 | 62.96 | 62.96 | 62.96 | 62.96 | 44.74 | Instruction Naive | | 197.1 | 202.31 | 193.36 | 147.29 | 88.12 | 143.83 | 126.84 | 54.62 | 44.59 | 34.3 | 23.35 | 11.93 | limits Simulated | | 338.72 | 338.72 | 338.72 | 338.72 | 338.72 | 203.54 | 203.54 | 175.52 | 154.28 | 130.58 | 113.19 | 80.88 | Bandwidth L1 | | 231.51 | 231.51 | 231.51 | 231.51 | 231.51 | 116.84 | 116.84 | 97.51 | 83.68 | 68.99 | 58.69 | 40.54 | limits streaming | | 197.1 | 202.31 | 193.36 | 147.29 | 88.12 | 143.83 | 126.84 | 54.62 | 44.59 | 34.3 | 23.35 | 11.93 | In-L1 Predicted | | 199.05 | 195.83 | 184.84 | 142.04 | 81.33 | 132.1 | 124.43 | 54.17 | 44.53 | 34.23 | 23.39 | 11.92 | Observed | | 197.1 | 202.31 | 193.36 | 147.29 | 88.12 | 116.84 | 116.84 | 54.62 | 44.59 | 34.3 | 23.35 | 11.93 | Streaming Predicted | | 83.9 | 86.62 | 96.23 | 119.99 | 67.44 | 74.21 | 59.69 | 42.64 | 38.37 | 28.26 | 22.56 | 12.37 | Observed | discretization increases register pressure because of the larger number of inputs that combine in each output. Opportunities for reuse of input elements expand with the stencil width up to the ability to keep them resident. In a P -point regular stencil (regardless of number of spatial dimensions) each input element is operated upon with a pre-stored weight P times: once in the stencil centered upon it, and once in each neighboring stencil with which its own stencil overlaps. Floating point arithmetic intensity increases in proportion to P . That is, if there are N elements in the input or output array, there are PN floating point multiply-adds per N floating reads and N floating writes. Explicit methods for nonlinear systems, especially with high-order discretization techniques such as Weighted Essentially Non-Oscillatory (WENO) or discontinuous Galerkin [36], have similar properties, including a larger number of input streams, but with much higher arithmetic intensity. S -stage chaining (as in the simultaneous accumulation of Ax , A 2 x , A 3 x , . . . A s x ) allows the output vector to be fed back as input before being written. Per output vector of N floating point writes, there are N/S reads and PN floating point multiply-adds. Therefore, up to the ability to keep the additional operands cached, both higher-order operators and chained operations improve the potential for the transformations described here. Irregular stencils require integer reads, in addition to floating point reads to determine which elements of the input vector go with each row of the matrix. This further dilutes advantages that lead to the great breakthroughs in stencils per second described here. Stencil operations with constant coefficients and sparse matrix-vector multiplies with general coefficients are similar when counting floating operations, but very different when it comes to data volume. Spatially varying coefficients require the loading of additional weights, each of which is used only once, each of which is of the same floating precision of the input and output vectors, P of them in the production of each output vector element, with each input vector element being combined with P different weights. While each input and output element can still be reused up to P times in the execution of one pass through the overall vector-to-vector map, the dominant array in the workspace is the coefficient matrix of weights so the benefits of reusing the vectors are minimal. This situation is typical for nonlinear problems when using Newton-Krylov and linear multigrid methods. However, when 'free flops' are available, the weights can also be recomputed on the fly as a nonlinear function of a given state and/or scalar coefficients. In this case, the number of input streams is similar to the linear constant coefficient case (perhaps larger by a factor of 2 or 3), but the number of floating point results is several times higher and involves the problem-specific 'physics.' Putting the physics inside the kernel like this suggests that there will be an emphasis on the ability to quickly develop high-performance kernels for new physics. ## Acknowledgments We are grateful to Andy Ray Terrel for his helpful commentary on an early draft of this paper. We are also indebted to Andrew Winfer for his support in conducting our numerical experiments on the Shaheen Blue Gene/P system at the KAUST Supercomputing Laboratory, and to Vernon Austel for his assistance in running experiments on the Blue Gene/P system at IBM's Watson Research Center. ## References - [1] Rajagopal Ananthanarayanan, Steven K Esser, Horst D Simon, and Dharmendra S Modha. The cat is out of the bag: cortical simulations with 109 neurons, 1013 synapses. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis , SC '09, pages 63:1--63:12, New York, NY, USA, 2009. ACM. - [2] Mauricio Araya-Polo, F´ elix Rubio, Ra´ ul De, Mauricio Hanzich, and Jos´ e Mar´ ıa. 3D seismic imaging through reverse-time migration on homogeneous and heterogeneous multi-core processors. Scientific Programming , 17:185-198, 2009. - [3] D.H. Bailey, E. Barszcz, J.T. Barton, D.S. Browning, R.L. Carter, L. Dagum, R.A. Fatoohi, P.O. Frederickson, T.A. Lasinski, R.S. Schreiber, H.D. Simon, V. Venkatakrishnan, and S.K. Weeratunga. The NAS parallel benchmarks. International Journal of High Performance Computing Applications , 5(3):63, 1991. - [4] MJ Berger and J Oliger. Adaptive Mesh Refinement for Hyperbolic Partial Differential Equations. Journal of Computational Physics , 53(3):484-512, 1984. - [5] D. Callahan, J. Cocke, and K. Kennedy. Estimating interlock and improving balance for pipelined architectures 1. Journal of Parallel and Distributed Computing , 5(4):334-358, 1988. - [6] S. Carr and K. Kennedy. Improving the ratio of memory operations to floatingpoint operations in loops. ACM Transactions on Programming Languages and Systems (TOPLAS) , 16(6):1768-1810, 1994. - [7] C.M. Chang, C.M. Chen, and C.T. King. Using integer linear programming for instruction scheduling and register allocation in multi-issue processors 1. Computers & Mathematics with Applications , 34(9):1-14, 1997. - [8] Matthias Christen, Olaf Schenk, and Helmar Burkhart. Automatic code generation and tuning for stencil kernels on modern shared memory architectures. Computer Science Research and Development , 26(3-4):205-210, April 2011. - [9] Matthias Christen, Olaf Schenk, Esra Neufeld, Peter Messmer, and Helmar Burkhart. Parallel data-locality aware stencil computations on modern micro-architectures. In 2009 IEEE International Symposium on Parallel & Distributed Processing , pages 1-10. IEEE, May 2009. - [10] Dac C. Pham, Tony Aipperspach, David Boerstler, Mark Bolliger, Rajat Chaudhry, Dennis Cox, Paul Harvey, Paul M. Harvey, H. Peter Hofstee, Charles Johns, Jim Kahle, Atsushi Kameyama, John Keaty, Yoshio Masubuchi, Mydung Pham, J¨ urgen Pille, Stephen Posluszn, Kazuaki Yazawa. Overview of the architecture, circuit design, and physical implementation of a first-generation cell processor. Solid-State Circuits, IEEE Journal of , 41(1):179-196, 2006. - [11] K. Datta, S. Williams, V. Volkov, J. Carter, L. Oliker, J. Shalf, and K. Yelick. Autotuning the 27-point Stencil for Multicore. In Proc. iWAPT2009: The Fourth International Workshop on Automatic Performance Tuning , 2009. - [12] Kaushik Datta. Auto-tuning Stencil Codes for Cache-Based Multicore Platforms . PhD thesis, EECS Department, University of California, Berkeley, December 2009. - [13] Jack Dongarra, Pete Beckman, Terry Moore, Patrick Aerts, Giovanni Aloisio, JeanClaude Andre, David Barkai, Jean-Yves Berthou, Taisuke Boku, Bertrand Braunschweig, Franck Cappello, Barbara Chapman, Chi Xuebin, Alok Choudhary, Sudip Dosanjh, Thom Dunning, Sandro Fiore, Al Geist, William Gropp, Robert Harrison, Mark Hereld, Michael Heroux, Adolfy Hoisie, Koh Hotta, Jin Zhong, Yutaka Ishikawa, Fred Johnson, Sanjay Kale, Richard Kenway, David Keyes, Bill Kramer, Jesus Labarta, Alain Lichnewsky, Thomas Lippert, Bob Lucas, Barney Maccabe, Satoshi Matsuoka, Paul Messina, Peter Michielse, Bernd Mohr, Matthias S. Mueller, Wolfgang E. Nagel, Hiroshi Nakashima, Michael E Papka, Dan Reed, Mitsuhisa Sato, Ed Seidel, John Shalf, David Skinner, Marc Snir, Thomas Sterling, Rick Stevens, Fred Streitz, Bob Sugar, Shinji Sumimoto, William Tang, John Taylor, Rajeev Thakur, Anne Trefethen, Mateo Valero, Aad van der Steen, Jeffrey Vetter, Peg Williams, Robert Wisniewski, and Kathy Yelick. The international exascale software project roadmap. International Journal of High Performance Computing Applications , 25(1):3-60, 2011. - [14] Hikmet Dursun, Ken-Ichi Nomura, Liu Peng, Richard Seymour, Weiqiang Wang, Rajiv K Kalia, Aiichiro Nakano, and Priya Vashishta. A Multilevel Parallelization Framework for High-Order Stencil Computations. In Proceedings of the 15th International Euro-Par Conference on Parallel Processing , Euro-Par '09, pages 642-653, Berlin, Heidelberg, 2009. Springer-Verlag. - [15] A.E. Eichenberger, P. Wu, and K. O'Brien. Vectorization for SIMD architectures with alignment constraints. ACM SIGPLAN Notices , 39(6):82-93, 2004. - [16] W. Feng and K. Cameron. The Green500 List: Encouraging Sustainable Supercomputing. COMPUTER , pages 50-55, 2007. - [17] A. Ganapathi, K. Datta, A. Fox, and D. Patterson. A case for machine learning to optimize multicore performance. In Proceedings of the First USENIX conference on Hot topics in parallelism , pages 1-1. USENIX Association, 2009. - [18] Amol Ghoting and Konstantin Makarychev. Indexing genomic sequences on the IBM Blue Gene. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis , SC '09, pages 61:1--61:11, New York, NY, USA, 2009. ACM. - [19] John L Hennessy and Thomas Gross. Postpass Code Optimization of Pipeline Constraints. ACM Trans Program Lang Syst , 5(3):422-448, 1983. - [20] T. Henretty, K. Stock, L.N. Pouchet, F. Franchetti, J. Ramanujam, and P. Sadayappan. Data Layout Transformation for Stencil Computations on Short-Vector SIMD Architectures. In Compiler Construction , pages 225-245. Springer, 2011. - [21] IBM Blue Gene Team. Overview of the IBM Blue Gene/P project. IBM J. Res. Dev. , 52(1/2):199--220, 2008. - [22] S. Kamil, C. Chan, L. Oliker, J. Shalf, and S. Williams. An auto-tuning framework for parallel multicore stencil computations. In Parallel & Distributed Processing (IPDPS), 2010 IEEE International Symposium on , pages 1-12. IEEE, 2010. - [23] Shoaib Kamil, Kaushik Datta, Samuel Williams, Leonid Oliker, John Shalf, and Katherine Yelick. Implicit and explicit optimizations for stencil computations. Proceedings of the 2006 workshop on Memory system performance and correctness - MSPC '06 , page 51, 2006. - [24] Shoaib Kamil, Parry Husbands, Leonid Oliker, John Shalf, and Katherine Yelick. Impact of modern memory subsystems on cache optimizations for stencil computations. Memory System Performance , pages 36-43, 2005. - [25] D.E. Keyes. Exaflop/s: The why and the how. Comptes Rendus M´ ecanique , 339(23):70-77, 2011. - [26] Sriram Krishnamoorthy, Muthu Baskaran, Uday Bondhugula, J Ramanujam, Atanas Rountev, and P Sadayappan. Effective automatic parallelization of stencil computations. ACM Sigplan Notices , 42(6):235, 2007. - [27] Pat Hanrahan Larry Seiler, Doug Carmean, Eric Sprangle, Tom Forsyth, Michael Abrash, Pradeep Dubey, Stephen Junkins, Adam Lake, Jeremy Sugerman, Robert Cavin, Roger Espasa, Ed Grochowski, Toni Juan. Larrabee: a many-core x86 architecture for visual computing. In ACM SIGGRAPH 2008 papers , pages 1-15. ACM, 2008. - [28] Z. Li and Y. Song. Automatic tiling of iterative stencil loops. ACM Transactions on Programming Languages and Systems (TOPLAS) , 26(6):975-1028, 2004. - [29] E. Lindholm, J. Nickolls, S. Oberman, and J. Montrym. NVIDIA Tesla: A unified graphics and computing architecture. Micro, IEEE , 28(2):39-55, 2008. - [30] C Mueller and B Martin. CorePy: High-Productivity Cell/BE Programming. Applications for the Cell/BE , 2007. - [31] C.J. Newburn, B. So, Z. Liu, M. McCool, A. Ghuloum, S. Du Toit, Z.G. Wang, Z.H. Du, Y. Chen, G. Wu, Peng Guo, Zhanglin Liu, and Dan Zhang. Intel's Array Building Blocks: A Retargetable, Dynamic Compiler and Embedded Language. Proceedings of Code Generation and Optimization , 2011. - [32] A. Nguyen, N. Satish, J. Chhugani, C. Kim, and P. Dubey. 3.5-D Blocking Optimization for Stencil Computations on Modern CPUs and GPUs. sc , pages 1-13, 2010. - [33] Liu Peng, Richard Seymour, Ken-Ichi Nomura, R K Kalia, Aiichiro Nakano, P Vashishta, Alexander Loddoch, Michael Netzband, W R Volz, and C C Wong. Highorder stencil computations on multicore clusters. IPDPS , pages 1-11, 2009. - [34] D F Richards, J N Glosli, B Chan, M R Dorr, E W Draeger, J.-L. Fattebert, W D Krauss, T Spelce, F H Streitz, M P Surh, and J A Gunnels. Beyond homogeneous decomposition: scaling long-range forces on Massively Parallel Systems. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis , SC '09, pages 60:1--60:12, New York, NY, USA, 2009. ACM. - [35] Gabriel Rivera and C W Tseng. Tiling optimizations for 3D scientific computations. In Proceedings of SC00 . IEEE Computer Society, 2000. - [36] C.W. Shu. High-order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD. International Journal of Computational Fluid Dynamics , 17(2):107-118, 2003. - [37] A. Solar-Lezama, G. Arnold, L. Tancau, R. Bodik, V. Saraswat, and S. Seshia. Sketching stencils. In Proceedings of the 2007 ACM SIGPLAN conference on Programming language design and implementation , pages 167-178. ACM, 2007. - [38] Sosa, C. and International Business Machines Corporation Organization. International Technical Support. IBM system Blue Gene solution: Blue Gene/P application development . IBM International Technical Support Organization, 2008. - [39] Yuan Tang, Rezaul Alam Chowdhury, Bradley C. Kuszmaul, Chi-Keung Luk, and Charles E. Leiserson. The pochoir stencil compiler. In Proceedings of the 23rd ACM symposium on Parallelism in algorithms and architectures - SPAA '11 , page 117, New York, New York, USA, June 2011. ACM Press. - [40] Gerhard Wellein, Georg Hager, Thomas Zeiser, Markus Wittmann, and Holger Fehske. Efficient Temporal Blocking for Stencil Computations by Multicore-Aware Wavefront Parallelization. 2009 33rd Annual IEEE International Computer Software and Applications Conference , pages 579-586, 2009. - [41] K Wilken. Optimal instruction scheduling using integer programming. ACM SIGPLAN Notices , 2000. - [42] Samuel Williams, Jonathan Carter, Leonid Oliker, John Shalf, and Katherine Yelick. Lattice Boltzmann simulation optimization on leading multicore platforms. 2008 IEEE International Symposium on Parallel and Distributed Processing , 69(9):1-14, 2008.