Frédéric Dadeau
email: frederic.dadeau@lifc.univ-fcomte.fr

Fabien Peureux
email: fabien.peureux@lifc.univ-fcomte.fr

Grey-box Testing and Verification of Java/JML

Keywords: symbolic execution, constraint solving, Java/JML, grey-box testing, set-theory

We present in this paper the application of constraint solving techniques to the validation and automated test cases generation for Java programs, annotated with JML specifications. The Java/JML code is translated into a constraint representation based on a subset of the set-theory, which is well-suited for modelling object-oriented programs. Symbolic code execution techniques can then be applied to produce test cases, using classical structural test selection criteria, or to detect possible runtime errors, and non-conformances between the Java code and its embedded JML model.

I. INTRODUCTION

Design-by-Contract (DbC) [START_REF] Meyer | Object-oriented software construction[END_REF] is a lightweight technique that aims at introducing formal methods for non-specialists. It consists of adding contracts to the code of a program for basic modelling purposes, namely, specifying invariants on classes, and pre-and post-conditions on methods. The annotations are written directly inside the code, which simplifies its use by developers, who do not need to learn a whole new modelling language, beside some simple keywords.

Over the past few years, several annotation languages have arisen, such as JML [START_REF] Leavens | Preliminary Design of JML: a Behavioral Interface Specification Language for Java[END_REF] for Java, Spec# [START_REF] Barnett | The Spec# Programming System: An Overview[END_REF] for C#, ACSL [START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF] for C, etc. Annotation languages are very convenient at two levels [START_REF] Leavens | How the Design of JML accommodates both Runtime Assertion Checking and Formal Verification[END_REF]. At the verification level, it is possible to take the annotations into account to ease the proof of the considered program, and strengthen the properties one wants to specify on the system. At the validation level, the contracts can be used for generating black-box test data, using method preconditions, and to assign the verdict of the test, using the postcondition by runtime assertion checking techniques.

Various approaches use constraint solving techniques together with annotated programs either to produce test cases or to verify program correctness. In most cases, the program is flattened into Single Static Assignment form and a set of constraints is computed as a path condition that represents the conjunction of the decision predicates on the path.

We propose here to reinvest constraint solving techniques that have been originally designed for the symbolic animation of B machines [START_REF] Bouquet | CLPS-B: A Constraint Solver to animate a B Specification[END_REF], in order to provide a flexible and expressive means for both verifying and validating a Java program annotated with JML constructs. This technique is based on a graph representation of the program structure/code. The symbolic evaluation is then done on-the-fly, when exploring the graph to cover a specific path, that can be extracted using classical structural coverage strategies (all-nodes, all-decisions, etc.). public class A { //@ invariant x > 0; short x; static byte y; A z; / * @ requires a != null && y > x && a != this; * @ assignable \everything; * @ ensures \result > 0; * @ signals (BarException e) \old(x) == x; * / public int foo(A a) throws BarException { ... } } Fig. 1. General Form of the Contracts in JML

We show in this paper how the underlying set-theory of B, tool-supported by a customized constraint solver, can be used to efficiently represent object-oriented programs/models. The paper is organised as follows. Section II presents an overview of the Java Modeling Language syntax and semantics. Then, Section III illustrates how Java programs (i.e., data model and method codes) can be represented in the set-theory handled by our constraint solver. Section IV presents our grey-box test generation technique, that takes into account both the program and its specification. Then, we present in Sect. V a verification technique, that aims at finding three types of non-conformances between the program and its specification. Finally, Section VI presents the related works, before concluding in Sect. VII.

II. OVERVIEW OF THE JML SYNTAX AND SEMANTICS

The Java Modeling Language is an annotation language based on Design-By-Contract paradigm, introduced by Eiffel [START_REF] Meyer | Object-oriented software construction[END_REF]. The contracts are written inside the Java code, which facilitates the use of JML by developers. In this section, we briefly introduce the syntax and the semantics of JML.

Contracts can be defined as illustrated in Fig. 1. Keyword requires indicates the precondition, namely the contract that has to be fulfilled by the caller of the method. The assignable clause provides the frame condition, namely the set of attributes that can be assigned inside the method body. The ensures clause designates the normal postcondition that is supposed to be established when the method terminates normally, i.e. without throwing an exception. Finally, the signals clause designates the postcondition that is supposed to be established if the method terminates abnormally, i.e. by throwing an exception. In this case, the specification indicates the exception and its possible postcondition. In the example, the presence of an exception of any type (specified by the generic Exception class) implies that attribute x is not modified (its after value, designated by x, has to be equal to its before value, designated by \old(x)).

Besides pre-and postconditions, JML contracts are completed by class invariants, that apply on class attributes.

The JML Runtime Assertion Checker (RAC) instruments the Java code (without modifying its behaviour) to performed on-the-fly verification that will raise a specific set of exceptions in case of verification failure. The verifications performed are the following: (i) at the beginning of the method, the class invariant is checked, along with the precondition. All \old expressions used in the postconditions are evaluated. (ii) at the end of the method body execution, the normal postcondition is checked, replacing the \old expressions with their previously-stored values. (iii) if the method body throws an exception, then there are two cases: either it is a JML-specific exception, raised inside the body of the method and it is thus propagated, or it corresponds to an exception of the program, and thus the corresponding exceptional postcondition (signals clause) is checked (if it exists), before propagating the exception. When executing tests, the absence of JML-specific exceptions means that the test passed, otherwise it failed and the corresponding assertion can be retrieved to help the debugging.

III. EXPRESSING JAVA PROGRAMS IN SET-THEORY

In order to be able to apply a constraint-based verification and validation techniques on Java programs, it is mandatory to provide a representation of the target language in constraints. We start this section by the description of our underlying constraint solver, namely CLPS-BZ. We then explain how to use a set-theoretical representation of object-oriented programs, inspired by the one used in a previous work [START_REF]SICStus Prolog[END_REF]. We then introduce a graph representation of the program code, that can be explored by our symbolic animation engine in order to perform symbolic execution of the code.

A. Using Constraint Solving to evaluate Java/JML Programs

The evaluation of Java/JML programs is based on the constraint logic programming technology called CLPS-BZ [START_REF] Bouquet | CLPS-B: A Constraint Solver to animate a B Specification[END_REF]. CLPS-BZ is a constraint solver implemented in SICStus Prolog [START_REF] Bouquet | Automated Boundary Test Generation from JML Specifications[END_REF], and originally designed to evaluate and to generate functional test cases from B and Z formal models [START_REF] Legeard | Automated Boundary Testing from Z and B[END_REF]. Because of the set oriented features of both B and Z notations, CLPS-BZ is mainly a set constraints solver.

Its constraints domain is the hereditarily finite sets with nested pairs that makes it possible to manage relations and functions, but CLPS-BZ also supports constraints over integers, atoms, pairs, sets, relations and functions, including equality, inequality, membership, subset, range, domain constraints, etc.

The constraints solving techniques are based, on the one hand, on set interval propagation [START_REF] Gervet | Interval Propagation to reason about Sets: Definition and Implementation of a Practical Language[END_REF], and, on the other hand, on specific propagation rules in order to compute variable instantiations. For basic arithmetic operations and constraints solving on set cardinality, there is a cooperation between CLPS-BZ constraints propagation and the integer finite domains solver CLP(FD) of SICStus Prolog. CLPS-BZ involves partial consistency techniques [START_REF] Jaffar | Constraint Logic Programming : A survey[END_REF] during constraints propagation; the completeness is ensured by classical labelling based on AC-3 algorithm [START_REF] Tsang | Foundations of Constraint Satisfaction[END_REF].

More precisely, CLPS-BZ is composed of four subsystems: an executer, a reducer, a constraint solver and a constraint store. The executer manages the symbolic execution of the Java/JML programs: it translates this source code into a labelled transitions system representing its control flow graph, and passes this system to the reducer. The reducer then achieves the effect of the execution by translating each transition label expression into lower-level constraints that can be recursively acquired on-the-fly by the constraint solver. Finally, the constraint store module allows recording the (possibly symbolic) current state during all the execution. After all constraints have been added, if a unique solution is not obtained, the executer applies a labelling procedure to calculate a solution to the constraint system, and providing one instantiation to cover the original path of the Java/JML program given as input.

We now describe the rewriting rules applied to translate Java/JML expressions into CLPS-BZ constraints, and define the approach, based on CLPS-BZ, to generate test data for Java/JML methods using a control flow graph representation.

B. Expressing Data Model with Sets and Functions a) Instances:

