Boubacar Diouf

Albert Cohen

Fabrice Rastello

A Polynomial Spilling Heuristic: Layered Allocation

Keywords: compilers, algorithms, performance, optimizations, register allocation compilateurs, algorithmes, les uns avec les autres optimisation, allocation de registres

Register allocation is one of the most important, and one of the oldest compiler optimizations. Its purpose is to map temporary variables to either machine registers or main memory locations and explicit load/store instructions. The latter option is referred to as spilling.

This paper addresses the minimization of the spill code overhead, one of the dicult problems in register allocation. We devised a heuristic approach called layered. It is rooted in the recent advances in SSA-based register allocation. As opposed to the conventional incremental spilling approaches, our method incrementally allocates clusters of variables. We describe a new polynomial method, the layered-optimal allocator, and demonstrate its quasi-optimiality on standard benchmarks and on two architectures.

Introduction

Register allocation is an important compiler optimization. Its goal is to map temporary variables in a program to either machine registers or memory locations. Register allocation is subdivided into two sub-problems: rst, the allocation selects the set of variables that will reside in registers at each point of the program; then, the assignment or coloring picks a specic register where a variable will reside. Usually, all the variables of code cannot reside in registers. The variables not held in registers should reside in memory, these variables are called spilled variables. The spilling problem [START_REF] Farach-Colton | On local register allocation[END_REF][START_REF] Bouchez | On the complexity of spill everywhere under ssa form[END_REF] decides which variables should be stored in memory to make the assignment possible; it aims at minimizing the overhead of loads and stores. The coalescing [START_REF] Bouchez | On the complexity of register coalescing[END_REF] and alienation (when repairing is enabled [START_REF] Colombet | Graph-coloring and treescan register allocation using repairing[END_REF]) problems aims at minimizing the overhead of moves between registers. Spilling and coalescing are correlated problems that are, in classical approaches, done in the same framework. Live-range splitting (i.e., adding register-to-register moves) to reduce register pressure is sometimes considered in such a framework [START_REF] Cooper | Live range splitting in a graph coloring register allocator[END_REF], but it is very hard to control the interplay between spilling and splitting or coalescing.

Building on the properties of the static single assignment form (SSA), it is now possible to decouple the allocation from the assignment. Indeed, the interference graph of a program in SSA form is a chordal graph [START_REF] Hack | Register allocation for programs in SSA-form[END_REF]. Since coloring a chordal graph is easy, it follows that the assignment problem is also easy. Finding a valid coloring whenever it exists can thus be solved optimally with a greedy, linear algorithm on chordal graphs, called tree-scan [START_REF] Colombet | Graph-coloring and treescan register allocation using repairing[END_REF]. It follows the spirit of the linear-scan [START_REF] Poletto | Linear scan register allocation[END_REF], but applied to the dominance tree instead [START_REF] Rong | Tree register allocation[END_REF]. Thus, performing register allocation under SSA has led to new approaches where the remaining dicult problems, spilling and coalescing, are treated separately. When spilling, MaxLive, the maximal number of variables simultaneously live at a program point, is used as a criterion to guarantee that the forthcoming assignment will be performed without any spill. If MaxLive is lower or equal to R, the number of available registers, then all the variables will be assigned without any spill. This decoupled approach is advocated by Fabri [START_REF] Fabri | Automatic storage optimization[END_REF], Appel and George [START_REF] Appel | Optimal spilling for CISC machines with few registers[END_REF], and Hack [START_REF] Hack | Register allocation for programs in SSA-form[END_REF].

Apart from allowing the design of more ecient coalescing heuristics [START_REF] Bouchez | Advanced conservative and optimistic register coalescing[END_REF], the main advantage of this decoupled approach concerns the spilling problem: checking if the register pressure, MaxLive, is low enough is much simpler than checking the colorability of a general graph.

Because of this, existing graph based heuristics use node degree to guide the spilling decision.

This can lead to spilling a variable because its corresponding node has a high degree, while it is not live at any point of high register pressure. In other words, spilling this variable is useless in helping the assignment problem anyhow. This point can be illustrated by Figure 1: variable a 2 has 6 neighbors of high weight (spilling each variable h i is very costly), so spilling it looks like a good idea for the graph coloring. But in terms of register pressure, there are no more than three variables simultaneously live inside the loop: liveness set on the control ow graph, provides this information very naturally and a decoupled approach does not require spilling neither a 2 nor any h i .

This observation has led several researchers to design program-based heuristics to lower register pressure, opposing the new decoupled approach to the old graph-based spilling heuristics.

Remember that the decoupling approach eases the question of whether spilling a variable is useful or not, but did not make the optimization of allocation problem polynomial (nding a set of variables to be spilled of minimum cost). A simple heuristic consists in considering program points one after another, and when the register pressure at the current point is too high, incrementally spilling some variables to lower it. This incremental scheme needs a notion of protability to choose which variable to spill among the set of all live variables at a given point. To this end, Belady's furthest rst strategy works very well on an interval grapha single basic block in SSA: the idea is to consider spilling the variable which live-range goes furthest to be the most prof-

RR n°8007

itable. Consider our running example again where a 1 and a 2 have been coalesced into variable a.

The generalization of the notion of furthest use to a general control ow graph would consider a to be more protable than d. The furthest use is not the notion we want. Instead, protability of spilling a variable should be related to the number of program points of high register pressure within the live-range of the variable. Here a would cover three high register pressure program points, while d would cover ve. Going a little bit further, we understand that protability is not just about coverage of high register pressure program points: spilling a variable is protable because it avoids spilling some other variables, and because spilling it is less costly than spilling those other interfering variables. But what if we have to spill these interfering variables anyhow?

As an example, spilling a avoids spilling d and e, but as it is very protable to spill d (more than e), e should impact more than d the protability of a, leading to an inductive denition of protability... To break out of this loop, let us recall that under the SSA form, there exists a perfect mapping between maximal cliques and live variables at a given program point [START_REF] Hack | Register allocation for programs in SSA-form[END_REF]. In other words, a maximal clique cover (which is polynomial for a chordal graph) allows to express the notion of register pressure, exactly for chordal graphs and reasonably accurately for non SSA graphs occurring in the real world [START_REF] Magno | Register allocation via coloring of chordal graphs[END_REF]. On our running example, d would be part of two maximum cliques of size higher than three, while a would be in only one. This debunks the main motivation for not using a graph based approach. 1. The pseudo-polynomial complexity of allocation in the number of registers [START_REF] Bouchez | On the complexity of spill everywhere under ssa form[END_REF] suggests a heuristic that solves (optimally) roughly R over step allocation problems on step registers each. The nal allocation being the layered of the stepwise allocations, we call it the layered-optimal heuristic.

… = c + 1 e = b + 1 f = e + 1 … = c + 1 e = b + 1 f = e + 1 a 2 = a 1 h 1 = a 2 + 1 h 2 = h 1 + 1 a 2 = a 1 h 1 = a 2 + 1 h 2 = h 1 + 1 g = d + e … = d
h 3 = h 1 + 1 h 4 = h 2 + 1 h 5 = h 3 + 1 h 6 = h 4 + 1 h 1 = h 5 + 1 h 2 = h 6 + 1 h 3 = h 1 + 1 h 4 = h 2 + 1 h 5 = h 3 + 1 h 6 = h 4 + 1 h 1 = h 5 + 1 h 2 = h 6 + 1 … … livein = { a 2 } liveout = { Ø} liveout = { Ø} a 2 h 1 h 2 h 3 h 4 h 5 h 6
2. Stepwise optimality does not guarantee an overall optimal allocation, but we will show that it comes very close to optimal, even with step = 1. Intuition for this comes from recent work by Diouf et al. [START_REF] Diouf | Split register allocation: Linear complexity without the performance penalty[END_REF], observing that allocation decisions tend to be a monotonic function of the number of registers.

This approach is still incremental, but it allocates layers of variables instead of spilling one variable at a time. Thanks to the pseudo-polynomial property mentioned above, the choice of this set constituting a layer can be decided optimally in polynomial time.

To make a long story short, this paper addresses the spill-everywhere problem in a decoupled context. We introduce layered allocation a new strategy that incrementally allocates variables instead of incrementally spilling variables. We evaluate our approach in the context of decoupled register allocation.

The outline of the paper is as follows. Section 2 presents the rationale for our new approach in detail. Section 3 surveys the important concepts and results upon which our approach is built.

Section 4 presents our layered-optimal allocator for SSA programs. Section 5 adapts this scheme into a non-optimal heuristic for general, non-chordal interference graphs. Section 6 evaluates the algorithm and compares it with state of the art methods. Section 7 discusses related work and Section 8 concludes the paper.

A Graph-Based, Incremental Allocation Approach

Our approach is motivated by the facts that register lowering is pseudo-polynomial in the number of registers and stepwise allocation is quasi-optimal. We explain here in details these two observations, but before that, we rst explain why we think that the spill everywhere problem is relevant.

2.1

Why Spill Everywhere?

