Marc Sergent

David Goudin

Samuel Thibault

Olivier Aumage

Controlling the Memory Subscription of Distributed Applications with a Task-Based Runtime System

Keywords: memory control, task-based run-time systems, compressed linear algebra, distributed computing

The ever-increasing supercomputer architectural complexity emphasizes the need for high-level parallel programming paradigms. Among such paradigms, task-based programming manages to abstract away much of the architecture complexity while efficiently meeting the performance challenge, even at large scale. Dynamic run-time systems are typically used to execute task-based applications, to schedule computation resource usage and memory allocations. While computation scheduling has been well studied, the dynamic management of memory resource subscription inside such runtimes has however been little explored. This paper studies the cooperation between a task-based distributed application code and a run-time system engine to control the memory subscription levels throughout the execution. We show that the task paradigm allows to control the memory footprint of the application by throttling the task submission flow rate, striking a compromise between the performance benefits of anticipative task submission and the resulting memory consumption. We illustrate the benefits of our contribution on a compressed dense linear algebra distributed application.

I. INTRODUCTION

Efficiently parallelizing an application involves the two fundamental steps of exposing the application parallelism, and of mapping the exposed parallelism onto the available computing resources. Both steps require large amounts of expertise from programmers, in all but the simplest cases. Task-based parallel programming models and their associated runtime systems are becoming popular among such programmers, due to their ability to transparently handle the mapping step, while reducing the burden of the parallelism exposition step. Such programming environments rely on the programmer to define elementary pieces of work (the tasks), and to detail how such tasks relate to each other (the dependencies). Still, expressing task dependencies by hand can be error prone. The Sequential Task Flow (STF) task parallelism model [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF], [START_REF] Bueno | Implementing ompss support for regions of data in architectures with multiple address spaces[END_REF], [START_REF] Gautier | Locality-aware work stealing on multi-cpu and multi-gpu architectures[END_REF] therefore simplifies the programmer job even further by automatically building a graph of task nodes and dependence edges from the sequential algorithmic flow of task submissions to the runtime system, by detecting when a task accesses pieces of data modified by previously submitted tasks, in which case the submission ordering between those tasks should be enforced.

The resulting graph of tasks is then mapped on the platform to schedule the execution of the application. We have shown in previous work [START_REF] Augonnet | StarPU-MPI: Task Programming over Clusters of Machines Enhanced with Accelerators[END_REF], [START_REF] Agullo | Harnessing Supercomputers with a Sequential Task-based Runtime System[END_REF] that this can be efficiently extended to distributed-memory by using MPI for instance.