We consider the set of existing objects to be a set of atoms, that are named constants. Each value designates an address of the object in the memory heap. Class instances are sets of atoms, subsets of this superset. Class instances are denoted by I C where C is the class name. Inheritance relationships make it possible to partition the sets of objects, such that: if class B and class C both directly extend class A, then

I B ⊆ I A ∧ I C ⊆ I A ∧ I B ∩ I C = ∅
Due to finite domains constraints solving, the set of existing objects has to be included in a superset of all possible instances, (similar to a set of possible object addresses), that are (or not) assigned to objects, and denoted I.

b) Object attributes: Instance attributes are represented as total functions mapping existing instances to their values. Static attributes (whose values are shared by all instances) are directly valued. Values can be either a range of integers (for primitive integer-like types: byte, short, int1 , characters, booleans), or a class. In this case, the domain is defined by the set of existing instances of the corresponding class augmented by a special value null, designating a null reference. Each aforementioned attributes value is named by a unique identifier, concatenation of its originating class and attribute.

Example 1 -Attributes naming and typing. Considering the example given in Fig. 1, attributes x, y and z will be named and typed as follows:

A.x ∈ I A → -32768..32767, A.y ∈ -128..127 and A.z ∈ I A → I A ∪ {null} where → designates a total function. c) Arrays: In Java, arrays are represented as objects of a specific type. We thus represent existing arrays similarly to existing objects using a set of dedicated instances, depending on the data type contained in the array. For example, I ashort represents the set of arrays containing short integers. Array contents are described in dedicated variables, depending on the possible content type. Arrays are typed as a set of triplets (a, i, v) where v is the content of array a at index i. For example, the set of arrays containing short integers is typed as: arrayContent ashort ∈ I ashort × 0..M AXIN T × -32768..32767. Notice that, when executing Java code, we assume that the sizes of the considered arrays are always known.

d) System States: Since we will consider the symbolic execution of the code, we need to be able to represent system states and their evolutions. Systems states are stored by the CLPS-BZ constraint store module in a map associating to the attribute name its value. This map is updated when an attribute is assigned, and referred when a value is needed. For example, knowing instance a of class A, the value of a.x is given by the image of a by the function representing the attribute x, i.e. A.x(a).

C. Expressing Method Code with Graphs

Inspired of our previous work on the symbolic animation of B machines, we have decided to represent the method code as labelled transitions systems encoding the control flow graph of the operation.

Symbolic execution of a method consists in simulating the execution of the code on the constraint representation of the program. To achieve that, we follow the execution paths of the code, guided by the control-flow graph. For each label on the path, we perform the following actions:

• if the label designates a condition, it is evaluated w.r.t.

the current values of the state variables, namely, the references to the data are replaced by the (symbolic) values provided by the system states and a constraint is added to the constraint store. When this latter becomes inconsistent, the considered path prefix is declared infeasible.

• if the label designates a statement, it is executed and the system state is updated. In a large majority of cases, the statements are in fact assignments, for which the assigned expression is first evaluated, and the corresponding data in the system state is then updated with the (possibly symbolic) value resulting of this evaluation. for computing the greatest common divisor of two positive integers arguments x and y. We provide its associated control flow graph representation.

Considering path [START_REF] Meyer | Object-oriented software construction[END_REF][START_REF] Leavens | Preliminary Design of JML: a Behavioral Interface Specification Language for Java[END_REF][START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF][START_REF] Leavens | How the Design of JML accommodates both Runtime Assertion Checking and Formal Verification[END_REF][START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF][START_REF] Bouquet | CLPS-B: A Constraint Solver to animate a B Specification[END_REF][START_REF]SICStus Prolog[END_REF]. We start with an initial state: x = X 1 , y = Y 1 , with associated constraints X 1 ∈ -32768..32767, Y 1 ∈ -32768..32767. We then perform the interpretation of the labels that change the system state:

• after 1→2: x = X 1 , y = Y 1 , r = undef • after 2→4: x = X 1 , y = Y 1 , r = undef and constraint X 1 ≥ Y 1 is added to the constraint store • after 4→5: x = X 1 , y = Y 1 , r = undef and constraint X 1 mod Y 1 = 0 is added to the constraint store • after 5→4: x = Y 1 , y = R 1 , r = R 1 and constraint R 1 = X 1 mod Y 1 is added to the constraint store • after 4→6: x = Y 1 , y = R 1 , r = R 1 and constraint Y 1 mod R 1 =
0 is added to the constraint store • after 6→7: no more constraint is added. Finally, an instantiation of the constraint store provides initial values for the parameters; here, x= -16384 and y= -32768.

Notice that the first part of the path contains conditions that represent explicit typing informations on the parameters, and on the instance on which the method is invoked.

Considered subset of Java/JML. As explained, our representation makes it possible to express the classes, their instances, and arrays. At the moment, integers are restricted to those supported by SICStus Prolog (i.e. 24-bits integers), and reals (float and double in Java) are not supported. In addition, arrays are (preferably) of one dimension.

IV. GREY-BOX TEST GENERATION

Our grey-box testing approach considers both the Java code for extracting structural test cases, and the JML specification, in order to provide the preconditions and the class invariant, that represent the invocation context of the methods.

We propose to apply classical structural coverage criteria on the graph in order to extract paths. These paths are provided as an input to our process. If the path is infeasible, it will be discarded. Otherwise, the process produces test data for input parameters and user-defined values.

The code coverage criteria that we are able to apply are the classical ones from the literature (as listed by Myers [START_REF] Myers | Art of Software Testing[END_REF]): all-statements, all-decisions, all-k-paths. In addition, condition coverage criteria can be applied to cover the disjunctive predicates in the decisions (e.g. MC/DC, MCC, etc.) [START_REF] Vilkomir | Formalization of software testing criteria using the Z notation[END_REF].

Example 3 -Grey-box test generation. Consider again the gcd example introduced previously. We add the following contract.

/ * @ requires x >= 0; @ assignable \everything; @ ensures \old(x) % \result == 0 && @ \old(y) % \result == 0; @ * / short gcd(short x, short y) { ... } Using the test data produced initially (x= -16384, y= -32768), and running the tests on an JML-instrumented program raises a PreconditionFailedException indicating an invocation of the method without satisfying the precondition contract. In JML, a direct invocation of the method out of its precondition is not considered as a failed test (but it is when the method is called out of the precondition from another method, raising in this case an InternalPreconditionFailedException). If we take the precondition into account, we produce the following test data for path [START_REF] Meyer | Object-oriented software construction[END_REF][START_REF] Leavens | Preliminary Design of JML: a Behavioral Interface Specification Language for Java[END_REF][START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF][START_REF] Leavens | How the Design of JML accommodates both Runtime Assertion Checking and Formal Verification[END_REF][START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF][START_REF] Bouquet | CLPS-B: A Constraint Solver to animate a B Specification[END_REF][START_REF]SICStus Prolog[END_REF]:

