Jens Gustedt

Maxime Mogé

Memory access classification for vertical task parallelism

Keywords: automatic parallelization, iterative algorithms, parallel control flow, memory access classification

We present a paradigm and implementation of a parallel control flow model for algorithmic patterns of two nested loops; an outer iteration loop and an inner data traversal loop. It is centered around memory access patterns. Other than dataflow programming it emphasizes on upholding the sequential modification order of each data object. As a consequence the visible side effects on any object can be guaranteed to be identical to a sequential execution. Thus the set of optimizations that are performed are compatible with C's abstract state machine and compilers could perform them, in principle, automatically and unobserved. We present two separate implementations of this model. The first in C++ uses overloading of the operator[] to instrument the memory accesses. The second in Modular C uses annotations and code transformations for the two nested loops. Thereby the code inside the loops may stay as close as possible to the original code such that optimization of that code is not impacted unnecessarily. These implementations show promising results for appropriate benchmarks from polybench and rodinia.

Classification d'accès mémoire pour la parallélisation verticale en tâches

Résumé : Nous présentons un paradigme et une implémentation d'un modèle de flot de contrôle parallèle pour des motifs algorithmiques de deux nids de boucles imbriqués; une boucle externe d'itération et une boucle interne de parcours de données. Ce paradigme est basé sur les motifs d'accès mémoire. Contrairement à la programmation par flot de données, il assure l'ordre séquentiel des modifications de chaque objet. Par conséquent, les effets de bords visibles sur tous les objets sont garantis identiques à une exécution séquentielle. Ainsi l'ensemble des optimisations effectuées sont compatibles avec la machine abstraite d'état de C, et les compilateurs peuvent les effectuer, en principe, de façon automatique et non-observée. Nous présentons deux implémentations distinctes de ce modèle. La premiere, effectuée en C++, utilise la surcharge de l'opérateur [] pour instrumenter les accès mémoire. La seconde, effectuée en Modular C, utilise des annotations et transformations de code pour les deux boucles imbriquées. Ainsi le code à l'intérieur des boucles peut rester aussi proche que possible de l'original pour ne pas impacter inutilement son optimisation. Ces implémentations montrent des résultats prometteurs pour des codes de référence tirés de polybench et rodinia.

Mots-clés : parallélisation automatique, algorithmes itératifs, flot de contrôle parallèle, classification d'accès mémoire

Introduction and overview

Race conditions are one of the principle challenges of automatic parallelization, and over last decades a lot of effort has been put into understanding and mastering them, in particular into their detection in faulty executions, see e.g [START_REF] Banerjee | A theory of data race tection[END_REF][START_REF] Matsubara | Model checking with program slicing based on variable dependence graph[END_REF][START_REF] Wang | Data race detection based on dependence analysis[END_REF]. Another direction of research has been to attempt to provide a race free task schedule once that a dependency graph is given [START_REF] Sánchez | On efficient distributed deadlock avoidance for real-time and embedded systems[END_REF]. But generally, the need for a programmer to manually specify a dependency graph, even in the form a data flow graph is a burden that defeats automatic parallelization.

In our opinion, to be closer to current programming practice, detection of dependencies and thus possible parallelism should be based on the features that are present in every day's imperative programming languages. In these, dependency between program statements are usually not explicit but they are forced indirectly via data accesses when one statement that reads a data has to be scheduled after another statement that wrote it.

Therefore parallelization of loops with irregular access patterns requires a runtime data dependency analysis that is followed by a tool that ensures that dependencies are enforced. There are currently mainly two types of tools: speculative and scheduling based.

Speculative techniques such as LRPD [START_REF] Dang | The r-lrpd test: speculative parallelization of partially parallel loops[END_REF] identify the dependencies during the execution of the loop. The loop is run in parallel as a doall and followed by a test if the computations are correct. If not, a rollback to the previous correct state is performed and a sequential execution of the loop is carried out. For certain types of data accesses, namely affine or nearly affine functions over loop counts, polyhedral models can be used to provide race-free static task parallelism [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF], and these approaches have also recently be used successfully to build runtime systems that are able to rollback execution if access violations occured [START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF]. For programs with low effective parallelism, this can lead to a significant overhead due to rollback followed by sequential execution when many conflicts occur between threads.

Scheduling based methods preprocess the loop to compute a dependency graph at runtime and use this information to do a parallel scheduling. The Inspector-Executor model [START_REF] Saltz | Run-time parallelization and scheduling of loops[END_REF][START_REF] Venkat | Automating wavefront parallelization for sparse matrix computations[END_REF] consists of three phases: inspection, scheduling and execution. First the program is instrumented to explicitely compute a dependency graph. Then a parallel scheduling of the iterations is derived, and then the iterations are run in parallel in wavefronts, using synchronization between consecutive wavefronts. Huang et al. [START_REF] August | Automatically exploiting cross-invocation parallelism using runtime information[END_REF] propose a scheduling based method that does not explicitely compute the dependency graph, but instead overlaps the inspection and scheduling. A dedicated scheduling thread dynamically ensures that there is no conflict between threads, and allows cross invocation parallelism. However this method still needs an inspector that is run at the beginning of each iteration. The limitations of these methods are inherent to the Inspector-Executor model: it needs one inspector per iteration. In the case of a cyclic dependency between data and computation of shared addresses the inspector is basically the whole loop body.

In contrast to that, we will not restrict ourselves to a specific fine grained access projection model, but instead focus on certain type of applications, namely those that repeat the same data access pattern (coined data traversal) over a set of iterations, Section 2. This is a pattern that is found widely in the field, e.g many iterative algorithms that traverse matrices or geometric objects fall into that category. Using the assumption of a constant data access pattern, we detect dependencies and derive an implicit scheduling at runtime using Ordered Read-Write Locks (ORWL) [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF] during some initial iterations. After this instrumentation phase, we do not need to check any condition or recompute dependencies, thus eliminating the drawbacks of the Inspector-Executor model.

As a basis for our argumentation, here, we will use C's model of side effects on data [START_REF]Programming languages -C[END_REF]. C has an abstract state machine that can be used to describe the effects of any valid program and compilers are allowed to perform all optimizations that respect that machine model. Our present work emphasizes on the fact that our parallelization fully respects the computation in that machine model. A fine grained (theoretical) model is presented that guarantees that all computational results are exactly as they would have been produced by the originating sequential program, Section 3. Based on [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF] we are able to prove that our parallel execution model is correct, fair and deadlock free.

Being much too fine-grained for practical utilization, we have to coarsen our model by grouping programming steps into "meta-steps" and by classifying objects into "meta-objects", Section 4. We are able to prove that by doing so the good properties of our model are maintained, and present two different strategies that can be used for a memory access classification.