Most of the research work on dependent tasks has until now been dedicated to the execution stage. However, the task submission stage impacts the application performance as well: if the application submits tasks as short, interspersed sequences, the runtime system does not get sufficiently far look-ahead to perform good task scheduling and to anticipate on data transfer requirements to overlap data transfers time with computations. Conversely, if the application submits tasks as long, bulk sequences, the runtime system may itself waste valuable resources especially in terms of memory, notably in the case of distributed memory, for which reception buffers have to be allocated to receive contributions from other nodes, resulting in sub-optimal performances, and even reaching fatal out-of-memory condition. This paper hence studies the possible cooperation between the application and the runtime system to control the task submission flow from a memory consumption point of view, to allow for sufficient parallelism while keeping the memory footprint under control. The main idea is to integrate a throttling mechanism on the task submission flow, and to drive it by monitoring the memory consumption level. When the memory usage exceeds a predefined upper threshold, the throttling temporarily blocks the task submission routine. Once the memory usage falls below a lower threshold as the result of previously submitted tasks running to completion, the task submission routine is unblocked again. This control flow on task submission is made possible thanks to the property of the STF model to be immune to dependencerelated deadlocks by design. This does not prevent memoryrelated deadlocks (as typically prevented by the Banker's algorithm), but as discussed in section V-C, in our MPI distributed-memory context, this is not a problem in practice because temporary buffers used to store contributions from different MPI nodes easily provide memory to release, and the real concern is to control the submission rate to avoid out of memory conditions.

We experiment our proposal by considering a Boundary Element Methods (BEM) [START_REF] Sauter | Boundary element methods[END_REF] code designed at CEA, which uses a compressed dense linear algebra solver where the compression ratio of each matrix tile cannot be estimated a priori. This application emphasizes the need for a control mechanism between the task submission flow and the memory consumption level, to be able to run large test cases to completion without wasting prohibitive amounts of memory into grossly over-estimated safety margins.

The contributions of this paper are the following:

• A mechanism to monitor the level of memory subscription induced by a flow of submitted tasks. • A mechanism to throttle the flow of tasks submission to smooth out the memory subscription peaks induced by burst task submission sequences, drived by the memory subscription monitor. • The evaluation of these combined mechanisms on a compressed dense linear algebra application, running successfully to completion with a much narrower tolerance gap in terms of memory consumption than without these mechanisms. This paper is structured as follows: Section 2 presents the STF model and its properties that make such a task submission flow control possible, as well as the STF-based StarPU runtime system we used for our experiments. Section 3 exposes the state of the art. Section 4 introduces our contribution about the control of the memory subscription of StarPU. Section 5 describes the distributed, compressed dense linear algebra application we used as the target for the study and the consequences on memory control. Section 6 presents our performance results. Section 7 concludes the paper and gives some prospects of this work.

II. SEQUENTIAL TASK FLOW

The Sequential Task Flow class of parallel programming models characterizes task parallelism environments involving two stages running concurrently: a sequential tasks submission stage, and a parallel tasks execution stage. The submission stage is performed by the application core part, which unfolds its "taskified" algorithm sequentially, submitting the tasks constituting the algorithm one at a time, for subsequent asynchronous accomplishment by the execution stage.

For each submitted task, the application also specifies the pieces of data accessed by the task (R: Read-Only), as well as those possibly modified by the task (R/W: Read-Write or W: Write-Only). Taking into account the pieces of data, their so-called access modes, and the sequential task submission order, the execution stage implemented as a runtime system can build a graph of task dependence edges (usually a Directed-Acyclic Graph (DAG)) and use this graph to decide which sets of tasks may safely run in parallel, and which sets of tasks should strictly be executed following the submission order, to preserve the expected application semantics. For instance, Figure 1 shows a StarPU version of the tiled Cholesky factorization algorithm. The for for (j = 0; j < N; j++) { POTRF(RW,A[j][j]); A key property of STF is that by construction, no submitted task may depend on the result of an as-yet unsubmitted task. Consequently, the submission and execution stages are loosely coupled: The submission stage may temporarily stop submitting tasks at any point in the task flow without any risk of causing deadlocks; Also, the submission stage may wait for the completion of some batch of submitted tasks at any time, but is not conceptually required to do so, except at the very end of the application execution.

for (i = j+1; i < N; i++) TRSM(RW,A[i][j],R,A[j][j]); for (i = j+1; i < N; i++) { SYRK(RW,A[i][i],R,A[i][j]); for (k = j+1; k < i; k++) GEMM(RW,A[i][k],R,A[i][j],R,A[k][j]); } } starpu_task_wait_for_all();

A. The StarPU Task-Based Programming Environment

StarPU is a task-based parallel programming environment and associated runtime system developed by the STORM research team in Bordeaux, France. It targets multicore and heterogeneous, accelerated, computing platforms. StarPU exposes a sequential task flow API to the application. It dynamically schedules and executes submitted tasks on the available cores, following the programmer selected stock scheduling algorithm or a specifically tailored one. When running on an accelerated node, a distributed shared memory (DSM) embedded within StarPU transparently handles data transfers between the main memory space and the accelerators (or GPUs) memory spaces, and manages data replica consistency. The scheduling and DSM engines of StarPU also cooperate together to perform data prefetching and overlap transfers with computations.

Beyond intra-node parallelism, StarPU also provides two supports for interacting with MPI in a distributed context [START_REF] Augonnet | StarPU-MPI: Task Programming over Clusters of Machines Enhanced with Accelerators[END_REF], [START_REF] Agullo | Harnessing Supercomputers with a Sequential Task-based Runtime System[END_REF]. The explicit MPI support defines a set of starpu_mpi_ * wrappers on top of common MPI routines such as MPI_Isend or MPI_Irecv that insert proper data dependence edges between computation tasks and MPI requests. When using StarPU's implicit MPI support instead, the application has to provide an initial data distribution over the computing nodes. Then, each node submits the whole sequence of application tasks identically, and every node is responsible for effectively executing a subset of the global set of tasks in accordance with data ownership; for instance, a node executes a given task if and only if it owns the piece of data the task writes to. Since the whole set of tasks A node really only needs to submit its own tasks, as well as the neighbor tasks at the other end of its own tasks' dependence edges. We have shown in previous work [START_REF] Agullo | Harnessing Supercomputers with a Sequential Task-based Runtime System[END_REF] that thanks to these optimizations, while exposing a programming model which is simple, our approach exhibits performance that is competitive with state-of-the-art distributed environments such as ScaLAPACK or the DPLASMA framework [START_REF] Bosilca | Distributed dense numerical linear algebra algorithms on massively parallel architectures: DPLASMA[END_REF], [START_REF] Wu | Hierarchical dag scheduling for hybrid distributed systems[END_REF] over clusters of more than a hundred hybrid nodes, and scalability perspectives seem to reach at least thousands of nodes.

In order to support various data types (vectors, matrix blocks, sparse matrix blocks, etc.), StarPU uses one data interface per type of data. This interface defines methods for dealing with the data. Figure 2 shows the methods which will be of interest for this article. The allocate and free interfaces are used to allocate and release a piece of data. In them, the application can allocate data with the underlying memory allocator of his choice (malloc, hoard or any other third-party allocator). Some private data also allows the interface to store information for its own use, such as the tile size, leading dimension, number of vector elements, etc. For the most simple cases such as vector or matrix blocks, sending and receiving a piece of data is implemented natively inside StarPU with direct MPI_Isend and MPI_Irecv calls. For more complex cases such as sparse matrix blocks or the BLR-compressed blocks we will discuss in Section V, the pack_data and unpack_data methods can be defined to describe respectively how to pack the piece of data into a contiguous block of bytes (before sending over MPI), and how to unpack such block of bytes (after reception from MPI) into the original data content.

Such kind of interfaces are usually defined once for good: StarPU provides interfaces for most basic data types, and knowledgeable users can define their own interfaces, which can then easily be reused as such by other users in various applications without additional development. In Section V, we describe how we implemented an interface for BLR-compressed matrix tiles, which can now be easily used by various BLR-based applications.

III. STATE OF THE ART OF DISTRIBUTED TASK-BASED

RUNTIMES

Beside StarPU, numerous task-based runtime systems have been proposed to drive applications on top of distributed HPC platforms. Some runtimes such as OmpSs/S-tarSS [START_REF] Bueno | Implementing ompss support for regions of data in architectures with multiple address spaces[END_REF], Kaapi [START_REF] Gautier | Locality-aware work stealing on multi-cpu and multi-gpu architectures[END_REF], Qilin [START_REF] Ravi | Scheduling concurrent applications on a cluster of cpugpu nodes[END_REF] or others [START_REF] Beri | A scheduling and runtime framework for a cluster of heterogeneous machines with multiple accelerators[END_REF] also follow the principle of full run-time task discovery of the sequential task-flow paradigm. Some environments such as Charm++ [START_REF] Mei | Optimizing a parallel runtime system for multicore clusters: A case study[END_REF] are based on the message-driven parallelism instead. Others, such as the DPLASMA framework [START_REF] Bosilca | Distributed dense numerical linear algebra algorithms on massively parallel architectures: DPLASMA[END_REF], [START_REF] Wu | Hierarchical dag scheduling for hybrid distributed systems[END_REF], associate a compilation stage to build a parametric representation of the task graph with an execution stage to map the resulting parametric graph on computation units.

The problem of memory-aware task scheduling or workflow mapping has been the subject of theoretical research works [START_REF] Eyraud-Dubois | Parallel scheduling of task trees with limited memory[END_REF] as well as some implementations such as Charm++ [START_REF] Dooley | A study of memory-aware scheduling in message driven prallel programs[END_REF], or the Bounded Memory Scheduling (BMS) model of dynamic task graphs recently implemented in the concurrent-collection based Qt-CnC runtime system [START_REF] Sbîrlea | Bounded memory scheduling of dynamic task graphs[END_REF]. A prototype for a more general resource-aware task scheduling was developed in OmpSs [START_REF] Tillenius | Resource-aware task scheduling[END_REF]. The topic was also studied with keen interest in the domain of embedded computing, where memory resources may drastically be limited [START_REF] Arras | List scheduling in embedded systems under memory constraints[END_REF], [START_REF] Desnos | Memory analysis and optimized allocation of dataflow applications on shared-memory mpsocs[END_REF]. We refer to Section 4.3 of [START_REF] Agullo | Implementing multifrontal sparse solvers for multicore architectures with sequential task flow runtime systems[END_REF] for a survey of theoretical and practical memory-aware algorithms.

Involving the application in the memory footprint control was explored on the algorithmic side [START_REF] Solomonik | Communication-optimal parallel 2.5d matrix multiplicatoin and lu factorization algorithms[END_REF], and as a collaborative approach between the application and a runtime system [START_REF] Agullo | Multifrontal qr factorization for multicore architectures over runtime systems[END_REF]. However, to the best of our knowledge, no support fully handled on the runtime side, in a transparent manner, has yet been proposed. This is the purpose of our contribution.

Scope of our contribution

A well-known issue of task-based runtime systems when scaling up on distributed supercomputers is to regulate the memory resource subscription without a prohibitive performance penalty.

An attempt of collaborative approach between an application and a task-based runtime system to control the memory subscription has been studied by Agullo et al. in [START_REF] Agullo | Implementing multifrontal sparse solvers for multicore architectures with sequential task flow runtime systems[END_REF]. However, that work is focused on an application which does not fit in the memory of a single computer. It proposes an application-level memory subscription control coupled with the StarPU runtime system, which provides memory consumption feedback to the application thanks to its memory manager engine. Instead, our work aims to control the memory subscription of applications at the runtime system's level and to extend the memory control feature to distributed applications thanks to the STF model.

We discuss in the following section a method based on the STF paradigm to tackle that issue and how we implemented it into the StarPU runtime system.

IV. CONTROLLING THE MEMORY SUBSCRIPTION

With the property of the STF model that decouples the submission step and the execution step, and guarantees that no submitted task depends on the result of unsubmitted tasks, it is possible to temporarily put the submission step on hold without risking dependence-related deadlocks: a submitted task will always be able to complete. The principle of addressing the memory subscription issue is to anticipatively compute, at the task submission time, the memory space that needs to be booked to guarantee its successful completion when its execution time comes later. If this memory space requirements cannot be known exactly at submission time, an upper bound overestimation should be used instead. The runtime system then adds up this amount to the memory subscribed for all the already submitted and not yet completed tasks, and may temporarily block the task submission flow when the estimation of memory booked by the runtime system reaches a maximum subscription threshold. Once the runtime system detects that enough memory space has been freed by tasks completions, it may resume the task submission flow. We implemented our contribution in the StarPU runtime system which implements the STF model, but any runtime system which can hold the task submission flow without introducing deadlocks (e.g. by respecting the sequential order of the task submission) and can track the memory consumption of each submitted task can implement the mechanism we introduce in this paper.

A. Implementation in StarPU with the memory manager

Our implementation extends the use of the memory manager engine of StarPU for the control of the task submission by introducing the starpu_memory_allocate and starpu_memory_deallocate functions, to be called by the data interface (as defined in Section II-A) in its allocate and free methods, i.e. along with actually allocating or deallocating any piece of data. These calls enables StarPU's accounting of the amount of data allocated by the interface. That separation of allocation concern vs accounting concern allows the user to allocate data with his memory allocator, and declare the allocation to StarPU. We then introduce some SLEEP and WAKEUP thresholds for stopping and restarting the task submission. When the submission thread of the application submits a task, an estimation of the memory subscription needed for its execution is computed from the data interfaces of each piece of data involved with that task. The memory manager is then advised of this memory booking with a call to the starpu_memory_allocate routine before actually submitting the task. If the memory manager detects that the memory subscription of StarPU would exceed the SLEEP threshold, the submission thread blocks, thus avoiding to overrun the threshold.

While the submission thread is blocked, already-submitted tasks still proceed progressively to execution and eventually free some memory as the result of their completion. The runtime system is notified of this by calls to starpu_memory_deallocate from the free method of the data interface. When the memory manager detects that the current memory subscription of StarPU has reduced below the WAKEUP threshold, it wakes up the submission thread, which thus resumes submitting tasks. We use different thresholds for SLEEP and WAKEUP to reduce the number of wake-ups and thus the overhead of this mechanism. The SLEEP threshold can be defined by default to the size of available memory on the machine minus a couple of GiB for the operating system, and the WAKEUP threshold be defined to 90% of the SLEEP threshold.

The starpu_memory_allocate/deallocate functions have been used by [START_REF] Agullo | Implementing multifrontal sparse solvers for multicore architectures with sequential task flow runtime systems[END_REF] to ensure that enough memory space is available for the execution of a given set of tasks that need to allocate temporary pieces of data of known size. That was however integrated into the application loop, while we here propose to encapsulate the implementation inside the data interface, thus relieving the application from code changes. Moreover, in the case of our targeted compressed dense linear algebra application, there is no temporary data when executing on a single node: data is assembled completely before starting submitting tasks, and computation is performed in place. However, when employing several distributed nodes instead, tasks do indeed expect contributions from other nodes, which we will call remote data in this paper. The resulting incoming messages have to be stored in temporary buffers, whose memory footprint should be accounted for as well, as discussed in the next section.

B. The distributed case: memory footprint of incoming remote contributions

To store remote data on the receiving node which will execute the task, StarPU automatically allocates a buffer. Since the allocate method of the data interface calls starpu_memory_allocate along with allocating the buffer, the StarPU memory manager is made aware of this additional consumption of memory, and will thus eventually block the submission thread once a lot of buffers have been allocated this way.

By default, StarPU will cache the received remote data [START_REF] Augonnet | StarPU-MPI: Task Programming over Clusters of Machines Enhanced with Accelerators[END_REF], so that when several tasks need the same remote data, only one communication is performed. As a consequence, the corresponding buffers will be kept allocated. The application can, while submitting tasks, make a call to starpu_mpi_cache_flush to specify when a given piece of data will probably not be used by the application algorithm in the close future. For instance, Figure 3 shows how such calls can be quite naturally added to the tiled Cholesky factorization. In that case the buffer containing the cached remote data will be released once all the previouslysubmitted tasks referencing this data have completed. In the Cholesky example of Figure 3, diagonal A[j][j] tiles can for instance be flushed once the corresponding TRSM tasks have completed. The free method of the data interface calls starpu_memory_deallocate along with the actual release so that the memory manager can then possibly wake the submission thread since memory has been released.

for (j = 0; j < N; j++) { POTRF(RW,A[j][j]); for (i = j+1; i < N; i++) TRSM(RW,A[i][j],R,A[j][j]); starpu_mpi_cache_flush(A[j][j]); for (i = j+1; i < N; i++) { SYRK(RW,A[i][i],R,A[i][j]); for (k = j+1; k < i; k++) GEMM(RW,A[i][k],R,A[i][j],R,A[k][j]); starpu_mpi_cache_flush(A[i][j]); } } starpu_task_wait_for_all();
To summarize, this mechanism allows to accurately control the flow of task submission from the application, according to how much memory gets used for the buffers required for receiving remote data from MPI. In the end, StarPU lets the application submit as many tasks as possible (and thus achieve most parallelism), while taking into account the entailed memory consumption.

C. Dealing with remote data of unknown size

In our compressed dense linear algebra case, the size of the remote data received via MPI may be unknown to the receiving side, since the receiving node does not know the compression ratio achieved by the sending node for a given tile. In that situation, the data interface can announce an overestimation to StarPU: for instance, a possible worst case upper bound estimation may be obtained from considering that a compressed tile cannot get larger than an uncompressed tile. The data interface will thus pass the size of an uncompressed tile to starpu_memory_allocate at submission time, to remain on the safe subscription side. When MPI_Irecv completes, the actual size of the compressed tile gets known, and in the unpack method the data interface calls starpu_memory_deallocate with the difference, to account for the amount of memory which did not get used in the end.

As a result, since initially a lot of MPI requests are pending on the receiver side, we first get an overestimation of the memory that will be required to receive them. That overestimation however gets progressively fixed as MPI requests complete, thus releasing memory subscription, and thus allowing more tasks to be submitted by the application for better parallelism, which entail more MPI requests with over-estimations of the memory required to receive them, and so on.

D. Fragmentation issues

Allocating data with varying sizes would often lead to allocation fragmentation, which would waste memory, and thus make our memory accounting underestimate the actual memory usage. In our study, fragmentation was very limited and has thus not posed problems. If fragmentation becomes a concern, a way to deal with it would be to periodically request statistics from the allocator, more precisely the ratio between the amount of allocations vs the corresponding memory usage, and multiply the StarPU memory accounting by this ratio, thus fixing the estimation of memory use for the allocated buffers.

E. Discussion

To summarize, using the starpu_memory_allocate and starpu_memory_deallocate memory accounting functions allows StarPU to guarantee a safe execution of the application, by blocking and restarting its submission loop thanks to thresholds.

For all applications for which the size of the data does not change during the execution (the dense Cholesky factorization falls within this scope for example), StarPU simply needs to deal with the remote data received via MPI, where the size of each piece of data is completely known. Since the behaviour of the memory consumption is similar to the one presented on Figure 4 (and will be discussed in Section VI-B), we do not present those experimental results due to lack of space.

However, some applications may manipulate data which grows in size over the execution, e.g. matrix fill-in. In some cases, this fill-in can be easily estimated, and a static overestimation exposed to StarPU instead of the growing real size, to keep on the safe side.

In other cases, the fill-in can not be predicted accurately, and overestimation would have to be very gross, leading to a poor use of available memory. This is the case we discuss in the next section, taking a tiled Cholesky factorization with BLR compression as a case study, to show how our contribution allows to strike a compromise between performance and memory consumption for such applications.

V. TASK-BASED BLR SOLVER

A. BLR: dense linear algebra and compression techniques

The need for scalable direct solvers in dense linear algebra faces the memory footprint issue. Nowadays, the necessity of storing the full dense matrix in memory limits the maximum size of problem that can be computed at least as much as the computing time.

An algorithmic solution has been developed by Bebendorf to approximate boundary element matrices in [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF]. The idea is that dense blocks of matrices can be represented as a product of two tall and skinny matrices. Some of the numerical precision is lost in the process, but the loss is controlled by a cut factor called , used during the compression step. A dedicated arithmetic for such compressed matrices has also been introduced by Bebendorf on this paper. This arithmetic relies on re-compression steps during the update of the trailing matrix to keep it as compressed as possible during the execution. However, after the re-compression step, the size of each updated piece of data data can change, and thus the pair of matrices that represents each dense block must be reallocated to match that new size.

An implementation of this type of solver exists at CEA for their applications. For coupling this solver with the StarPU runtime system, we use the Chameleon software, which is a collection of solvers on top of task-based runtime systems developed by the HiePACS research team at Bordeaux, France. We developed in Chameleon a StarPU-compliant data interface (as defined in Section II-A) for compressible matrices. This allowed an easy coupling of CEA's solver with the StarPU runtime system, which is one of the runtime systems supported by Chameleon. We then simply kept the same tiled Cholesky algorithm source code, as shown on Figure 3, but extended the underlying kernel into using the CEA-provided BLAS kernels which can operate on various combinations of uncompressed and BLR-compressed data.

B. Issues entailed by using BLR

The size of the compressed tiles will typically tend to grow in an unpredictable way while the application algorithm progresses, due to data contributions from other tiles, notably from the tiles which could not be compressed. This means that starpu_memory_allocate will have to be called by the data interface to account for the additional required space, to delay the submission of tasks which would have used that space. In compressed linear algebra algorithms, the contributions coming from other MPI nodes are not modified, so only the local tiles will grow. We may still end up overflowing the memory, if the inflation of local tiles as tasks complete is bigger than the amount of the cached MPI remote data that can be released when tasks complete.

This issue can be avoided altogether safely by overestimating the amount of memory needed for local tiles, by assuming for instance that they are all uncompressed and making the data interface announce the corresponding size instead of the actual size of the compressed block. That will however make a poor use of the available memory, and StarPU will have to let the application submit fewer tasks at the same time, thus reducing the available parallelism.

On the other hand, since tiles have various compression ratios, tasks will have varying completion time: a GEMM (matrix-matrix product) task with two dense tiles as input and a compressed tile as output will take roughly one order more time to execute than a GEMM task with only compressed tiles. It is thus crucial to let the application submit enough tasks in advance for the runtime scheduler to be able to reorder them to compensate the irregularity.

The user could increase the amount of available memory by running the computation on more MPI nodes, but in practice this is not reasonable because the goals of using compression is not only to reduce the amount of computation, but also to reduce the number of nodes needed to fit the computation in memory, as will be discussed in Section VI-B. Assuming that tiles are uncompressed would for instance lead typically to a 10× overestimation of the required number of MPI nodes.

We are thus aiming at a compromise between achieving good performance by submitting enough tasks in advance, while using as few MPI nodes as possible, and guaranteeing the termination of the execution as much as possible.

C. A compromise approach

Generally speaking, it is not possible to know how much fill-in will happen overall, because matrices can be arbitrarily complex and it is possible that the computation ends up filling the matrix in completely. Users can thus prefer to decide to play on the safe side by assuming that local tiles are uncompressed when making memory subscription, even if that will probably lead to unused memory. They can then choose between executing on few MPI nodes (and thus little available memory, thus little available parallelism as well as lower achieved performance), or spending more MPI nodes on the computation to get the result faster.

With the matrices used in practice, the computation however does not behave so badly, and application users may have a rough idea of how much fill-in will happen overall with the matrices being worked on, and assume e.g. that local tiles will not grow more than a 2× factor. This is for instance true in the experimentation conducted in Section VI (ratio between initial memory and current memory). We can then simply apply this factor on the size initially announced by the data interface for the compressed tiles, it is a much better overestimation than using the uncompressed size mentioned in the previous paragraph. If there can be a guarantee on this inflation factor, execution is then guaranteed to complete without overflowing the memory. The application may even refine the estimation during the execution: when it knows a better estimation for the maximum future size of a tile, the data interface can be made to call starpu_memory_deallocate to release the difference, thus allowing for more parallelism.

In the case of compressed linear algebra, we however do not have any strong guarantee since the inflation of tiles over the whole application execution is not formally bound. In practice, though, the average inflation factor per task is very small. There may be tiles which inflate noticeably, but on average the tile expansion is smaller than the amount of data released when flushing the cached remote data from other MPI nodes once tasks complete.

In that case, we can get bolder, and not overestimate the size of local compressed tiles any more, and when their size increases, the data interface calls starpu_memory_allocate to declare the additional usage. This will make StarPU prevent the application from submitting more tasks for yet more time, until that increase in local tile memory consumption gets compensated by the release of cached remote data from other MPI nodes. We just have to lower the SLEEP/WAKEUP thresholds slightly to have enough memory margin in between. In practice this works really well, as will be seen in the experimental section, e.g. on Figure 4: the amount of subscribed memory nicely evolves between the two SLEEP/WAKEUP thresholds. StarPU can then make fairly good use of all the available memory for caching MPI remote data and allowing the application to submit a lot of tasks at a time.

If users underestimated the inflation factor, i.e. they did not use enough MPI nodes for the computation, the inflation of the local tiles may consume more and more memory until filling the total memory. In that case StarPU will warn the user, and still try to continue the execution, but a single task at a time, to minimize the memory footprint. Resorting to out-of-core mechanisms would be the only way for the runtime system to support this case, and while supported in StarPU, it is not an option on the platform used for the studied test case and is outside the scope of this paper.

In the end, since the memory subscription mechanism is encapsulated once for good inside the data interface, application code can remain as simple as Figure 3, and the user can decide between full safety, better performance, and the number of MPI nodes that will have to be reserved for the computation.

VI. EXPERIMENTS

To validate the compromise approach described in the previous section, we analyze the memory behaviour of a tiled Cholesky factorization on BLR-compressed matrices, typically used at CEA for simulations, on a cluster of several dozens of nodes. We show that the amount of memory consumed by the application during the execution is successfully constrained inside the memory bounds given to the runtime system by the user. We show that the throttling of the task submission flow does not alter the performance of the application compared to other runs of the same application without any limitations on the task submission flow. Finally, we show that our memory subscription control mechanism allows a better memory filling of computing nodes without risking out-of-memory scenarios, which allows to use fewer nodes to perform the computation. We have also applied our approach on various linear algebra applications which use another compression strategy, H-matrices, and a similar approach on sparse linear algebra [START_REF] Agullo | Implementing multifrontal sparse solvers for multicore architectures with sequential task flow runtime systems[END_REF] for the QR factorization with the same conclusions; we here only present results for the BLR-compressed Cholesky case which is much simpler to analyze.

A. Experimental context

We conducted our experiments on the TERA100 Hybrid cluster [START_REF]TERA-100 Hybrid[END_REF] of the CEA. Each computing node is composed of 2 Quad-Core Intel Xeon E5620 CPUs running at 2.4 GHz, and we used between 9 and 81 nodes for our experimental campaign.

For our experiments, we coupled the distributed multicore implementation of the BLR-based dense double-precision complex tiled Cholesky direct solver of the CEA with the Chameleon software.

The application's parameters have been tuned as follows : the size of a block is 512, the chosen MPI data layout is the two-dimensional block cyclic distribution, and we only use square grids of processes. The block size choice is a trade-off between small blocks to compress blocks as much as possible and big blocks to fully exploit the computing units. Square grids of processes with the two-dimensional block cyclic distribution allow to minimize the number and volume of communications and to balance them over all the nodes for dense linear algebra [START_REF] Dongarra | A look at scalable dense linear algebra libraries[END_REF]. However, the BLR-based Cholesky factorization holds most of its dense blocks on the diagonal of the matrix. This means that extra-diagonal nodes will receive more dense blocks, and thus as much amount of data, than diagonal nodes. In summary, square grids of processes with the two-dimensional block cyclic distribution is a rather challenging case in terms of volume of communication for the BLR-based Cholesky factorization.

We used the prio scheduler of StarPU which allows the application to give priorities to tasks at submission time. We designed a set of priorities for each task of the BLR-based Cholesky factorization to enhance the performance of the application, but how we designed them is out of the scope of this paper.

We define two different setups that use different matrix sizes and memory thresholds in Table I. The first setup is used to show the behaviour of the memory control of StarPU and prove its feasibility, while the second is used to study the scalability of our contribution. We see on the top graph of Figure 4 that without memory control, the overestimation curve greatly exceeds the total memory of the machine. Thus, there is no guarantee that the memory will not overflow until the end of the execution. We can also observe that the RSS curve of the top graph stops around 9.5 GB of memory consumption while the one of the bottom graph only goes up to 7 GB. In the top (non-controlled) case, all the memory allocations needed for remote data are performed as soon as possible, which is why the RSS curve's gradient is steep at the beginning of the execution. After that, the gradient of the curve follows the growth of local data, which is caused by the filling phenomenon of the local blocks of the compressed matrix.

B. Results

Instead, the bottom graph shows that with the memory subscription control and well-chosen thresholds (as discussed in IV-A), the memory footprint of the process is kept in check. The growth of the memory footprint of the process is smoother because the memory allocations for remote data are spread throughout the execution. This allows the memory allocator to re-use the memory space used by previously received pieces of data for newly received ones.

Finally, we note that the execution times are almost the same, which lets us assume that the memory control does not affect the performance of the execution.

This assumption is confirmed on Figure 5, which presents the performance results in terms of CPU.hours from 16 to 81 computing nodes on Setup 2.

On the first hand, we can observe that the performance results with and without memory control nicely scale similarly. The amount of consumed CPU.hours for an execution on 81 nodes is only about 6% higher than for an execution on 16 nodes, and the effect on performance of using memory control, which would be thought to harm parallelism, is negligible. This shows that the throttling of the task submission is sufficiently relaxed not to hinder the performance of the computation.

On the other hand, we can see in this plot that the curve without memory control has no points for 25 and 16 nodes. This means that the runs reached an out-ofmemory condition, which will be discussed in the next paragraph. Furthermore, those runs with memory control required slightly less CPU.hours than those with 36 nodes or more. This shows that our contribution allows the user to run their solvers on fewer nodes than before while consuming less CPU.hours, which is really interesting in an industrial context where access to computational resources is constrained, in particular large cluster reservations. Figure 6, which is formatted the same way as Figure 4, shows what happens to the memory subscription when the execution faces an out-of-memory condition without memory control, and how the execution successfully goes to its end with memory control. These runs have been made with 25 computing nodes on Setup 2. On the top plot, StarPU submits all the tasks as soon as possible, and allocates the memory space needed to receive all the remote data that those tasks will need in the future, which causes a quick growth of the RSS curve. When all the data is finally allocated, the growth of local data caused by the filling of the compressed matrix slowly fills up the remaining memory space until the out-of-memory killer is called when the RSS reaches up to 20.8 GB. With the memory control, the memory allocations for receiving remote pieces of data are diluted throughout the execution. Thus, the allocator can re-use the memory space of previously allocated pieces of data for newly allocated ones in a more effective way than without memory control. The bottom plot shows that, indeed, the execution with memory control goes to its end, with a maximum memory footprint of 10GB.

C. Experimental conclusion

These experimental results enlightens three important points of our contribution. 1) Memory control efficiently bounds the memory consumption and limits the memory footprint of the application, as shown in Figures 4 and6.