x= 3, y= 2.
The instantiation of the constraint system not only provides the values of the parameters but also an invocation context that matches the contract associated to the considered method under test. From the abstract system state representation provided with instantiated sets and functions, one can easily build JMLUnit test cases that will be run on the concrete system under test. Indeed, the Java reflection API makes it possible to access (and directly assign) object attributes, as shown in the example hereafter.

Example 4 -Test case concretization. Suppose one wants to invoke method foo(A) whose contract is shown in Fig. 1. The set of constraints associated to the object hierarchy and the contract (precondition + class invariant) is as follows. ). This set of constraints can be instantiated by the following (minimal) solution:

I = {a 0 , a 1 , a 2 ,
I = I A = {a 0 , a 1 }, A.x = {a 0 → 1, a 1 → 1}, A.y = {a 0 → 2, a 1 → -128}, A.z = {a 0 → null, a 1 → null}, this = a 0 , a = a 1
This state is then used to produce the following J(ML)Unit test case:

public void testFoo() { //

V. DETECTION OF ERRORS AND NON-CONFORMANCES

We present in this section the verifications we aim at performing using our approach. These verifications are classified into two categories: runtime errors that only focus on the Java code, and non-conformances w.r.t. the JML specification.

A. Runtime Errors

The runtime errors represent classical Java errors that may happen at runtime. Usually, such errors raise an exception that is a subclass of RuntimeException. We will focus on the most frequent ones. The presence of these runtime errors can be checked by the satisfiability of an additional constraint involving that represents a condition that has to be fulfilled for the exception to happen.

We give here the list of considered runtime errors, for which we provide the context of the error detection, and the considered additional constraint that represents the condition under which the exception will be raised.

Error

Context Constraint Null pointer deref.

x.y is evaluated x == null Index out of bounds t[i] is accessed i > t.length Division by zero

x/y or x%y y == 0

When the additional constraint is acquired, a valuation is immediately asked to the solver to check the overall consistency of the constraint store. If the valuation succeeds then it provides a counter-example that shows the potential presence of a runtime error.

Example 5 -Runtime error detection. Consider the example of the gcd method introduced before, when performing the exploration of the control flow graph. When encountering the condition labeling transition 4→5, the additional constraint y == 0 is stored and a labeling is performed. In this case, the constraint system produces an example of test data, x= 0 and y= 0, that leads to an error of division by zero.

Notice that, as a side effect, our methodology makes it possible to detect infeasible paths. However, in a large majority of cases, the paths are extracted blindly from the graph structure and are, by definition, not necessarily feasible. Nevertheless, detecting loop-free prefixes that are not feasible indicates dead code.

Example 6 -Dead code detection. Suppose we enrich the precondition of the gcd example to add the following information:

/ * @ requires y > 0 && x > y; @ ... @ * / short gcd(short x, short y) { ... }

In this case, we are able to detect that all paths beginning with [START_REF] Meyer | Object-oriented software construction[END_REF][START_REF] Leavens | Preliminary Design of JML: a Behavioral Interface Specification Language for Java[END_REF][START_REF] Barnett | The Spec# Programming System: An Overview[END_REF] are unfeasible since the condition on the label of transition 2→3 is contradictory with the JML precondition.

B. Non-Conformances w.r.t. the JML Specification

The second kind of verifications we consider involves the verification of the contracts that makes it possible to detect non-conformances between the program and its specification. Namely, we will consider three kinds of non-conformances that can be checked at runtime. We briefly detail these nonconformances and informally indicate how these verifications can be performed.

• Provider non-conformance: when the method terminates, it does not establish its contractual postcondition. Such a non-conformance is detected by negating the postcondition at the end of the method's execution, and verifying the satisfiability of this negation.

Client non-conformance: when the method invokes another internal method, the precondition of this latter is not fulfilled. Such a non-conformance is detected by negating the precondition of the invoked method, and verifying the satisfiability of this negation. • Integrity non-conformance: after the execution of the method, the class invariant is not satisfied. Such a nonconformance is detected by negating the invariant at the end of the method's execution. Since these verifications can not be done for all the possible paths of the methods code, we propose to apply them on demand, on a relevant subset of the paths that is computed using the classical code coverage criteria (all-statements, alldecisions, etc.).

Example 7 -Provider Non-Conformance. Consider again the example of the GCD algorithm presented in Fig 2 . Suppose one wants to check that path [START_REF] Meyer | Object-oriented software construction[END_REF][START_REF] Leavens | Preliminary Design of JML: a Behavioral Interface Specification Language for Java[END_REF][START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF][START_REF] Leavens | How the Design of JML accommodates both Runtime Assertion Checking and Formal Verification[END_REF][START_REF] Baudin | ACSL: ANSI C Specification Language[END_REF][START_REF] Bouquet | CLPS-B: A Constraint Solver to animate a B Specification[END_REF][START_REF]SICStus Prolog[END_REF] does not present a provider non-conformance. Considering the state resulting of the execution of this path:

x = Y 1 , y = R 1 with associated constraints X 1 ∈ -32768..32767, Y 1 ∈ -32768..32767, Y 1 > 0, X 1 > Y 1 , X 1 mod Y 1 = 0, R 1 = X 1 mod Y 1 , Y 1 mod R 1 = 0
The negation of the postcondition provides the additional constraint:

X 1 mod R 1 = 0 or Y 1 mod R 1 = 0
Each of these literals is added independently to the previous store and the satisfiability of the resulting constraint system is checked, to look for counter-example. In both cases, the solver concludes that each of these new constraints breaks the consistency of the constraint store. It is thus possible to conclude that the considered path does contain a nonconformance.

VI. RELATED WORKS

The present work is an extension of an original work on test generation from JML specifications [START_REF]SICStus Prolog[END_REF], for which the same data model representation was used. The main difference is the use of the Java code, that was not considered previously.

Many works deal with the verification and validation of Java programs possibly using the JML specifications. JMLUnit [START_REF] Cheon | A Simple and Practical Approach to Unit Testing: The JML and JUnit Way[END_REF] and Jartege [START_REF] Oriat | Jartege: A Tool for Random Generation of Unit Tests for Java Classes[END_REF] use random test data generators to respectively build unit test cases and test sequences. In both cases, the JML specification is used to filter irrelevant test data that may be produced, and to assign the test verdict by runtime assertion checking [START_REF] Leavens | How the Design of JML accommodates both Runtime Assertion Checking and Formal Verification[END_REF].

(j)Cute [START_REF] Sen | Cute and jcute: Concolic unit testing and explicit path model-checking tools[END_REF], DART [START_REF] Godefroid | DART: Directed Automated Random Testing[END_REF], EXE [START_REF] Cadar | EXE: Automatically Generating Inputs of Death[END_REF], Pathcrawler [START_REF] Williams | PathCrawler: Automatic Generation of Path Tests by Combining Static and Dynamic Analysis[END_REF] use symbolic execution for the automated structural test generation using branch coverage strategies. Some of these tool performs both symbolic and concrete execution of the code (also known as concolic) so as to generate test data and check the coverage of the corresponding control flow graph. Recently Pathcrawler and Frama-C2 have been coupled by the SANTE tool [START_REF] Chebaro | Combining Static Analysis and Test Generation for C Program Debugging[END_REF] to combine test generation and verification of C programs.

In the context of verification of JML/Java programs, ESC/Java [START_REF] Flanagan | Extended Static Checking for Java[END_REF] uses theorem provers to check verification conditions extracted from Java/JML program. The KeY project [START_REF]Verification of Object-Oriented Software: The KeY Approach[END_REF] aims at similar purposes. Bogor [START_REF] Robby | Checking JML Specifications using an Extensible Software Model Checking Framework[END_REF] uses a model checking approach to explore Java system states and check the JML assertions.

Our approach is inspired by the InKA approach [START_REF] Gotlieb | Automatic Test Data Generation using Constraint Solving Techniques[END_REF], a goal oriented structural test generator using constraint solving. More recently, Euclide [START_REF] Gotlieb | Euclide: A Constraint-Based Testing Framework for Critical C Programs[END_REF] acts in a similar way, taking into account C programs annotated with ACSL contracts. This work represents a continuation of the work done in the DANOCOPS project 3 . A similar approach has been proposed by Collavizza et al. [START_REF] Collavizza | A Dynamic Constraint-Based BMC Strategy For Generating Counterexamples[END_REF], using a bounded model-checking approach to generate counter-examples of program properties using the Comet platform.

VII. CONCLUSION AND FUTURE WORKS

We have presented in this paper the use of a set-theoretical representation to express Java programs annotated with JML contracts. This representation makes it possible to perform symbolic execution of the Java code of the methods, under the assumptions of the JML preconditions, which can be used for both generating structural test data, or performing code-based verification of assertions, by looking for non-conformances between the program and its specification.

Our future works will consider the improvement of our constraint solving techniques in two directions. First, we have noticed that, in the absence of errors, the AC-3 consistency algorithm used to store the constraints requires the constraints system to be eventually instantiated. In these cases, it is important to notice that a given number of symmetric solutions are enumerated, slowing down the verification process, without changing the result of the computation. The second lead to investigate is to remove the finite domains assumption concerning the set of existing objects. On one hand, the overall set of objects has to be finite to guarantee the termination of the consistency checks, that requires an enumeration. On the other hand, the selected bound might be insufficient to dispose of enough objects for falsifying a verification condition. We are now investigating a means to address these two issues together, in order to improve the execution time and the correctness of the verification and test generation processes. Another option to solve these problems would be to consider another constraint solver that would also handle reals, such as Choco in Java [START_REF] Choco Team | choco: an Open Source Java Constraint Programming Library[END_REF]. Nevertheless, it would be necessary to also handle functions and relations, so as to be able to conserve this convenient representation of object-oriented data model.

Finally, we are currently implementing a fully automated tool-chain of this test generation process (especially to automate the translation of Java/JML programs into labelled transition systems representation). This tool-chain will take Java/JML code as input, and automatically derive executable JMLUnit test cases. This achievement will make it possible to more intensively carry out various experiments and get relevant feedback on the scalability of our approach.

Example 2 -Fig. 2 .

 22 Fig. 2. The GCD example

  a 3 , ...}, // bounded set of existing objects I A ⊆ I, // structural constraint for class hierarchy A.x ∈ I A → -32768..32767, // typing of attribute x A.y ∈ -128..127, // typing of static attribute y A.z ∈ I A → I A ∪ {null}, // typing of attribute z ∀p ∈ I A ⇒ A.x(p) > 0 // invariant of class A this ∈ I A , this = null, // existence of this a ∈ I A ∪ {null}, // typing of parameter a a = null, // precondition of foo A.y(this) > A.x(this), // precondition of foo a = this // precondition of foo Notice that this example uses native operators of the CLPS-BZ solver (∈, ∀, ⊆, etc.

with a restriction on this data type since SICStus Prolog does not represent integers on 32 bits, but on

bits

http://frama-c.cea.fr

http://lifc.univ-fcomte.fr/DANOCOPS


Editorial: Viscoplastic Fluids: From theory to application 2013

Xavier Chateau a , Anthony Wachs b, * a Laboratoire Navier (UMR 8205), IFSTTAR, University Paris-Est, Marne-la-Vallée, France. b IFP Energies nouvelles, Fluid Mechanics Department, Rond-point de l'Echangeur de Solaize, BP 3, 69360 Solaize, France.

This issue of the Journal of Non-Newtonian Fluid Mechanics includes a series of papers based on work presented at the international workshop on Viscoplastic uids: from theory to application, held Nov. [18][19][20][21]2013 in Rueil Malmaison, France. A list of participants is provided in Table 1. This was the fth biannual meeting on this subject. The previous meetings were held in Ban (Alberta, Canada), Monte Verita (Ascona, Switzerland), Limassol (Cyprus) and Rio de Janeiro (Brasil).

Like previous editions, the aim of the workshop was to bring together leading researchers in the eld of viscoplastic uids across several disciplines to foster the awareness and the transfer of ideas, both from academic research and industry. The program consisted in a single technical session and three invited keynote lectures. A total of 49 talks and 9 posters were presented, spanning fascinating topics from the coating of viscoplastic uid on a plate to the numerical simulation of the transition of viscoplastic uid ows to turbulence.

The workshop enjoyed an unprecedented number of 80 participants, a popularity that emphasizes that viscoplastic uids are a "hot" topic in the eld of non-Newtonian Fluid Mechanics.

Invited speakers were John Tsamopoulos (University of Patras, Greece), Guillaume Olvarez (Laboratoire Navier, Université Paris-Est, France) and Fabrice Toussaint (Lafarge Centre de Recherche, France). John Tsamopoulos opened the meeting with an invited talk on the numerical simulation of yield stress uid ows and its applications to the problem of the rising of a single bubble. He both elaborated on the technical details of the numerical tools available nowadays for the simulation of this class of ows and provided insight on the eect of elasticity and pressure oscillations. On the second day of the workshop, Guillaume Olvarez gave an overview on the rheological behavior of suspensions of particles and bubbles in yield stress uid, illustrating his ndings with a large number of experimental results.

On the third day, Fabrice Toussaint bridged academic reseaach and industrial concerns in the eld of concrete rheometry, ranging from the fundamental behavior of concrete slurry and their characterization in rheometers to the fascinating last i-Phone application designed to tell the truck driver how quick the concrete slurry properties are evolving in the truck's rotating drum while driving to the construction site. Standard talks ranged from experiments, theory to numerical simulations, with an emphasis on the rheological behavior of carbopol gels and the more general question of how to experimentally characterize yield stress materials with thixotropy, a question that has been mobilizing the attention of the yield stress community for a few decades. The realistic and accurate modelling of the behavior of viscoplastic and thixotropic still remains an unsolved question in the eld, as already underlined two years ago in Rio de Janeiro, but also progress has been made, with new measurement techniques as, e.g., low amplitude oscillatory shear (LAOS). In particular, the workshop featured intense discussion on thixotropy and its mathematical modelling. Eorts in designing new numerical approaches with enhanced accuracy and fast convergence has seemed to slow down and the workshop was an occasion to collectively acknowledge that this research direction should be revived.

The workshop took place at IFPEN-Rueil Malmaison, in the suburb of Paris, France.

IFPEN is a large national research center in the eld of energy, with nice facilities and well designed services. The social program started with an ice-breaking cocktail in a Frenchstyle coee shop named "Café Lee" in the heart of Rueil Malmaison, which was a casual occasion to chat with long-term colleagues and enables new participants to VPF to immerse into the yield (no) stress community and yield to the pleasure of enjoying French wines and the warm and friendly atmosphere of the place. The grand dinner was organized in the magnicent "Chateau de la Malmaison", a genuine castle from the Napoleon era (early 18th century) now home to the Tuck Foundation since 1992. This year the organizers had decided not to deliver any "best talk" or "best poster" prices, in an attempt not to discriminate the contributions, and also because all were of equal amazing quality ! However, the grand dinner was the opportunity to give the "Bingham chocolate medal" to a long-term and high-level contributor to the eld of viscoplastic uids. The spirit of this prize is towfold: serious and casual, serious means that the organizing committee strongly believes that the scientic contribution of the laureate is outstanding, while the "chocolate" adjective reminds everyone that giving a prize to someone in particular is a highly subjective choice and should be considered with a certain degree of detachment. The prize was given to Prof. Ian Wilson for his work of the viscoplasticity of uids in forming and extrusion process over the grand dinner and the audience had the privilege to listen to an amazing and hilarious speech from Prof. Ian Frigaard, which tted the spirit of the prize, i.e., a combination of sincere credit to Ian Wilson's scientic contribution to the eld as well as both genuine and imaginary stories about his career. As in Rio two years before, participants could compete in a rheology contest which consisted in guessing the plastic viscosity and yield stress of four materials including chocolate paste, mustard, cancoillotte

(

  Figure 1 Group photo in the park at IFPEN-Rueil Malmaison

A. Ahuja G. Georgiou L. Muravleva A. Alexandrou B. Geraud M. Naccache D. Andrade S. Gonzalez-andrade C. Negrao C. Anglade E. Gouillard N. Nernabeu C. Barentin I. Hénaut A. Nir F. Bertails-descoubes B. Herzhaft G. Ovarlez M. Briceno R. Holm T. Palermo M. Bryan S. Hormozi Y. Peysson T. Burghelea R. Huilgol R. Phillips C. Castelain I. Ionescu T. Pringuey X. Chateau L. Jossic Quignon G. Chambon L. Jorgensen F. Rasschaert T. Chevalier I. Karimfazli S. Rosales anzola R. Chhabra Z. Kebiche M. Rukhadze M.-F. Chirac E. Lac E. Sanchez P. Coussot M. Lambert P. Saramito Y. Damianou F. Link V. Shelukhin A. Davaille Y. Liu A. Syrakos G. Daviet S. Lovett L. Talon J. De bruyn L.-H. Luu F. Toussaint M. Denn M. Maillard J. Tsamopoulos P. De souza mendes R. Mendes O. Thual B. Dollet G. Moises G. Vinay J. Ferec M. Moyers gonzalez A. Wachs S. Frey M. Mouzouri I. Wilson I. Frigaard E. Muravleva N. Zeraibi Table 1 Participants to the 5th workshop on Viscoplastic uids: from theory to application


Jean-Charles Papin

Christophe Denoual

Laurent Colombet

Raymond Namyst

SPAWN: An Iterative, Potentials-Based, Dynamic Scheduling and Partitioning Tool

Keywords: Simulation, dynamic load-balancing, tasks, manycore, pair potential, graph partitioning

Many applications of physics modeling use regular meshes on which computations of highly variable cost can occur. Distributing the underlying cells over manycore architectures is a critical load balancing step that should increase the period until another step is required. Graph partitioning tools are known to be very effective for such problems, but they exhibit scalability problems as the number of cores and the number of cells increases. We introduce a dynamic task scheduling approach inspired by physical particles interactions. Our method allows cores to virtually move over a 2D/3D mesh of tasks and uses a Voronoi domain decomposition to balance workload among cores. Displacements of cores are the result of force computations using a carefully chosen pair potential. We evaluate our method against graph partitioning tools and existing task schedulers with a representative physical application, and demonstrate the relevance of our approach.

INTRODUCTION

Many physics simulation applications rely on large meshes where each cell contains a few elementary computing elements (e.g. particles, finite elements, or finite difference cells) and is linked to its neighboring cells. Figure 1 illustrates the concept of a regular mesh. The computing cost of each cell varies with the modeled material and the applied action (e.g. shock wave or distortion). Thus, distributing cells among computing units must preserve locality of neighbors and balance the computing load. Mesh partitioning can either be achieved with common task scheduling tools or graph partitioning tools.

Task-based parallelism is now a well-known concept. Along with specific libraries [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF][START_REF] Kukanov | The foundations for scalable multi-core software in intel threading building blocks[END_REF]21] and languages extensions [16], many languages now integrate tasks support in their standard [START_REF]The rust programming language[END_REF], such as C++11 with its future variables concept. This allows programmers to easily exploit multiple computing units. Tasks are generally associated with a data set and can exchange data with other tasks, allowing to define an affinity criteria for a given task. Therefore, a particular attention shall be paid to the place (i.e. the computing unit) where the task will be executed. However, the still increasing cores-on-chip and chips-in-system number [START_REF] Hemmert | Green hpc: From nice to necessity[END_REF], forces applications to spawn a large number of tasks so as to prevent cores from reaching a starving state. This eventually leads to an important complexity: since manifold scheduling policies try to distribute tasks among available cores by satisfying these affinity criteria while minimising the overall application computing time, the number of combinations increases with the number of available cores and tasks.

Exploiting the fact that tasks have strong affinities with their neighborhood (direct neighbors and neighbors of 2nd-3rd rank, see fig. 1), such meshes can direcly be represented by graphs where each node is a task, and where a link is a neighboring connection. Graph partitioning tools [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF][START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] can thus be used to distribute tasks over cores in this context. Partitioning tools provide very efficient algorithms, for partitioning very large regular and non-regular graphs, which can take into account affinity characteristics. However, such tools do not perform well when dynamic scheduling is needed. Slight task weight variations will usually lead to the generation of a completely different graph partition, which will incur numerous data transfers to re-assign tasks. Support for dynamic scheduling was thus added to recent and distributed versions of partitioning libraries [START_REF] Chevalier | PT-Scotch: A tool for efficient parallel graph ordering[END_REF][START_REF] Karypis | Parallel multilevel k-way partitioning scheme for irregular graphs[END_REF][START_REF] Devine | Design of dynamic load-balancing tools for parallel applications[END_REF]: they offer refinement functionality by using the previous computed distribution to compute the next one. On modern architecture (i.e. the Intel R Xeon Phi), the number of cores is significant (typically 61 cores for Intel R Xeon Phi 7100 series) and the number of available threads unit doubles or quadruples this number. The amount of memory per core is consequently a limitating factor (less than 300 Megabytes which must be shared among threads) and low memory footprint load balancing mechanisms mush thus be designed.

We present Spawn, a physical interaction inspired scheduler that produces compact and optimal voronoi domains. In our case, Voronoi diagrams maximize per-core data locality, by providing numerous advantages: cache usage improvements, more efficient NUMA-aware memory allocation/accesses and less Point-to-Point communications. This scheduler has the advantage of being efficient to compute and offers an automatic refinement. This paper is organized as follows: section 2 presents the rationale behind using physical interactions for dynamic scheduling while section 3 focuses on physical background and implementation details. Section 4 presents a set of experiments that compare our approach with graph partitioning tools and dynamic task schedulers. Concluding remarks and future work are discussed in the last section.

TASK SCHEDULING WITH POTENTIALS

We represent the set of cells by a multi-dimensional grid (see fig. 2). Each cell has its computing cost of its own, which can evolve over time. We insert a virtual representation of a computing unit on the grid, which we call a vCore. We then gather tasks around vCores, with a Voronoi tessellation [START_REF] Aurenhammer | Voronoi diagrams[END_REF], by using each vCore as a Voronoi site. Voronoi diagrams insure compact sets and provide geometric stability [START_REF] Reem | The geometric stability of voronoi diagrams with respect to small changes of the sites[END_REF]. They are used in many scientific applications (biology, chemistry, computer graphics, . . . ). These vCores are able to move on the grid. When moving, vCores modify Voronoi domains, which re-equilibrates in turn vCores load. Their displacements result of forces computation (with a pair potential): attractive or repulsive forces between vCores are balanced by the vCores computing load (total cost of Voronoi cell). Overloaded cores are more attractive than underloaded cores, which are then more repulsive. Table 1: Analogy summarizes between molecular dynamic simulation and cell scheduling.

PHYSICAL CONSIDERATIONS

Molecular Dynamic (MD) describes the evolution of a set of particles (a system) over time by computing forces between particles, velocities and displacements. The behavior of the whole system is defined by a potential that characterizes interactions between particles and, once derived, to forces between particles. In the following, we make an analogy between MD and task scheduling/graph partitioning (a cell can be either a task or a graph node). Table 1 summarizes this analogy.

The computing load Qi of a vCore is the sum of the cost q k of all the cells in its Voronoi domain (see eq. 1). Hence, a slight cell cost variation has a direct influence on the corresponding vCore load.

Qi =

k∈V oronoi(i)

q k . ( 1
)
Based on this load information, we define forces between vCores (see fig. 3 and the next section), and then, we move vCores. This point is described in section 3.2, while the next section gives physical background used to compute vCores displacements.

Evolution of a set of charged particles

Let us consider the Coulomb force between two particles i, j separated by a vector rij = xixj with electrical charge Ci, Cj:

Fij = λ rij |rij| 3 . with λ = CiCj (2)
Different charge signs lead to attractive forces whereas charges with the same sign produce repulsive ones. Pairs interactions of this N-Body system are calculated by exploiting the symmetry of interactions, i.e., Fij = -Fji. In order to get the relaxed state only, we minimize the potential by using a steepest descent algorithm: dx/dt = -αFij (with α a positive scalar), and by lumping α and time increment into a simple scalar k:

x i (t) = x i (t -1) -k j Fij (3)
We rescale k for every step so that the distance x(t) -x(t -1) is a fraction of the cell box dimension, which ensures convergence to stable or metastable states.

Force

Our potential should minimize the overall computing time.

For this purpose, we define the optimal vCore load m as the sum of all vCores load divided by their number (this choice is discussed in section 5):

m = 1 N i=1,N
Qi with N , the number of VCores.

The difference m -Qi thus represents the distance between the load of a particle i and the average (optimal) load m. For two interacting VCores i and j, we take a force proportional to the sum of two distances (m -Qi and m -Qj) instead of the product of charges CiCj.

Fij = λ rij |rij| 3 + rij |rij| 5 with λ = 1 - Qi + Qj 2m . ( 5
)
This potential works in this way: an underloaded vCore repulses neighboring vCores through repulsive forces. Thus, its own Voronoi domain will grow by collecting cells. In the same way, when a vCore is overloaded, it has an attractive force, leading other vCores to grab cells from it. The last term of equation 5 is used as a short-range repulsive term preventing dipole formation (see details in [START_REF] Papin | Dynamic load balancing with pair potentials[END_REF]).

Iterative model

Since we need to compute particle motions, we face an iterative model. Thus, we have introduced two input parameters, P and I, that respectively characterize the requested precision and the maximum number of iterations used to reach this precision. The precision is a percentage that represents the maximum allowed imbalance. This interface allows a fine-grained usage of the scheduler which permits us to compute a task scheduling only when VCores are unbalanced, during the idle time of computing units (see 4.3 for details).

Figure 4 shows that after 100 iterations, the benefit per iteration is close to zero. Therefore, trying to reach a better task distribution (i.e., trying to reduce P ) is pointless. Thus a trade-off between P and I must be found. Section 4.3 shows the way we solve this problem.

Oscillations

In the following, we analyse the relevance of our optimisations and implementation by considering the variance of the vCore load, an optimal vCores configuration being characterized by a null variance.

Slight oscillations near optimal configuration are visible (see figure 4). These oscillations are due to cells at Voronoi boundaries that can throw out of boundary zone for small variation of vCore location. In order to avoid this, we weight the cost of cells by a scalar that decreases with the distance from the boundary of the Voronoi domain (see eq. for the plotted curve). This implies that cells at boundaries will have a lower weight than cells close to the Voronoi domain center, leading to a more stable Voronoi site. For a given cell i and its two closets Voronoi site A and B, we define the weighted cost, f (Q), defined as:

f (Q) = Q • 1 2 1 -tanh(k ÃB • Ãi - 1 2 | ÃB|) (6)
In this equation, the scalar k defines the width of the transition area, i.e, the width of the boundary in which cells are smoothed. For our needs, we simply set k = 1, so that only one row of cells is smoothed. Figure 4 shows the affect on oscillations. Therefore, equation 5 becomes:

Fij = λ rij |rij| 3 + rij |rij| 5 with λ = 1- i k=0 f (Qi) + j k=0 f (Qj) 2m .
(7)

Voronoi diagrams

The most time-consuming part of our implementation lies in the Voronoi diagram construction. Voronoi diagrams are computed frequently, since every partitioning triggers a new diagram construction. Since our algorithm actually computes Voronoi sites which compose the whole diagram, we do not need to rely on an external voronoi diagram construction library [START_REF]Boost. Boost c++ libraries[END_REF][START_REF] Cgal | Cgal -voronoi diagrams[END_REF][START_REF]Voro++. Voro++ -a 3d voronoi cell software library[END_REF]. We only need to attach each cell to a site in order to create Voronoi domains. Hence, we need to determine, for each cell, the closest vCore. A naive implementation is in O(n•m) (n, the number of cells, m, the number of vCores), but several more efficient nearest-neighborsearch algorithms obviously exist [START_REF] Louis | Multidimensional binary search trees used for associative searching[END_REF][START_REF] Kuan | Fast k nearest neighbour search for r-tree family[END_REF][START_REF]Boost. Boost c++ libraries[END_REF]. The Vantage-Point tree [START_REF] Peter | Data structures and algorithms for nearest neighbor search in general metric spaces[END_REF], which has good construction (O(n • log(n)) and search performance (O(log(n)) is considered in the following.

RESULTS

To evaluate our scheduler, we have developed two representative stencil-like model applications. The first one, uses a virtual cells domain, where every cell has a given cost generated by a random Perlin noise [START_REF] Perlin | Perlin noise is a computer-generated visual effect developed by ken perlin, who won an academy award for technical achievement for inventing it[END_REF]. We don't do any computation with cells. This allows us to compare the quality and the time needed for the partition of the domain (see 4.1.1). The second one uses a real cells domain and provides an intensive memory usage with data exchanges between cells. Every cell has 4 direct neighbors (four cardinal directions), and needs data (a matrix of double) from its neighbors to compute its own matrix. These exchanged matrices are the result of the previous iteration. In a distributed memory (see section4.1.2), this allows us to evaluate the number of MPI communications, induced after the new domain decomposition and, during computation (data exchange between meshes). We evaluate our scheduler with sequential partitioning tools and parallel partitioning tools (with partitioning refinement enabled). Section 4.2 evaluate the same application, but in a shared memory, by comparing our scheduler against common task scheduling strategies. Thus, our analysis focuses on cache miss rates.

Comparison with Graph Partitioners

Scheduling quality & efficiency

This section focuses on the time needed to compute one partition, and the quality of the resulting partition. For this specific purpose, we compare our scheduler against Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF] (LaBRI/INRIA Bordeaux Sud-Ouest) and Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] (University of Minnesota). With Spawn, we compute the domain of each vCores, and with Scotch or Metis, we compute a partition with [number of vCores] parts. For every single configuration (dimension of the grid, number of vCores), we run thousands of sequential tests on a E5-2650 Intel R Xeon (2.0 GHz) processor. Figure 6 shows an example of generated Perlin noise along with different generated partitions, in terms of shape, after a call to each of these tools. Figure 7 shows results for a grid of 256x256 cells, with a variable number of vCores. For our algorithm, we use P = 20 and I = 20. These two charts represent the time (x-axis) needed to compute a partition of a given quality (y-axis). The quality of the partition is the distance (in percentage) to the optimal configuration in the interval [optimal, 2 × optimal]. The left chart shows the time needed to compute a partition for a given initial configuration, while the right chart represents the time needed to compute a partition after a load variation on the task domain. Scotch and Metis perform at least twice faster than Spawn when computing the initial partition. This is related to our iterative model and to the initial vCores position that can be far from the optimal position. On the other hand (right chart), Spawn shows better results after a load variation. This is a consequence of the model, we use the previous decomposition to compute the new one. Thus, if the previous decomposition is relatively close to the new optimal, after a load variation, we need a few iterations to reach the new optimal task distribution. If we decrease the P parameter (in order to reach a better task distribution), the benefit vanishes and requires a higher computing time (see bottom chart of figure 7).

We can notice an important computing time and a far-fromoptimal partitioning when using 512 vCore. We explain this behavior by a over demanding P parameter, 20% in this case, with respect to the configuration. Since we use a Voronoi tessellation, when vCores move, they exchange more than one task, and possibly a full boundary. Consequently, the 20% imbalance parameter represents a cost that is weaker than the cost vCores may exchange. Figure 8 confirms this observation: we compare the needed time (y-axis) to reach a given requested imbalance degree (x-axis), within at most 1,000 iterations (I parameter). We can see that for the same number of vCores, but for a variable domain size, we are unable to go further a given imbalance. In conclusion, this approach produces a task distribution with a quality comparable to the one produced by Scotch or Metis, with a slight overhead. However, due to the geometric shape of cells, leading the optimal task distribution might be difficult in regard to the system configuration (size and number of vCores) and to input parameters P and I. Nevertheless, our way of using our scheduler (see sec. [START_REF] Louis | Multidimensional binary search trees used for associative searching[END_REF] background and as long as an imbalance exists.

Performance

This section focuses on performance of Spawn, in a dynamic task scheduling context, by using an iterative model application (as described in the beginning of section 4). For both partitioning tools and Spawn usage, a new partition is computed between each new iteration, with newest task costs information. Since such tools use graphs as input, we use tasks as graph vertices, and communication pattern between tasks as edges. For both tools, we use a criterion that produces compact sets of vertices, with an equilibrium between the load of sets.

Shared Memory. Figure 9 compares performance of Spawn in relation to Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF] and Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] libraries. Left figures refer to the static load case while right charts refer to the dynamic load case, and, upper charts refer to the speedup while the bottom ones refer to the average number of L1 data cache misses per threads for one iteration.