Our approach has been implemented with two different approaches, Section 5. The first uses C++ and its ability to overload the operator[]. By that we are able to dynamically instrument the access pattern of complex code without having to rely on the programmer. A second implementation in Modular C [START_REF] Gustedt | Modular C[END_REF] provides #pragma annotations for a code transformation that shows to be more efficient, but that needs manual annotations of the data accesses. Both implementations are tested with a set of benchmarks from the polybench and the rodinia collections, Section 6, and show very satisfying speedups.

Iterations over data transversals

In this work we restrict ourselves to a specific framework, namely programs that are dominated by two levels of loops, see Fig. 1.

1. An outer loop that we call iteration, e.g, an iteration over a simulated time or an iterative approximation of some numerical quantity. 2. One or several inner loops that we call data traversal tasks or just tasks, e.g loops that visit all elements of a matrix or that iterate over a geometric domain such as the facets of a polyhedral object description.

Our main assumption for data accesses is the following:

The data access pattern does not depend on the iteration. That is, we assume that each outer iteration visits the same data in the same order as the previous ones, e.g a traversal of a matrix or of geometry elements would be done in the same order in each iteration. That does not mean that a task can't use the iteration variable for its computation, only that its data access should not depend on it.

The inner part of a data traversal task t is called a step, denoted by T t (i, p), characterized by the triple (i, t, p), where i is the actual iteration, t is the task and p is the actual position in the traversal. That is, we assume that the program identifies code segments that perform these steps, and that the information about the current iteration, task and position in the traversal is available.

Steps are sequenced, that is there is a linear order in which the steps are executed by the program. This order corresponds to the lexicographic ordering of the indices (i, t, p).

Listing 1 illustrates what we have in mind here. It shows an idealized code excerpt of rodinia's hotspot3D benchmark, where we use a split directive (pragma) and a special task ∶∶ for construct to identify a data traversal task. A step is then the entire code inside the {}. It has access to the iteration variable i and the traversal position z.

Another feature of that implementation is also shown, namely data accesses are instrumented explicitly by means of directives that proclaim that the step may modify (lvalue) or just read (rvalue) a specific data.

Execution model: loop exchange and parallelism

In general, we suppose that a sequential program P has the following form:

for (all iterations i) for (all tasks t) for (all positions p in T t data traversal) T t (i, p);

Our goal is to improve the performance of P by parallelizing it. The main idea is to automatically perform a loop exchangetransformation as exemplified in Fig. 1. After the exchange the outer loops are in fact over the positions p of the data traversals and (at least conceptually) an independent thread is launched for each such position. Inside such a thread, the iterations are performed over the individual steps that are all associated to the same position p.

As remarked in the comments, data races could occur when two inner steps would access the same data concurrently. Therefore we will have to develop a mechanism that helps us to detect when one step has to wait for another one because of such a concurrent access.

As already implicit in the hotspot3D example, our notion of dependency of steps is not an abstract model of data or control flow, but directly deduced from access to objects of the program. It directly follows C's notion of side effects [START_REF]Programming languages -C[END_REF], that is, of modifications effected to data.

Step S 1 directly depends on step S 0 if -S 1 is sequenced after S 0 and -S 1 reads or writes a data object o that S 0 has written.

By that, direct dependency between two steps can never be induced by only read accesses to a common data.

Step S n (indirectly) depends on step S 0 if there is a direct dependency chain between them, that is, if there is a, possibly empty, sequence S n-1 , . . . , S 1 such that S j directly depends on S j-1 for all 0 < j ≤ n.

Dependence in that sense establishes a partial ordering ℘ the dependency relation on the steps that are executed by a program. Two steps S 0 and S 1 are independent if neither of them depends on the other and if so, they can be executed concurrently without introducing a race condition. Observe that two independent steps may be from the same or different tasks, iterations and traversal positions.

The notions of dependence and independence easily extend to sets of steps S 0 and S 1 . In fact, S 0 depends on S 1 if there is a S 0 ∈ S 0 that depends on S 1 ∈ S 1 ; S 0 and S 1 are independent if all pairs S 0 ∈ S 0 and S 1 ∈ S 1 are independent. Observe Inria that here the dependency relation between sets of steps may have cycles if we are not careful enough in choosing such sets.

Parallelizations

Fig. 2 shows an example of the initial section of a direct dependency graph and different scheduling strategies that have been applied. In all four examples, lines corresponds to concurrent executions of steps at a given time. Columns visualize all steps that deal with the same position in their respective data traversal. For the sake of the example, we make the simplified assumption that all steps have the same execution time. Sequential execution. A sequential scheduling of the graph is shown in Fig. 2(a). Here, all steps are strictly executed in iteration and traversal order. E.g first we schedule all steps T 0 (0, .) of task T 0 . Among these we respect the traversal order, that is we execute first T 0 (0, 0), then T 0 (0, 1), T 0 (0, 2) and T 0 (0, 3). The rectangle in the figure shows the pattern of the steady state, that is the pattern of the graph that is repeated over and over again until the final iteration. Horizontal parallel execution is similar to what OpenMp would do with a for loop that is prefixed with a parallel for directive. But without keeping track of dependencies, OpenMp would allow all steps from the same task to run in parallel, and thus lead to incorrect results.

In fact, introducing horizontal parallelism has several problems that are not easily solved with static scheduling how about horizontal data dependencies? -how about cache accesses? -how about execution time variations? Fig. 2(b) shows that our example only allows for a limited horizontal parallelization that respects dependencies: the small rectangles inside the steady state visualize the steps from the same task that can be parallelized. We see that we can place our 8 steps in 6 time units, so on average we have 8 6 = 1.3 parallel threads in the steady state. Vertical parallel execution is illustrated in Fig. 2(c). Here, we start the steps as early as possible under the constraint to still execute the steps of the same task and iteration in order. In our example, we see that at time 2, step T 1 (0, 0) already starts in parallel to step T 0 (0, 2). Generally, for all i > 0, the four types of sets of steps

T 0 (i, 0), T 1 (i -1, 2) T 0 (i, 1), T 1 (i -1, 3) T 0 (i, 2), T 1 (i, 0) T 0 (i, 3), T 1 (i, 1)
are always scheduled in parallel. Thus there are 2 parallel threads running at all times in the steady state. Vertical and horizontal parallel execution as a combination of both parallel modes, shows the best available parallelism, see Fig. 2(d). Here it not only allows to run steps of the same iteration in parallel (e.g T 1 (0, 0) and T 1 (0, 1)) but also steps of the same task but from different iterations, e.g T 0 (0, 3) and T 0 (1, 0). Or time 0 T 0 (0, 0)