2) Memory control does not cause any performance decay at high scales, as shown in Figure 5. 3) Memory control allows the successful execution of bigger test cases on a small number of nodes while slightly reducing the required CPU.hours to perform the run, as shown in Figure 5.

VII. CONCLUSION AND PROSPECTS

This paper proposed a runtime-level control of the memory subscription growth during application execution, based on throttling the task submission flow. The runtime system monitors the memory subscription of the application by maintaining an account of the memory space required for each piece of data of each submitted task. When the exact memory footprint of a piece of data is not yet known at the time of submission, the accounting is temporarily given an overestimated quantity, later revised when the exact data size gets settled. With that information, the runtime system can decide to stop the task submission flow when the memory subscription reaches a maximum threshold, and to restart submitting tasks when enough memory space has been freed by the application, without risking dependencerelated deadlocks thanks to the STF model which imposes a sequential ordering for task submission. This mechanism guarantees the successful execution of the application within the memory bounds given by the user.

To deal with a compressed dense linear algebra application for which the size of most of the pieces of data tend to grow in ways unpredictable prior to the execution, we proposed a compromise approach between performance and successful execution within the memory bounds. We validated this approach by coupling a BLR-based tiled Cholesky direct solver of the CEA with the StarPU runtime system. We demonstrated on a real-life compressed matrix of the CEA that the memory consumption growth is constrained thanks to the runtime memory control instead of risking out-of-memory. Our results also showed that throttling the task submission flow has no performance impact when using dozens of nodes. Furthermore, the memory control allows to execute this application on fewer nodes than before, while slightly decreasing the amount of consumed CPU.hours. This approach has also been successfully used for dealing with other compressed linear algebra applications (H-matrices).