Here, all partitioning strategies produce compact sets, minimizing data transfers between tasks, but Spawn has better performance in both static and dynamic load variations. As explained in section 1, Scotch and Metis produce a new and totally different task distribution after a task costs variation. Due to hardware and Operating System behavior, a task cost may vary over time, even if it's the same task and it does the same things. Hence, such partitioning tools produce different task distributions after each iteration. This involves data displacements between sockets and higher cache miss rates (see bottom charts in fig: 9).

Furthermore, in section 4.1.1, we saw that Spawn has a higher computing cost. Figure 9 puts in evidence that this cost is fully hidden by the quality and benefits of the resulting scheduling. Distributed Memory. Figures 10 and 11 present the same evaluation but in a distributed memory environment (using MPI libraries). Measurements focus on speedup and data transfer rate between each iteration. We use matrices of 128x128 doubles in a domain of 512x512 tasks (262,144 tasks). Tests run on a cluster of Intel R Xeon X7650 nodes, interconnected with an Infiniband QDR network. Figure 10 shows the average number of exchanged tasks on the MPI network (in percentage of the entire domain) and the resulting application speedup. Scotch and Metis involve an important data exchange since they aren't designed for dynamic partitioning: nearly the entire task domain is exchanged for each iteration while less than 5% of the domain is exchanged with Spawn. Parallel and Refinement-Capable Graph Partitioners.