5 5 1 T 0 (0, 1) Ø Ø 5 5 2 T 0 (0, 2) 5 5 3 T 0 (0, 3) Ø Ø 4 T 1 (0, 0) 5 T 1 (0, 1) Ø Ø 5 5 6 T 1 (0, 2) 7 T 1 (0, 3) Ø Ø 8 T 0 (1, 0) 5 5 9 T 0 (1, 1) Ø Ø 5 5 10 T 0 (1, 2) 5 5 11 T 0 (1, 3) Ø Ø 12 T 1 (1, 0) 13 T 1 (1, 1) Ù Ù 5 5 14 T 1 (1, 2) 15 T 1 (1, 3) (a) sequential execution time 0 T 0 (0, 0) 5 5 1 T 0 (0, 1) Ø Ø 5 5 2 T 0 (0, 2) 5 5 3 T 0 (0, 3) Õ Õ 4 T 1 (0, 0) T 1 (0, 1) Õ Õ 5 5 5 T 1 (0, 2) T 1 (0, 3) Ø Ø 6 T 0 (1, 0) 5 5 7 T 0 (1, 1) Ø Ø 5 5 8 T 0 (1, 2) 5 5 9 T 0 (1, 3) Õ Õ 10 T 1 (1, 0) T 1 (1, 1) Ö Ö 5 5 11 T 1 (1, 2) T 1 (1, 3)
(b) horizontal parallelism: if possible steps of the same task are executed in parallel. On average, 8 6 = 1.3 parallel threads are in the steady state. time 0 T 0 (0, 0)