The contribution described here considered applications for which the amount of data allocated on each node does not exceed the available memory on the machine. Using out-of-core techniques to overcome this limitation is part of ongoing work in StarPU. Another way we intend to explore to address this issue would be to combine memory-aware task scheduling heuristics with the runtime-level memory control we proposed.

In this paper, we suppose that the sequential submission order of the tasks by the application guarantees that the memory consumption of the application cannot reach an out-of-memory condition if all the tasks are executed sequentially. This is not always the case, as many applications divide their work in assembling (allocating) and computation (deallocating) phases. To address this issue, an idea would be that, instead of blocking the task submission flow when there is no memory space available, to force allocation tasks to wait for memory and continue to submit tasks until deallocation tasks are submitted and executed, which then unlocks the waiting allocation tasks.

Figure 1 .

 1 Figure 1. Source for the tiled Cholesky factorization in the STF model

Figure 3 .

 3 Figure 3. Source code of tiled Cholesky with flush notifications.

Figure 4

 4 Figure4shows the behaviour of the memory subscription of StarPU and the resulting memory footprint of the process (which is shown by the RSS curve) with and without the

Figure 4 .

 4 Figure 4. Memory state of node 1 during the execution on 9 nodes on Setup 1, without (top) and with (bottom) memory control

Figure 5 .

 5 Figure 5. CPU Hours from 16 to 81 nodes on Setup 2