Tools like PTScotch [START_REF] Chevalier | PT-Scotch: A tool for efficient parallel graph ordering[END_REF] and ParMETIS [START_REF] Karypis | Parallel multilevel k-way partitioning scheme for irregular graphs[END_REF] can perform parallel graph partitioning. ParMETIS supports refinement graph partitioning, i.e., it can take into account the actual partition to produce a new partition with the newest task costs.

Zoltan [START_REF] Devine | Design of dynamic load-balancing tools for parallel applications[END_REF] is another partitioning tool that can use not only its internal algorithms but also Scotch/PTScotch or ParMETIS to compute a graph partitioning.

One limitation of these tools is that they need to know, for every local tasks and for every neighboring local or remote task, the identifier of the node that owns it. This implies to maintain on each MPI node, and with every call to a partitioning routine, the list of neighboring nodes, and the list of tasks they share between them.

Next evaluation focuses on Zoltan and ParMETIS (with Zoltan) tools. We evaluate Zoltan with PHG and RCB partitioning algorithms. PHG is the internal Parallel Hypegraph and Graph partitioning method. It supports initial partitioning, refinement and re-partitioning operations. Both repartitioning and refinement operations reuse the current partition, but refinement operation is stricter. The RCB method is a geometric cutting algorithm that recursively divide the computational domain in two sub-domains, of equal work load, until the number of sub-domains is equal to the number of requested parts. We use one internal algorithm option, RCB_REUSE, that indicates whether previous cuts should be used as initial guesses for the current RCB partition.