The spilling problem can be considered at dierent granularity levels: the highest, so called spill everywhere, corresponds to considering the live range of each variable entirely. A spilled variable will then lead to a store after the denition and a load before each use. Of course, in practice, if the variable can stay in a register between two consecutive uses, a load is saved. The nest granularity, so called load-store optimization, corresponds to optimize each load and store separately. The latter, also known as paging with write back, is NP-complete [START_REF] Farach-Colton | On local register allocation[END_REF] on a basic block, even under SSA form. The spill-everywhere problem is much simpler, applicable to justin-time compilation, and many instances are polynomial under SSA form [START_REF] Bouchez | On the complexity of spill everywhere under ssa form[END_REF]. The algorithms we propose can be applied to both spill everywhere and load-store optimization problems. We focus here on the former for its simplicity, because our past experience summarized in the 4 following points tends to conrm the practical eectiveness of the spill everywhere problem: RR n°8007

1. The complexity of the load-store optimization problem comes from the asymmetry between loads and stores. Also, most SSA variables have only one or two uses in practice, and the cost of the store favors spilling the entire live range instead of two sub-ranges of dierent variables.

2. The queuing mechanism present in most architectures behave like a small, extremely fast cache. But it is highly sensitive to the number of simultaneously spilled variables.

3. In the other extreme situation where stores have no cost, a variable can be considered to be either in memory or in register but not in both. Such a formulation [START_REF] Appel | Optimal spilling for CISC machines with few registers[END_REF] is strictly equivalent to a spill everywhere formulation where live ranges are split at every use.

4. Last, a solution to the spill-everywhere problem gives to a load-store optimization problem the global view lengthily discussed so-far that existing heuristics lack. In other words a spill-everywhere solution can play the role of an oracle.

Allocation Instead Of Spilling

After giving the reasons that support our work on the spill everywhere problem, let us stress the dierence we want to make here between spilling and allocation. Spilling aims at nding which variable to evict from registers while allocation aims at nding which variable to keep in registers. Of course, one is the dual of the other, so conceptually spilling and allocation are the same. Now suppose you have a set of variables and you want to evict (spill) a minimum amount of them such that MaxLive is lowered by just one. As shown in [START_REF] Bouchez | On the complexity of spill everywhere under ssa form[END_REF] this problem is NP-complete even for the simplest SSA program instance. On the other-hand, consider you have already a set of allocated variables and you aim at allocating a maximum number of additional ones such that at every program point the register pressure Live is increased by at most 1.

Then as outlined before, this problem is, under SSA, polynomial with a complexity of O(Ωn). Ω being the maximum simultaneously live variables that remains to be allocated; n being the size of the program. Hence, in a way allocation is simpler than spilling. Our approach pushes this distinction further: Conceptually, every variable is initially in memory, and we evaluate the gain of allocating a given one instead of considering every variables to be initially in a virtually unbounded register le and evaluate the cost of evicting it. As we will see in this paper, this allows to be much more accurate concerning the modeling of gain/cost that accounts for ABI and register constraints.

Stepwise Allocation Is Close To Optimal

In a recent paper, Diouf et al. [START_REF] Diouf | Split register allocation: Linear complexity without the performance penalty[END_REF] studied the question to know whether or not the variables spilled on an optimal allocation with R registers are included in the set of variables spilled on an optimal allocation with R -1 registers (R > 0). Conceptually, this is equivalent to telling that the variables allocated on an optimal allocation with R -1 registers are included in the set of variables allocated when R registers are available. The answer to the question is no and this is illustrated in Figure 2. In this gure, we give the graph version of the example used by Diouf et al. [START_REF] Diouf | Split register allocation: Linear complexity without the performance penalty[END_REF]. The cost of each variable is represented by the number close to its corresponding node.

Dashed black circles correspond to spilled variables. Figure 2(a) depicts the optimal allocation performed when R = 1 and Figure 2(b) shows the optimal allocation performed when R = 2.

When R = 1, to perform the allocation of lowest cost, we need to spill b and d, which form the optimal spill set. When R = 2, we need to spill c, which is the optimal spill set. We clearly see that the optimal spill set when R = 2 is not included in the optimal spill set when R = 1. Even if, theoretically, the answer to the question of spill set inclusion is no, Diouf et al.

experimentally validated that when varying the number of registers from R min , the minimum number of registers to enable code generation, to the number of registers allowing to allocate all the variables, the inclusion property holds for 99.83% of the SPEC JVM98's methods. This experimental evaluation has been done with JikesRVM [START_REF] Alpern | The Jikes RVM project: Building an open source research community[END_REF], the research virtual machine of IBM.

This also proves, empirically, that the stepwise allocation is close to optimal.

3 Baseground

We now summarizes some denitions and results on graphs and chordal graphs upon which our approach is based.

In the rest of this paper, we assume that an estimated spill cost has been computed for each variable. A spill cost represents the access frequency of a variable, it is high when the variable is frequently accessed and low when it is not. We denote R the number of available registers.

Programs are usually represented as graphs, within graph coloring frameworks, and live sets within linear scan frameworks. Thus the spilling problem is naturally solved over these two representations. Our approach is compatible for both representations, but in the rest of this section we will focus on the graph representation.

Graphs and Weighted graphs

A graph G = (V, E) consists of two sets, V the set of vertices or nodes, and E the set of edges.

Every edge (v 1 , v 2 ) of E has two end points v 1 ∈ V and v 2 ∈ V . We say that v 1 and v 2 are adjacent(s) or are neighbor(s) if (v 1 , v 2 ) ∈ E.
The number of neighbors of a vertex v is called the degree of v. Here, We only consider undirected graphs, i.e., we do not make dierence between the edges (v 1 , v 2 ) and (v 2 , v 1 ). Figure 3(a) shows an arbitrary graph. In contrast to a clique, a stable set or an independent set is a subset S ⊆ V that does not contain two vertices that are adjacent.

A sequence of vertices [v 0 , v 1 , v 2 , . . . , v l , v 0 ] is called a cycle of length l + 1 if (v i-1 , v i ) ∈ E for i = 1, 2, l and (v l , v 0 ) ∈ E. A subset A ⊆ V is called a clique of G if every two distinct vertices of A are adjacent. A clique A is maximal if it is not properly contained in any other clique of G. A clique is maximum if there is no clique of G of larger cardinality. A vertex v of a graph G is simplicial if its neighbors RR n°8007
Assuming each vertex v of G = (V, E) is associated with a non-negative number w(v), the weight of a subset S ⊂ V is expressed as:

w(S) = v∈S w(v)
The graph G associated with the function w is called a weighted graph and denoted G w . The maximum weighted stable set is the stable set of maximal weight.

From a graph representation of a program, if we associate to each vertex a weight corresponding to its cost (it is assumed that a cost has been computed for each variable), the spilling problem becomes equivalent to the problem of choosing the set of vertices of minimal weight to remove from a weighted graph to make a coloring/assignment possible. This problem is much more complicated on arbitrary graphs, since the coloring problem is already NP-complete on these category of graphs. This coloring problem becomes easy on chordal graphs which are discussed below.

Chordal Graphs

The static single assignment (SSA) form is an intermediate representation with very interesting properties. A code is in SSA form when every scalar variable has only one textual denition in the program code. Most compilers use a particular SSA form, the strict SSA form, with the additional so-called dominance property: given a use of a variable, the denition occurs before any uses on any path going from the beginning of the program (the root) to a use. One of the useful properties of such a form is that the dominance graph is a tree and the live ranges of the variables (delimited by the denition and the uses of a variable) can be viewed as subtrees of this dominance tree. The intersection graph of these subtrees of the dominance tree represents the interference graph. An important result of graph theory states that the intersection graph of a family of subtrees of a tree is a chordal graph [START_REF] Charles | Algorithmic Graph Theory and Perfect Graphs[END_REF]. It follows that the interference graph of a program in SSA form is a chordal graph.

A graph G is chordal, triangulated or rigid-circuit if every cycle of length four or more has a chord, a chord being an edge joining two vertices of the cycle, that are not consecutive. The

Algorithm 1 MaximumWeightedStableSet

Require: σ: a perfect elimination order Require: w: a map associating to each vertex its weight Require: adj: a map associating to each vertex the list of its neighbors Var: w': a map associating to each vertex its current weight during computation Var: n: the number of vertices Var: marked_red: a (last in rst out) list keeping track of vertices marked red Var: marked_blue: a list keeping track of vertices marked blue 1: for i = 1 → n do 2:

v ← σ(i) An interesting property that we are going to use below in the paper is that it is easy to compute the maximum weighted stable set of a chordal graph [START_REF] Frank | Some polynomial algorithms for certain graphs and hypergraphs[END_REF] and this with a complexity of

O(|E| + |V |).
Before explaining the Frank's algorithm which computes the maximum stable set of a weighted graph, we need to explain the notion of perfect elimination order. An ordering v 1 , v 2 , . . . , v n of the vertices of a graph G is a perfect elimination order (PEO) if each v i is a simplicial vertex in G {vi,vi+1,...,vn} , the graph remaining from G when all the vertices preceding v i in the ordering have been removed. It has been proven in graph theory that a graph is chordal if and only if it has a perfect elimination order [START_REF] Charles | Algorithmic Graph Theory and Perfect Graphs[END_REF]. For instance [a, f, d, e, b, g, c] is a PEO of the chordal graph given in Figure 4.