Figure 6 .

 6 Figure 6. Memory state of node 1 during the execution on 25 nodes on Setup 2, without (top) and with (bottom) memory control


Joan Sol Roo
email: joan-sol.roo@inria.fr

Martin Hachet
email: martin.hachet@inria.fr

Interacting with Spatial Augmented Reality

Keywords: Mixed reality, Spatial Augmented Reality, Interaction techniques, Tangible User Interfaces. ACM Classification H.5.1. Information Interfaces and Presentation (e.g. HCI): Multimedia Information Systems

Moving from flat screens into the physical world is an ongoing trend of HCI. Spatial augmented reality (SAR) augments the real world using projectors. It is used in multiple disciplines, especially in design and creative environments such as prototyping and artistic installations. Still, interaction with SAR is little explored. The objective of this thesis is to study SAR and its rich space for interaction, and find novel techniques to take advantage of it. First, previous work on SAR must be studied, then experimental prototypes have to be implemented and evaluated. This paper presents a brief introduction to technical aspects that need to be addressed, in addition to a study of the design spaces related to SAR. A first proof of concept is provided: using augmented objects in front of the screen, as part of a standard desktop environment workflow. Future work will involve iteratively creating prototypes exploring the identified strengths and weaknesses of SAR, analyzing the viability of such prototypes via dedicated user studies.

INTRODUCTION

Augmented Reality (AR) [START_REF] Bimber | Modern approaches to augmented reality[END_REF] combines virtual and real information, in a cohesive manner. In contrast with traditional See-through Augmented Reality (STAR), which requires devices such as glasses or mobile phones, Spatial Augmented Reality (SAR) [START_REF] Raskar | Interacting with spatially augmented reality[END_REF] generates the augmentation directly onto the physical environment. This opens a new realm of possibilities.

Moving the augmentation from a fixed window in front of the user (i.e. glasses or screens) to the environment creates the illusion that the a real object has different physical characteristics (material, behavior, and to a lesser extent, geometry). The augmentation can also be shared by multiple users, and over a wide range of sizes: from a single bean of rice 1 , small surfaces and objects [START_REF] Bandyopadhyay | Dynamic shader lamps: Painting on movable objects[END_REF], rooms [START_REF] Benko | Dyadic projected spatial augmented reality[END_REF][START_REF] Jones | Magical experiences enabled by scalable, adaptive projector-camera units[END_REF] to whole buildings (e.g. [START_REF] Huber | Lightbeam: interacting with augmented real-world objects in pico projections[END_REF].000 square meters on the iMapp Bucharest 555 project 2 ).

The restrictions imposed by the fixed nature of projectors have been overcome on the previous years, using either self contained devices such as PlayAnywhere [START_REF] Wilson | Playanywhere: a compact interactive tabletop projection-vision system[END_REF], or using moving projectors [START_REF] Wilson | Steerable augmented reality with the beamatron[END_REF] or mirrors [START_REF] Pinhanez | The everywhere displays projector: A device to create ubiquitous graphical interfaces[END_REF]. With the creation of small portable laser projectors, called pico projectors, new possibilities arise [START_REF] Willis | Motionbeam: a metaphor for character interaction with handheld projectors[END_REF][START_REF] Willis | Sidebyside: ad-hoc multi-user interaction with handheld projectors[END_REF][START_REF] Huber | Lightbeam: interacting with augmented real-world objects in pico projections[END_REF][START_REF] Molyneaux | Interactive environment-aware handheld projectors for pervasive computing spaces[END_REF][START_REF] Willis | Hideout: mobile projector interaction with tangible objects and surfaces[END_REF]. Cellphone prototypes with embedded projectors3 provide on-board computational power, and extended interaction spaces (i.e. where and how the human-computer interaction takes place). This allows to envision wereable, self contained, smart devices that can project information into the surrounding space (e.g. SixthSense [START_REF] Mistry | Sixthsense: A wearable gestural interface[END_REF]).

Closely related with Spatial augmented reality are Tangible user interfaces (TUIs), which give physicality to virtual information and operations, taking advantage of the human natural skills for interaction with physical objects. Another feature of TUI is reducing the distance between input and output devices, enhancing the embodiment [START_REF] Fishkin | A taxonomy for and analysis of tangible interfaces[END_REF]. Moving to physical objects opens new possibilities, in particular the use of materials with different characteristics, providing not just passive haptic feedback, but also rich expressiveness (for example, non-rigid interfaces). As a result, it provides more intuitive interaction techniques, which facilitate collaboration and learning. Also related with SAR and TUI are smart objects, as envisioned by ubiquitous computing [START_REF] Weiser | The computer for the 21st century[END_REF], organic user interfaces [START_REF] Holman | Organic user interfaces: designing computers in any way, shape, or form[END_REF] and radical atoms [START_REF] Ishii | Radical atoms: beyond tangible bits, toward transformable materials[END_REF]. They envision a world where objects have advanced capabilities for displaying, holding information, communicating with each other and providing interactive functionality [START_REF] Heun | Smarter objects: using ar technology to program physical objects and their interactions[END_REF] and dynamic shape [START_REF] Follmer | inform: dynamic physical affordances and constraints through shape and object actuation[END_REF]. Such technologies can enhance the human communication and thinking process [START_REF] Victor | Humane representation of thought: A trail map for the 21st century[END_REF]. SAR allows us to prototype such objects with technology that is already available. SAR is also called "projected augmented reality", to differentiate from spatial augmentations that focus other senses, such as aural or haptic feedback. The combination of multisensorial spatial augmentations generates rich and immersive experiences, similar as the ones created on Virtual Reality (VR) CAVEs [START_REF] Cruz-Neira | Surround-screen projection-based virtual reality: the design and implementation of the cave[END_REF], but anywhere.

