Paraiso : An Automated Tuning Framework for Explicit Solvers of Partial Differential Equations

We propose Paraiso, a domain specific language embedded in functional programming language Haskell, for automated tuning of explicit solvers of partial differential equations (PDEs) on GPUs as well as multicore CPUs. In Paraiso, one can describe PDE solving algorithms succinctly using tensor equations notation. Hydrodynamic properties, interpolation methods and other building blocks are described in abstract, modular, re-usable and combinable forms, which lets us generate versatile solvers from little set of Paraiso source codes. We demonstrate Paraiso by implementing a compressive hydrodynamics solver. A single source code less than 500 lines can be used to generate solvers of arbitrary dimensions, for both multicore CPUs and GPUs. We demonstrate both manual annotation based tuning and evolutionary computing based automated tuning of the program.


Introduction
Today, computer architectures are becoming more and more complex, making it more and more difficult to predict the performance of a program before running it. Parallel architectures, GPUs being one example and distributed memory machines another, forces us to program in different ways and our programs tend to become longer. Moreover, to optimize them we are asked to write many different implementations of such programs, and the coding task becomes time-consuming. However, if we describe the problem we want to solve in a domain-specific language (DSL) from which a lot of possible implementations are generated and benchmarked, we can automate the tuning processes.
Examples of such automated-tuning DSL approaches are found in fast Fourier transformation library FFTW [1], linear algebra library ATLAS [2], digital signal processing library SPIRAL [3]. In these works the authors regarded automated tuning not just as a tool to avoid manual tuning, but as a necessary tool to have "portable performance" -the practical way of optimizing domain-specific codes for various complicated architecture we have today and in the future. Now, our domain is explicit solvers of partial differential equations (PDEs). The main portions of such solvers consist of stencil computations, and a few global reduction operations needed to calculate Courant-Friedrichs-Lewy (CFL) conditions. Stencil computations are computations that update mesh structures, and the next state of a mesh depends only on the states of its neighbor meshes. Due to its inherent and coherent parallelism, stencil computation DSLs have been actively studied. For example, [4] have demonstrated ×1.5 − ×5.6 speedup on various hardwares via hardware and memoryhierarchy aware automated tuning. [5] have demonstrated generating multi-GPU codes that weakly scales up to 256 GPUs.
Still, reports are limited to simple equations such as diffusion equations and Jacobi solvers of Poisson equations, and to implement solvers of more complicated equations such as compressive hydrodynamics, magneto-hydrodynamics or general relativity, and their higher-order versions, we are forced to manually decompose the solving algorithms to imperative instructions that are often tens and hundred thousands of lines.
It is a problem common to all the languages that is designed to be compatible with C or Fortran, such as OpenMP, CUDA, PGI Fortran and OpenACC -that it is hard, if not impossible at all, to avoid repeating yourself in these languages. Despite all the efforts made so far to extend these languages, the programs in those languages fundamentally lack the ability to manipulate programs themselves and other abstract concepts. Every time a new generation of language appear, we are forced to make painful choices of porting a huge amounts of legacy codes to the new language. It is also a pain that various numerical techniques are mixed in one code, and can hardly be reused. If we want to make computers automatically combine such numerical techniques, compose a variety of implementations, and search for the better ones, the abstraction power is necessary. A complementary approach is needed here, that works together with the parallel languages.
Our contribution, Paraiso (PARallel Automated Integration Scheme Organizer ), enables us describe PDE-solving algorithms succinctly using algebraic concepts such as tensors. Paraiso also enables us to describe various manipulations on algorithms such as introducing higher-order interpolations, in composable and reusable manner. That is, once we define a certain interpolation method in our language as a transformation of basic solver to a higher-order one, we can reuse the transformation for any equations. In this way Paraiso reduces the cost of rewriting when we search for better discretizations or interpolations; enables us to port the collections of algorithms to new parallel languages, without changing the Paraiso source codes but by updating the common code generator; allows yet another layer of automated tuning at translating Paraiso codes to other parallel languages.
In writing Paraiso, we wanted to define arithmetic operations between tensors and code generator fragments. We want our tensors to be polymorphic, but we don't want to allow addition between tensors of different dimensions. We wanted to make generalized functional applications. We wanted to traverse over data structures. We need to manage various contexts, like context of code generation and the context of serialization from a program to a genome.
Programming language Haskell supports all of these. Moreover, they are supported not as built-in features that users cannot change; they are libraries that Haskell users can freely combine or create their own. This flexibility of Haskell is based on its strong, static, higher-order type system with type classes and type inference. Thus Haskell essentially allows us to develop our own type-systems within it, which gives it a unique advantage as a platform for developing embedded DSLs. Many parallel and distributed programming languages has been implemented using Haskell [6]. Nepal [7] and Data Parallel Haskell [8] are implementations of NESL, a language for operating nested arrays in Haskell. Accelerate [9] and Nikola [10] are languages to manipulate arrays on GPUs written in Haskell. Finally, Liszt [11] is a DSL for solving mesh-based PDEs based on functional programming language Scala.
Our contributions are the following: • A domain-specific language (DSL) embedded in Haskell, with which one can describe explicit solvers of partial differential equations (PDEs) in a succinct and organized manner, using tensor notations.
• A code-generation mechanism that takes the DSL and generates OpenMP and CUDA programs.
• A compressible Euler equations solver implemented in the DSL, and tuning experiments using the solver. Written in tensor notations, the solver can be applied to problems of arbitrary dimension without changing the source code at all but a single type declaration that sets the dimension of the solver.
• An annotation mechanism with which one can give hints to code generators, which makes it drastically easy to search for better implementations manually. By adding just one line of hint to the solver causes overall refactoring, adding 4 subroutines to the generated code, making it consume 1.3 times more memory but 6.42 times faster.
• An automated benchmarking and tuning mechanism based on parallel simulated annealing and genetic algorithms, which generates a lot of different implementations of the PDE solvers and search for faster implementations. It speeds up the unannotated solver by factor of 5.26, and the annotated solver further by factor of 1.78.
Just for a comparison, the Paraiso framework is about 5'000 lines of code in Haskell, and the compressible Euler equations solver implemented in Paraiso is less than 500 lines. From that, we have generated more than 500'000 instances of the solver, each being 3'000 -10'000 lines of code in CUDA. The automatically tuned codes are faster than the manually tuned codes reported by other groups. Paraiso is ready to optimize solvers of other equations, and all the solvers written in Paraiso can migrate to new parallel languages and new hardwares once the Paraiso code generator supports them.
Our work, Paraiso is the first consistent system that combine those previously studied techniques of symbolic computations, DSLs, GPU computations and automated tuning. We demonstrate the utility of such a system in the domain of explicit solvers of PDEs.
This paper is organized as follows. In section 2, we describe the overall design of Paraiso as well as its components. In section 3, we describe our automated tuning mechanism, which is a combination of genetic algorithm and simulated annealing that is designed to utilize the varying number of available nodes in a shared computer system. In section 4, we introduce the compressible hydrodynamics solver we choose as the tuning target, and describe the manual and automated tuning experiments. In section 5, we analyze the automated tuning history. In section 6 are concluding remarks and discussions.