Algorithm 1 computes the maximum weighted stable of a weighted graph G w . It receives σ a perfect elimination order of the weighted chordal graph, adj a map that associates to each vertex the list of its neighbors and w the weight function that associates to each vertex of G w its weight. Algorithm 1 goes through the list of vertices according to the order of σ. At the step i, it tests if the current weight of the vertex v that occupies the i -th position in σ is positive. If it is negative Algorithm 1 goes to the next step. Otherwise, v is marked red and each neighbor u of v has its current weight w (u) reduced by w (v). Any weight that becomes negative is altered to 0 and w (v) is set to 0. At the end of this process, the vertices marked red are visited in the reverse order of their insertion in the list. A vertex v is marked blue if it is not a neighbor of all the vertices previously marked blue. Finally, Algorithm 1 returns the set of vertices marked blue, that is stable set of maximum weight. The column iteration presents the iterations of the second for-loop of Algorithm 1 that modies the set of vertices marked red. The second column keeps track of the values of w . The vertices are ordered according to the perfect elimination order. The last column shows how the set of marked red evolves. The rst row shows, before the beginning of the loop, the values of w for each vertex and the set of vertices marked red which is empty. At the rst iteration, the weight of a is 1, thus a is marked red and the weights of its neighbors d and f are decreased by 1. The weight of a is then set to 0, which is underdrew. At the second iteration, the weight of f is 5. Thus, f is marked red and the weights of its neighbors a, d and e are decreased by 5. The weight of f is then set to 0, which is underdrew. Finally, we obtain the set of vertices marked red which are composed of b, f , a. blue. This is performed with the while-loop of Algorithm 1. At the rst iteration, the vertex b is chosen and is inserted in the set of vertices marked blue. At the second iteration the vertex f is chosen and inserted in the set of vertices marked blue. The vertex a is adjacent to f and cannot be added to the set of vertices marked red. Thus, a is removed from the list of vertices marked red. We then end up with a set of vertices marked blue composed of f and b of weight 8.

Layered-Optimal Register Allocation

We will focus here on the spilling problem for SSA programs. In the next section, we will give an extension of our approach which works on general graphs.

Based on the two observations explained on Section 2, we present here our solution which solves the spill minimization problem for R registers by layered optimal solutions to simpler problems on few registers. Each of this simpler problem is considered to have step, which is a small number lower or equal to R, available registers. Stepwise optimality does not guarantee an overall optimal allocation, but we will show that it comes very close to optimal, even with step = 1.

Algorithm 2 LayeredOptimalAllocation

Var: candidates: the list of vertices that are candidate to an allocation Var: allocated_list: the list of so far allocated variables 1: count ← 0 2: while candidates = ⊥ ∧ count < R do 3:

result ← OptimalAllocation(candidates)

4:

add every vertex of result to allocated_list 5:

remove every vertex of result from candidates 6:

count ← count + 1

7: end while 8: return allocated_list Algorithm 2 implements the layered-optimal heuristic. It takes as input candidates, the list of variables that are candidates to register allocation. It then returns as result allocated_list, the list of variables that have been allocated with R registers. Algorithm 2 calls OptimalAllocation which returns the optimal allocation set minimizing the spill cost among the variables that have not yet been allocated (currently in candidates). This set is added to allocated_list and removed from candidates. In its last step, Algorithm 2 nds the set of variables that minimizes the spill cost among the variables remaining in candidates.

RR n°8007

The function OptimalAllocation used in Algorithm 2 solves the allocation problem on a subset of the candidates variables when a unique register is available. On a chordal interference graph, this problem is equivalent to the problem of the maximum weighted stable set. Thus, the OptimalAllocation can be implemented with Algorithm 1. When assuming that step ≥ 2, OptimalAllocation can be implemented through dynamic programming [START_REF] Bouchez | On the complexity of spill everywhere under ssa form[END_REF].

In the following, we restrict ourselves to a step of one. The complexity of the layered-optimal

allocator is O(R(|V | + |E|)).
Algorithm 2 is a solid basis for an incremental allocation, but we have found two ways to improve it: biasing the cost of the variables, and iterating further on the set of allocated variables until we reach a point where we could not allocate more variables.

4.1

Biasing the weights Before we explain how we bias the costs/weights of variables/weights, let us rst have a look on Figure 6. We assume we have two registers, the step is set to one, and we are looking for the set of variables to allocate for the weighted graph given in Figure 5(a). We call this graph G w . When called on G w , Algorithm 2 will rst look for the maximum weighted stable set on the graph and consider the variables within this set as allocated. It will then look for the maximum weighted stable set on the graph remaining when the allocated variables are removed. G w has two maximum weighted stable sets with a weight of 8. The rst one shown in Figure 6(a) is composed of vertices b and f , in dashed lines. The second one composed of c and f , also in dashed lines, is shown in Figure 6(b). If b and f are chosen, at the next step, Algorithm 2 will look for the maximum weighted stable set on the graph remaining when b and f are removed from G w . This graph is represented by the black nodes and edges in Figure 6(a). The maximum weighted stable set for this graph is composed of d and g and has a cost of 6. This leads to spilling variables a, c and e with a spill cost of 4. In contrast, if we choose c and f , at the next step the maximum weighted stable of the graph, shown in black nodes and edges in Figure 6(b), will be composed of b and d with a cost of 7. This lead to a spill cost of 3.

This example shows that the choice among dierent maximum weighted stable sets has an impact on the next iterations of Algorithm 2. Arbitrarily Choosing a maximum weighted stable set can deteriorate the global register allocation. Our intuition to ameliorate the choice of the maximum weighted stable set is that: it is almost always better to choose the maximum weighted stable set that removes the most interferences in the graph on non-allocated variables.

Our approach to achieve this is to bias the cost/weight of variables/vertices with the number of neighbors of a vertex and we dene, to this purpose, the new weight function w as:

w (v) = w(v) × |V | + |adj(v)|
Inria where |V | is the number of vertices of the graph and adj(v) is the number of neighbors of the vertex v.