The design space for spatial augmented reality is broad and hold great potential, but it is also heterogeneous and complex. Studying the technical requirements and also the possible interaction techniques, while identifying the ones suited for each case, is the objective of this work. This will enable the creation of new techniques and systems to interact with SAR.

TECHNICAL ASPECTS

The generated illusions by SAR hold great potential, when everyday surfaces gain new aspect and functionality on an almost magical way: a clear example of this is the amount of artistic performances involving Spatial augmented reality (on this context, SAR is called projection mapping). But the illusions are fragile. In order to generate the augmentation is required to align a virtual scene with the real world, and in some cases take into account the user(s) head position. The basic issues to address are:

• Calibration: Scene elements, such as cameras and projectors, have different positions and orientation, and an unifying coordinate system is required. Additionally, optical devices are not perfect and possess lens distortion, which needs to be accounted for [START_REF] Audet | A user-friendly method to geometrically calibrate projector-camera systems[END_REF].

• Geometry: Knowing the surface geometry is indispensable to align the virtual and real information. This requires a 3d model, which can be generated either manually or by scanning [START_REF] Curless | From range scans to 3d models[END_REF][START_REF] Rocchini | A low cost 3d scanner based on structured light[END_REF].

• Position: knowing the geometry is not enough, it is also required to know where the geometry lies in relationship with the real scene. This can be measured manually, using sensors [START_REF] Lee | Automatic projector calibration with embedded light sensors[END_REF][START_REF] Raskar | Rfig lamps: interacting with a self-describing world via photosensing wireless tags and projectors[END_REF] or tracking, either with [START_REF] Rekimoto | Matrix: A realtime object identification and registration method for augmented reality[END_REF] or without markers. Geometry and position can be obtained simultaneously using SLAM [START_REF] Klein | Parallel tracking and mapping for small ar workspaces[END_REF][START_REF] Izadi | Kinectfusion: real-time 3d reconstruction and interaction using a moving depth camera[END_REF].

• Material: the color and reflectance of the surface where a pixel is projected will affect the observed color.

A nice example of how to overcome this is Illumi-Room [START_REF] Jones | Illumiroom: peripheral projected illusions for interactive experiences[END_REF].

• Light condition: Similar to the material, ambient light would affect the final result, and it will compete with the projector luminosity.

• Shadow casting: when a pixel path from the projector reaches an object, it stops traveling. This produces shadows. More than one projector can be used [START_REF] Sukthankar | Dynamic shadow elimination for multi-projector displays[END_REF], which requires additional considerations.

• Dynamic environments: any of the previous factors can change over time, and the system requires to take this into account. On highly dynamic environments, it can be required to predict the change, taking into account any computation delays.

The selected strategy for each item will depend on the complexity of the required solution. Basic video mapping installations use controlled light conditions and static objects [START_REF] Dalsgaard | 3d projection on physical objects: design insights from five real life cases[END_REF], while highly interactive spaces require more robust approaches [START_REF] Jones | Magical experiences enabled by scalable, adaptive projector-camera units[END_REF]. Either way, the process has not yet reached the point where it can be done completely automatically. These are key research areas to explore on the future.

INTERACTION TECHNIQUES

Related work has provided both classifications of the application space, and detailed technological solutions for interaction with augmented spaces. A brief overview is presented below.

Application space

In order to know more about a system it is required to answer some basic questions. Here are the questions on the context of SAR:

How is it done?

Ens et al. [START_REF] Ens | Ethereal planes: A design framework for 2d information space in 3d mixed reality environments[END_REF] presented a thorough study of augmented spaces with planar interfaces, finding seven key dimensions (Perspective, movability, proximity, input mode, tangibility, visibility and discretization) grouped in three categories (frame of reference, spatial manipulation and spatial composition). A factor not taken into account is the scale of the augmentation. Objects small enough to be handled can benefit of direct interaction. Reaching the room size, whole body interaction could be preferable. Larger augmented spaces might require indirect approaches, such as maquetes or pointing techniques.

Dubois and Nigay [START_REF] Dubois | Augmented reality: Which augmentation for which reality?[END_REF] study where the focus of the task is (either the virtual or real object) and what the nature of the augmentation (evaluation or execution). They show the interaction technique and augmentation mode will depend on these factors. Tangibility is closely related with execution, while the graphical quality is more important for evaluation.

What does it mean?

Fishkin et al. [START_REF] Fishkin | A taxonomy for and analysis of tangible interfaces[END_REF] presented a semantic space for TUI, taking into account metaphor and embodiment of the interface. Abstraction is also relevant when talking about interfaces. Burner [START_REF] Bruner | Toward a theory of instruction[END_REF] presents three learning levels: enactive, iconic and symbolic, each one with increasing abstraction, which are closely related with Fishkin's taxonomy. TUI involves more natural and intuitive interaction [START_REF] Jacob | Reality-based interaction: a framework for post-wimp interfaces[END_REF], less abstract than traditional I/O systems. TUIs specificness simplifies the interface, but it also reduces its flexibility. Conversely, highly abstract systems such as desktop computers, while harder to learn, have proven to be flexible and empower the users for precision tasks.

How long will it last?

The lifespan of the interface is also a factor to take into account. Doring et al. [START_REF] Döring | A design space for ephemeral user interfaces[END_REF] study the design space for ephemeral interfaces, that is, interfaces where at least one component is created to last a limited amount of time. This kind of interfaces have strong emotional impact.

Working with SAR implies working with light, which is itself an ephemeral element, and it can be combined with other ephemeral materials, such as soap, fog [START_REF] Nakamura | Bubble cosmos[END_REF] or ice [START_REF] Virolainen | Cool interaction with calm technologies: experimenting with ice as a multitouch surface[END_REF]. Interfaces can also be created by a temporary arrange of otherwise persistent objects [START_REF] Walsh | Ephemeral interaction using everyday objects[END_REF].

Input techniques

The interaction space is rich, and several alternatives are available, with complementary characteristics.

Manipulation

This input technique is the one used on TUIs, but is not reduced to simply manipulate rigid bodies. The possibility to directly touch anywhere, and to do so on an expressive way, is an ongoing area of research using both on-object sensors [START_REF] Sato | Touché: enhancing touch interaction on humans, screens, liquids, and everyday objects[END_REF] and vision [START_REF] Harrison | Worldkit: rapid and easy creation of ad-hoc interactive applications on everyday surfaces[END_REF]. Also, the interaction with non-rigid materials is currently being studied, with great advances both on technology (tracking deformable objects with [START_REF] Zollhöfer | Real-time non-rigid reconstruction using an rgb-d camera[END_REF] and without [START_REF] Newcombe | Dynamicfusion: Reconstruction and tracking of non-rigid scenes in real-time[END_REF] a rigid scan) and the application space [START_REF] Troiano | User-defined gestures for elastic, deformable displays[END_REF][START_REF] Holman | Organic user interfaces: designing computers in any way, shape, or form[END_REF][START_REF] Ramakers | Paddle: highly deformable mobile devices with physical controls[END_REF][START_REF] Troiano | Deformable interfaces for performing music[END_REF].

Augmented tools provide augmented behavior. They can be either specific [START_REF] Marner | Physical-virtual tools for spatial augmented reality user interfaces[END_REF] ("workbench metaphor") or generic [START_REF] Walsh | Ephemeral interaction using everyday objects[END_REF], and provide a range of abstraction according to the degree of embodiment.

The main drawback of this technique in SAR is that shadows are harder to avoid when two objects are close together, but it has not prevented this technique from providing applications, such as Illuminating clay [START_REF] Piper | Illuminating clay: a 3-d tangible interface for landscape analysis[END_REF].

Body and speech

Elepfandt et al. [START_REF] Elepfandt | Multimodal, touchless interaction in spatial augmented reality environments[END_REF] studied the interaction space around the user, and concluded that gaze, gestures and speech are best suited for SAR. Gestural interaction has become more popular since the release of commodity depth sensors, such as Microsoft Kinect [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF]. Hand based interaction can be implemented using either sparse approaches or dense position estimation using kinetic models (Leap Motion and kinect2 [START_REF] Sharp | Accurate, robust, and flexible real-time hand tracking[END_REF]). Such natural interaction could be incorporated easily on our everyday interaction with agents and artistic performances, but might fall short for intense or precise activities.

Pointing