5 5 1 T 0 (0, 1) { { 5 5 2 T 1 (0, 0) T 0 (0, 2) 5 5 3 T 1 (0, 1) { { 5 5 T 0 (0, 3) { { 4 T 0 (1, 0) 5 5 T 1 (0, 2) 5 T 0 (1, 1) { { 5 5 T 1 (0, 3) { { 6 T 1 (1, 0) T 0 (1, 2) 5 5 7 T 1 (1, 1) 1 1
T 0 (1, 3) (c) vertical parallelism: task starts executing as soon conflicts with preceding tasks are resolved. There are 2 parallel threads in steady state.

time 0 T 0 (0, 0) 5 5 1 T 0 (0, 1) { { 5 5 2 T 1 (0, 0) T 1 (0, 1) { { % % T 0 (0, 2) 5 5 3 T 0 (1, 0) 5 5
T 0 (0, 3)

{ { 4 T 0 (1, 1) { { 5 5 T 1 (0, 2) T 1 (0, 3) { { 5 T 1 (1, 0) T 1 (1, 1) 1 1 T 0 (1, 2) 1 1
(d) vertical and horizontal parallelism: steps are scheduled with first-fit. On average, 8 3 = 2.6 parallel threads are in steady state.

A fine grained execution model

We aim to develop an execution model that is compatible with C's abstract machine in a way that the parallelizations that we propose do not change the observable state of a program.

Observable changes in the states of that abstract machine occur through socalled side effects. By neglecting IO and similar external events for the moment, the interesting side effects in our context are modifications of the data objects that a program manipulates. Such manipulations occur through operations (such as assignment = or increment ++) or C library functions (such as memcpy).

For each object o in an execution of a program the C standard requires that accesses to o are properly sequenced, that is, that it can be deduced from the program structure if a write operation A on o provides the value for a read operation B of o; -two write operation A and B on o are not separated by another write operation.

Read-only operations in turn are not necessarily expected to be properly sequenced. E.g the two evaluations of i in a function call f (i, i) may be executed in any order. For our context here, we can assume that an access A to object o is characterized by a triple (i, t, p) for step T t (i, p), and, if it is a write access, by the new value that is stored in o at the end of step T t (i, p). In particular, we don't distinguish if a step T t (i, p) does several accesses to the same object o and whether read and write access occur inside the same step.

The modification order M (o) = M 0 , M 1 , . . . is the precise list of changes that are applied to an object o. The access order L(o) is M (o) together with the read accesses to o, that is where we insert the read accesses in sequence order between write accesses:

L(o) = W 0 , R 0 0 , R 1 0 , . . . , W 1 , R 0 0 , R 1 0 , . . .
Since we are assuming that memory accesses do not depend on the iteration, L(o) is of the form

L(o) = L(o, 0), L(o, 1), . . .
where L(o, x) is a predetermined list (x, t 0 , p 0 ), (x, t 1 , p 1 ), . . . where we textually substitute x by i for each iteration i.

By changing the perspective from data objects to steps, for each step T t (i, p) we can now identify the objects o = o 0 , o 1 , . . . that it accesses as well as the list position (i, t, p, o) in L(o, i) where this access occurs.

Fifo execution

The lists that we defined in the previous section can now be used to establish a parallel execution model for the program steps, that takes care of the coherence requirements of Fig. 1.

For each object o we establish a FIFO data structure F(o) that follows L(o) and for which we define the following operations:

init(F(o)) initializes the FIFO position to the start of L(o).

release(F(o), i, t, p) announces that iteration i of task t at position p has ended.

acquire(F(o), i, t, p) blocks step T t (i, p)

• if T t (i, p) is a write access until all previous accesses to position (i, t, p) in L(o) have been released; • if T t (i, p) is a read-only access until the last write access previous to position (i, t, p) in L(o) has been released.

Observe that the conditions for acquire ensure that each write access can be performed exclusively, and that all read-only accesses that follow the same write access can be honored concurrently. Now we can augment each step of the program such that it is protected from race conditions by acquiring all FIFOs for all accessed objects and by releasing them thereafter. A protected step Prot⟪T t ⟫(i, p) derived from T t (i, p) is:

Prot⟪T t ⟫(i,p) { for (all objects o accessed by T t (i, p)) acquire(F(o), i, p); T t (i,p); for (all objects o accessed by T t (i, p)) release(F(o), i, p); } Clearly, substituting Prot⟪T t ⟫ for T t for all t in the whole program does not change its visible side effects: all acquire operation would only be issued after all previous positions in L(o) are released, and all operations T t themselves would still appear in sequence order.

We can now proceed to our final goal, the loop exchange. A step iteration

StepIt⟪T t ⟫ performs the iteration loop

StepIt⟪T t ⟫(p) { for (i in iteration) Prot⟪T t ⟫(i,p); } Now our parallelized program can be formulated as for (all tasks t) do parallel { for (all positions p in T t data traversal) do parallel StepIt⟪T t ⟫(p); } Theorem 1. The parallelized program is well defined and produces the same side effects as the original program.

Inria

A proof of this theorem follows from the fact that the dependency relation ℘ is acyclic and from the theory developed in [START_REF] Clauss | Iterative Computations with Ordered Read-Write Locks[END_REF] for Ordered Read-Write Locks (ORWL). This formalizes a model for autonomous execution of tasks or processes that manage their shared resources by means of cyclic FIFOs that are attached to data resources. The ORWL model here can guarantee correctness, equity and liveness of execution. It has a reference implementation that has first been described in [START_REF] Gustedt | Relaxed Synchronization with Ordered Read-Write Locks[END_REF], and that has proven to be very efficient for coarse grained parallelism.

Expected efficiency: Generally for the framework that we described so far, we are unfortunately not in a coarse grained setting. We have made no restriction about the code that is found inside a step, and for many sensible programs that we want to cover a step could just be several instructions. In particular, we have to resolve the following problems:

-Data traversals can have thousands, millions or billions of steps. Thus our parallel program would have to run as many threads or processes, exceeding the possible number of threads for current architectures. -A protected step Prot⟪T t ⟫ introduces runtime overhead for each memory access, namely calls to acquire and release the corresponding FIFO. Therefore we'd have to expect the cost of a step at least to triple by that procedure. -Additionally, we also would have to maintain a FIFO for each data object.

So the memory footprint of the application would at least double. -There is an important overhead to initiate FIFOs for all objects.

Coarsening

To be useful as such, our step-model of is too fine grained. We have to coarsen it with respect to two aspects. First, we have to ensure that the number of steps (and thus the number of threads) does not exceed the usual limits of existing platforms. Second, we have to avoid that each data object gives rise to its own FIFO, but must regroup objects into classes such that the three types of overhead that we describe above are limited.

Groups of steps

To reduce the number of steps, we transform the original program by grouping steps:

for (i in iteration) for (t in tasks) for (group g of positions in T t traversal) for (p in g) T t (iteration, p);

The easiest strategy to form such groups of positions is to divide the positions p 0 , p 1 , . . . , p r in the traversal into a fixed a amount m of intervals of successive positions.

g k = (p a k , . . . , p b k ) , with a k = kr m , b k = (k + 1)r m .
This allows to rewrite a group g of steps into a "meta-step" MetaStep m ⟪T t ⟫: As a consequence this rewrite of the program fits exactly the same model as the original program; it consists of an iteration loop that is nested with data traversal loops such that the data access pattern does not change with the iterations. In particular we have Lemma 1. The induced dependency relation ℘ m between meta-steps is acyclic.

With that grouping strategy the total amount of threads can be limited as a function of m, the number of groups that we fix per task.

Classes of objects

To bundle the management of FIFOs for several data objects at a time we partition the set of data elements

D into m ′ meta-objects C = {c k ⊂ D}. D = ⋃ k<m ′ c k
Such a classification of objects into meta-objects gives rise to an induced dependency relation ℘ m C between meta-steps in a obvious way. We have Lemma 2. ℘ m C is an partial order extension of ℘ m .

and in particular

Corollary 1. ℘ m C is an acyclic, irreflexive and transitive relation.

It is also easy to see that the access list L(c) of a meta-object c ∈ C is the lexicographic merge of all lists L(o), with o ∈ c. That is, we consider that an meta-object is accessed if any of its members is accessed. We then attribute one FIFO F(c) to each meta-object c ∈ C and define the protected meta-step analogously as before.

MetaProt⟪ MetaStep⟪T t ⟫⟫(i,g) { for (c ∈ C touched by MetaStep⟪T t ⟫(i,g)) acquire(F(c), i, g); MetaStepMetaStep⟪T t ⟫(i,p); for (c ∈ C touched by MetaStep⟪T t ⟫(i,g)) release(F(c), i, g); }
Clearly, using the resulting scheduling is more restrictive, but, because of Corollary 1 its validity is not impacted.

Theorem 2. The parallelized program that is coarsened with respect to a step grouping of m and an object classifier C is well defined and produces the same side effects as the original program.

Inria

To assemble the meta-objects we investigated two strategies. Objects such as vectors or matrices for which sub-objects occupy contiguous memory can be classified in a similar way as we grouped steps by using access intervals. If this is not the case or if the access pattern is too irregular, we have to keep track of individual accesses.

Classify groups of contiguous memory accesses in ranges. Suppose the main data access is through a vector V[j] for j = 0, ..., N -1 and that the access to that vector is regular such that the progression of the data traversal corresponds to a linear access to the vector.

To classify the access to V we can then fix a chunk size M and define the meta-objects as

contig(V[k]) = ⌊k M ⌋ Contig(j) = V[k] contig(V[k]) = j
The advantage of Contig(.) is that the membership for any data object to its class is easily computed. Therefore the overhead that is needed to establish the FIFOs for the meta-steps is relatively small. Classify random memory access by ownership. Vector classification does not work if data access is through lists, e.g if we do a data traversal over a complicated data structure or if our access pattern follows some random order. In these cases for each object o we can designate the owner (owning step) and use that owner to define a partition of the set of objects: owner(o) = first mega-step g that writes to o in the sequence ordering

Owner(g) = o o object, owner(o) = g
We may prepare the computation of Owner(g) by first running several sequential iterations:

-In a first counting phase we determine the number of objects that will be accessed by all mega-steps and dimensions a hash table with the appropriate size. -A second owning phase notes the first write access for all objects in the hash table.

-A third attribution phase determines the owner of each object. Now the data objects are partitionned by their access in the meta-step. We have at most as many meta-objects as there are meta-steps.

After this initial setup of the meta-objects, the remaining iterations can be done in parallel.

Compared to Contig, Owner classification has a startup that is is more expensive. On the other hand, no assumption is made about access order or about a specific organization of the data.

Implementations

We have implemented the Owner classifier with two different approaches. The first uses C++ and overloading of the operator[] and the second uses Modular C [START_REF] Gustedt | Modular C[END_REF] and #pragma directives for a source-to-source transformation.

Both build on the EiLck library for the FIFO functionality. EiLck is a standalone library written in Modular C for Linux. It implements cyclic access FIFOs on the basis of Linux' futex system call that efficiently combines atomic access with a scheduling queue in the system kernel, [START_REF] Franke | Fuss, futexes and furwocks: Fast userlevel locking in Linux[END_REF]. All scheduling of tasks is done implicitly via these FIFOs; EiLck executes the application without dedicated management threads. EiLck is compatible with different thread libraries (C11, C++11, POSIX, OpenMp) and is interfaced to C and C++. It efficiently avoids overhead for busy waiting if a FIFO access is congested and provides a minimal run time overhead in case it is not.

C++: Operator overloading and explicit thread creation

Besides the use of Eilck, our C++ implementation is based on standard tools, mainly operator[] and the explicit creation of threads for the execution of the tasks.

operator[] is used to instrument every single memory access to a set of vectors that the programmer designates as being shared between steps. In addition, the programmer also identifies the iteration loop and registers the data traversal tasks as functions that are executed in separate C++11 threads.

Then, during execution, data traversals are grouped and protected automatically. For the grouping of the steps into meta-steps, each task counts its number of calls to operator[] and uses this count as an indication of the traversal position. For the classification of the memory accesses we use a dynamic std::map to register the owners and then classify each access according to that owner. A set of FIFOs is constructed from this, and used in subsequent iterations.

This approach has the advantage that we don't even have to identify steps. Any code that supplies sufficient amount of work can automatically be split by this operator[] approach into steps and meta-steps, respectively. The programmer only needs to designate the iteration loop, a set of tasks (data traversal) and the data vectors they share. The modifications needed for this parallelization are simple, as shown in Listing 2 for the hotspot3D benchmark. The approach has several disadvantages, though:

-All meta-steps of a task are executed in the same thread. Parallelism is only achieved between meta-steps that are part of different tasks. -Control operations by operator[] interleave the computation. This breaks local optimization opportunities, e.g vectorization. -The operator[] needs a thread local variable for the position in the traversal.Access to such variables is generally slow.

Inria

Listing 2. Parallelization of rodinia's hotspot3D benchmark using operator overloading -the code was simplified for clarity (e.g. omission of template parameters). is for a construct that precedes some statement or block and potentially changes its interpretation. The second, amend, #pragma CMOD amend mctask keyword = arguments ... code snippet to be rewritten ... #pragma CMOD done can be used to rewrite an entire code snippet as necessary. Here, both directives refer to an external script, mctask, that implements our parallelization. The main features shown in that example are the following -An iterate directive replaces the outer iteration loop.

-A split directive prefixes the data traversal.

-A task ∶∶ for (or task ∶∶ while) loop specifies the data traversal.

-Data access directives lvalue and rvalue indicate all memory accesses that could be subject to race conditions.

The iterate directive uses amend to transform the code of the iteration. It creates two modified copies of the code, one for an initial instrumentation phase (similar to the C++ version) and one for the steady state iterations. For the first phase, lvalue and rvalue are replaced by calls to appropriate access classification functions, for the second they are removed such that the steady-state iterations can use the original code directly.

The split directive together with the task ∶∶ for ensures that the loop is interpreted as a single task. It groups the steps into meta-steps and ensures that the protection against race conditions is inserted between them. But in contrast to a steps directive, see below, it still only launches one thread per task, and we only will see vertical parallelism if we only use split.

The transformed parallel program can be quite effective because the inner data traversal is exactly as programmed originally. In particular, the program maintains the same optimization opportunities, e.g for vectorization. Also, compared to our C++ version, there is no need for thread local variables, the sourceto-source transformation has enough possibilities to provide the necessary information . On the other hand, this implementation here requires that the structure of the program makes the iteration and data traversal loops apparent. E.g hiding the data traversal inside a function would not be possible.

Listing 3 also features two other directives that are be quite useful. With the innermost amend, the alternate directive implements an automatic duplication that alternates the roles of tOut and tIn. The next level of amend with a duplicate directive duplicates its inner part, again, such that at the end we have four copies of the inner part, giving rise to four different tasks.

In total we have 4×8 = 32 meta-steps. Here, because of the split directive, we have one thread for each task, so 4 tasks in total that are parallelized vertically. To gain also horizontal parallelism we can replace The steps directive creates a separate thread for each meta-step, so here this results in 32 threads in total.

Benchmarks

We tested our implementations with those benchmarks from the polybench and the rodinia series that fulfill our requirements, and inside a real life application, the SOFA framework.

The test platform is a Linux machine with an Intel Xeon E5-2650 v3 processor, 2 cpu sockets, 10 real cores (and 20 hyperthreaded cores) each. Although we did not expect much parallelism beyond a speedup factor of 2 or 3, we chose such a platform with relatively many cores to avoid an artificial bound on the parallelism. For the C++ implementation, we adapted the code using C++ vectors.

Fig. 3 shows the speedups that we obtained with these relatively simple tests, once we have reached the steady state. We compare the original sequential implementation seq to our implementations using Modular C with split and steps directives, and C++. For some tests we unroll the iteration loop (2x or 4x factor) to increase the number of tasks. With Modular C, we see that all but one of the benchmarks achieve worthy speedup factors that range from 1.2 or 1.3 (srad2D), over 1.5 (fdtd), 1.7 (adi and heat3d), 2.2 (hotspot3D), to up to 4 for seidel_2d. Indeed, the latter parallelizes perfectly with the number of duplicates of the data traversal that we produce.

Most of these parallelizations are straight forward code annotations as we have seen them for the hotspot3D benchmark above. There is no clear winner between the split and steps directives; adding more threads is not always an advantage. The parallelization of srad2D is a little bit more involved than the others, because the original sequential code concentrates most of the computation in one of the traversals. So here a simple split even slows the execution down. But when we divide that traversal loop into several such that we obtain 4, 5 or 7 tasks in total, we observe a speedup of about 1.3.

Only the nearest neighbor (nn) benchmark from the rodinia suite detaches from this. It is completely dominated by the IO for reading the input data. The computation itself is negligible, so we can't gain much by parallelization. Nevertheless, we also don't lose anything, which shows the validity of our approach, where we are able to integrate the IO seamlessly into a data flow parallelism.

With the C++ version, most of the tested benchmarks achieve worthy speedup, too, although lower than with Modular C. The overhead induced by the overloading and the inhibited compiler optimizations can seriously impact the performances (e.g. only 1.1 speedup for heat3d and even a slowdown for adi). However, a simple unrolling of the iteration allows us to create more tasks and gives a much better speedup (1.7 for heat3d with duplicated tasks, and up to 2.9 with 4x unrolling).

The cost of the startup steps is highly dependent on the data objects considered and the application. With default settings, the Modular C version uses fine grained data objects (cache lines) with a cost ranging from 2 to 120 sequential iterations. This granularity can be manually adapted to reduce this cost. The C++ version, uses coarser data objects (rows of 2D matrices) with a cost of 3 sequential iterations for seidel2D, 12 for heat3d, 16 for hotspot3D.

We also tested our C++ implementation on a large scale application, the SOFA framework, with interesting performance: we were able to maintain the runtime overhead induced by our parallelization mechanism low enough to get a significant speedup on a specific test simulation.

Conclusion and outlook

We have shown that a semi-automatic parallelization tool for a special class of iterative algorithm can be implemented. Provided that the data access pattern does not change in the iteration loop, execution is guaranteed to follow the sequential modification order for each data. Thereby we can e.g guarantee that numerical iterations lead to exactly the same results as sequential execution and that all proofs for convergence translate directly to the parallel version. We showed that our implementations leads to satisfying speedups for a variety of benchmarks.

For the moment our implementations both need intervention from the programmer. Both need an identification of the iterations and data traversals. Also, they are limited in their ability to instrument the memory access. The C++ version has a certain runtime overhead for the overloaded operator[]. The Modular C version needs manual annotation of the accesses and has a high startup cost because of its use of fine grained data objects. We think that all these shortcomings Inria can be overcome by implementing the memory access instrumentation directly in the compiler.

The need of identification could be circumvented by using a framework such as Apollo [START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF] that allows to speculate on a particular pattern, and then to rollback if the expected pattern was not respected.

Listing 1 .

 1 annotation of a data traversal loop of rodinia's hotspot3D benchmark for (size_t i = 0; i < numiter; ++i) { # pragma CMOD insert mctask split task∶∶for(size_t z = 0; z < nz; z++) { size_t z0 = ... ; size_t z1 = ... ; / * insert some miraculous data coherence enforcement * / # pragma CMOD insert mctask lvalue = tOut[z] # pragma CMOD insert mctask rvalue = tIn[z], tIn[z0], tIn[z1] for(size_t y = 0; y < ny; y++) { size_t y0 = ... ; size_t y1 = ... ; for(size_t x = 0; x < nx; x++) { size_t x0 = ... ; size_t x1 = ... ; / * do the computation * / tOut[z][y][x] = tIn[z][y][x]*cc + tIn[z][y0][x]*cn + ... ; } } } // same data traversal, but inverting roles of tOut and tIn ... }

Fig. 2 .

 2 Fig. 2. Schedulings of the direct dependency graphs of steps. Arrows represent data dependencies, numbers on the left represent execution times of the steps in their line. Rectangles illustrate patterns of the steady state. Inria

#

  pragma CMOD insert mctask groups = 8 #pragma CMOD amend mctask iterate = numiter #pragma CMOD amend mctask duplicate = 2 // unroll for more parallelism #pragma CMOD amend mctask alternate tOut tIn // unroll to alternate buffers # pragma CMOD insert mctask split task∶∶for(size_t z = 0; z < nz; z++) { size_t z0 = ... ; size_t z1 = ... ; / * insert some miraculous data coherence enforcement * / # pragma CMOD insert mctask lvalue = tOut[z] # pragma CMOD insert mctask rvalue = tIn[z], tIn[z0], tIn[z1] for(size_t y = 0; y < ny; y++) { size_t y0 = ... ; size_t y1 = ... ; for(size_t x = 0; x < nx; x++) { size_t x0 = ... ; size_t x1 = ... ; / * do the computation * / tOut[z][y][x] = tIn[z][y][x]*cc + tIn[z][y0][x]*cn + ... ; } } } #pragma CMOD done #pragma CMOD done #pragma CMOD done 5.2 Modular C: pragmas and source-to-source transformation Listing 3 shows a complete parallelization of the hotspot3D benchmark using our Modular C implementation. Here we use two forms of Modular C directives. The first, insert, #pragma CMOD insert mctask keyword = arguments ... prefixed statement or block ...

#

  pragma CMOD insert mctask split task∶∶for(size_t z = 0; z < nz; z++) { ... } by Inria #pragma CMOD insert mctask steps = nz { size_t z = task∶∶step(); ... }

Fig. 3 .

 3 Fig. 3. Speedups of steady state iterations achieved for benchmarks from polybench and rodinia with the Modular C and the C++ implementations.

  Listing 3. The complete paralellization of rodinia's hotspot3D benchmark. alternate and duplicate unroll the iteration loop four times such that four tasks can be generated.

	void computeTask(VectorWrapper tOut_w, VectorWrapper tIn_w, ...) { // encapsulate computation in a separate function
	for(size_t z = 0; z < nz; z++) {
	size_t z0 = ... ; size_t z1 = ... ;
	for(size_t y = 0; y < ny; y++) {
	size_t y0 = ... ; size_t y1 = ... ;
	for(size_t x = 0; x < nx; x++) {
	size_t x0 = ... ; size_t x1 = ... ;
	tOut_w[z][y][x] = tIn_w[z][y][x]*cc + tIn_w[z][y0][x]*cn + ... ;
	} } }
	}
	VectorWrapper tOut_w(tOut), tIn_w(tIn); // overload operator[] using a vector wrapper
	IterationLoopIdentifier loopId = ...; // identify an iteration loop
	for (size_t i = 0; i < numiter; ++i) {
	beginLoopIteration(loopId);
	threadPool->submit(new TraversalTask(computeTask,tOut_w,tIn_w)); // start task 1
	threadPool->submit(new TraversalTask(computeTask,tIn_w,tOut_w)); // start task 2
	endLoopIteration(loopId);
	}

RR n°9182

rue du Jardin Botanique CS20101 54603 Villers-lès-Nancy Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Ide-Flore Kenmogne
email: ide-flore.kenmogne@inria.frvincent.drevelle

Vincent Drevelle

Eric Marchand
email: eric.marchand@irisa.fr

Cooperative Localization of Drones by using Interval Methods

Keywords:

Introduction

In this paper, we address the problem of cooperative pose estimation [START_REF] Spletzer | A Bounded Uncertainty Approach to Multi-Robot Localization[END_REF] in a group of N unmanned aerial vehicles (UAV), each equipped with a camera that sees landmarks with known positions. The UAVs communicate, exchange poses and measure distances with neighbours and a base station. Our aim is to compute the pose domain of each robot assuming the errors on measurements are bounded.

Single robot pose estimation

Each robot first estimates its pose domain r = (x, y, z, φ, θ, ψ) using camera and base distance constraints. To get the camera constraints, the perspective projection equation (Eq. ( 1), pinhole camera model) of a 3D world point w X in the camera frame represented in normalized coordinates x = ( c x, c y) is used (see [START_REF] Kenmogne | Imagebased UAV localization using Interval Methods[END_REF] in camera only case).

x = Π c T r r T w (r) w X (1)
with r T w the unknown transformation matrix between the world reference frame and a frame attached to the robot and c T r is the known rigid transformation between the camera and the robot frames. Equation (1) is applied for each visible landmark w X i (i ∈ 1..m) and the following constraints can be derived

C i :    ( c X i , c Y i , c Z i ) = c T r r T w (r) w X i c x i = c X i c Z i , c y i = c Y i c Z i , c x i ∈ [ c x i ], c y i ∈ [ c y i ], c Z i > 0.
(2)

The image/range-based pose estimation problem is then defined as a constraint satisfaction problem (CSP)

H :   r ∈ [r], {C i , i ∈ 1...m} C dist  
Where C dist is the additional distance constraint between the robot position p = (x, y, z) and the base station B used to get a tighter pose estimate

C dist : d = p -b 2 , d ∈ [d]
with b the known position of the base station. A robot R k computes a domain in a form of an outer subpaving S + r k , that contains all the feasible poses, using SIVIA [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF] to solve H.

Robots cooperation: data exchange

At each time step, once the pose domain S + r k is computed, the robot exchanges the bounding box of its position domain [p k ] = proj p S + r k , where is the bounding box operator, and proj p is the projection onto the position space. The position [p k ] is transmitted to all neighboring robots R j , j ∈ N (k), and the distances d k,j between R k and R j are simultaneously measured (with N (k) the neighbours of R k ).

At reception of information (position boxes [p j ] and bounded-error distances measurements [d k,j ]) from neighboring robots, R k tries to refine its actual pose domain, by propagating the new distance constraints between R k and each of its neighbours. A CSP is also built and SIVIA is used to refine the pose domain.

Experimental results

The proposed method has been tested with data acquired on Parrot AR-Drone2 UAV, with 5 landmarks represented by AprilTag markers. The image measurement error bounds are set to ±0.5 px and the range measurement error is assumed to be within ±5 cm.

Figure 1 :

 1 Figure 1: Pose domain for 4 robots.


Ide-Flore Kenmogne
email: ide-flore.kenmogne@inria.frvincent.drevelle

Vincent Drevelle

Eric Marchand
email: eric.marchand@irisa.fr

Quantifier elimination for robot positioning with landmarks of uncertain position

Keywords: Interval Analysis, Pose estimation, Quantified set-inversion

We consider a mobile robot equipped with a camera that moves in a 2D environment. We aim to compute a domain for the pose (position and orientation) r = (x, y, θ) in which we are sure the robot is situated [START_REF] Spletzer | A Bounded Uncertainty Approach to Multi-Robot Localization[END_REF][START_REF] Kenmogne | Image-based UAV localization using Interval Methods[END_REF]. The robot knows landmarks positions X w (x w , y w ) in the world frame, and measures their projection u in the image (represented in normalized coordinates, assuming known camera calibration parameters). The perspective projection of a world point in the camera frame is given as follows:

u = f (r, X w ) = 1 zc ((x w -x) cos(θ) -(y w -y) sin(θ)) z c = (x w -x) sin(θ) + (y w -y) cos(θ) > 0 (1)
Problem The problem amounts in determining the set of all poses that are consistent with 1-D image measurement and its 2-D world corresponding point which is known with a bounded error. Assuming the image measurement uncertainty is represented by the interval [u] and the landmarks positions are known to be inside the box [X w ], the pose domain is the solution of the following quantified set inversion [START_REF] Jaulin | Resolution of nonlinear interval problems using symbolic interval arithmetic[END_REF] problem:

S = {r ∈ R 2 × [0, 2π], ∃X w ∈ [X w ], f (r, X w ) ∈ [u], z c > 0} (2)
Since the expression of f in Eq. 1 is not well posed when z c gets close to 0, the constraints f (r, X w ) ∈ [u], z c > 0 can be re-written as two inequalities:

g : g 1 : (x w -x)(cos(θ) -u sin(θ)) -(y w -y)(sin(θ) + u cos(θ)) > 0 g 2 : (x w -x)(cos(θ) -u sin(θ)) -(y w -y)(sin(θ) + u cos(θ)) < 0
where u and u denote the lower bounds and the upper bound of [u]. Quantified set inversion with a branch and bound algorithm s.a. SIVIA requires bisecting also on the quantified variables (otherwise, parts of the domain may remain undetermined). We therefore need to eliminate the quantifier (∃) to avoid bisecting on the landmark position box [X w ].

Quantifier free equations We adopt a geometric approach to perform the elimination. Let C = {(x w , y w ), (x w , y w ), (x w , y w ), (x w , y w )}, the set of corners formed by the box [X w ]. The pose is found consistent if g 1 is verified by at least one corner from C, g 2 is also verified by at least one corner, and the box [X w ] is not behind the camera.

S = r ∈ R 2 × [0, 2π], Xw∈C g 1 (r, X w , [u]), Xw∈C g 2 (r, X w , [u]), ¬behind(r, X w )
The first two conditions can also be written as: max Xw∈C ((x w -x)(cos(θ) -u sin(θ)) -(y w -y)(sin(θ) + u cos(θ))) > 0 min Xw∈C ((x w -x)(cos(θ) -u sin(θ)) -(y w -y)(sin(θ) + u cos(θ))) < 0

The ¬behind constraint expressing the fact that the camera is front-looking is given by:

     cos α ≥ 0 =⇒ x < x w cos α ≤ 0 =⇒ x > x w sin α ≥ 0 =⇒ y < y w sin α ≤ 0 =⇒ y > y w ⇔      cos θ -u sin θ < 0 ∨ x < x w cos θ -u sin θ > 0 ∨ x > x w
sin θ + u cos θ > 0 ∨ y > y w sin θ + u cos θ < 0 ∨ y < y w with α the observed ray angle in the world frame:

Results An example with 3 uncertain landmarks is depicted in Fig. 1.

Figure 1 :

 1 Figure 1: Pose domain in X-Y plane without and with quantifier elimination.


Ide-Flore Kenmogne
email: ide-flore.kenmogne@inria.fr

Vincent Drevelle
email: vincent.drevelle@irisa.fr

Eric Marchand
email: eric.marchand@irisa.fr

Using constraint propagation for cooperative UAV localization from vision and ranging

Keywords: Intervals, Cooperative localization, Constraints propagation

. Then, through position boxes exchange, positions are cooperatively refined by constraint propagation in the group. Results are presented with real robot data, and show position accuracy improvement thanks to cooperation.

Introduction

In this paper, we consider the problem of cooperative localization [START_REF] Roumeliotis | Distributed multirobot localization[END_REF] in a group of N unmanned aerial vehicles (UAV). The robots are equipped with cameras, able to see landmarks of known positions. A communication and ranging system provides to each robot R k a means of exchanging data and measuring distances with its neighbours and a base station B (Fig. 1). The goal for each robot is to compute a domain for its pose (position and orientation), assuming bounded error measurements.

The paper is organized as follows: we first present how each robot is able to independently compute a domain for its pose using constraints from camera

R 1 R 2 B d 1 d 2 d 1,3 Landmarks R 3 d 3 d 1,2 d 2,3
Fig. 1: Cooperative localization with camera and range measurements measurements and distance to the base station. Then, in a second part, a cooperative localization method is introduced, in which neighbours positions and distances are introduced as additional constraints to tighten the pose domain of each robot, thanks to data exchange. We finally provide experimental results obtained with quadcopter mini-drones.

Vision-based pose computation

This section addresses pose estimation [START_REF] Marchand | Pose estimation for augmented reality: a hands-on survey[END_REF] from camera measurements. To compute the pose r = (x, y, z, φ, θ, ψ) of a UAV amount in estimating the transformation r T w between the world reference frame and a reference frame attached to the robot. We assume that the rigid transformation c T r between the camera and the robot frame is known from calibration [START_REF] Tsai | A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf tv cameras and lenses[END_REF], and the camera is calibrated. For a known landmark of 3D coordinates w X in the world reference frame, the normalized coordinates x = ( c x, c y) of its projection in the camera frame are given by the pinhole model [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

x = Π c T r r T w (r) w X ( 1
)
where Π is the perspective projection operator.

For each visible landmark w X i (i ∈ 1..m), we can derive the following constraints:

C i :            ( c X i , c Y i , c Z i ) = c T r r T w (r) w X i c x i = c Xi c Zi , c y i = c Yi c Zi , c x i ∈ [ c x i ], c y i ∈ [ c y i ] bounded error camera measurement c Z i > 0 front looking camera (2)
We then define the image-based pose estimation problem as a constraint satisfaction problem (CSP) as:

H img : r ∈ [r], {C i , i ∈ 1...m} (3)
This CSP is solved with Contractor Programming [START_REF] Chabert | Contractor programming[END_REF] and Set Inversion via Interval Analysis (SIVIA) [START_REF] Jaulin | Set inversion via interval analysis for nonlinear boundederror estimation[END_REF]. The altitude [z], pitch [θ] and roll [φ] components of the initial domain [r] are set thanks to onboard sensors (altimeter and inertial measurement unit). An outer approximation of the feasible domain for the pose r is obtained in the form of a subpaving. An example of solution set for one robot is shown in Fig. 2a (see [START_REF] Kenmogne | Image-based uav localization using interval methods[END_REF]).

Assuming a bounded error measurement [d] of the distance d between the robot and the base station B is available, it can be used as an additional constraint on the robot position p = (x, y, z), to get a tighter pose estimate:

d = p -b 2 , d ∈ [d] ( 4
)
where b the known position of the base station. We obtain a tighter domain as presented in Fig. 2b.

a. using camera only b. using camera and range to base Dealing with camera tracking failures The subpaving computed with SIVIA can be empty. This situation corresponds to inconsistencies in the measurements, due to failure of the visual tracking of one or several landmarks (see [START_REF] Kenmogne | Image-based uav localization using interval methods[END_REF] for more details). In this case, our approach put aside the image measurements and use only the distance to the base station B to compute S + r k . The solution-set is in this case a ring centered on B.

Using range measurements for cooperative localization

Assuming each UAV is equipped with a communication system with ranging capabilities, each robot therefore measures range to its neighbours and cooperate with them by exchanging its position.

Sharing positions Lets consider one robot R k (k ∈ 1..N ) from the group. N (k) denotes the neighbours of R k , i.e. the robots within communication range.

At each time step, R k first computes an outer subpaving S + r k , that contains all the feasible poses, considering the camera and base distance bounded-error measurements (as presented in Section 2).

Once the pose domain S + r k is computed, the robot computes the bounding box of its position domain [p k ] = proj p S + r k , where is the bounding box operator, and proj p is the projection on the position space. [p k ] is transmitted to all neighboring robots R j , j ∈ N (k), and the distances d k,j between R k and R j are simultaneously measured.

Pose contraction At reception of information (position boxes [p j ] and boundederror distances measurements [d k,j ]) from neighboring robots, each robot R k tries to refine its actual pose domain, by propagating the new distance constraints between R k and each of its neighbours. Recalling that p k = (x k , y k , z k ) is the position of R k , each robot R k contracts a local CSP H k defined as follows:

H k :     p k ∈ proj p (S + r k ), p j ∈ [p j ], j ∈ N (k) d k,j ∈ [d k,j ], j ∈ N (k) d k,j = p k -p j 2 , j ∈ N (k)     (5)
We use interval constraint propagation to solve H k , in order to reduce the pose domain S + r k .

Fixed-point

The constraint network formed by the group of robots contains cycles spanning several robots. The contraction of the local CSPs H k has to be propagated again through the network to improve the pose domain reduction of each robot. If after solving H k , the robot position bounding box [p k ] is reduced, then the robot R k retransmits its updated [p k ] to its neighbourhood. This process is iterated until a fixed-point is reached (no more significant improvement of the robots positions bounding boxes).

Experimental results

The proposed method has been tested with data acquired on Parrot AR-Drone2 UAV, with 5 landmarks made with AprilTag markers (Fig. 3). Image measurement error bounds are set to ±0.5 px and range measurement error is assumed to be within ±5 cm. Subpavings obtained with 4 robots in cooperative localization are presented in Fig. 4. Figure 4b clearly shows how cooperative localization reduces the feasible pose domain when some robots cannot clearly see the landmarks, by propagating position information of the neighbours.

Table 1 shows how making more robots cooperate in the fleet improves localization, first by reducing the width of the computed position domain (Table 1a), and also by improving the precision when using the center of the domain as a point estimate (Table 1b). With 4 UAVs, the average horizontal position error is less than 5 cm for all the drones. In this paper, we proposed a method to solve cooperative localization in a group of UAVs. Computations rely on interval constraint propagation, assuming bounded error image and distance measurements. Each UAV first independently compute a pose domain from its camera measurements, and then exchange position information with its neighbours to further reduce the position domains of the robots in the group. The method has been applied to real data, and enables to improve the positioning precision of the AUVs thanks to position information propagation. The experiments show that increasing the number of robots in the group provides additional constraints on position, and yields smaller uncertainty of the computed robots poses.

Fig. 2 :

 2 Fig. 2: Pose domain computation for a single drone (horizontal projection)

Table 1 :Fig. 3 :Fig. 4 :

 134 Fig. 3: Onboard cameras views at t = 8 s. Landmarks are boxes with a printed pattern.