Figure 11 shows data transfer rates and speedups of Zoltan (with PHG, RCB and ParMETIS) in comparison with Spawn.

We can see that Spawn performance are similar to graph partitioning and geometric algorithms when the network isn't used (i.e. one computing node used). For a few number of MPI used nodes, Spawn is still better than Zoltan or ParMETIS, while, for a higher number of MPI nodes, RCB algorithm outperforms others algorithms. In a more detailed point of view, if we compare the RCB algorithm to Spawn, we can see (left chart), that Spawn produces less data transfers (1-3% of the entire domain versus 5-35% for RCB). This explains that our algorithm involves a better speedup when the network isn't used (time spent in internal communication is regained by the time not spent in exchanging data).

On the other hand, with a higher number of nodes, central behavior of Spawn is penalizing and collective MPI operations considerably slow down the global computation.

This observation is, however, incorrect if we compare Spawn to the PHG method or to ParMETIS. In comparison with Spawn, data transfer rates are more important for PHG, and really similar for ParMETIS. Thus, we could expect better results for, at least, ParMETIS. By considering that partitioning tools have scalability limitations (because of internal synchronizations, exacerbate by the network), we can consider that data transfers are less penalizing than scalability limitations, which explain PHG and ParMETIS results shown in figure 11.

Comparison with Common Scheduling Strategies