Interacting on a higher level of abstraction can be done by controlling a cursor. Direct pointing is frequently used in VR [START_REF] Hinckley | A survey of design issues in spatial input[END_REF] and have been used in SAR, either using tools [START_REF] Bandyopadhyay | Dynamic shader lamps: Painting on movable objects[END_REF][START_REF] Thomas | Spatial augmented reality-a tool for 3d data visualization[END_REF] or fingers [START_REF] Ridel | The revealing flashlight: Interactive spatial augmented reality for detail exploration of cultural heritage artifacts[END_REF]. Relative pointing in space have been studied for flat [START_REF] Nacenta | Perspective cursor: Perspective-based interaction for multi-display environments[END_REF] and arbitrary surfaces [START_REF] Gervais | Pointing in Spatial Augmented Reality from 2D Pointing Devices[END_REF], proving to be a valid option to remote and/or precise tasks leveraging standard 2D/3D input devices. Flexible switching between absolute and relative pointing has also being studied, but only for flat surfaces [START_REF] Forlines | Hybridpointing: fluid switching between absolute and relative pointing with a direct input device[END_REF].

Desktop-based augmented spaces such as "the office of the future" [START_REF] Raskar | The office of the future: A unified approach to image-based modeling and spatially immersive displays[END_REF] and augmented surfaces [START_REF] Rekimoto | Augmented surfaces: a spatially continuous work space for hybrid computing environments[END_REF] extended the "desktop metaphor" to the environment, while preserving the potential of the WIMP paradigm (windows, icons, menus, pointer).

The presented interaction techniques allow to interact with augmented reality on different levels of abstraction. The selected interaction technique will strongly depend on the task nature. Allowing the user to switch between interaction techniques could provide an unifying interaction framework.

INTEGRATED AUGMENTED DESKTOP

As a first proof of concept of unified interaction with SAR, we created an augmented working space in front of a standard desktop computer 4 . Using perspective cursor [START_REF] Nacenta | Perspective cursor: Perspective-based interaction for multi-display environments[END_REF] allows us to point on real objects. By restricting the cursor influence to a window in front of the screen, a seamless space is created, where a tangible objects can be reached by a cursor, as if they were virtual 3D objects. Such objects are not virtual, and can be handled outside of the screen space, enabling direct and gestural interaction. 4 The mentioned work was submitted to TEI '16 This system enables to iterate back and forth between precise and powerful desktop applications and natural interaction, in order to create rich interactive tangible objects.

CONCLUSIONS

On this work, the interaction space with SAR was studied. The key taxonomies and the available input techniques were selected from the bibliography, in order to understand the exploration space. Also, the technical issues to tackle were considered, and the currently available solutions were presented.

Once the context of interaction with SAR is defined, the next step is to study the presented techniques in action, by constructing prototypes, and to explore novel interaction techniques and uses cases for the rich possibilities of SAR. A first prototype is presented. In the future, we will continue exploring novel and integrated ways to interact with SAR.

https://vimeo.com/130165596

https://vimeo.com/107203878

http://www.cnbc.com/id/102713692


A E Koudou
email: efoevi.koudou@iecn.u-nancy.fr

A Matsumoto-Yor property for Kummer and Wishart random matrices

Keywords: Wishart distribution, Matsumoto-Yor property, Matrix Kummer distribution, Matrix Beta distribution

Introduction