For two vertices u and v, the two following properties will always be veried with the new weight function: We introduce here the second improvement we want to perform on Algorithm 2. We assume we have two registers, the step is set to one, and we are looking for the set of variables to allocate for the weighted graph given in Figure 7(a). This graph is called G w and has four maximal cliques which are: {a, d, f {b, c, e}, {c, d, e}, and {d, e, f }.

if w(u) < w(v) then w (u) < w (v) if w(u) = w(v) then w (u) ≤ w (v) if adj(u) ≤ adj(v)

RR n°8007

11: end for

The example given in Figure 7 shows that we could miss interesting allocations when naively using the Algorithm 2. Algorithm 3 rst performs the layered allocation iterating at most R times (lines 1 to 7). It then calls Algorithm 4 which increments, for each freshly allocated vertex, the number of allocated vertices of each clique which contains it. If a clique has R of its vertices allocated, all the vertices of this clique are removed from candidates and this clique is removed from the list of allowed_cliques, which are the cliques that can have one of their non-allocated vertices allocated at next rounds. After Algorithm 4 nishes, Algorithm 3 calls Inria the function OptimalAllocation and Algorithm 4 iteratively until it reaches a x point, that is, an allocation that cannot be improved by subsequent calls of Algorithm 3.

Spilled variables

When a variable is spilled, it does not completely disappear from the interference graph. It is replaced by a set of short-lived variables which must be taken into account. For instance on RISC architectures, memory can only be accessed through load and store instructions. For example, before using a variable v spilled at address a, the value of v must be loaded from address a into a register. The extra instructions inserted to reload spilled variables form the spill code.

Some approacheslike the JikesRVM implementation of the linear scanspill locally an allocated variable when there is no free register to assign a reloaded variable. On CISC architectures like the x86, we also can take advantage of complex addressing modes to get operands directly from memory (at most one such operand on x86). On the other hand, graph coloring heuristics iteratively rebuild interferences after spill. Symmetrically, we can iteratively update the interferences after allocation.

Layered-Heuristic Allocator

Although the spill minimization problem is only pseudo-polynomial on SSA programs, the method also applies to general programs. The layered approach remains applicable, but the Frank's algorithm is not applicable as the graphs are not chordal.

The layered-heuristic algorithm clusters the nodes of an interference graph as stablesor independent setsusing a greedy heuristic. It clusters the variables according to their interference and spill costs, then allocates registers into layered cluster-based allocations.

We cluster variables that do not interfere and that corresponds to a stable set in the interference graph. Since we cannot compute the maximum weighted stable set in an arbitrary graph in a polynomial time, we approximate it. The performance of the allocator will depend on the quality of the approximated weighted stable set of maximum weight, which is called a cluster.

The clusters are computed incrementally, that is, we rst compute a cluster and then another cluster which does not contain any variable of the rst cluster, and so on, until we put all the variables into clusters. To approximate a cluster, we incrementally merge high-weights nodes which do not interfere with the variables already present in the stable set.

The clustering is performed by Algorithm 5 which transforms candidates, a list of vertices of a graph sorted by decreasing weight, into cluster_list, a list of clusters. It constructs a new cluster at each iteration of the outer while-loop. In order to compute cluster, the new cluster, all the variables still in candidates are added to potentials. The list potentials keeps tracks, at each round of the inner while-loop, the vertices that do not interfere with the vertices already into cluster. Every time a vertex v is added to cluster, all the neighbors of v are removed from potentials. At the end of the inner while-loop, the computed cluster is added to the cluster_list, and the next round of the outer while-loop starts. Finally, Algorithm 5 ends when every variable is in a cluster.

After the clusters have been computed, Algorithm 6 decides which clusters should be allocated to registers. The R clusters that maximize the sum of their weights are allocated.

The complexity of the layered-heuristic allocation is O(R × (|V | + |E|)). Indeed the algorithm iterates at most R times for R registers, the clustering step visits every neighbor of a node only once.

RR n°8007

Algorithm 5 ClusterVertices Require: candidates: a list of vertices, ordered by decreasing weight, that are candidate to an allocation Require: adj: a map associating to each vertex the list of its neighbors Var: cluster_list: a list of clusters 1: while candidates = ∅ do 2:

cluster ← ⊥ Our approach is very well suited to SSA programs, but we show that, it also yields excellent results on arbitrary interference graphs from non-SSA programs.

6.1 Chordal Graphs: SSA Programs

Methodology

We evaluated our approach on chordal interference graphs resulting from programs compiled with the Open64 compiler for the ST231 VLIW processor and for the ARM Cortex A8 (ARMv7). For the former, we generated the interference graphs for the SPEC CPU 2000int, the lao-kernels (an internal suite from STMicroelectronics) and the eembc benchmarks. We only used the lao-kernels for the ARMv7 processor.

For each of the considered benchmarks, we computed the spill costs based on the basic blocks's frequency and on the number of accesses to the variables within the basic blocks. We studied the impact of the register count, ranging from 1 to 32. For each instance of the register allocation problem and for each conguration, we compared the following algorithms:

GC The Chaitin-Briggs, optimistic graph coloring algorithm.

Inria

Optimal An optimal ILP-based allocator.

NL The layered allocation method implemented without the two improvements presented in Section 4.

FPL The layered allocation method with the xed point improvement that incrementally allocates variables until no variable could be allocated any more, with the given register count.

BL The layered allocation method that biases the spilling cost of nodes in order to choose between two stable sets of same (un-biased) cost the one that has more interferences.

BFPL The layered allocation method which biases the cost and iterates until a xed point in the allocation.

Results and discussion

The results obtained from the evaluation of chordal graphs generated from SPEC CPU 2000int, lao-kernels and eembc benchmarks are very similar. For the sake of exposition, we reported here the results for conguration with a register count of 1, 2, 4, 8, 16 and 32 registers. For all the conguration, BL, FPL, BFPL are close to optimal on average and are better than GC. On conguration with register counts up to 8, BL is RR n°8007 also quasi-optimal, but for conguration with 16 and 32 registers, we notice a performance degradation. This is reinforced by Figure 9 and on Figure 10 we also notice a performance degradation, when the register count is 32, of the FPL approach; it suggests that the biased improvement is very helpful on the lao-kernels benchmark suite, which is made of small benchmarks and thus can be more impacted by a bad allocation choice. to the optimal allocation for the specic benchmark. This gure depicts the distribution of these normalized allocation costs. GC, and to a lesser extent NL, show a high variability. This indicates that some benchmarks yield poor allocations for these allocators. On the contrary, BL, FPL and BFPL are consistently successful at computing close-to-optimal allocations. This is conrmed by Figures 12 and13 on the other benchmark suites. Notice a slight variability for FPL and registers on the lao-kernels targetting the ARMv7 (Figure 13).

Inria

Extension To Non-Chordal Graphs

We evaluate our approach on general, non-SSA progams, studying the SPEC JVM 98 benchmark suite (a benchmark set to measure the performance of Java virtual machines). We use the JikesRVM just-in-time compiler; its intermediate representation is not in SSA, and the interference graphs are not chordal in general.

We considered dierent congurations of register count going from 2 to 16. For each instance of the register allocation problem and for each conguration, we compared the following algorithms:

LS The original linear scan algorithm as implemented in JikesRVM.

RR n°8007

BLS A variant of the linear scan relying on Belady's furthest-rst strategy to make spilling decisions if their costs are close enough according to a chosen threshold.

GC The Chaitin-Briggs, optimistic graph coloring algorithm.

Optimal The globally optimal allocation implementing an ILP model proposed by Diouf et al. [START_REF] Diouf | Split register allocation: Linear complexity without the performance penalty[END_REF].

LH Our layered-heuristic method. Figure 14 shows the allocation costs for all SPEC JVM together, normalized over the cost of the optimal allocation's cost. Congurations with dierent register counts going from 2 to 16 registers. For almost all the register counts, the layered-heuristic allocator is close to optimal, except for the congurations with 14 and 16 registers. This can be explained by the accumulation of approximations in the incremental construction of maximal weighted stables, a consequence of the non-chordality of the interference graphs. Figure 15 reports for each individual benchmark the normalized allocation costs when we have a register count of 6 registers. We see that here again, the layered-heuristic allocator performs close to optimal allocations, and outperforms all the other allocation heuristics. For check, jess, javac, and jack, the overhead can reach 60% of the optimal, but the cost is still better than the conventional heuristics.

Inria 7 Related Work

Register allocation algorithms often rely on spilling algorithms to perform spill minimization.

In static compilation the dominant approach to register allocation is the graph coloring in which the spilling and coloring (assignment) algorithms are interleaved. During the simplify phase, whenever all the remaining nodes have at least R degrees, a node needs to be marked as spilled or pushed onto the stack (optimistic coloring) and removed from the graph. A natural intuition is to choose a node that has a low spilling cost and which interferes a lot. Many of the graph coloring variant are based on this intuition and use the quantity cost(v)/deg(v) to choose the variables to spill [START_REF] Chaitin | Register allocation & spilling via graph coloring[END_REF]. Thus, the spilling algorithm uses a global information over the whole program that combines the interference degree and the spilling cost.

In the context of just-in-time (JIT) compilation, register allocation time is part of the global execution time and (quasi-)linear complexity remains a driving force in the design of optimization algorithms. Moreover, when embedded systems are addressed, the limited memory resources is also an important issue. The linear scan which is one of the most used register allocation algorithm on JIT compilers has a worst case complexity of O(n × R), where n is the number of variables in the program and R is the number of available registers on the target architecture.

The original spilling heuristic used in linear scan [START_REF] Poletto | Linear scan register allocation[END_REF] is based on the Belady's furthest rst algorithm [START_REF] Belady | A study of replacement algorithms for virtual storage computers[END_REF]. This algorithm relies on local information to perform spilling: At a point p where registers are not enough to hold all the live variables, spill the variables whose live ranges go farther in the future. Recent versions of linear scan use more elaborate algorithms based on variables' spill cost estimation and sharing some of the global spilling decisions of graph coloring [START_REF] Sarkar | Extended linear scan: An alternate foundation for global register allocation[END_REF].

The idea of improving the spill minimization in a decoupled approach, where the allocation is decoupled from the assignment, has been explored by Proebsting and Fischer [START_REF] Proebsting | Probabilistic register allocation[END_REF], and by Braun and Hack [START_REF] Braun | Register spilling and live-range splitting for ssa-form programs[END_REF]. Braun and Hack generalized the Belady's furthest rst algorithm which works very well on straight-line code to control-ow graphs. Their approach, while being applicable as a pre-spill phase in any compiler, is more adapted to SSA-based register allocation.

They reported a reduction in the number of reload instructions by 54.5% compared to the linear scan and by 58.2% compared to the graph coloring. An other approach approach by Pereira and Palsberg rely on maximal cliques to drive spilling decisions with similar goals on chordal graphs, and generalizability to general graphs [START_REF] Magno | Register allocation via coloring of chordal graphs[END_REF]. Like the latter approaches, layered allocation is fast and can be used in a non-decoupled context for general programs, in a decoupled context for SSA programs, and as a pre-spill phase in any compiler. Unlike Braun and Hack, we experimentally show how that our layered-optimal algorithm performs close-to-optimal allocations.

Conclusion

Combining key observations in SSA-based, decoupled register allocation, we designed a new, polynomial approach to the spill-cost minimization problem: layered allocation. Our method contrasts with decades of work on register allocation by incrementally allocating clusters of variables to registers, while conventional heuristics incrementally spill variables. The criterion to form these clusters, rooted in the maximal clique problem (polynomial on chordal graphs), is also original. Our algorithm produces allocations that are very close to optimal on SSA programs, outperforming higher complexity heuristics such as the graph coloring methods. We also adapt our method to design an allocation heuristic for general, non-SSA programs.

These fundamental results pave the way to a simpler and very eective register allocation framework. Several steps remain to be taken to integrate it in a production compiler: studying the interactions with the register coalescing and other downstream optimizations, studying load/store RR n°8007 optimization variants (with transparent, ne-grain live range splitting), and reducing the number of incremental allocations to compete with the slightly faster linear scan allocators.

Figure 1 :

 1 Figure 1: Program based, versus graph based spilling approaches on a program with three registers.

Figure 2 :

 2 Figure 2: Counter example to spill set inclusion.

Figure 3 :

 3 Figure 3: (a) An arbitrary graph. (b) A weighted graph.

Figure 3 (

 3 Figure 3(b) shows a weighted graph based on the arbitrary graph presented in Figure 3(a). Each vertex has its weight close to it. For instance the vertex a has a weight of 2.

Figure 4 :

 4 Figure 4: An example of chordal graph.

Figure 5 :

 5 Figure 5: Looking for the maximum weighted stable set with Algorithm 1.

Figure 5 (

 5 Figure 5(b) and Figure 5(c) depicts the general steps of Algorithm 1 when applied to the

Figure 5 (

 5 Figure 5(c) explains how from the vertices marked red we compute the set of vertices marked

Figure 6 :

 6 Figure 6: Example showing the benet of biasing the weight.

Figure 7 :

 7 Figure 7: Example showing the benet of iterating until a xed point.

Figure 7 (

 7 Figure 7(c) without adding a clique of size 3 and thus making a coloring with 2 colors impossible. Unlike f , the vertices c and e are not contained in a maximal clique that has 2 vertices already allocated. It follows that either c or e can be added to the graph of allocated vertices shown in Figure 7(c) without making a coloring, with two colors, of the graph impossible since the size of the maximum clique of the resulting graph, shown in Figure 7(d), will be 2.

Figure 9 :

 9 Figure 9: Allocation cost for the EEMBC benchmark suite on ST231.

Figure 8

 8 Figure 8 presents the average of the allocation cost of all the application of the SPEC CPU 2000int. For the sake of exposition, we reported here the results for conguration with a register

Figure 10 :

 10 Figure 10: Allocation cost for the LAO-KERNELS benchmark suite on ARMv7.

Figure 11 :

 11 Figure 11: Distribution of the allocation costs over individual programs of the SPEC CPU 2000int benchmark suite on ST231.

Figure 11 studies

 11 Figure 11 studies how the allocation results vary across individual interference graph for all the benchmark programs in the SPEC CPU 2000int suite. Each allocation result is normalized

Figure 13 :

 13 Figure 13: Distribution of the allocation costs over individual programs of the lao-kernels benchmark suite on ARMv7.

Figure 14 :

 14 Figure 14: Layered-heuristic allocator compared to other algorithms for dierent register counts.

Figure 15 :

 15 Figure 15: Layered-heuristic compared to other allocators when the register count is 6.

  : candidates: the list of vertices that are candidate to an allocation Var: allocated_list: the list of so far allocated variables Var: allowed_cliques: the list of cliques that do not have more than R allocated vertices, it is

	Algorithm 3 FixedPointLayered (FPL)
	initialized to the list of maximal cliques Var: allocated_Per_clique: a map associating to each maximal clique the number of its allocated
	vertices 1: count ← 0 2: while candidates = ⊥ ∧ count < R do 3: result ← OptimalAllocation(candidates) 4: add every vertex of result to allocated_list 5: remove every vertex of result from candidates 6: count ← count + 1
	7: end while 8: Update(candidates, allocated_list, allowed_cliques, allocated_Per_clique) 9: while candidates = ⊥ do 10: result ← OptimalAllocation(candidates)) 11: remove every vertex of result from candidates 12: Update(candidates, result, allowed_cliques, allocated_Per_clique)
	13: end while 14: return allocated_list
	Algorithm 4 UPDATE
	vertices 1: for all v ∈ allocated_list do 2: for all clique ∈ allowed_cliques do 3: if v ∈ clique then 4: increment allocated_Per_clique(clique) by one 5: if allocated_Per_clique(clique) ≥ R then 6: remove all the vertices of clique from candidates 7: mark clique as to be removed at the end of the loop
	8:	end if
	9:	end if
	10:	end for

Require

Require: candidates: the list of vertices that are candidate to an allocation Require: allocated_list: a list of allocated variables Var: allowed_cliques: the list of cliques that do not have more than R allocated vertices Var: allocated_Per_clique: a map associating to each maximal clique the number of its allocated

RR n°8007RESEARCH CENTRE PARIS -ROCQUENCOURT Domaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


John Douglas

Philippe Gueguen

Emmanuel Chaljub

Fabrice Cotton

Peter Suhadolc

Giovanni Costa

Donat Faeh

Eva Spühler

Andrej Gosar

Enrico Priolo

Carla Barnaba

Roberto Paolucci

Carlo Cauzzi

Claudio Eva

DISSEMINATION OF ALPINE ACCELEROMETRIC DATA

.

Introduction

Due to tourism and industrial development, urbanized areas located in Alpine valleys have undergone massive development during recent years. These infrastructures are located on deep-filled sediment (post-glacial alluvium) alpine valleys, close to faults where moderate to large earthquakes have occurred in the past. Recent studies have shown that seismic motion is strongly amplified by local soil conditions increasing seismic hazard. Improving seismic hazard evaluation requires better instrumentation to record earthquakes for a better knowledge of the hard rock site response, the non-linear site response of alpine alluviums deposits and the complex 2D-3D site effects in Alpine valleys. Moreover, a share of the available information and common database is useful for elaboration of alpine geotechnical profiles which could be used in calculations and the 'alpine records' can be used to define specific seismic motions.

In the framework of the SISMOVALP project , the compilation of the available accelerometric data was planned, as well as the collection of new data from temporarily seismic experiments planed in selected alpine valleys. This data collection, in the form of a CD ROM has just been released. This short article summarises how the data was selected, processed, the contents of the CD ROM and also introduces the browser software provided alongside the data.

Data selection

The purpose of this new CD ROM is to present earthquake engineers and seismologists with a large set of accelerometric data for which the parameters (source, path and site) required for most analyses are complete. This CD ROM only contains data from the Alpine region (France, Italy, Slovenia and Switzerland) collected by permanent and/or national networks or by temporarily networks deployed for the SISMOVALP project. Additional European data from other regions can be downloaded from ISESD and imported into SMDN for analysis, if required.

The same strict selection criteria imposed by [START_REF] Ambraseys | Dissemination of European Strong-Motion Data[END_REF] for selecting records for the CD ROM of strong-motion data from the European and Middle East was not applied here since data from the Alpine region is sparse and their selection criteria would significantly reduce the quantity of data on the CD ROM.

The following criteria have been adopted to select records for this CD ROM:

1. M L ≥4.0 if record from RAF. 2. M L ≥3.0 if record from RAP.

3. Uncorrected acceleration time-history is available. 4. Natural frequency and damping are available if analogue instrument. 5. All three (two horizontal and one vertical) components are available. 6. Acceleration time-histories are of high-quality.

Figure 1 shows the distribution of records in the selected set in terms of M w, if available, or local magnitude (M L ), if M w not available, and distance to surface projection of rupture, if available, or epicentral distance otherwise. In total, the CD ROM contains 434 triaxial strong-motion records from 140 earthquakes and 77 stations in Alpine valleys in Europe.

Correction technique

In order to provide corrected acceleration, velocity and displacement time-histories that are reliable a correction procedure was adopted that used individual filter cut-off frequencies for each record. The correction technique implemented in the Basic Strong-Motion Accelerogram Processing (BAP) software [START_REF] Converse | BAP Basic Strong-Motion Accelerogram Processing Software[END_REF] was used for the correction of all of time-histories contained on this CD ROM. This method consists of: interpolation to 200 samples a second; a correction for the instrument response, when the required characteristics are known, and high-cut filtering, with a cosine transition from the roll-off frequency to the cut-off frequency; followed by low-cut bidirectional Butterworth filtering of the acceleration after padding the time-history with zeros. The main problem with filtering strong-motion records is the selection of appropriate cut-off frequencies for the high-cut and, particularly, low-cut filters. For this CD ROM a method based on the estimated signal-to-noise ratio of each record was used. Time-histories from digital instruments often include long-enough pre-event portions to use as an estimate of the noise. Therefore for those records with pre-event portions the Fourier amplitude spectrum of this noise estimate was computed as was the Fourier amplitude spectrum of the rest of the signal and the ratio was calculated. The low cut-off frequency was chosen from this spectrum as the frequency at which the signal-to-noise ratio becomes less than two. The instrument-corrected and filtered displacement timehistory was then plotted and the cut-off frequencies altered slightly if the displacement trace did not look realistic. The cut-off frequencies often did not need changing. The vast majority of time-histories on this CD ROM are from digital instruments.

It is more difficult to choose cut-offs for records from analogue instruments because often estimates of the noise do not exist. Most films from analogue instruments feature a fixed trace, which record the position of light beams reflected from mirrors attached to the instrument case [START_REF] Hudson | Reading and interpreting strong motion accelerograms[END_REF]. Since these traces do not record the ground motions, if they were digitized in the same way as the ground motion traces they would provide the best estimate of the recording and digitization noise. Unfortunately, however, they are not often digitized or disseminated. For those time-histories that have such a digitized fixed trace they were used to select the low cut-off frequency in the same way as was done with the records with pre-event portions.

For those records from analogue instruments with no available fixed trace, the upper edge of the average hand digitization noise estimate from Figure 8 of Trifunac & Todorovska (2001) was used as an estimate of the noise. This estimated noise spectrum may be too high for many records that were digitized on automatic or semi-automatic digitizing tables but because the digitization procedure is unknown for most records it was thought to be a conservative estimate. This signal-to-noise ratio spectrum and the signal's Fourier amplitude spectrum were used to estimate the low cut-off frequency. The Fourier amplitude spectrum was examined following the suggestion of [START_REF] Zarè | Strong motion dataset of Turkey: Data processing and site classification[END_REF] that Fourier amplitudes that do not tend to zero at low-and high-frequencies are evidence for noise. The estimated cut-off frequency were often varied if it was found that the displacement traces were not realistic or if it was found that a less strict cut-off frequency could be used and still obtain a realistic displacement trace. Even with the techniques adopted here for the selection of cut-off frequencies for records from analogue instruments with no digitized fixed traces there is still some subjectivity in the process.

After choosing the appropriate cut-off frequency for each component (two horizontal and one vertical) of a strong-motion record a single cut-off frequency was chosen for all three components for consistency. The choice of this single cut-off was difficult for some records because the generally lower amplitude vertical ground motions imply that often the appropriate low-frequency cut-off for vertical components is higher than for the corresponding horizontal components. The high-frequency filtering was accomplished using the commonly-chosen roll-off frequency of 23 Hz and a cut-off of 25 Hz for records from analogue instruments and a roll-off of 50 Hz and a cut-off of 100 Hz for records from digital instruments (e.g. [START_REF] Converse | BAP Basic Strong-Motion Accelerogram Processing Software[END_REF]. Since most digital instruments have natural frequencies of about 50 Hz, and some of those with lower natural frequencies correct for the instrument response automatically, the effect of instrument correction is not large and therefore the requirement to apply a high-cut filter is less than for records from analogue instruments.

This procedure is the same as followed by [START_REF] Ambraseys | Dissemination of European Strong-Motion Data[END_REF].

Strong-Motion Datascape Navigator version 1.1

For their strong-motion data dissemination project [START_REF] Ambraseys | Dissemination of European Strong-Motion Data[END_REF] developed a new piece of software called Strong-Motion Datascape Navigator (SMDN) to select, view and analyze the data contained on the CD ROM produced. In addition, it allowed the user to import and analyze data from the Internet Site for European Strong-motion Data (ISESD, http://www.isesd.cv.ic.ac.uk). The program was written using Matlab 6.5.1 from Mathworks and then compiled (using Matlab Compiler 3.0.1) to create a stand-along program that can be run on PCs with or without Matlab installed.

To create the CD ROM of the SISMOVALP project it was decided to use the same browser software although slightly modified to allow the user to access extra information on the site conditions at strong-motion stations within Alpine valleys. To complete the information on ground motion in Alpine valleys, the browser also allows the user to view a list of references (with an electronic copy, when available) dealing with the geotechnical description of the valley and the seismic response. The new software is called Strong-Motion Datascape Navigator version 1.1.

Figure 2 shows a screenshot of the main window in SMDN v1.1. This window lists all the currently selected records grouped by earthquake in ascending date order. Under a header line for each earthquake giving the main parameters for the earthquake (name, date, time, location, magnitude and mechanism), the associated waveforms for this earthquake are listed along with their main parameters (station, source-to-site distances, local site category and a selection of strong-motion parameters). Earthquakes and records can be selected by clicking on their entries in this list. SMDN calculates a large number of strong-motion parameters, many of which can be used as search criteria (see below). Also it allows the user to plot the following graphs that are useful for strong-motion analysis:

1. acceleration, velocity and displacement time-histories; 2. linear elastic response spectra; 3. normalized linear elastic response spectra; 4. linear elastic response spectra using the tripartite representation;

5. pseudo-spectral acceleration versus spectral displacement; 6. inelastic ductility demand spectra; 7. inelastic strength demand and inelastic spectral velocity spectra; 8. strength demand versus ductility demand for a given natural period and damping; 9. Husid plot, i.e. Arias intensity versus time; 10. normalized energy density versus time; 11. absolute Fourier amplitude spectra (these can be smoothed); 12. horizontal to vertical ratio of the Fourier amplitude spectra (including user-defined smoothing).

All of these graphs (except number 1) can be plotted for a number of records simultaneously, thereby aiding comparison of different records. Graphs displaying linear elastic response spectra can be plotted to also display a chosen design spectra.

Maps showing the locations of the strong-motion stations that recorded accelerograms during a particular earthquake can be drawn using SMDN v1.1. These maps also show the epicentre and if known the location of the causative fault. On these maps the recorded acceleration, velocity and displacement time-histories can be plotted to aid understanding of the spatial variability of the strong ground motions.

Like in the original SMDN, if there is a photograph of a recording station and/or a shear-wave velocity profile then these are displayed when the station information is viewed. In SMDN additional information on site conditions is also displayed, if available. This information includes the horizontal-to-vertical spectrum from microtremor measurements and the estimated natural frequency of the valley. Figure 3 shows an example of the window showing station information. Like in the original SMDN, SMDN v1.1 allows the user to import data from ISESD and use the data in a similar way. Since ISESD does not provide all the associated parameters as provided for the inbuilt data, nor does it provide constant strength or constant ductility inelastic spectra some features are disabled when SMDN is used with ISESD data.

Searching

In previous CD ROM strong-motion data collections and on Internet sites providing strong-motion data (e.g. ISESD and COSMOS, http://db.cosmos-eq.org/) the search criteria that can be used are limited to geophysical parameters such as earthquake magnitude, source-to-site distance and local site category and possibly a handful of simple strong-motion parameters, such as peak ground acceleration. Such criteria, although useful, do not allow the user to find records that match criteria defined in terms of strong-motion parameters, which are often the only parameters researchers or designers know. For example, if a user wants a selection of time-histories to use to test a structure with a given natural period and damping subjected to shaking specified by a seismic design code then the spectral ordinate for that given period and damping is known and not the magnitude and distance of the target event. Therefore, SMDN allows the geophysical and strong-motion parameters specified in Table I to be used as search criteria. The thresholds and limits, for example those used for the computation of the absolute bracketed duration or the spectral intensity, used for the different strong-motion parameters can be altered by the user.

Brief analysis of strong-motion data from the Vallorcine earthquake

On the 8 th September 2005, as this CD ROM was being finalized, an earthquake of M w 4.5 (SED) occurred close to Vallorcine (Haute Savoie) in France. This earthquake had an almost pure strike-slip mechanism and a centroid depth of 8km (SED). According to the macroseismic intensity map reported by the Bureau Central Sismologique Français this earthquake was associated with an epicentral EMS-98 intensity of V. This earthquake clearly shows the importance and benefit of collection and dissemination of data from a variety of networks in the Alpine region since it occurred close to the France-Italy-Switzerland borders and valuable data was recorded by at least three accelerometric networks in these three countries. This section briefly discusses the available data from this earthquake.

Comparison between observed and predicted peak ground acceleration

Figure 4 shows a comparison between the peak ground accelerations recorded during the Vallorcine earthquake at all stations that triggered (and not just those included on the SISMOVALP CD ROM) and the values predicted by the empirical ground motion estimation equation of [START_REF] Ambraseys | Equations for the estimation of strong ground motions from shallow crustal earthquakes using data from Europe and the Middle East : Horizontal peak ground acceleration and spectral acceleration[END_REF]. It can be seen that, on average, PGA is well predicted by the model of [START_REF] Ambraseys | Equations for the estimation of strong ground motions from shallow crustal earthquakes using data from Europe and the Middle East : Horizontal peak ground acceleration and spectral acceleration[END_REF] up to a distance of about 150km even though [START_REF] Ambraseys | Equations for the estimation of strong ground motions from shallow crustal earthquakes using data from Europe and the Middle East : Horizontal peak ground acceleration and spectral acceleration[END_REF] used only data from earthquakes larger than Mw 5.0 (the magnitude of the Vallorcine earthquake is 4.5) and also all data used by [START_REF] Ambraseys | Equations for the estimation of strong ground motions from shallow crustal earthquakes using data from Europe and the Middle East : Horizontal peak ground acceleration and spectral acceleration[END_REF] came from more seismically active parts of Europe and the Middle East than the western Alps. For distances greater than about 200km the observed ground motions are of lower amplitude than predicted. This could be due to anelastic attenuation, which is not explicitly included in the model of [START_REF] Ambraseys | Equations for the estimation of strong ground motions from shallow crustal earthquakes using data from Europe and the Middle East : Horizontal peak ground acceleration and spectral acceleration[END_REF]. Also note, as expected, ground motions at sites in different countries at similar distances are of comparable amplitudes.

Figure 1 .

 1 Figure 1. Distribution of selected records in terms of moment magnitude (Mw), if available, or local magnitude (ML), if Mw not available, and distance to surface projection of rupture (df ), if available, or epicentral distance (de) otherwise.

Figure 2 .

 2 Figure 2. Screenshot of main screen in SMDN v1.1 showing main characteristics of available records.

Figure 3 .

 3 Figure 3. Screenshot of the window showing station information.

Figure 4 .

 4 Figure 4. Comparison between the peak ground accelerations recorded during the Vallorcine earthquake (Mw 4.5, strike-slip mechanism) and the values predicted by the ground-motion model of Ambraseys et al. (2005).

Table I .

 I Allowable search criteria in SMDN

	Geophysical Source	Path	Site
		Name	Epicentral distance	Name
		Country	Hypocentral distance	Country
		Date	Fault distance	Location (Lat. &
		Location (Lat. & Long.)	Rupture distance	Long.)
		Focal depth		Site class
		Magnitude (M w , M s , m b		Vs,30
		& M L )
		Mechanism
	Waveform	Amplitude	Duration	Frequency
		PGA	Absolute bracketed	Spectral intensity
		Effective PGA	Relative bracketed	Spectral shape
		PGV	Absolute uniform
		PGD	Relative uniform
		Arias intensity	Absolute significant
		Spectral intensity	Relative significant

Acknowledgements

The CD ROM was created with the support of the Alpine Space Program of the European Community. It is a deliverable of the project: 'Seismic hazard and alpine valley response analysis -SISMOVALP' (F/1-2/3.3/25) supported by the CIP INTERREG III B Alpine Space Program. For details of this project see http://www-lgit.obs.ujfgrenoble.fr/sismovalp. We would like to thank the organizations and agencies members of this project, which generously contributed strong-motion records to the SISMOVALP project: Observatoire de Grenoble, Grenoble, France; Université Joseph Fourier, LGIT, France; Centro Ricerche Casaccia, Ente per le Nuove Tecnologie, l'Energia e l'Ambiente (ENEA), Rome, Italy; Dipartimento di Scienze della Terra, Universita degli Studi di Trieste, Trieste, Italy; Dipartimento per lo studio del territorio e delle sue risorse, Universita degli Studi di Genova, Italy; Instituto Nazionale di Oceanografia e Geofisica Sperimental, Trieste, Italy; Politecnico di Milano, Italy; Servizio Sismico Nazionale (SSN), Ufficio Sistemi di Monitoraggio, Rome, Italy; Societa Gestione Impianti Nucleari (SOGIN), Rome, Italy; Environmental Agency of the Republic of Slovenia, Slovenia; Swiss Seismological Service, Swiss Federal Institute of Technology, Zurich, Switzerland.

The partners of this project are interested in using the browser written for this project for the dissemination of other strong-motion data. Because the data and program are kept separate the use of the browser for the analysis of other strong-motion data is straightforward. Therefore any organisations that are keen to use this software should contact John Douglas (j.douglas@brgm.fr).


N Mellado

P Reuter

C Schlick

Semi-automatic geometry-driven reassembly of fractured archeological objects

Keywords: (according to ACM CCS): I, 3, 5 [Computer Graphics]: Computational Geometry and Object Modeling-Geometric algorithms I, 3, 6 [Computer Graphics]: Methodology and Techniques -Interaction techniques

3D laser scanning of broken cultural heritage content is becoming increasingly popular, resulting in large collections of detailed fractured archeological 3D objects that have to be reassembled virtually. In this paper, we present a new semi-automatic reassembly approach for pairwise matching of the fragments, that makes it possible to take into account both the archeologist's expertise, as well as the power of automatic geometry-driven matching algorithms. Our semi-automatic reassembly approach is based on a real-time interaction loop: an expert user steadily specifies approximate initial relative positions and orientations between two fragments by means of a bimanual tangible user interface. These initial poses are continuously corrected and validated in real-time by an algorithm based on the Iterative Closest Point (ICP): the potential contact surface of the two fragments is identified by efficiently pruning insignificant areas of a pair of two bounding sphere hierarchies, that is combined with a k-d tree for closest vertex queries. The locally optimal relative pose for the best match is robustly estimated by taking into account the distance of the closest vertices as well as their normals. We provide feedback to the user by a visual representation of the locally optimal best match and its associated error. Our first results on a concrete dataset show that our system is capable of assisting an expert user in real-time during the pairwise matching of downsampled 3D fragments.

Motivation

Archeological objects are often broken and fractured into a large amount of fragments, and the real-world reassembly of the fragments is sometimes impossible due to their size, weight, fragility, or inaccessibility. With the increasing availability of 3D laser scanners, the scanning of fragments is becoming more and more popular, resulting in large collections of detailed fractured archeological 3D objects (see for example [TFK * 09,Lev00]). Fragments can then be reassembled virtually, and gained information used as a blueprint to reconstruct the real-world object.

Let us first recall some previous work about virtual reassembly. On the one hand, there is a variety of techniques to fully automatically reassemble fractured objects. Most of these techniques are based on a pairwise matching of ge- tunately, purely automatic reassembly methods are highly data-dependent and generally fail when fragments are missing or strongly deteriorated.

On the other hand, fragments can be reassembled manually on the computer by archeologists, through the help of an efficient computer-human interface, thus taking into account human expertise and knowledge. To overcome the difficult problem of positioning and orienting 3D models relative to each other, a bimanual tangible user interface for the efficient reassembly of fractured archeological objects has been presented before. It allows the user to reassemble two fragments as if they were in his hands [RRC * 07]. Unfortunately, purely manual reassembly lacks precision, and does not benefit from a computer assistance, although the geometry is available.

Our contribution consists in a semi-automatic pairwise fragment matching technique that takes into account both the archeologist's expertise, as well as the power of automatic geometry-driven matching algorithms. Previous work on semi-automatic reassembly has not taken into account the user-specified relative positions and orientations during the interaction. For example, it either requires the user to specify beforehand corresponding landmarks [RJ04] or annotations [KTNL05] for guiding the matching, or afterwards to make the final choice between multiple assembly propositions [PKT01, TFK * 09].

In contrast to these methods, our new semi-automatic reassembly approach integrates the user in a real-time interaction loop: an expert user specifies approximate initial relative positions and orientations between two fragments that are used by the system to "snap" them to the locally best match. Once aligned, the two fragments are virtually linked and considered as a single new fragment. The complete reassembly of a complex broken object can thus be done by successive pairwise matchings.

The contributions of this paper are threefold:

• A semi-automatic reassembly approach that integrates the user in the real-time interaction loop and steadily finds the best local match according to the user-specified initial pose,

• a bounding sphere hierarchy combined with a k-d tree for the identification of the potential contact surface between the two fragments to match as well as for efficient closest point queries, and

• an ICP-based matching method that takes into account both the influence of the positions and normals of the contact surface by using a bi-factorial weighting function, and that is robust to slight changes of the initial pose.

This paper is organized as follows. In section 2, we present the idea of an interaction loop that ensures the integration of expert user knowledge. In section 3, we show how we ensure the robust geometric matching to operate in realtime. In section 4, we present qualitative and quantitative results, before we conclude with directions to future work in section 5.

Semi-automatic interaction loop

The central idea of this paper is to increase the efficiency of the reassembly process by integrating real-time geometrydriven matching algorithms into the user interaction loop, as can be seen in figure 2.

User

Initial pose

Locally optimal pose Geometric matching

Interaction technique Feedback

Figure 2: Semi-automatic interaction loop.

There are three principal prerequisites for this semiautomatic interaction loop: first, we need an efficient interaction technique for the specification of approximate initial relative positions and orientations between two fragments. Second, we need a real-time geometric matching algorithm that finds the locally optimal match with respect to this relative initial pose. Third, we have to provide a feedback that makes it easy for a user to validate the match, or refine the input. In the following, we show the choices that we have made for each of these prerequisites.

Interaction technique

The role of an expert in pairwise semi-automatic reassembly is to specify relative positions and orientations of the fragments that are used as input for the matching algorithm. In order to keep concentrated on the actual archeological previous work dealing with bimanual tangible user interfaces [RRC * 07]: as can be seen in figure 3, in each hand, the user manipulates an electromagnetically tracked prop, and the translations and rotations are directly mapped to the corresponding virtual fragments on the display. For each hand, a corresponding foot pedal is used to activate props tracking (or "clutch").

Geometric matching

We rely on an algorithm based on the iterative corresponding point (ICP) [START_REF] Besl | A method for registration of 3-d shapes[END_REF] for real-time matching, as it was already made in previous work on automatic virtual reassembly [HFG * 06, BTFN * 08]. Whereas the ICP is mostly used for registration to transform the surface sheets from different scans of a single object in a common coordinate frame, we use it for the alignment of contact surfaces for scans coming from two different objects. Consequently, the ICP has to be modified in order to reject all the surface parts that do not belong to the contact surface. As the name suggests, the ICP algorithm calculates the registration iteratively, and it directly works on the scanned 3D points that define the surface (in the remainder of this paper, we will refer to the scanned points as vertices). In each iteration, the algorithm selects the best corresponding vertices in order to calculate the locally optimal rotation and translation for surface alignment, by minimizing an objective function. In section 3, we explain what we understand by "best corresponding vertices" for the alignment of the contact surface, and how we select and weight them in real-time.

Feedback

Once the locally optimal match with respect to the userspecified initial pose has been determined, we have to provide feedback to the user about the locally best match. We use a visual transparent 3D representation of the matching result, as can be seen in figure 1. Based on this visual feedback, the user can evaluate if the matching result is geometrically plausible and coherent with his intent and validate the proposition or specify new initial positions and orientations. As can be seen in figure 1, we also provide a graphical indicator of the local error of the ICP algorithm, i.e. the root mean square (RMS) error of the matching of the contact surface.

Note that with this visual feedback, the reassembly "snaps" to the locally best corresponding match, like a magnet that sticks the fragments together, as long as the provided pose is within the distance threshold, and as long as it is not closer to a different local minimum. A similar "snap" metaphor has proven to be efficient in 2D vector graphics applications, where imaginary grid lines at a coarse spacing help to precisely align 2D objects despite a roughly aligned input.

Optimizing the geometric matching

3D models of fragments are acquired with 3D scanning devices; they are thus often noisy and consist of millions of vertices. However, the interaction loop involved in our semi-automatic matching algorithm (figure 2) requires our ICP-based matching algorithm to operate in real-time, and to robustly align the two surface sheets even in the presence of noise. According to Rusinkiewicz's fast ICP variants [RL01], an iteration is composed of six steps that are likely to be optimized: data selection, pairwise vertex matching, weighting pairs, rejecting pairs, computing an error, and minimizing the error. In the following, we show how we optimize the first four steps for efficient reassembly. More precisely, we will focus on the one hand on the vertex pair selection in order to reduce the amount of data and to fulfill the real-time constraint, and on the other hand on the weighting and rejection of the pairs for efficient and robust matching of noisy data.

Efficient vertex pair selection

For an efficient alignment of the two surface sheets, in a preprocess, we first downsample the vertices of the fragments to some tens of thousands of vertices [START_REF] Turk | Zippered polygon meshes from range images[END_REF]RL01]. For a further optimization, we only need to select the vertices that are present on the potential contact surface, while rejecting all the others. Note that in our semi-automatic reassembly method, the user specifies the relative initial positions and orientations of the two fragments so that the potential contact surface is already roughly aligned. As a consequence, we can consider that the vertices of both contact surface sheets are close to each other: their distance is less than a distance threshold d T that can be adjusted by the user.

For the efficient detection of contact surfaces and closest vertex queries, we construct a k-d tree for each of the two fragments in a preprocess. Furthermore, inspired by classical collision detection algorithms [Lin93], we mix the k-d tree with a bounding sphere hierarchy: for every node of the k-d tree, we determine a sphere that encloses all vertices of the descendant nodes.

Consequently, in the real-time interaction loop, for every user-specified initial pose, we can efficiently determine all vertices of the contact surface by intersecting the two hierarchies recursively (see the bold spheres bs and b ′ s in figure 4). This pruning is slightly different compared to collision detection since we do not only want to obtain the sphere intersections, but also all vertices where the distance is less than the distance threshold d T . Our solution is conservative: we increase the sphere radii by dT 2 to capture all necessary vertices, but also vertices at a distance greater than d T (see the dotted spheres in figure 4). To obtain the closest vertex pairs, for every vertex of the potential contact surface of the first fragment, we efficiently determine the closest vertex of the second fragment by using the k-d tree.

N. Mellado & P. Reuter & C. Schlick / Semi-automatic geometry-driven reassembly

× × A × B × × bs ×C × × × × × bs ′ d BC d T 2 d T d AC × × × × × × × × × × × × × × × × × × × × × × bs bs ′ bs ′ × × × × × bs Figure 4: (left)
The two bounding sphere hierarchies. (right) Determining the contact surface by conservative pruning of the two hierarchies: even though d BC < d T , the sphere intersection does not detect it. By increasing the sphere radii by dT 2 , the intersection is successfully determined, but too distant vertices are detected as well: d AC > d T .

Vertex pair weighting and rejection

Recall that once all the closest vertex pairs are conservatively determined, we have to reject all the pairs with a distance d greater than d T in order to respect the contact surface. Moreover, we have to determine the influence of each pair. Following Pulli [Pul99], we take into account normal coherence in addition to vertex distance: pairs with normals of opposite direction should have a greater influence (marked in green in figure 5) compared to other neighboring pairs of the contact surface (marked in yellow in figure 5). Pairs where the dot product of the normals η is greater than a threshold η T should not be taken into account at all. The result is a better convergence of the ICP algorithm. Due to the interaction loop in our semi-automatic reassembly, we must ensure that the locally optimal match evolves coherently over time. Indeed, when the user slightly changes the relative pose of the two fragments, the matching should only vary gradually in order to prevent from popping artifacts when some vertex pairs oscillate around the distance threshold d T or the normal coherence threshold η T . Therefore, we propose to smooth both the distance of the vertex pairs and their associated normals by the decay function f (x) = (x 2 -1) 2 . This results in the weighting functions f d for the vertex distance d (equation 1) and fn for the normal coherence (equation 2). A zero weight corresponds to a pair rejection.

   f d (d) = d d T 2 -1 2 if 0 ≤ d ≤ d T f d (d) = 0 otherwise (1)    fn(η) = 1+η 1-η T 2 -1 2 if -1 ≤ η ≤ -η T fn(η) = 0 otherwise (2)
We combine both weighting criteria in the following bifactorial weighting equation, where α ∈ [0; 1] is a parameter that can be tuned by the user.

B f (d, η, α) = e - α f d (d) * e - 1 -α fn (η) e -1
(3)

Results

We now show results for each of the choices we have made for the interaction technique, the matching algorithm and the feedback provided to the user. All the results are based on an informal user study, as well as on our own experiments on a variety of fractured objects that we have tried to reassemble pairwisely.

Interaction technique

Concerning the bimanual tangible user interface [RRC * 07], we have observed a good acceptance and an efficient manipulation after a short learning process. However, the lack of haptic feedback and the absence of a proper collision detection is sometimes irritating during fragment manipulation. Moreover, the ICP algorithm has difficulties to converge when objects interpenetrate at the initial pose due to the high number of closest vertices that are not at the contact surface, as explained in the following section.

Geometric matching

We analyzed the results of our ICP-based geometric matching technique quantitatively and qualitatively.

Concerning the qualitative results, we first analyzed the c The Eurographics Association 2010. influence of the distance threshold d T , the normal threshold η T , and the parameter α that adjusts the influence of the positions and normals of the vertices of the contact surface (equation 3).

When choosing a too small distance threshold d T or normal threshold η T , the ICP algorithm does not have enough input points to converge to the same local optimum over time (i.e. with slight changes of the initial pose). On the contrary, when choosing too large thresholds, there are too many selected vertex pairs that do not belong to the contact surface, resulting in an inefficient and wrong matching. As a general rule, the distance threshold d T has to be defined specifically for each fragment pair, since for the alignment of the two subparts of the fragments we cannot know in advance the size of the contact surface with respect to the size of the entire objects. Inversely, the normal threshold η T is somehow easier to adjust: in our experiments, we have generally used η T = 0.5 corresponding to a normal cone with a radius of 45 degrees. This threshold can further be reduced when there is less noise in the datasets.

As said above, the parameter α adjusts the influence of the positions and normals of the vertices of the contact surface. Whereas using only the distances of the vertex pairs (α = 1) results in an inefficient matching by taking into account the interpenetrating parts of the surfaces, using only the normals (α = 0) does not filter out surface parts that are too far away and do not spatially belong to the contact surface. As expected, taking into account both positions and normals produces better results, and we have tried several values for different datasets. In general, we experienced the best results for α ≈ 0.2 corresponding to a higher influence of the normal coherence compared to the distance.

Concerning the quantitative results, we analyzed the number of vertex pairs located in the contact zone that are feasible to process in real-time. Of course, the required time for the ICP-based matching depends on the number of iterations. Usually, the ICP converges quite rapidly (3-4 iterations) on well-defined thresholds d T and η T combined with a good initial pose, and we are able to deal with almost 10, 000 vertex pairs in the contact surface per second on a single core Intel Pentium 4 at 3.0Ghz running Linux. When the ICP converges too slowly or diverges, our real-time interaction loop does not wait for the result and treats the next initial pose. In order to fulfill the real-time constraint, the desired frame rate can be user-adjusted to either take into account more vertex pairs, or to ensure a faster processing and thus feedback. Recall that for scanned fractured objects with millions of vertices, we first have to downsample the objects to some tens of thousands of vertices in order to satisfy the real-time constraint. For example, both fragments of the fractured head of figure 1 were downsampled to some 15, 000 vertices.

Feedback

Concerning the feedback provided to the user, the informal user study shows that the visual feedback helps to analyze the position of the best match and to reason about its plausibility. Moreover, the graphical indicator allows the user to rapidly detect whether the local matching of the contact surfaces converges well. By an interpretation of both visual feedbacks, the user has a complete understanding of the global and local coherence of the matching.

Conclusion and future work

In this paper, we have presented a semi-automatic reassembly approach based on a real-time interaction loop. This interaction loop consists of an efficient bimanual interaction technique, a real-time matching algorithm, and a way to provide a visual feedback to the user about the best match and its associated error. As a consequence, we consider the user as the key of your approach: his knowledge and his capacity to integrate semantic knowledge in the reassembly process are used to increase the performance of the matching.

Our first results of an informal user study show that our system is capable of assisting an expert user in real-time during the pairwise matching of downsampled 3D fragments. Although our algorithm is optimized with spatial data structures, it could further be accelerated by a better exploitation of the system resources (e.g. multithreading or GPU computation). Note also that for every pair of fragments, there are three parameters that have to be adjusted beforehand which confers flexibility to the technique. However, it can sometimes be tedious to correctly choose these parameters, especially for noisy and eroded objects. This is a direct consequence of using the ICP algorithm for matching, as it only takes into account positions and normals.

In future work, we plan to integrate higher order derivatives in the semi-automatic reassembly process, as for example curvatures that we can represent by second order polynomials. We also believe that an a priori analysis of the entire fragment for salient feature detection at different scales could be used to identify potential matching candidates and their initial relative poses that can then be validated by expert users by means of the visual feedback. This first selection could reduce the number of potential matching candidates and thus makes sequentiel fragment matching more efficient.

Figure 1 :

 1 Figure 1: Pairwise reassembly

Figure 3 :

 3 Figure 3: The bimanual tangible user interface task and the involved high-level knowledge without being distracted with complex interaction techniques, we rely on c The Eurographics Association 2010.

Figure 5 :

 5 Figure 5: The normal coherence.

c The Eurographics Association 2010.

Acknowledgments

3D models are courtesy by Vienna University of Technology. This work was supported by the ANR SeARCH project, grant ANR-09-CORD-019 of the French Agence Nationale de la Recherche.