Overall Design
Paraiso is to tackle the ambitious problem of generating fast and massively-parallel codes from human-friendly notations of algorithms. To divide the problem into major components, which may be conquered by different set of people, we set the overall design of Paraiso, as illustrated in Fig. 1.
The users of Paraiso must manually invent a discretized algorithm for the partial differential equations they want to solve. The solving algorithms are described in Paraiso using Builder Monads. Builder Monads generate programs for Orthotope Machine (OM), a virtual parallel machine designed to denote parallel computations on multidimensional arrays. Then the back end OM compiler generates the native machine source codes such as C++ and CUDA. Finally, native compilers translate them to executable codes.
When someone translates a simulation algorithm in his mind, to different languages such as Fortran, C and CUDA, he starts from building mathematical notations in his mind, and then gradually decomposes it to machine level. There should be the last common level that is independent of the detail of the target languages or the target hardwares. The level consist of very primitive operations on arrays. Orthotope Machine (OM) is designed to capture this level: The instruction set of the OM is kept as compact as possible, while maintaining all the parallelisms found in the solving algorithms. In forthcoming technical report we plan to provide the formal definition of the Orthotope Machine.
In order to generate the OM data-flow graphs from tensor expressions, and to translate them as native programs and further to apply manual and automated tuning over them, Paraiso introduces a number of abstract concepts centered around the OM. These components are quite orthogonal to each other, and some may even be useful outside the context of Paraiso. The data-flow graph for the OM. c.f. §2.2, Fig. 2, Fig. 3.
type Annotation in Language.Paraiso.Annotation The collection of annotations that are added to each OM data-flow graph node. c.f. §2.6, §3.1.
data Plan in Language.Paraiso.Generator.Plan The fixed detail of the code to be generated such as amount of memories and what to do in each subroutine. c.f. Fig. 5. To see how a plan is fixed c.f. Fig. 7.
data Program in Language.Paraiso.Generator.Claris The subset of C++ and CUDA syntax tree which is sufficient in generating codes in scope of Paraiso. c.f. Fig. 5.
newtype Genome in Language.Paraiso.Tuning.Genetic The set of annotations that belongs to an individual encoded as a string of letters, with which one can mutate, cross and triangulate. The evolution algorithm is in §3.2. interface, we can combine various stencil computation applications with state-of-theart techniques developed so far, such as cache-friendly data structures [4], overlapping communication with computation [12,5] or heterogeneous utilization of CPU/GPU [13]. On the other hand, various concepts of numerical simulations has been decomposed to elemental calculations before the OM level, so the problem of how to build a parallel computation from components such as new spatial interpolations, time marching methods and approximate Riemann solvers, can be addressed and developed separately from the detail of the hardware. To achieve such orthogonality is one of the aims of the Paraiso project.