If X, Y are two real positive and independent random variables and if f : ]0, ∞[→]0, ∞[ is a decreasing and bijective function, define U := f (X + Y ) and V := f (X)f (X + Y ). For f (x) = 1/x, if X follows a generalized inversed Gaussian (GIG) distribution and if Y has a gamma distribution with a suitable matching of the parameters, then U and V are independent. This fact, partially observed in [START_REF] Matsumoto | An analogue of Pitman's 2M -X theorem for exponential Wiener functional, Part II: the role of the generalized inverse Gaussian laws[END_REF] and completely proved in [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF], is known in the recent literature as the Matsumoto-Yor property. Furthermore, it characterizes the product of GIG and gamma laws, as shown in [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF].

In [START_REF] Koudou | Independence properties of the Matsumoto-Yor type[END_REF], other independence properties of the Matsumoto-Yor type have been proved, by observing that there are essentially four possible functions f , and characterizations of the corresponding laws have been given under the assumptions of existence and smoothness of densities (see also [START_REF] Koudou | Which distributions have the Matsumoto-Yor property? To appear in Electronic Communications in Probability[END_REF], where these assumptions have been slightly relaxed). For f (x) = log(1 + x) -log x the corresponding independence property states the following: if X, Y are independent random variables such that the law of is a Kummer law

K (2) (a, b, c)(dx) := Cx a-1 (1 + x) -a-b e -cx 1 (0,∞) (x)dx, a, c > 0, b ∈ R (C is the normalizing constant) and such that the law of Y is the gamma distribution γ(b, c)(dx) = c b Γ(b)
x b-1 e -cx 1 (0,∞) (x)dx, then the random variables

U := 1 + 1 X+Y 1 + 1 X , V := X + Y (1.1) are independent. Moreover, U ∼ Beta(a, b) and V ∼ K (2) (a+b, -b, c) where Beta(a, b)(dx) = Γ(a + b) Γ(a)Γ(b) x a-1 (1 -x) b-1 1 {0<x<1} dx.
The Kummer and gamma distributions exist also as distributions on the set of positive definite matrices (the Wishart distribution in the gamma case), and matrix versions exist for the GIG-gamma Matsumoto-Yor property (see Letac and Weso lowski, 2000 and Massam and Weso lowski, 2006). Therefore, a natural question is whether the independence of the random variables U and V defined by (1.1) still holds if X and Y are r × r random matrices (of course, one has to formulate U and V as matrices, see (2.4)). We prove that the answer is affirmative. Although the proof only relies on a Jacobian calculation, the framework of matrices makes it non-trivial.

In Section 2 we state the result after a recall of the definition of matrix Kummer, Wishart and beta distributions. Section 3 is devoted to the proof and Section 4 to a short discussion.

The result

Let r ≥ 1 be an integer. Denote by M r the Euclidean space of r × r real symmetric matrices, and by I the unit element of M r . Consider the inner product A, B = tr(AB) on M r . Let M + r be the cone of positive definite matrices in M r . For Σ ∈ M + r and b ∈ {0, 1 2 , 1, 3 2 , . . . , r-1 2 } ∪ ( r-1 2 , ∞), the Wishart distribution γ(b, Σ) (see for instance Letac and Weso lowski, 2000 and Massam and Weso lowski, 2006) is defined as the law of a random matrix Y valued in the closure of M + r , with Laplace transform

IE(e σ,Y ) = det Σ det(Σ -σ) b , Σ -σ ∈ M + r . If b > r-1
where Γ r is the multivariate gamma function defined for any complex number z with Re(z)

> (r -1)/2 by Γ r (z) = π r(r-1)/4 r j=1 Γ z - j -1 2 .
For α, β > r-1 2 , the beta distribution Beta(α, β) (also called matrix variate beta of type I, see [START_REF] Nagar | Matrix-variate Kummer-beta distributions[END_REF] on

M + r is Beta(α, β)(du) = Γ r (α + β) Γ r (α)Γ r (β) (det u) α-r+1 2 (det(I -u)) β-r+1 2 1 U (u) du, (2.2)
where U is the set of matrices u in M + r such that

I -u ∈ M + r . For a > r-1 2 , b ∈ R, Σ ∈ M + r , the matrix Kummer distribution on M + r is defined by K(a, b, Σ)(dx) = 1 Γ r (a)ψ(a, a -b + r+1 2 ; Σ) (det x) a-r+1 2 (det(I+x)) b exp(-Σ, x )1 M + r (x) dx, (2.
3) where ψ is the confluent hypergeometric function of second kind with matrix argument (see Joshi and Joshi, 1985, formula (2)). This distribution is called in the literature (see e. g. [START_REF] Gupta | Matrix-variate Kummer-Dirichlet distributions[END_REF] the Kummer-Gamma distribution or the Kummer distribution of type II. In this paper we concisely call it the Kummer distribution.

We now state the result of the present paper.

Theorem 2.1 Consider a, b, Σ such that a > r-1 2 , b -a > r-1 2 and Σ ∈ M + r .
Let X and Y be two independent random matrices valued in M + r . Assume that X follows the Kummer distribution K(a, b, Σ) and Y the Wishart distribution γ(ba, Σ). Then, the random matrices

U := [I + (X + Y ) -1 ] 1/2 [I + X -1 ] -1 [I + (X + Y ) -1 ] 1/2 , V := X + Y (2.4)
are independent. Furthermore, U ∼ Beta(a, ba) and V ∼ K(b, a, Σ).

3 Proof

For x, y ∈ M + r define u := [I + (x + y) -1 ] 1/2 [I + x -1 ] -1 [I + (x + y) -1 ] 1/2 , v := x + y. (3.1)
Observe that, since x and y are positive definite matrices, it is also the case for u, v and Iu. (It is obvious for u and v. To check it for Iu, one can write

I -u = [I + (x + y) -1 ] 1/2 [I + (x + y) -1 ] -1 -[I + x -1 ] -1 [I + (x + y) -1 ] 1/2
and use the fact that if A and B are positive definite matrices such that A -B is positive definite, then B -1 -A -1 is positive definite). Therefore, the transformation

T : (x, y) → (u, v) maps M + r × M + r to U × M + r .
It is clearly bijective. We now compute the Jacobian of T -1 .

The Jacobian

It is convenient to write T = T 2 • T 1 , where T 1 : M + r × M + r → U × M + r and T 2 : U × M + r → U × M + r are defined by T 1 (x, y) = (w, z) := [I + x -1 ] -1 , x + y and T 2 (w, z) = (u, v) := [I + z -1 ] 1/2 w[I + z -1 ] 1/2 , z . We have (x, y) = T -1 1 (w, z) = [w -1 -I] -1 , z -[w -1 -I] -1 (3.2) and (w, z) = T -1 2 (u, v) = [I + v -1 ] -1/2 u[I + v -1 ] -1/2 , v . (3.3) Lemma 3.1 The Jacobian of the transformation T -1 is J(u, v) = det(I + v) r+1 2 det(I + v -uv) -r-1 det v r+1 2 .
(3.4) PROOF: We compute the Jacobian of T -1 via the Jacobians J (u,v)→(w,z) and J (w,z)→(x,y) . J (u,v)→(w,z) is the determinant of the differential of T -1 2 at (u, v). This differential is a linear endomorphism of M r × M r which, by (3.3), can be writen as D E 0 Id Mr where E does not need to be computed and D is the differential of the function u → [I + v -1 ] -1/2 u[I +v -1 ] -1/2 for fixed v. Thus, J (u,v)→(w,z) equals det D (we use the same notation for the determinant of a matrix and the determinant of a linear operator). According to [START_REF] Muirhead | Aspects of Multivariate Statistical Theory[END_REF], Theorem 2.1.7, we have

det D = det[I + v -1 ] -1/2 r+1 = det[I + v -1 ] -r+1 i.e. J (u,v)→(w,z) = det[I + v -1 ] -r+1 2 . (3.5)
Let us now compute J (w,z)→(x,y) . By (3.2), we have x = [w -1 -I] -1 so that the differential G(w) of x with respect to w is the linear map

h → [w -1 -I] -1 w -1 hw -1 [w -1 -I] -1 on M r . Its determinant is, again by Muirhead (1982), equal to det(G(w) = det [w -1 -I] -1 w -1 r+1 = det(I -w) -(r+1) .
Thus, by (3.2), we have

J (w,z)→(x,y) = det G(w) 0 H Id Mr = det(G(w) = det(I -w) -(r+1) (3.6)
where no detail is needed about H. It follows from (3.5)and (3.6) that

J(u, v) = J (u,v)→(w,z) J (w,z)→(x,y) = det[I + v -1 ] -r+1 2 det(I -w) -(r+1) . (3.7) We have det[I + v -1 ] = det v -1 det(I + v). (3.8)
Since det(I + AB) = det(I + BA) for any symmetric matrices A and B, we write

det(I -w) = det I -[I + v -1 ] -1/2 u[I + v -1 ] -1/2 = det I -u[I + v -1 ] -1 = det I -uv[I + v] -1 = det(I + v) -1 det(I + v -uv).
(3.9) Substituting ( 3.8) and (3.9) in (3.7) we obtain (3.4).

Density of (U, V )

Denote by f X and f Y the densities of X and Y respectively, and by f (U,V ) the density of (U, V ). Since X and Y are independent,

f (U,V ) (u, v) = f X (x)f Y (y)J(u, v)1 U (u)1 M + r (v) (3.10)
where (x, y) = T -1 (u, v) and J(u, v) is the Jacobian given by (3.4). Recall that X ∼ K(a, b, Σ) and Y ∼ γ(ba, Σ). By (2.2) and (2.3), denoting by C the product of all constants, (3.10) can be written, for any (u, v)

∈ U × M + r , f (U,V ) (u, v) = C (det x) a-r+1 2 (det(I + x)) -b exp(-Σ, x ) ×(det y) b-a-r+1 2 exp(-Σ, y ) J(u, v) = C (det x) a-r+1 2 (det(I + x)) -b (det y) b-a-r+1 2 × exp(-Σ, x + y ) J(u, v). (3.11)
In the following lemma we express det x, det y, det(I + x) in terms of u and v.

Lemma 3.2 Consider x, y, u, v such that T (x, y) = (u, v) as in (3.1). Then,

det x = det (I + v -uv) -1 det u det v, (3.12
)

det(I + x) = det (I + v -uv) -1 det(I + v), (3.13)
det y = det v det (I + v -uv) -1 det(I + v) det(I -u). (3.14)
PROOF: By (3.1) we have

x -1 = I + v -1 1/2 u -1 I + v -1 1/2 -I = I + v -1 1/2 u -1 -[I + v -1 ] -1 I + v -1 1/2 ,
which gives

x = I + v -1 -1/2 u -1 -[I + v -1 ] -1 -1 I + v -1 -1/2 . (3.15) Since I + v -1 = v -1 (I + v) we have u -1 -[I + v -1 ] -1 = u -1 (I + v -uv)(I + v) -1 , so that ( 3.15) gives x = v 1/2 (I + v) -1/2 (I + v)(I + v -uv) -1 uv 1/2 (I + v) -1/2 = v 1/2 (I + v) 1/2 (I + v -uv) -1 uv 1/2 (I + v) -1/2 . (3.16)
Taking the determinant we obtain (3.12). The definitions of u and v imply

det u = det(I + v -1 ) det(I + x -1 ) -1 = (det v) -1 det(I + v) det x det(I + x) -1
from which we get det(I + x) = (det u) -1 (det v) -1 det(I + v) det x. Then, (3.13) follows from (3.12).

To obtain (3.14), we use (3.16) to write To complete the proof of the theorem, we plug (3.12), (3.13), (3.14) and (3.4) into (3.11), replace x + y with v and the result follows. Here are the details: for any

y = v -x = v 1/2 Iv 1/2 -v 1/2 (I + v) 1/2 (I + v -uv) -1 u(I + v) -1/2 v 1/2 = v 1/2 I -(I + v) 1/2 (I + v -uv) -1 u(I + v) -1/2 v 1/2 .
(u, v) ∈ U × M + r , f (U,V ) (u, v) = C det(I + v -uv) r+1 2 -a det u a-r+1 2 det v a-r+1 2 × det(I + v -uv) b det(I + v) -b (det v) b-a-r+1 2 det(I + v -uv) a-b+ r+1 2 det(I + v) b-a-r+1 2 det(I -u) b-a-r+1 2 det(I + v) r+1 2 det(I + v -uv) -r-1 det v r+1 2 exp(-Σ, v ) = C det u a-r+1 2 det(I -u) b-a-r+1 2 det v b-r+1 2 det(I + v) -a
from which we see that U and V are independent and follow the laws Beta(a, ba) and K(b, a, Σ) respectively (see (2.2) and (2.3).

Discussion

1. As said in the introduction, the independence property established in the present paper has been proved in the case r = 1 by [START_REF] Koudou | Independence properties of the Matsumoto-Yor type[END_REF], where the converse has been proved too, thus providing a characterization of the product of Kummer and gamma laws under the assumption of existence and smoothness of densities. It is highly likely, although not easy to prove, that this characterization holds also in the case of matrices. Namely, one can make the following conjecture: Let X and Y be two independent random matrices valued in M + r with smooth densities. The random matrices U and V defined by (2.4) are independent if and only if there exist a, b, Σ such that a > r-1

2 , ba > r-1 2 , Σ ∈ M + r and such that X follows the Kummer distribution K(a, b, Σ) and Y the Wishart distribution γ(ba, Σ). 2. The Matsumoto-Yor property has been discovered in a particular case by Matsumoto and Yor (2001) while studying certain exponential functionals of the Brownian motion. An interesting question is whether there exists a stochastic process whose behaviour could be linked to the independence property shown in the present paper in the case r = 1 as well as in the matrix framework.

  As a consequence, det y = det v det I -(I + vuv) -1 u = det v det (I + vuv) -1 (I + vuvu) and (3.14) follows since I + vuvu = (Iu)(I + v).

(we refer to this case as the non-degenerate case), the Wishart distribution γ(b, Σ) has a density:γ(b, Σ)(dy) = (det Σ) b Γ r (b) (det y) b-(r+1)/2 exp(-Σ, y )1 M + r (y) dy,(2.1)