We have extended the StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] runtime (the development version, see the next section for details) with our scheduler to evaluate Spawn performance against common task scheduler. Data block of tasks are managed by StarPU through the starpu_block_data_register routine, which enables automatic and backgrounded data transfers. We evaluate our scheduler on a two-sockets Intel R Xeon E5-2680 (Ivy-Bridge, 2x10 cores, 2.8GHz) machine. We did a hundred of runs with a domain of 100x100 tasks, where each task contains a matrix of 8x8 doubles. We evaluate performance in both static and dynamic load variations over time.

Integration

StarPU comes with a variety of scheduling policies, but also offers the capability to define new ones. This allows us to manage StarPU workers (abstract representation of computing units) ourselves. Spawn is used by StarPU through an intermediate meta-scheduler that implements interfaces required by StarPU (the starpu_sched_policy C-structure).

When tasks are submitted to StarPU, they are forwarded to the scheduling layer, and so, to our meta-scheduler. This meta-scheduler asks our library (Spawn), to which worker this task should be assigned to (more precisely, to which Voronoi domain this task is attached to).

If it is the first time that our library has been called, then it hasn't any information on tasks, therefore, one initial distribution is used. This initial distribution depends on the initial position of vCores, which can be, for example, a regular grid or a random position. Once tasks are distributed among StarPU workers, their execution (i.e. the computing time) is monitored and then sent to our library.