Outline of The Orthotope Machine(OM)
The Orthotope Machine (OM) is a virtual machine much like vector computers. Each register of OM is multidimensional array of infinite size. Arithmetic operations of OM work in parallel on each mesh, or loads from neighbor cells. We have no intention of building a real hardware: OM is a thought object to capture parallel algorithms to data-flow graphs without losing parallelism.
The instruction set of OM resembles those of historical parallel machines such as PAX computer [14], and is subset of partitioned global address space (PGAS) languages such as XcalableMP [15].
Each instance of OM have a specific dimension (e.g. a two-dimensional OM of size N 0 × N 1). The variables of OM are either arrays of that dimension, or a scalar value. All of the arrays must have a common size. The actual numbers (N 0, N 1) are fixed at native code generation phase. We say that arrays are variable with Local Realm , while the scalars are variable with Global Realm.
OM has a set of Static variables which denotes the current state of the simulation. Each Static variable has a string name. OM has a set of Kernels -they are subroutines for updating the Static variables. Inside each Kernel, you can generate Temporal values in static single assignment (SSA) manner.
To summarize, the lifetime of an OM variable is either Static or Temporal, and the Realm of an OM variable is either Local or Global. Static variables survive multiple Kernel calls, while Temporal variables are limited to one Kernel. Local variables are arrays. Global variables are scalar values; in other words, they are arrays whose elements are globally the same.
OM cannot handle array of structures; Local variables may contain only simple objects such as Bool or Double, but not composite ones such as complex numbers or vectors. Nevertheless we can easily use such composite concepts in Paraiso programs at Builder Monad level and translate them to OM instructions, for example by using the applicative programming style [16].
An OM Kernel is a directed bipartite graph consisting of NInst nodes and NValue nodes, as illustrated in Fig. 2. Each node has arity :: a -> (Int, Int), number of incoming and outgoing edges. NValue nodes have arity of (1,n), where n is the number of NInst nodes that uses the value and can be an arbitrary integer. Arities of NInst nodes are inherited from the instructions they carry.
OM has nine instructions: Here, we limit ourselves to C pseudo-code of OM instruction semantics. The formal definition of the OM will be in forthcoming technical reports. Note that we do not translate the OM instructions one by one to C codes listed here. Instead, the instructions are merged into sub-graphs and then translated to C loops (or CUDA kernels) with much larger bodies, to make efficient use of computation resource and memory bandwidth (c.f. §3.1).
Imm arity (0,  Arith perform various arithmetic operations. The arity of this NInst node is inherited from its operator. The realms of the inputs and outputs must match, and can be either of Global or Local. If Local, array elements at matching index are operated in parallel (i.e. zipWith).
For example,

An Example of Orthotope Machine Data-Flow Graph
Here we show the use case of the Orthotope Machine operators introduced in §2.2 within a simple PDE solver. We solve the following linear wave equation: where t is time, x is one-dimensional space coordinate and c is the signal speed. By introducing the time derivative g = ∂f /∂t, Eq. (1) is rewritten as following system of first-order PDEs: and satisfies the following conservation law: We will write a Paraiso code that simulate the discrete system of Eqs. (2,3) and tests the conservation law. Let f n [i] denote the discrete field where n and i are the discrete time and space coordinate, respectively. We choose the following 2nd-order, Lax-Wendorf and Leap Frog scheme to solve Eqs. (2,3) : and measure the following discrete form of the conserved quantity: Table 2 shows the Paraiso implementation for the algorithm Eq. (6,7). The corresponding OM data-flow graph is visualized in Fig. 3. See how the data-flows  Table 2. For brevity, NValue nodes are represented as little circles and only NInst nodes are shown; the Arith instruction nodes are replaced by arithmetic operators inside them, and the instruction arguments are not shown but for Load, Store and Reduce. proceed :: Builder Vec1 Int Annotation () proceed = do c <-bind $ imm 3.43 n <-bind $ loadSize TLocal (0::Double) $ Axis 0 f0 <-bind $ load TLocal (0::Double) $ fieldF g0 <-bind $ load TLocal (0::Double) $ fieldG dx <-bind $ 2 * pi / n dt <-bind $ dx / c f1 <-bind $ f0 + dt * g0 fR <-bind $ shift (Vec:~-1) f1 fL <-bind $ shift (Vec:~1) f1 g1 <-bind $ g0 + dt * c^2 / dx^2 * (fL + fR -2 * f1) store fieldF f1 store fieldG g1 dfdx <-bind $ (fR -fL) / (2*dx) store energy $ reduce Reduce.Sum $ 0.5 * (c^2 * dfdx^2 + ((g0+g1)/2)^2) * dx from Load node to Store node. Here, we have c = 3.43 and 0 ≤ x < 2π. We have performed numerical simulations using this program and the following initial condition: g(x) = cos(3x), with resolution varying from N = 8 to N = 3072, and have confirmed that the fluctuation of discrete conserved quantity, Eq. (8) was smaller than of the order of 10 −13 .

Typelevel Tensor
We introduce typelevel-tensor library, to abstract over the dimensions; the use of tensor notations allow us to describe the algorithms for different dimensions in a single source code. The benefit of having information on tensor dimensions and ranks at the type-level is that we can detect erroneous operations like adding two tensors of different dimensions at compile time. Also we have much less need for explicitly mentioning the tensor dimensions in our programs, because the dimensions will be type-inferred. Our approach is similar to that in [17], where two constructors Z and :. are used to inductively define multi-dimensional tuples. Their n-dimensional vectors are actually n-tuples, so it is possible that the set of n elements consist of different types. In contrast, we needed that all the n elements of a vector are of the same type, because we wanted to operate on them, for example by applying a same function to all of them.
So instead of [17]'s approach: infixl 3 :. data Z = Z data tail :. head = tail :. head We have these: Here, Vec is the type constructor for 0-dimensional vector, and :~is a type-level function that takes the type constructor for n-dimensional vector as an argument and returns the type constructor for n + 1-dimensional vector. We define type synonyms for successively higher vector: Vec0, Vec1, Vec2, and Vec3 are all types of kind * -> *, meaning that it takes one type (the element type) and returns another (the vector type). For example, the three-dimensional double precision vector type is Vec3 Double. Higher-rank tensors are defined as nested vectors; for example Vec3 (Vec3 Int) is a 3 × 3 matrix of integers.
Since we know that all the elements of our tensor are of the same type a, we can make our tensors instances of Traversable type class [16,18]: The benefit of making our tensors instances of Traversable is that we can traverse on them: traverse :: Here, suppose t is our tensor type-constructor and f is some context -for example, a code generation context. a and b are elements of our tensor. Then the type of the traverse function means that if we have code generators for the computation of one element (a -> f b), and we have a t a, a tensor whose elements are of type a, then we can deduce the code generator for computation of the entire tensor f (t b).

Builder Monad
Builder monad is a State monad whose state is the half-built data-flow graph of the Orthotope Machine. To represent the data-flow graph, we use Functional Graph Library (FGL) [19,20]. Authors learned from the Q monad [21] how to encapsulate the construction process. The graph carried by the State monad has the following type: Graph takes three type arguments. vector :: *->* is a type constructor that denotes the dimension of the OM. gauge is the type for the indices of the arrays, which is usually an Int. anot is the type the nodes of the graph are annotated with. Such annotations are used to analyze and optimize the data-flow graph. The three types vector, gauge, and anot are passed to the Nodes of the graph. Nodes are either NValue or NInst. Two types vector and gauge are further passed to the instruction type Inst, because instructions such as shift requires the information on the array dimension and indices. Every graph nodes are also annotated by type anot.
On the other hand, the edges contain none of the three types. They are just unordered edges EUnord or edges ordered by an integer EOrd Int. For example, we can exchange two edges going into addition instruction, but cannot exchange those into subtraction.
We define various mathematical operations between Builder Monad in a consistent manner (c.f. Fig. 4). For any operator ⊕, where Value i is the value computed by Program i which is generated by Builder i. For example, a helper function that takes an operator symbol op, two builders builder1 and builder2, and create a binary operator for builder, is as follows: mkOp2 :: (TRealm r, Typeable c) => A.Operator --^The operator symbol -> (Builder v g a (Value r c)) --^Input 1 -> (Builder v g a (Value r c)) --^Input 2 -> (Builder v g a (Value r c)) --^Output mkOp2 op builder1 builder2 = do We first extract the graph nodes from left hand side and right hand side builders in Builder context, add an NInst node that contains op symbol, add an NValue node after that, and return the node index.
In Haskell, defining mathematical operators between a data type is done by declaring the data type as an instance of the type class that manages the operator. In this sense type classes in Haskell are the parallels of algebraic structures such as group, ring and field, that manages addition, multiplication, and division, respectively. We use a Haskell package numeric-prelude that provides such algebraic structures [22].
For example, an Additive instance declaration of Builder is as follows: instance (TRealm r, Typeable c, Additive.C c) => Additive.C (Builder v g a (Value r c)) where zero = return $ FromImm unitTRealm Additive.zero These type class instances, together with typelevel-tensor library, allows us to write tensor equations used in Paraiso application programs. Moreover, such equations are good for arbitrary instances of the type class. For example, here are the definition of momentum and momentum flux in our Euler equations solver: These functions can be used to directly calculate momentum vector and momentum flux tensor whose components are of type Double. The very same functions are used to generate the solvers for CPUs and GPUs. At that time, their components are inferred to be Builder types. In addition to that, these functions can handle tensors of arbitrary dimensions.

Backend
The back end converts the data-flow graph of Kernels native codes, and create a C++ class corresponding to an OM. The code generation processes of Paraiso is illustrated in Fig. 5. First, various analysis and optimizations are applied. Analysis and optimizations are functions that takes an OM and returns an OM, so we can combine them in arbitrary ways. Though some analysis are mandatory for code generation.
Paraiso has an omnibus interface for analysis and optimization using dynamic programming library Data.Dynamic in Haskell [23] : Here, analyzers as well as human beings can add annotations of arbitrary type a to the graph. On the other hand, optimizers can read out annotations of what type they recognize and perform transformations on the data-flow graph. The set of annotations can be serialized to, and deserialized from genomes, which are binary strings used in automated tuning phase.
Examples of annotations are the choices of whether to store a value on memory and reuse or not to store and recompute it as is needed (c.f. Fig. 6); the boundary analysis result used for automatically adding ghost cells; dependency analysis; and labels used for dead code eliminations.
Once the analysis and optimizations are done, an OM is translated to a code generation Plan. Here, decisions are made for how many memory are used, what portion of computation goes into a same subroutine, and so on. The nodes in data-flow graph are greedily merged as long as they have no dependence and can be calculated in the (a) Some nodes need to be Manifest.
(b) For each of other nodes, the genome specifies whether it is Manifest or Delayed.
(c) The sets of nodes that can be calculated in the same loop are greedily merged into write groups.
(d) The Kernel is divided to multiple subKernels, each of which representing one write group. Each subKernel is in turn translated to a function in C or a global function in CUDA. same loop (Fig. 7). The Plan is further translated to CLARIS (C++-Like Abstract Representation of Intermediate Syntax). CLARIS is subset of C++ and CUDA syntax which is sufficient in generating codes in scope of Paraiso.
Finally, CLARIS is translated to native C++ or CUDA codes.

Tuning Targets
The objects we want to optimize are the implementations of a partial differential equations solving algorithm. We call them individuals, adopting the genetic algorithm terminology. Each individual has a genome that encodes how it have chosen to implement the algorithm. The choices are (1) how much computation speed and memory bandwidth to use, (2) where to synchronize the computation, and (3) the CUDA kernel execution configuration. The fitness of the genome is the benchmark score of the generated code, measured in cups (the number of fluid cells updated per second) which we want to maximize. A simple example program Fig. 6 indicates that we have implementation choices for intermediate variables: whether to store its entire contents on the memory (Manifest) or not to store them and recompute them as they're needed (Delayed). The terms Manifest and Delayed are inherited from REPA [17]. If we increase the number of Manifest nodes, we consume less arithmetic units but more memory and its bandwidth; decreasing the number of Manifest nodes have the opposite effect. There are two extreme configurations, one is making as many nodes Manifest as possible and the other is making as few nodes Manifest as possible. In most cases both of them result in poor performance and moderate configurations are faster. Paraiso generate codes for many possible combinations of Manifest / Delayed choices (Fig. 7) and searches for such configurations by automated tuning.
Another tuning done by Paraiso is the choice of synchronization points. In CUDA, inserting syncthreads(), especially before load/store instructions cause the next instruction to coalesce and increase the speed of the program. Inserting too much synchronization, on the other hand, is a waste of time. Again, Paraiso searches for better configurations by automated tuning.
The last tuning done by Paraiso is to find the optimal CUDA kernel execution configuration, i.e. how many CUDA threads and thread blocks to be launched simultaneously. These tuning items summed, the genome size is approximately 6000 bits for our hydrodynamics solver, which means that there are 2 6000 possible implementations. Brute-force searching for the fastest implementation from this space is inutile and in the first place impossible. Instead, the goal of Paraiso is to stochastically solve such a global optimization problem in this genome space, where the function to be optimized is the benchmark score of the PDE solver generated from the genome. proceed = do x <-bind $ loadLD density y <-bind $ x * x z <-bind $ y + y store density z proceed = do x <-bind $ Anot.add Sync.Pre <?> loadLD density y <-bind $ Anot.add Alloc.Manifest <?> x * x z <-bind $ y + y store density z The generated .cu program file, abbreviated:  Table 3. An example of changes in the genome and the implemented algorithm caused by annotations. Table 3 shows a simple Paraiso kernel and how its genome and implementation is altered by adding annotations. This kernel named proceed performs the following calculations: It first loads from an array named density, performs a multiplication, then an addition, and then stores the result to density again. The header file on the right side of the Table 3, compared to the left one, allocates an additional device vector and declares an additional subkernel as results of a Manifest annotation. The right .cu file, compared to the left one, differs in two points: one is that it performs the multiplication y = x * x and the addition z = y + y in separate CUDA kernels and stores the intermediate result to the additional device vector, which is yet another result of the Manifest annotation. Another difference is that it calls syncthreads() just before the load instruction, which is the result of the synchronization annotation.

Parallel Asynchronous Genetic Simulated Annealing
Simulated annealing has been widely used to solve global optimization problems. However, a standard simulated annealing method is not suitable for parallelization, and if the annealing schedule (how we cool down the temperature as function of time) is too quick, it tends to fall into local minima. Replica-exchange Monte Carlo method [24] solves these drawbacks by introducing multiple replicas of heat baths with different temperatures, and by allowing replica exchange between adjacent heat baths. Thus, replicas can be computed in parallel, and the annealing schedule is spontaneously managed by replica migration.
We further extend this method to fit into benchmark based tuning on shared cluster computer systems. First, in shared systems it is hard to maintain a fixed number of replicas because the available amount of nodes changes with time. Second, it is not efficient to perform replica-exchange synchronously, since the wall clock time required to calculate the fitness function varies with replicas. For these two reasons, we modify replica-exchange Monte Carlo method to a master/worker model (c.f. Fig. 8), where a master asynchronously launches varying number of workers in parallel.
The master holds the genomes and scores of all the past individuals in a database (DB). The role of the master is to draw individuals from the DB, to create new individuals using their genomes, and to launch the workers when the computer resource is available. The role of the workers, on the other hand, is to generate codes from the given genome, to take the benchmark, and to write the results into the DB.
In this way we can launch any number of workers asynchronously as long as the DB is not a bottleneck. In addition to that, we eliminate a common weak point of simulated annealing and genetic algorithms -that individuals of older generation are overwritten and are inaccessible.
For each individual I, the generated code is benchmarked 30 times. We record the mean µ(I) and the sample standard deviation σ(I) of the score. It is important to record the deviation. When we benchmark each individual only once, some individuals receive over-evaluated score by chance, infest the system and tend to stall the evolution.
At the end of each benchmark, the individual is briefly tested if the state of the simulation has developed substantially from the initial condition and there is no NaN (not a number). The test is chosen because not evolving at all and generating NaN are two dominant modes of failure, and the individual test time needs to be kept smaller than benchmark itself. If an individual fails a test, its score is 0. At the end of a tuning experiment, the champion individual is extensively tested if it can reproduce various analytic solutions, which takes hours.
A draw is the operation to randomly choose an individual from the DB. Each draw has a temperature T . In a draw of temperature T , the probability the individual I is where I is the individual with the largest µ(I).
Low-temperature draws strictly prefer high-score individuals while hightemperature draws does not care the score too much. And the difference small compared to σ(I ) + σ(I) is gracefully ignored.
Every time master creates a new individual, it chooses the draw temperature T randomly, so that the probability density of log T is uniformly distributed between log(σ(I )) and log(µ(I )).
To summarize, our method is a master/worker variant of replica-exchange Monte Carlo method [24], using genetic algorithms as neighbor generators. Thus, our method • can utilize parallel and dynamically varying computer resource.
• can find global maximum without hand-adjusted annealing schedule.
• can combine independently-found improvements.
An extensive review on the applications of the evolutionary computation in astronomy and astrophysics is found in [25].

Three Methods of Birth for Generating New Individuals
We use three different methods of birth to generate new individuals (c.f. Fig. 9).  The first method is mutation. We draw one individual from the database (DB), take its genome, overwrite it randomly and create a new genome. The second method is crossover. We draw two individuals from the DB, split their genomes at several random points and exchange the segments. The third method is triangulation, or three-parent crossover.
Multi-parent crossovers are crossovers that takes more than two parents and create one or multiple children. The proposed multi-parent crossover methods include taking bit-wise majority vote [26,27] and exchanging segments between multiple parents [28]. Multi-parent crossovers for real-number coded genomes are also studied [29]. The novelty of our three-parent crossover (triangulation) is to take the scores of the parents into consideration. Triangulation is designed to efficiently combine independent improvement found in sub-problems.
FFTW [1] uses the divide and conquer approach to its target problems. It first recursively decomposes large FFT problem into smaller ones and solves the optimization problem from smaller part. In contrast, Paraiso deals with monolithic data-flow graphs. It is not obvious how to decompose them into subgraphs -decompositions are actually the target of optimization. Therefore, we optimize subgraphs in vivo; we optimize subgraphs keeping them embedded into the entire graph.
Aiming to merge two distinct optimized subgraphs, we draw three individuals from the DB. Then we sort them in ascending order of the score and name them Base, Secondary and Primary. We perform bit-wise operation as in Table 4 to create their child. For each bit, if at least one of Primary or Secondary has been changed from Base, then we adopt the change.
When creating a new individual, the master chooses mutation, crossover, or triangulation at equal probability of 1/3. Then the master draws the needed number of parents from the DB. If two of the parents are the same the master retries the draws. Otherwise, the master creates a child with the chosen method of birth. Then if the genome of the newly created individual is already in the DB, the master discards the individual and retries the mutation until it gets a genome not found in the DB. With each retry, draw temperature T is multiplied by 1.2, so that the master eventually succeed in creating a new genome.

The Target Program
We implemented a 2nd order compressible hydrodynamics solver [30] in Paraiso, and then optimize it. Here we detail numerical scheme. Let where ρ , v a , m a , p and E are the density, the velocity, the momentum, the pressure and the total energy, respectively, The primitive and conserved variables are related by m a = ρv a , E = E i (p) + 1 2 ρv 2 where E i is the internal energy of the gas. By assuming the adiabatic equation of state for perfect gas, E i (p) = (γ − 1) −1 p where γ is the ratio of specific heats.
We numerically solve the Euler equations ∂ t U = ∂ a F a , or more specifically: The numerical method we used is as follows, based on a function Flux a (U0) which calculates the fluxes across the cell surfaces using a Riemann solver. Flux a (U0) is defined as follows. First, V0 are the primitive variables calculated from U0: The interpolated primitive variables VL and VR are where we use piecewise linear interpolation with minbee flux limiter [31,30]: Then, fluxes across the boundaries are defined using the HLLC Riemann solver [32].
This is the flux defined by Flux a (U0). Then, linear addition of flux to conserved variable AddFlux(∆t, F a , U) is defined by where ∆r a is the mesh-size along the a-axis. Using these notations, we construct the second-order time marching as follows, where ∆t is the time step determined by the CFL-condition: U2 is the set of conserved variables for the next generation. Notice how we rely on human mind flexibility when we convey the algorithms in the forms like above. In languages like Fortran or C, we are often forced to decompose those expression into elemental expressions and the source codes tend to become longer. In Haskell, we can express these ideas in well-defined machine readable forms while keeping the compactness and flexibility as in above.
We declare both the primitive variables Eq. (13) and conserved variables Eq. (14) as instances of the type class Hydrable, the set of variables large enough for calculating any other hydrodynamic variables. We used the general form by default, and bind either the primitive or conserved variables when we need the specific form. By doing so requisite minimum number of conversion code between primitive and conserved variables are generated.
When we defined the minbee interpolation for a triplet of real numbers Eq. (21) and then applied it to primitive variables Eq. (20), we implicitly extended a function on real numbers to function on set of real numbers. In Haskell, we can define how any such generalized function application over the set of hydrodynamic variables should behave, just by making it an instance of Applicative type class.  Table 5.
The codes annotated by hand and their performances. First four columns are the ID of the individuals, the CUDA kernel execution configuration, and the two choices of annotation -whether to Manifest or Delay (1) the result of spatial interpolation (2) the result of the Riemann solver. Next three columns are the properties of the generated codes as results of these choices made. They are namely the size of the code (number of lines), the number of subKernels in the code, the memory consumption in proportion to fluid cell number N . The last two columns are the speed of the generated codes for single precision (SP) and double precision (DP). The speed of the codes are measured by Mcups (10 6 cell update per second).
When we calculate the space-derivatives of the fluxes in Eq. (23), the component of the flux we access (index a in F a ) and the direction in which we differentiate (indices a in [i − 1 2 e a ] and [i + 1 2 e a ]) should match. Although the flux F and the spatial array indices i have very different types, Haskell's type inference guarantees that they are both tensors of the same dimensions, and we can access both of them by the common tensor index a. And we can sum over a, because tensors are Foldable and their components are Additive. Tensors are Traversable as mentioned before, and any type constructors that are Traversable are also Foldable.

Annotating By Hand
Before we start the automated tuning experiment, we generate several individuals by hand. The initial individual Izanagi is the one with the least number of Manifest nodes as possible. Izanami is the same individual with the CUDA execution configuration suggested by the CUDA occupancy calculator. Based on it, we add several Manifest annotation and create new individuals (c.f. Table 5). Manual annotations are not blindsearch process; each annotation has clear motivation such as "let us store the result of the Riemann solvers because it is computationally heavy", "let us store the result of interpolations because it moves a lot of data," etc. For example, the sample code in Table 6 shows the implementation of the piecewise linear interpolation with the minbee The interpolation code for Izanami: interpolateSingle order x0 x1 x2 x3 | order == 1 = do return (x1, x2) | order == 2 = do d01 <-bind $ x1-x0 d12 <-bind $ x2-x1 d23 <-bind $ x3-x2 let absmaller a b = select ((a*b) 'le' 0) 0 $ select (abs a 'lt' abs b) a b d1 <-bind $ absmaller d01 d12 d2 <-bind $ absmaller d12 d23 l <-bind $ x1 + d1/2 r <-bind $ x2 -d2/2 return (l,r) For Iwatsuchibiko, the last line is modified as follows: return (Anot.add Alloc.Manifest <?> l, Anot.add Alloc.Manifest <?> r) Table 6. An annotation made in Paraiso source code that causes the results of the interpolations to be stored on memory and reused. flux limiter, Eq. (21) in Paraiso, and how to annotate the return values of the limiter as Manifest.
We use Izanami as the base line of the benchmark, and use Izanagi as the initial individual of some experiment to see if the automated tuning can find out the optimal CUDA execution configuration by itself.
The codes are benchmarked on TSUBAME 2.0 cluster at Tokyo Institute of Technology. Each individual was benchmarked on a TSUBAME node with two Intel Xeon X5670 CPU(2.93GHz, 6 Cores × HT = 12 processors) and three M2050 GPU(1.15GHz 14MP x 32Cores/MP=448Cores).
The abstraction power of Builder Monad lets us change the code drastically with little modification. For example, Paraiso source code of Shinatsuhiko and Izanami differs by just one line of annotation. This introduces the 16 bits of difference in their genome, causing Shinatsuhiko generate a code that has 1/3 lines, which contains four more subroutines, consume 1.31 times more memory and is 6.42 times faster.

Automated Tuning
Next, we performed several automated tuning experiments c.f. Table 7.
To distinguish the contribution of the three methods to generate new individuals, we also performed automated tuning experiments with either crossover or triangulation turned off. Table 8 shows the initial possibility of the master node attempting each RunID prec.
initial score wct best ID/total high score GA-1 DP  Table 7.
The statistics of GB-* experiment series. The columns are RunID, precision, the score of initial individual, the ID of the best individual and the number of individuals generated, the high score (in Mcups). Experiments started with Shinatsuhiko. The nine experiments were performed in parallel, and took 92252 minutes of wall-clock time. method of birth. Note that the actual frequencies of crossover and triangulation are smaller than these value because the master may default to mutation. Table 9 shows the result of experiments. The resolutions were 1024 2 for GA-1 and GA-S1, problems, 100 3 for GA-3D problems and 512 2 for other GA-* series. The resolutions were 512 2 in GB-* series. The automated tuning system can generate and benchmark approximately 10 000 individuals per day. 20 − 100 workers were running at the same time. It takes a few days to tune up Izanami to speed comparable to Shinatsuhiko, or speed up Shinatsuhiko by another factor of 2. The best speed obtained was 35.3Mcups for double precision, and 53.7Mcups for single precision. Our automated tuning experiments on 3D solvers mark 42.4Mcups SP. These are competitive performances to hand-tuned codes for single GPUs; e.g. Schive et. al. [33] reports 30Mcups per C2050 card (single precision, note that their code is 3D. Asuncióna et.al. [34] reports 6.8Mcups per GTX580 card (single precision, 2D). Fig. 10 shows the result of convergence tests for individuals GA-1.20756, GA-4.33991, GA-D.59841, GA-DE.41250, GA-S1.33958, Izanami and Shinatsuhiko. The resolution N is varied from 16 2 to 1024 2 . The codes are tested for entropy wave propagation, sound wave propagation and Sod's shock-tube problem. The detail of test initial conditions are as follows. For all tests, numerically solved domains are (0 < x < 1, 0 < y < 1), out of which the analytic solutions are continuously substituted as boundary conditions. The system was numerically developed until t = 1 for entropy wave and sound wave problems, and t = 0.125 for shock tube problem, and then the numerical solutions were compared to analytic solutions using the norm of L 1 error vector.
The initial conditions are, for entropy wave problem: for sound wave problem: where the amplitude A = 10 −5 , and for Sod's shock-tube problem: x > 0.5 While the "coaxial" tests used the above initial conditions as they are, the "oblique" tests used the initial conditions rotated about (x, y) = (0.5, 0.5) for 60 • .

Overview of The Simulated Evolution and Analyses
The functions of score against evolution progress exhibit self-similar structure of repeated cliff and plateau. For example, see the evolution history of experiment GA-D illustrated in Fig. 11 and its two enlarged views Fig. 12 and Fig. 13, with the family tree of the best individual superimposed. The evolution of the high scores in the experiments are shown in Fig. 14.
To design more efficient auto-tuning strategies, we investigate what have happened and which part of the evolution contributed to create better individuals in our automated tuning experiments. In section 5.2, we measure the contributions from three different tuning items. In section 5.3, we classify the individuals by their aspects such as their method of birth, their distance to the champion in the family tree, and their fitness relative to their parents. In section 5.4 we study how these classes contribute to the evolution by performing correlation analyses among these classes. In section 5.5, we study how the method of birth of parents affect their children. In section 5.6 we summarize and conclude the analyses. The ID and the scores (in cups) of all the individuals generated in automated tuning experiment GA-D. The individual with the highest score as well as its ancestors are marked by crosses, and each of them is connected to its parents with lines.   The average of GB-370-* was always higher than that of GB-307-*, and the average of GB-333-* was always higher than that of GB-370-* for the 89% of the time.

Contributions of The Three Genome Parts
We assign the symbols to genome parts with different functions as follows: (C): the CUDA kernel execution configuration, (M): which data to store on the memory (to make them Manifest), and (S): when to synchronize the computation. We first investigate how these three components contributed to the score by component-wise artificial crossover between the initial individual and the best scoring individual (Table  10). In the case of GA-S1.33958, introducing improvement only in C, M, S part increase the score by 14%, 30%, and 0%, respectively. In the case of GA-4.33991, the increase are 0%, 85%, 0%. Both cases exhibit synergy effect. Introducing several modifications simultaneously have more effect than the sum of the separate effects, they multiply. So addition in log-space explains 98% and 88% of the progress for GA-S1.33958 and GA-4.33991, respectively, but the score of the final individuals are still slightly higher than predicted.
The score of the individuals created by artificial crossover between the initial individual I 0 and the best scoring individual I . The second to fourth columns indicate which component was taken from which individual. Columns C,M,S correspond to CUDA kernel execution configuration, Manifest/Delay choice, synchronization timing, respectively. For each individual I the fifth column shows µ(I) ± σ(I), the sixth column shows µ(I) µ(I ) − µ(I 0 ) ± σ(I) µ(I ) − µ(I 0 ) , and the seventh column shows log µ(I) − log µ(I 0 ) log µ(I ) − log µ(I 0 ) ± σ(I) (log µ(I ) − log µ(I 0 ))µ(I) .
100%, and 14%, respectively. Removing C,M,S part from GA-S1.33958 decreases the score by 71%, 87%, and −3%, respectively. Again we see the synergy effect, except that GA-S1.33958 were faster if S part were removed. We conclude that the Manifest/Delayed trade-off plays the central part in improving the score. Fixing the Manifest/Delayed nodes, determines the decomposition of the data-flow graph. Then tuning the synchronization timing and CUDA kernel execution configuration help further improve the score.

Classifications of Individuals
To measure how each individual I contributed in generating one of the best individuals, we define contribution distance d(I), and classify the individuals to those who contributed to the evolution and those who didn't. To begin with, let P(I) be the set of individual I's parents. (we use p for probability.) For individual I born by mutation, crossover, and triangulation, the size of the parent set n(P(I)) is 1,2, and 3, respectively.
We define d(I) as follows: • d(I ) = 0 where I is the individual whose µ(I) was the largest in the history.
• d(I) = 0 if one of I's children I c satisfies d(I c ) = 0.
We say that I is [ ] are children significantly faster than any of their parents; [ ] are children whose scores are comparable to their fastest parents; [ ] are children significantly slower than any of their parents; and [≤] are children significantly slower than their best parent but not significantly slower than their slowest parents. Note that children of rank [≤] are never born by mutations since there is only one parent. Table A8 -A14 shows the classifications for experiments.

Chi-squared Tests of Correlations between Classes
To study which set of individuals are contributing to produce best species, we perform Pearson's chi-squared test on pairs of predicates on individuals; see Table A15 and A16. We use 10 −3 significance, or X 2 > 10.83 unless otherwise mentioned. Our observations are as follows: • (row 1.) Significant negative correlations are detected between being born of mutation (n(P(I)) = 1) and being d(I) = 0 for all experiments.
• (row 2.) Significant positive correlations are detected between being born of crossover (n(P(I)) = 2) and being d(I) = 0 for all experiments but GA-S1.
• (row 3.) Significant positive correlations are detected between being born of triangulation (n(P(I)) = 3) and being d(I) = 0 for all experiments.

Correlation along the Family Line
Next, we analyze the correlation along the family line by interpreting the family tree as Markov processes (c.f. Table 11). For each individual I such that d(I) = 0, we trace back its family tree in all possible ways for n steps and obtain an n-letter word. For example, "312" means that I is born of crossover, and one of its parent was born of mutation from an individual who was born of triangulation. With such bag of words obtained, we investigate how the last letter of a word is correlated with letters in front of them. The second columns show X 2 for null hypothesis "The bag of words is a 0th order Markov chain," i.e. in the two-letter words the second letter is not correlated to the first. Formally written, the null hypothesis is: where a, b, c, ... denote characters and ab, abc, ... denote words. The third columns show X 2 for null hypothesis "The bag of words is a 0th order Markov chain," i.e. the probability of the occurrences of the letters depend only on their immediate predecessors and the distribution of three-letter words can be determined  Table 11. Chi-squared test of the family tree being lower-order Markov processes. The each column of the table shows the X 2 statistics of the null hypothesis the family tree being a n-th order Markov process and having no longer correlation.
from the distribution of two-letter words. Formally written, the null hypothesis is: Table 11 shows the result of the Markov chain analysis. Note that, the degree of freedom k for 0th-order and 1st-order Markov chains are 4 and 16, respectively, and the chi-squared values for 10 −3 significance are X 2 > 18.47 and X 2 > 37.70, respectively. From the table, we can deny the 0th-order Markov chain model for all of the experiments, and deny the 1st-order Markov chain model for 8 out of the 10 experiments. We also studied whether an occurence of a character is correlated to its prefix word. "3" was significantly likely to be followed by another "3" in all of the experiments, while "2" was significantly unlikely to be followed by another "2" in 6 out of 10 experiments. On the other hand, three consecutive occurences of the same letter was not significant at least in 8 out of 10 experiments.

Summary of The Analyses
Mutation is not efficient in directly producing individuals with good scores. Nevertheless mutations are indispensable part of the evolution because it is the only way that introduces new genomes to the genome pools and that have chance of finding new improvements.
Triangulations and crossover are two different methods to combine independentlyfound improvements, and are efficient in directly producing good individuals. Triangulations are relatively more efficient in producing d(I) = 0 individuals. However, crossover is relatively good at making large jumps that belongs to [ ] . Triangulations are working in another way; the statistics indicate that one triangulation are likely to be followed by another triangulation in a family tree, and the sequence of Triangulations work to accumulate minor improvements that belongs to [ ] to make a larger one. Fig. 14 -although the number of samples are too few to make a statistically significant statement -suggests that having both mode of combiner is better than having only one.
The main source of performance improvement was the change in memory layout and subroutine re-organization caused by making different Manifest/Delayed trade-offs. Tuning the synchronization timing and CUDA kernel execution configuration provided additional improvements.
The case of GA-S1.33958 is an evidence of the current genetic algorithm not finding the optimal individual. Random mutations tend to increase the genome entropy even under the selection pressure. Adding the rule-based mutations such as "remove all synchronization," which are improbable under random mutations, may contribute to further optimization.

Conclusion and Discussion
We have designed and implemented Paraiso, a domain specific language for describing explicit solvers of partial differential equations, embedded in Haskell. Using the typelevel-tensor and algebraic type classes, we can explain our algorithms with simplicity of mathematical equations. Then the algorithm is translated to OpenMP code or CUDA code. The generated code can be optimized both by applying annotations by hand, and by automated tuning based on a genetic algorithm. Although we present just one example here -a second order Euler equations solver -the front end of Paraiso with typelevel-tensor and Builder Monad readily accept various other algorithms. The code generation and automated tuning can revolutionize the way we invent, implement and optimize various partial differential equation solvers.
The website of Paraiso is http://www.paraiso-lang.org/wiki/, where you can find codes, documentations, slides and videos.
Making more efficient searches for fast individuals is an important future extension of Paraiso. One way is to use machine learning, for example to make suggestions for genomes to introduce, or to predict the scores of the individuals before benchmarking them. Another approach is careful selection of the tuning items, importance of them given either by hand or by inference from the syntactic structure of Paraiso source code. The computing time invested in automated tuning is only justified if the optimized program is used repeatedly, or if optimized single-GPU program is used as a "core" of a multi-GPU program.
Making the automated tuning results more flexible and re-usable is another future challenge. For example, the size of the array e.g. (N 0, N 1) in §2.2 for the twodimensional OM, is fixed before native code generation under the current design of Paraiso. To change the size N 0 or N 1, the users need to re-run the code generation process. One solution is to improve the representation of the constant values in Paraiso, so that the users can set the values that are constant within a run but may vary across different runs. Possible examples of such constants include the OM size and the parameters of the algorithm.
On the other hand, the change in the dimension of the OM e.g. to three; (N 0, N 1, N 2), or the change in the algorithm described by the Builder Monad will change the number of nodes in the OM data flow graph and the bits in the genome. Therefore, the automatically tuned genomes before the change become useless. The genome must be redesigned to perform more re-usable, general automated tunings.
Generating codes for distributed computers / distributed and accelerated computers is another important future extension of Paraiso. To this end, instead of generating MPI codes by ourselves, we plan to generate parallel languages such as X10 [35], Chapel [36], and XcalableMP [15], or domain specific languages for stencil computations such as Physis [5].
Finally, the potential users and developers of Paraiso are obstructed by the programming language Haskell, which is a very computer-science language and is far from being mainstream in the fields of simulation science and high-performance computing. This is a disadvantage of Paraiso compared to other embedded DSL approaches based on more popular syntax such as Python (e.g. [37].) We make two remarks on this.
On the one hand, we show in this paper that abstract concepts developed by computer scientists and found in Haskell, is useful and was almost necessary in implementing a DSL that consistently handles the translation from mathematical notations to evolutionary computation. We hope that such powerful DSLs appear in many different fields, with help of the findings in computer science. On the other hand, we want many other researchers to become developers and users of Paraiso without requiring much knowledge in computer science.
For promotion of such development and use of Paraiso, we plan to provide interfaces to the code generation and the tuning stages in Paraiso, in terms of strings or in lightweight markup languages, along with their specifications. Those will allow many other programming languages can access those stages, such as OM data-flow graphs, their annotations, and the genomes. Adding a simulation scientist friendly scripting language layer on top of Builder Monad is another future work, so that people not familiar with Haskell can access Paraiso. We also plan to publish a series of tutorials and example programs in the Paraiso website. Such collaborations are fundamental to the realization of the future works listed in this section and computer-assisted programming in many fields of simulation science.  Table A1. Contribution distance analysis for experiment GA-1 : The numbers of individuals with specific d(I) that is born of mutation, crossover, and triangulation. The last column shows the total numbers of individuals with specific d(I). The decimal numbers inside the parentheses are the ratio of the number of individuals with specific d(I) of specific birth, divided by total number of individuals of specific birth. Note that the grand total is slightly smaller than that in Table 7; this is because some worker nodes failed to write into DB.   Table A8. Children relative fitness for Experiment GA-1: The individuals classified based upon how they were born and their score with respect to their parents. For each class, the first row shows the number of individuals of that class, and the third row shows the number of individuals of that class with d(I) = 0. The second and fourth rows are first and third row divided by the total number of individuals born in that way. The fifth row is the contribution ratio, the third row divided by the first row (or fourth divided by second).