The scheduling can then be achieved in an imperceptible way for the final user (see Fig. 12): it is computed by the first StarPU worker that has finished its own set of tasks, and does it as long as another worker is still working. Thus, we only compute a task distribution when an imbalance appears. This has also the advantage to allow us to set P = 1 and I = 1 as parameters of our algorithm (which can be interpreted as "as long as an imbalance exists, do one iteration to reach the most optimal distribution"). As explained, a new task distribution is based on the task information (i.e. computing cost) of the previous tasks execution. Therefore, the first two load balancing steps might not be optimal. data transfer time during task placement. More information on this can be found in the StarPU HandBook [START_REF]INRIA Runtime team. Starpu handbook[END_REF].

Results

In both static and dynamic load variation cases, Spawn is close to StarPU schedulers performance and outperforms it with the growth of the number of threads. Bottom parts of figures show that Spawn dramatically reduces the number of L1 cache misses with a factor between two and four, thus explaining the gain. This is due to the Spawn ability to provide compact sets of tasks, reducing data transfers between threads.

On the other side, if we compare task cache misses with Eager, DM and DMDA, we note that Eager doubles the number of cache misses that DM and DMDA induce, but has better performance (see the speedup in fig. 13). This is related to the behavior of these schedulers: by using a task queue, Eager naturally maximizes StarPU workers activity, in opposition to DM and DMDA that distribute tasks in advance as soon as they are ready in order to minimize overall execution time.

Since our application uses intensive memory accesses, the execution time is strongly related to memory accesses, and task execution time depends on the place where tasks are executed. Therefore, assigned tasks by DM or DMDA constitute a set of non-related tasks, like they would have been assigned with Eager. Hence, with DM or DMDA, StarPU workers are imbalanced, since the execution time of assigned tasks doesn't reflect their real execution time. A deeper analysis of StarPu execution traces (thanks to the StarPU support of FxT [START_REF]Fxt library for execution traces generation[END_REF] traces generation) (see fig 14) confirms this point. With this information, it appears that common tasks schedulers seem unable to perform well with such task configurations (neighboring and data dependencies).

DISCUSSION

Our force equation (see eq. 5) defines the optimal vCore load m as the sum of all vCores load divided by their number. This centralized characteristic implies distinct considerations depending on whether we are in a shared or in a distributed memory environment. In a shared memory environment, we need to store task costs. Depending on the size of the domain and the number of used vCores, it can imply an important memory cost. In a distributed memory environment, this implies to retrieve all the task costs information on one node. Beside the fact that this requires one ALL-to-ONE communication each time we do a partitioning, this also necessitates a lot of memory into one node to store all the information.

Despite these considerations, results presented in sections 4.1.2 and 4.2 show us interesting performance. An improved version of our algorithms, with less centralized properties, should scale well on multiple nodes.

CONCLUSION & FUTURE WORK

This Pair-Potential approach for task scheduling shows good scheduling properties. Data locality and dynamic load variations are well handled, for a small computing cost. In a shared memory, or with a few number of MPI nodes, experiments show that this method is really comparable and outperforms well known graph partitioners and common task schedulers. Nevertheless, the current implementation has some limitations. Central characteristics of the current algorithm involves a memory consumption and communications that increase with the number of computing units (cores, MPI nodes). Therefore, the current implementation isn't well scalable on a large number of MPI nodes. Our next work will thus focus on a parallel version of the algorithm, for both shared and distributed memory, with a distributed information, and an inter-nodes force computation.

Figure 1 :

 1 Figure 1: 2D Mesh example. The red cell has 8 direct neighbors, 16 2 nd rank neighbors and 32 3 rd rank neighbors. Every cell contains a set of elementary elements (atoms, finite elements or finite difference cells).

Figure 2 :

 2 Figure 2: 2D representation of a cells distribution across computing units with a Voronoi tessellation. Every cell represents a task or a graph node and every sphere represents a computing unit. The color of the cell shows the computing unit to which the task is attached to.

Figure 3 :

 3 Figure 3: Force interactions between vCores. To compute forces acting on the cell 0, Q charges of neighboring domains are needed. MD terminology Analogy Description Potential Potential Defines forces, and thus system behavior. Particle vCore A computing unit. Particle charge vCore load The load of a computing unit.

Figure 4 :

 4 Figure4: System convergence with oscillations and with oscillation corrections versus the number of iterations for a typical system (50x50 tasks with 64 vCores). This plot of the variance of the system shows instabilities when we are close to an optimal configuration (yellow curves). Red curve represents the same situation, but with an oscillation corrections.

Figure 5 :

 5 Figure 5: Voronoi boundary cell cost smoothing as a function of distance.

Figure 6 :

 6 Figure 6: Inital random Perlin noise (top-left), cells distribution with Metis (top-right), cells distribution with Scotch (bottom-left) and cells distribution with Spawn (bottom-right).

Figure 7 :

 7 Figure 7: Efficiency & quality domain partitioning comparison with a random load distribution (topleft), after a load variation (top-right) and after a load variation plus a modification of the P parameter (bottom-left). Two graph partitioning tools are used (Metis and Scotch) and compared with Spawn. These graphics show the time needed to compute the partition of a 256x256 domain of tasks, with a variable number of vCores (or parts).

Figure 8 :

 8 Figure 8: Time needed to reach a given vCores imbalance, within at most 1,000 iterations, for different domain sizes. 0% of imbalance means perfect task distribution while a 100% of imbalance means that we can have one cell that has at most twice the optimal load.

Figure 9 :

 9 Figure 9: Resulting speedups and L1 data cache misses with partitioning tools (Scotch, Metis) and with Spawn on an IvyBridge Xeon E5-2680 computing node.

Figure 10 :

 10 Figure 10: Resulting speedups and data transfer rates with partitioning tools (Scotch, Metis) and with Spawn.

Figure 11 :

 11 Figure 11: Resulting speedups and data transfer rates with Zoltan (RCB, PHG), ParMetis (through Zoltan) and with Spawn.

Figure 12 :

 12 Figure 12: Timeline representing StarPU calls (gray), task computation (green) and partitioning computation with Spawn (red). Partitioning computation is done only when an imbalance is present and continues as long as a StarPU worker has ongoing tasks.

Figure 13 comparesFigure 13 :

 1313 Figure 13 compares performance of Spawn in relation to internal StarPU schedulers. Left figures refer to the static load case while right figures refer to the dynamic load case.In the same way, upper charts refer to the speedup while the bottom ones refer to the ratio of L1 data cache misses per thread for one iteration in comparison with Spawn. We use three internal StarPU schedulers: Eager, DM and DMDA. Eager uses a simple FIFO-based greedy policy to schedule tasks over computing units. DM and DMDA use performance models allowing them to perform smart task placement over StarPU workers, with the objective to minimize the overall execution time. DMDA works like DM, but takes into account

Figure 14 :

 14 Figure 14: FxT traces with the DM scheduler. Green elements refer to tasks execution while red elements refer to inactivity of StarPU workers. One row symbolise the activity of one StarPU worker.