Manuel Moussallam

Gaël Richard

Laurent Daudet

AUDIO SOURCE SEPARATION INFORMED BY REDUNDANCY WITH GREEDY MULTISCALE DECOMPOSITIONS

Keywords: Simultaneous sparse approximation, audio source separation, greedy decompositions

This paper describes a greedy algorithm for audio source separation of repeated musical patterns. The problem is understood as retrieving from a set of mixtures the part that is redundant among them and the parts that are specific to only one mixture. The key assumption is the sparsity of all the sources in the same multiscale dictionary. Synthetic and real life examples of source separation of hand cut repeated musical patterns are exposed. Results shows that the proposed method succeeds in simultaneously providing a sparse approximant of the mixtures and a separation of the sources.

INTRODUCTION

There are at least two specific cases of audio source separation problems where redundancy plays a fundamental role. Common signal separation is a problem where, from a set of mixtures, one tries to recover a source that is shared among all of them. Practical applications range from film music extraction [START_REF] Liutkus | Separation of music+effects sound track from several international versions of the same movie[END_REF] to multichannel denoising [START_REF] Tropp | Simultaneous sparse approximation via greedy pursuit[END_REF][START_REF] Gribonval | Average case analysis of multichannel thresholding[END_REF]. Repeating pattern separation [START_REF] Rafii | A simple music/voice separation method based on the extraction of the repeating musical structure[END_REF], focuses on separating a varying component (e.g. the singing voice) from a repeating background (e.g. musical accompaniment).

These two separation problems can be linked because they share the same underlying source model. A mixture X i indexed by i, is understood as a combination of an individual source P i , specific to the mixture, and a source component X c that is shared among all the mixtures (though potentially distorted in a different manner in each mixture).

In the common signal separation problem, the individual sources are often considered as noise and the shared component is the signal of interest. Redundancy in this case is the result of a multisensor acquisition [START_REF] Gribonval | Average case analysis of multichannel thresholding[END_REF] or the existence of multiple versions [START_REF] Liutkus | Separation of music+effects sound track from several international versions of the same movie[END_REF]. In the repeating pattern separation framework, the shared source is the musical background (e.g This work was partly supported by the QUAERO Programme, funded by OSEO, French State agency for innovation.

LD is on a joint position between Univ. Paris Diderot and Institut Universitaire de France accompaniment) that remains stable while occurring several times in the music and the individual sources are the parts varying among these occurrences (e.g solo instrument, singing voice). Redundancy is here the consequence of musical repetitions.

State of the art methods addressing the problem of Repeating pattern separation (e.g the REPET algorithm [START_REF] Rafii | A simple music/voice separation method based on the extraction of the repeating musical structure[END_REF][START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF]) are based on element-wise classification of a Time-Frequency (TF) representation. A TF mask (usually based on the power spectral density of the mixtures) is constructed for the repeating musical background and the separation is performed by means of Wiener filtering relative to this mask. Often (e.g in [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF]) an assumption is made on the individual sources, namely that they are sparse in the TF domain. In the same spirit, in [START_REF] Huang | Singing-voice separation from monaural recordings using robust principal component analysis[END_REF], the authors also consider the individual sources to be sparse while the shared component is captured in a lowrank approximant of the spectrogram, a matrix factorization scheme known as Robust PCA [START_REF] Candès | Robust principal component analysis?[END_REF].

Interestingly, the same sparsity hypothesis is also at the core of methods addressing the Common signal separation problem. The basic assumption (e.g in Sparse Component Analysis (SCA) [START_REF] Gribonval | A survey of Sparse Component Analysis for blind source separation: principles, perspectives, and new challenges[END_REF], or in simultaneous approximation problems [START_REF] Tropp | Simultaneous sparse approximation via greedy pursuit[END_REF][START_REF] Gribonval | Average case analysis of multichannel thresholding[END_REF]) is that the shared component has a sparse expansion in a dictionary Φ of waveforms called atoms.

In this work, we address the Repeating pattern separation problem using a sparse decomposition of the mixtures in a redundant dictionary. However, we consider that the shared source and the individual ones are no different in nature, and thus may all be sparsely decomposed in the same dictionary. Section 2 details the problem formulation and the sparse source models adopted. Section 3 introduces the greedy algorithm proposed in this work. Section 4 presents a comparison of behavior with TF based methods on synthetic and real-life examples. Finally Section 5 exposes the proposed separation scheme as a byproduct of a more general compression system.

SIMULTANEOUS APPROXIMATION PROBLEM

Let us formulate the source separation problem as a simultaneous approximation paradigm. Indeed, the separation is obtained from jointly estimating both the shared source and the individual ones.

General formulation

Let us now denote X ∈ R I×N the matrix of I mixtures X i ∈ R N each being of dimension N , and Φ ∈ R D×N an overcomplete dictionary of D unit-normed waveforms called atoms. An approximant X of X on Φ is of the form: X = C X • Φ where C X ∈ R I×D is sparse, meaning a large part of its values are zeros. The simultaneous approximation problem consists in jointly minimizing the divergence between data and approximant, and the number of non-zero elements. One formulation is:

min C X 0 s.t f (X -C X • Φ) ≤ ǫ (1)
where f is a divergence measure of interest (e.g. a squared reconstruction error), . 0 is the l 0 pseudo-norm 1 and ǫ is a desired level of precision. Since a strict l 0 problem is NP-hard to solve, a commonly adopted reformulation is a penalized version:

C X = arg min f (X -C X • Φ) + λ • C X p,q (2
) where . p,q is a mixed norm2 . p and q can be chosen depending on the desired sparsity and λ is a parameter that controls the weight of the sparsity constraint. It has been shown [START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF] that mixed norms can enforce structured sparsity. Actually, a column of C X filled with non-zero elements indicates that the corresponding atom can be found in all the mixtures and thus belongs to the shared source.

While convex optimization algorithms have been proposed along with structured sparsity priors [START_REF] Gribonval | A survey of Sparse Component Analysis for blind source separation: principles, perspectives, and new challenges[END_REF][START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF], greedy methods solving this problem are variants of Simultaneous Orthogonal Matching Pursuit (SOMP) [START_REF] Tropp | Simultaneous sparse approximation via greedy pursuit[END_REF]. This formulation is adapted when one tries to recover a shared component that is sparse in the dictionary and implicitly makes the assumptions that components from the individual sources will not be selected. In this context, the separation can be explained as a denoising of a multichannel signal based on inter-channel redundancies.

Distinguishing two different sparsities

In some situations, including music source separation, the previous formulation is not fully satisfactory. While a separation of the background is still desirable, the assumption that the individual sources cannot be sparsely represented in the same dictionary as the shared one does not hold any more. Without any knowledge of the sources characteristics or production mechanism (e.g source/filter modeling for singing voice) there is no reason to consider the shared and the individual components to be of a different kind.

Actually, it has been shown [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF] that most musical signals are efficiently and sparsely decomposed in Fourier-based dictionaries (e.g. Gabor frames).

Although the shared source and the individual ones can be sparsely decomposed in the same Φ, atoms used to represent them will exhibit different kinds of sparsity. In a recent paper, [START_REF] Candès | Robust principal component analysis?[END_REF] surveillance video frames were modeled as the sum of a low-rank and a sparse matrix. In a similar fashion, we can decompose C X in a sum of two components: C X = B X + P X where B X is a structured sparse matrix and P X an unstructured sparse matrix. B X has a small number of columns of non zero elements. Each column indicates an atom that is spread among all mixtures. P X on the opposite, contains at least one zero per column. Its non-zero elements denote atoms that only belong to a subset of mixtures, hence to a subset of individual sources.

The interest of such model is obvious for source separation, the shared source can be modeled as X c = B X • Φ and the individual sources are the rows in the product P X • Φ. We can rewrite the problem so as to take this two-sparsities model into account:

C X = arg min f (X -C X • Φ) + λ • B X p,q + γ • P X p ′ ,q ′
(3) which allows to put different constraints (by means of λ, p, q and γ, p ′ , q ′ ) on the matrices according to the desired sparsities for B X and P X . We could have designed a pseudoconvex optimization algorithm to specifically solve this problem (in the spirit of the Principal Component Pursuit proposed in [START_REF] Candès | Robust principal component analysis?[END_REF]), however these algorithms are computationally intensive and memory consuming. In order to process real scale audio data, we propose a simple greedy algorithm.

JOINTLY ADAPTIVE MATCHING PURSUIT

We propose the use of a fast greedy algorithm of the Matching Pursuit [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF] family to find (potentially suboptimal) solutions to [START_REF] Gribonval | Average case analysis of multichannel thresholding[END_REF]. The separation could be addressed in a post-processing step, for example by clustering the selected atoms according to their projections across mixtures. However, we have found that much better results can be obtained when the separation process is integrated in the greedy algorithm. This integration takes the form of two modifications of the basic algorithm. These changes are: i) the atom selection criterion has been changed, and ii) after an atom is chosen, a decision mechanism is added, attributing it either to the shared source or to the (or multiple) individual sources.

Structure

A matrix of residuals is initialized from the matrix of mixtures R 0 = X. The algorithm iteratively builds the two matrices B X and P X by selecting an atom in Φ according to a criterion C(Φ, R n ). Then a decision is taken whether to attribute the selected atom to the shared source or to a subset of individual sources.

Algorithm 1 Jointly Adaptive Matching Pursuit (JAMP)

Input: X , Φ 1: R 0 := X , n = 0 2: repeat 3:

Step

1 : Select atom φ k ← C(Φ, R n ) 4:
Step 2 : Decide if φ k is background or not 5:

if φ k is background then 6: ∀i, B X [i, k] = φ k , R n i 7: else 8:
Find which channels J ⊂ I, φ k belongs to. Step 3 : Update residual :

R n = X -(B X • Φ + P X • Φ) n ← n + 1 12:
until a stopping condition is met Output: R n , B X and P X

STEP 1 : Atom Selection

For the sake of clarity, we denote r n i (φ) the squared absolute value of the projection of an atom φ onto R n i , the residual of the i-th mixture at the n-th iteration: i.e r n i (φ) = | R n i , φ | 2 . Four criteria have been investigated in this work:

C S (Φ, R n ) = arg max φ∈Φ I-1 i=0 r n i (φ) C M (Φ, R n ) = arg max φ∈Φ min i r n i (φ) C W (Φ, R n ) = arg max φ∈Φ w(φ, R n ) • I-1 i=0
r n i (φ)

C P (Φ, R n ) = arg max φ∈Φ I-1 i=0 r n i (φ) + i =j |r n i (φ) -r n j (φ)|
C S is simply an energetic criterion, it does not influence the choice of an atom from the background or the foreground. C M is a criterion that minimizes the risk of selecting an atom not belonging to the background. C W is a weighted variant, the weight being defined by the spectral flatness of the distribution of the atom projections on the I residuals. This flatness is the ratio of the geometric mean over the arithmetic mean. This criterion penalizes the selection of an atom if it does not belong to the background. Finally C P encourages the selection of atom from the individual sources by adding the inter-channel atom projection differences to the plain energetic criterion.

STEP 2 : Separating the sources

The decision making obviously depends on the chosen criterion. Using C M , no atoms from the individual sources should be selected (at least until the background has been approximated to a good precision), thus no specific mechanism is required. Using any other criterion, on the other hand, forces us to add an additional step. Let φ k be the chosen atom, the distribution of the r n i (φ k ) is informative. If φ k is efficient in representing the background, then this distribution will be flat (e.g. the r n i (φ k ) have small empirical variance). On the opposite, if there are great disparities in the values of r n i (φ k ), then one can assume that φ k should not be assigned to the background, but to a subset (potentially only one) of the individual sources.

Any statistical measure of the dispersion of the r n i (φ k ) values can thus be used. In this work we have used a simple relative standard deviation D = σ µ . This value is low when the dispersion is weak, thus a threshold τ can be defined so as to make a decision on the atom appertaining to the background or the individual sources. In this work, we have set τ = 0.5 and have not tried to optimize this parameter.

EXPERIMENTS

Comparing Selection criteria

To evaluate separation performances of the various selection criteria we have designed the following experiment. 4 short audio excerpt (5 seconds) were used to create 12 sets of 3 mixtures. For each set, one of the excerpts is used as the background source and is present in all the mixtures without distortion. Three different Individual-to-Shared source energy ratios were used, namely 5, 0 and -5 dB, so that a variety of mixing situations are tested. Performance is assessed by means of the widely adopted measures presented in [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. Table 1 gives the results in terms of Sources-to-distortion ratio (SDR) and Sources-to-Interferences ratio (SIR).

The C M criteria is used without any decision mechanism, the foreground sources are being estimated directly from the residual since only atoms from the background should be selected by the algorithm. We can see that this method gives substantially lower results. Interestingly the best technique appears to be using C P .

Comparing with Time-Frequency based Techniques

In [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF], the authors have presented a simple separation technique based on an estimation of a Time-Frequency mask for the background source as the median (respectively the minimum) of the mixture spectrogram. We have implemented this method, labeled REPET-Median (resp. REPET-Min) and [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF] with Joint Matching Pursuit with criterion Penalized for various offsets in background alignments compared it to our own. Since JAMP is iterative, we can follow the evolution of the separation scores through the decomposition process. Figure 1 presents such results for analyzing mean performances on the same set of signals as above. JAMP with a C P selection criterion can reach the same SDR level than REPET-Median in about 10000 iterations. In average, REPET-Min gives better SDR results, however the mean SIR values are much better using JAMP.

Additionally, the JAMP algorithm is designed to be more robust to the backgrounds being offset in the mixtures. Actually, a local optimization of atoms time localization for each mixture is performed. It effectively manages to reduce preecho artifacts [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]. Figure 2 presents the results of an experiment in which the background sources are offset in each mixture. Performances of the REPET-Min algorithm drops quite sharply for offsets of about 15ms, while they remain quite unchanged until 150ms for JAMP.

Real audio data

The real difficulty arises when the background is not perfectly identical (e.g when considering a repeated musical pattern). The experiment here consists of separating the singing voice from a repeated musical background. Due to variation in the execution, the background is not exactly the same nor perfectly aligned since tempo variation can occur, which makes As in [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF], we use audio material from the Beach Boys. We have cut musical excerpts in 5 songs by hand and constituted 5 sets of 4 mixtures. In each set, all the mixtures are occurrences of the same repeated pattern (usually from the verse) and last a few seconds (from 3 to 6). We have compared the results of the JAMP algorithm (run for 10000 iterations) with REPET using the min and median methods. We have also compared performances when using only I = 3 occurrences for the separation. Results are summarized in Table 2. Performances are quite comparable with a slight advantage to JAMP on the singing voice SDRs and SARs. For the background, REPET-Min gives the best scores except for SIRs where JAMP is clearly ahead.

These results are encouraging. While using a single generic dictionary, JAMP manages to sparsely decompose both the shared source and the individual components. Perceptively though, JAMP creates ringing artifacts, but those could be reduced (e.g. by pre-echo control methods [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]). Increasing the number of iteration leads to an increase of all JAMP scores but the Musical Background SIRs.

BEYOND SOURCE SEPARATION

The simultaneous approximation problem (i.e. finding good joint approximations of the signals) appears disjoint from the source separation problem addressed above. Actually, the global reconstruction error is not intuitively linked to the source separation performances.

With the synthetic dataset described in 4.2, we have found that the criteria C P gave the best separation performances. Figure 3 shows that it is also the one that minimizes the reconstruction error ǫ:

ǫ = 10 log X -C X Φ 2 F X 2 F
where . 2 F is the squared Frobenius matrix norm or the sum of the squares of the entries in the matrix. This comes as a surprise since source separation and error minimization could [START_REF] Sundman | Greedy pursuits of compressed sensing of jointly sparse signal[END_REF] support this prospective. The source separation would then be a nice additional feature of the compression.

CONCLUSION

The joint modeling of the shared source and the individual ones accounts to modeling the redundant and the nonredundant parts of the signal. Although further theoretical studies must be conducted on the matter, it is worth noticing that efficiently separating those parts enables the compression of the redundant parts, but also succeed in minimizing a global reconstruction error on the original mixtures. For musical signals, one cannot always make the assumption that those parts have sparse expansions on different dictionaries. The proposed method overcomes this limitation. JAMP is a simple, fast pursuit algorithm. Hence, it provides an interesting alternative to existing methods in the context of musical repeated pattern separation. Artifacts reduction should be the next matter of concern, and future work will try to embed the model on a broader signal structuration scheme.

9 :

 9 ∀j ∈ J, P X [j, k] = φ k , R

Fig. 1 .Fig. 2 .

 12 Fig. 1. BSS eval mean score for synthetic examples. Comparing TF masking technique from [5] with two Joint Matching Pursuit with criteria Weighted and Penalized

Fig. 3 .

 3 Fig.3. Mean normalized reconstruction error for various criteria. The criteria C P that maximizes the separation performances is also the one that minimizes the reconstruction error. it a difficult task.As in[START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF], we use audio material from the Beach Boys. We have cut musical excerpts in 5 songs by hand and constituted 5 sets of 4 mixtures. In each set, all the mixtures are occurrences of the same repeated pattern (usually from the verse) and last a few seconds (from 3 to 6). We have compared the results of the JAMP algorithm (run for 10000 iterations) with REPET using the min and median methods. We have also compared performances when using only I = 3 occurrences for the separation. Results are summarized in Table2. Performances are quite comparable with a slight advantage to JAMP on the singing voice SDRs and SARs. For the background, REPET-Min gives the best scores except for SIRs where JAMP is clearly ahead.These results are encouraging. While using a single generic dictionary, JAMP manages to sparsely decompose both the shared source and the individual components. Perceptively though, JAMP creates ringing artifacts, but those could be reduced (e.g. by pre-echo control methods[START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]). Increasing the number of iteration leads to an increase of all JAMP scores but the Musical Background SIRs.

Table 1 .

 1 Separation scores (mean and std) after 1000 iterations of JAMP with various selection criteria

		Background	Foreground
		SDR	SIR	SDR	SIR
	C S	6.0 (1.4) 35.8 (8.6) 6.2 (1.5) 17.6 (4.6)
	C M	1.2 (0.9) 16.9 (6.4) 1.2 (1.6) 2.8 (2.0)
	C W	6.0 (1.4) 34.2 (7.7) 6.6 (1.8) 17.3 (4.8)
	C P	7.8 (1.8) 35.0 (6.3) 7.1 (1.7) 20.4 (5.5)

Table 2 .

 2 Separation scores on repeating musical segments from the Beach Boys. JAMP stopped after 10000 iterations have been antagonistic optimization goals. The fact that the proposed method optimizes both objectives opens interesting perspectives. Minimizing ǫ is a desirable property in a compression context, hence, this work could be embedded in a broader distributed source coding scheme. Recent work on distributed compressive sampling

			3 Versions			4 Versions
	Method	SDR (dB)	SIR (dB)	SAR (dB)	SDR (dB)	SIR (dB)	SAR (dB)
				Musical Background
	REPET-Min 3.16 ± 1.7 3.41 ± 5.8 10.03 ± 1.9 3.47 ± 1.2 3.29 ± 4.7 11.23 ± 2.0
	REPET-Med 2.49 ± 0.6 8.08 ± 6.4	3.28 ± 1.5	2.62 ± 0.7 7.61 ± 6.3	4.23 ± 1.6
	JAMP -C P	1.96 ± 0.6 19.14 ± 7.2 -0.87 ± 2.2 2.06 ± 0.6 17.42 ± 6.0 -0.60 ± 2.3
				Singing Voice
	REPET-Min 1.67 ± 0.9 9.96 ± 3.2	0.25 ± 3.0	1.39 ± 0.7 11.17 ± 3.1 -0.55 ± 2.4
	REPET-Med 2.91 ± 0.6 5.47 ± 2.7	4.71 ± 1.8	2.92 ± 0.4 5.40 ± 2.2	4.96 ± 1.5
	JAMP -C P	3.62 ± 0.8 5.94 ± 2.5	5.21 ± 1.8	3.48 ± 1.0 6.03 ± 2.5	4.79 ± 2.3

20th European Signal Processing Conference (EUSIPCO 2012) Bucharest, Romania, August 27 -31, 2012 © EURASIP, 2012 -ISSN 2076-1465

The l 0 pseudo-norm X 0 counts the number of nonzero entries of X.

See[START_REF] Kowalski | Beyond the narrowband approximation: Wideband convex methods for under-determined reverberant audio source separation[END_REF] for proper definition


Valmi Dufour-Lussier

Florence Le Ber

Jean Lieber

Thomas Meilender

Emmanuel Nauer

Semi-automatic annotation process for procedural texts: An application on cooking recipes

Keywords: cooking, natural language processing, procedural texts, semantic annotation, semantic wiki

Taaable is a case-based reasoning system that adapts cooking recipes to user constraints. Within it, the preparation part of recipes is formalised as a graph. This graph is a semantic representation of the sequence of instructions composing the cooking process and is used to compute the procedure adaptation, conjointly with the textual adaptation. It is composed of cooking actions and ingredients, among others, represented as vertices, and semantic relations between those, shown as arcs, and is built automatically thanks to natural language processing.

The results of the automatic annotation process is often a disconnected graph, representing an incomplete annotation, or may contain errors. Therefore, a validating and correcting step is required. In this paper, we present an existing graphic tool named KCatoS, conceived for representing and editing decision trees, and show how it has been adapted and integrated in WikiTaaable, the semantic wiki in which the knowledge used by Taaable is stored. This interface provides the wiki users with a way to correct the case representation of the cooking process, improving at the same time the quality of the knowledge about cooking procedures stored in Wiki-Taaable.

Introduction

This paper presents how an automatic textual annotation process of procedural texts, like cooking recipes, can be improved, using a graphical interface plugin in a wiki and by involving wiki users for correcting and completing the result of the automatic annotation. This work has been done in the framework of Taaable, a case-based reasoning (CBR) system which adapts cooking recipes to user constraints. According to the user constraints, Taaable looks up, in the recipe base, whether some recipes satisfy these constraints. Such recipes, if they exist, are returned to the user; otherwise the system is able to retrieve similar recipes (i.e. recipes that match the target query partially), called source cases in the context of a CBR application, and to adapt these recipes, creating new ones. The knowledge required by Taaable is stored in semantic wiki named WikiTaaable. Adaptation consists first in substituting some ingredients of the source cases by the ones required by the user. Then, two additional adaptations are computed: the adaptation of ingredient quantities and the adaptation of the text of the preparation procedure [START_REF] Blansché | Taaable 3: adaptation of ingredient quantities and of textual preparations[END_REF]. This last adaptation requires to represent semantically the sequence of instructions composing the cooking process. Previous work, based on natural language processing (NLP), has been realised in order to transform the textual procedure into a semantic annotation automatically. This semantic annotation takes the form of a directed graph in which cooking actions and ingredients, among others, are represented as vertices, and the semantic relations between those are represented as arcs.

The main objective of this work is to make it possible for wiki users to edit the graph that is automatically generated for each recipe in WikiTaaable, in order to correct and eventually complete it. In this paper, we show how KCatoS, an existing graphic tool designed to display and edit decision trees, has been adapted and integrated in WikiTaaable.

The paper is organised as follows: Section 2 specifies the problem in its whole context and introduces Taaable and WikiTaaable. Section 3 introduces the graph representation used for recipes, the NLP methods used to extract those graphs, and textual adaptation. Section 4 explains our approach for editing and correcting the graphs using KCatoS. Section 5 concludes the paper.

Context and motivation

Taaable is a CBR system that has been designed to take part in the Computer Cooking Contest6 , an international contest which aims at comparing CBR system results on a common domain: cooking. Several challenges are proposed in this contest. Among them, two challenges (won by Taaable in 2010) are of specific interest to this work:

• the main challenge, in which CBR systems must return recipes satisfying a set of constraints given by the user, such as inclusion or rejection of ingredients, the type or the origin of the dish, or its compatibility with some diets (vegetarian, nut-free, etc.). For example: a user may ask for "a dessert, with rice and fig", as illustrated in Fig. 1. Systems have to search into a limited set of (ca. 1500) recipes for recipes satisfying the constraints, and if there is no recipe satisfying all the constraints, the systems have to adapt existing recipes into new ones. • the adaptation challenge, in which CBR systems must adapt a given recipe to specific constraints. For example, "adapt the 'My strawberry pie' recipe because I do not have strawberries". Like many CBR systems [START_REF] Riesbeck | Inside Case-Based Reasoning[END_REF], Taaable uses an ontology to retrieve the source cases that are the most similar to a target case (i.e. the query). Taaable retrieves and creates cooking recipes by adaptation. According to the user constraints, the system looks up, in the recipe base (which is a case base), whether some recipes satisfy these constraints. Recipes, if they exist, are returned to the user; otherwise the system is able to retrieve similar recipes (i.e. recipes that match the target query partially) and adapts these recipes, creating new ones. Searching similar recipes is guided by several ontologies, i.e. hierarchies of classes (ingredient hierarchy, dish type hierarchy, etc.), in order to relax constraints by generalising the user query. The goal is to find the most specific generalisation (with the minimal cost) for which recipes exist in the case base. Adaptation consists in substituting some ingredients of the source cases by the ones required by the user.

Principles of Taaable

To deal with the adaptation of a specific recipe (which is the adaptation challenge problem), Taaable uses the same hierarchy based generalisation/specialisation mechanism on a recipe base containing only the recipe that has to be adapted. For example, when adapting the "My Strawberry Pie" recipe (in which strawberries are required) to the constraint "no strawberry", Strawberry is generalised on Berry, which is then specialised in another berry (Raspberry, Blueberry, Blackberry, etc.). Substitutions (e.g substitute Strawberry by Raspberry) are proposed to the user.

WikiTaaable

WikiTaaable7 is a semantic wiki that uses Semantic Me-diaWiki [START_REF] Krötzsch | Semantic mediawiki[END_REF] as support for encoding knowledge associated to wiki pages. WikiTaaable contains the set of resources required by the Taaable reasoning system, in particular an ontology of the domain of cooking, and recipes.

The cooking ontology is composed of 6 hierarchies: a food hierarchy (related to ingredients used in recipes, e.g. Berry, Meat, etc.), a dish type hierarchy (related to the types of recipes, e.g. PieDish, Salad, etc.), a dish moment hierarchy (related to the time for eating a dish, e.g. Snack, Starter, Dessert, etc.), a location hierarchy (related to the origins of recipes, e.g. France, Asia, etc.), a diet hierarchy (related to food allowed or not for a specific diet, e.g Vegetarian, NutFree, etc.), an action hierarchy (related to cooking actions used for preparing ingredients, ToCut, ToPeel, etc.). In the semantic wiki, each concept of a hierarchy is encoded as a category page Category: <concept name>. Each concept is described by a short description, lexical variants (used by the annotation bot, for searching concepts in the full text of recipes), its sub-categories and super-categories. For berries, the wiki page indicates that Berry is a sub-concept of Fruit (corresponding to the Category:Fruit page), and sub-categories of Berry (e.g. Raspberry, Blueberry, etc.) are listed. Each action has properties, both syntactic and semantic, associated to it, which makes the automatic case acquisition process described hereafter possible.

The set of recipes contained in WikiTaaable are those provided by the Computer Cooking Contest, that have been semantically annotated according to the domain ontology. Each recipe, as the one given in the example of Fig. 2, is encoded as a wiki page, composed of several sections: a title, which is the name of the recipe, an "Ingredients" section containing the list of ingredients used in the recipe, each ingredient being linked to its corresponding Category page in the food hierarchy, a "Textual Preparation" section describing the preparation process, some possible "substitutions" which are adaptation knowledge, and "other information" like the dish type, for example.

Additionally, for each recipe, there is a graph representing formally the preparation process, which constitutes a highly structured case representation usable by a CBR engine. For instance, Taaable uses this to propose a fully adapted recipe text to the user. Previously, WikiTaaable used a simple tree representation as shown in Fig. 3. This representation is limited and has not been designed to be edited, which is why it is being replaced by the more complete yet easier to edit representation shown in Fig. 4. Initially, this graph is acquired from the ingredients and the textual preparation parts through an automatic process making use of NLP methods, which will be presented in the next section.

However, NLP never gives perfect results, and this case acquisition application is no exception. While efforts are being invested in making the acquisition more accurate, if higher quality cases are required in the short term, user intervention is required.

An interface is proposed to make it possible for users to edit the graphs. Graph edition is accomplished in interaction with the automatic acquisition application, such that for each change entered by a user, the system is able automatically to suggest further changes that seem to be required.

Increasing the quantity and quality of the worklfows in the wiki will improve the results of the reasoning system.

The next section details our approach for automatic case acquisition and illustrates how the case representations are being used currently in Taaable.

Textual CBR adaptation 3.1 Case representation

Two important methods that are used in CBR to represent processes, including recipes, are workflows [START_REF] Minor | Towards case-based adaptation of workflows[END_REF] and qualitative algebras [START_REF] Dufour-Lussier | Adapting spatial and temporal cases[END_REF]. A tree formalism presented earlier was also used in [START_REF] Dufour-Lussier | Text adaptation using formal concept analysis[END_REF]. It is important that the representation of a recipe in Wi-kiTaaable be general enough so that it could be translated to any of those formalisms.

To make semi-automatic annotation by users easier, the representation that is used is presented as a graph. However, each recipe graph is directly translatable to an equivalent knowledge base expressed for instance in description logics. In recipe graphs, each vertex represents the instantiation of a WikiTaaable category, for instance an action or a food item. An arc expresses a semantic link between two vertices. For one thing, each action vertex is related to a certain amount of food vertices through relations of type hasDOInput and/or hasPCInput, and to (generally) one food vertex through the relation hasOutput. hasDOInput refers to an input linguistically specified as a direct object, and hasPCInput as a prepositional complement: e.g. in "add milk to the batter", "milk" is a direct object and "batter" is a prepositional complement. Action vertices are also related to each others through relations such as isBefore or isDuring.

Normally, the preparation of a recipe comes to a point where all ingredients are mixed together to form a new whole. This indicates that a correct annotation of a recipe will normally be a connected graph.

NLP for procedural case acquisition

The case acquisition process proposed uses a combination of classical NLP tools which have been slightly adapted to give better results within the framework of procedural texts, as well as some new ideas that were implemented specifically for this type of texts. Other approaches are possible, e.g. [START_REF] Schumacher | Extraction of procedural knowledge from the web: A comparison of two workflow extraction approaches[END_REF] discusses two different methods based on information extraction to acquire cases as workflows.

A text is first tokenized (split in words), in order to make further linguistic analysis possible. Then, a part-of-speech tag are assigned to each word, indicating whether it is a verb, a noun, etc. This makes is possible to analyse sentences syntax to physically identify, for instance, actions and their arguments. For each action, an action vertex is added to the graph, as well as one (exceptionally, more than one) new food vertex corresponding to its output.

Those steps are not trivial, but effective tools exist, which are either based on stochastic machine learning (e.g. the Brill tagger [START_REF] Brill | A simple rule-based part of speech tagger[END_REF] used for part-of-speech tagging) or rule-based (e.g. the chunker [START_REF] Abney | Parsing by chunks[END_REF] used with hand-crafted grammar rules for syntactic analysis). But this is not sufficient, e.g. to identify properly which food items an action takes as input.

For instance, in a sentence such as "Peel the mangoes, slice lengthwise and remove the pits", it is easy as humans to understand that mangoes are being sliced and pitted, but some heuristic is needed is order for a computer to realise that. This is where the property in the action hierarchy of Wi-kiTaaable comes in handy. Each action has an arity which makes it possible to tell when an argument is missing. In the sentence above, for instance, "slice" requires a direct object which is missing, and "remove", a prepositional complement. Whenever this happens, it is taken that the missing argument is in fact the output of the last action, in this example, the mangoes. In that way, we are able to deal with anaphora, the phenomenon wherein a different word, or no word at all as it might be, is used to represent an object.

Other types of anaphora appear in recipes. The expression "seasonings ingredients" clearly refer to some set of food components, so the ingredient hierarchy is used to find all the nodes of the tree that fit under the "seasonings" category. A phrase such as "cover with sauce" is trickier because there is no obvious clue either in the text or in the ontology as to which food the word "sauce" may refer to. We built, from the analysis of thousands of recipes, "target sets" of ingredients that usually appear in the food components being referred to by word such as "sauce" or, say, "batter". This allows for probabilist association of the word with the proper food item with respect to the ingredients of this food.

Textual adaptation

One thing having a high structured case representation of recipes is useful for is to suggest to the user an adaptation at the formal representation level, or even at the text level.

For instance, in [START_REF] Dufour-Lussier | Text adaptation using formal concept analysis[END_REF], an adaptation method that consists in replacing a sequence of actions applicable to an ingredient α with a sequence more suited to a substitution ingredient β taken from another recipe was introduced. Compared with, for instance, complex and error-prone text generation, this method is economical in terms of textual adaptation, because it makes it possible to prune a piece of text and replace it with a piece from a different text, modulo minor linguistic adjustments (fixing capitalisation and punctuation).

Whatever the case representation formalism, the effect of any adaptation method will be to modify the formal representation of the retrieved case to make it suitable as a solution to the target problem. If the solution is to be presented to the user in text format, it is necessary that these modifications be applied to the text at the same time as to its formalised form.

This requires some annotation in the text to establish a mapping between the text and its formal representation. Minimally, clause or sentence segmentation must be marked in the text, so that each clause or sentence can be mapped to the action(s) it expresses. In order simply to replace whatever α-specific preparation steps in the recipe with β-specific preparation steps, this is sufficient.

Fig. 5 shows an example of this adaptation method: the tree branch corresponding to mangoes is being pruned at the arc shown in dotted style, and a fig branch from another recipe is being grated in its place. The corresponding adaptation is done in text at the same time. The result of the textual adaptation appears in Fig. 1.

In the new representation format used in WikiTaaable, the text-formal representation mapping is materialised in the graph by one vertex associated each clause from the text. Each action vertex is connected to exactly one clause vertex with an arc labelled with isRelatedToClause.

Editing graphs using KCatoS 4.1 KCatoS

Many processes or decision makings can be represented using decision trees. Unfortunately, Semantic Mediawiki does not provide any decision tree editor. That is why KCatoS has been created. KCatoS is a semantic decision tree editor, which provides a collaborative tool to simplify knowledge acquisition. Using a simple graphical language, KCatoS allows exporting decision trees to formalised knowledge, by proposing an original algorithm export to OWL [START_REF] Meilender | Semantic decision trees editing for decision support with KCATOS -Application in oncology[END_REF].

KCatoS is initially a part of a larger work about collaborative editing of clinical guidelines. KCatoS has been created as a SMW extension for Oncologik8 [START_REF] Meilender | From web 1.0 to social semantic web: Lessons learnt from a migration to a medical semantic wiki[END_REF], a semantic wiki that shares clinical practice guideline in oncology. Indeed, to represent decision making, guidelines use visual representations that can mostly be viewed as decision trees from which a meaning can be extracted.

KCatoS proposes various features. Among these features, a syntactic module can be used to check if the edited tree respects decision tree rules. Included in the interface, the module allows to validate trees step by step while drawing, by identifying shapes with mistakes. As an output, different formats are proposed: bitmap (PNG and JPG), vector graphics (SVG), and ontologies (OWL). Moreover, KCatoS includes its own versioning systems. As each tree is kept on a distant server, modifications are saved. Currently, only few functions dealing with history are available: previous versions of a tree can be viewed and restored with some information about authors and dates. However, some pieces of information are saved into XML files that will allow to add functionalities such as the comparison of versions and merging algorithms. Those improvements are planned to be integrated at a future time.

KCatoS decision tree editor is a web-based application using Google Web Toolkit9 (GWT) that allows to create complex Ajax applications. A few additional APIs dedicated to GWT are used to manage the interface. Drawing capabilities rely on SVG and JavaScript technologies while OWL export is done thanks to OWL API [START_REF] Horridge | The OWL API: A Java API for OWL ontologies[END_REF]. Thus, KCatoS is open to collaborative work and web services. Its framework can be integrated in most of content management systems.

KCatoS for recipes

KCatoS cannot be used without modifications in Wiki-Taaable. One reason is the size of the graphs. KCatoS was made to edit small decision trees, whereas a recipe with i ingredients and a actions will have at least 3a + i vertices (each action normally has one for vertex itself, one for its output, and one for its clause).

Is is therefore necessary to present the users with features designed to handle the size and complexity of the graphs, such as the "intelligent zoom" shown in Fig. 6. When the user selects an action for which the annotation needs to be corrected, the editor can automatically centre on this action and show only relevant vertices, such as the foods known to be available for use at the time when the action is executed. Another difficulty is that one simple error in the case acquisition process can result in numerous errors down the road. For instance, considering that ingredient names are not usually repeated more than once in the text, if one input for an action was missed, it will remain missing from the input of all the following actions. Correcting all these mistakes is bound to be a time-consuming and error-prone process. For this reason, KCatoS was modified to allow for semi-automatic annotation. The user is asked to obey to the text order when correcting mistakes. In exchange, the application is capable of using the human-validated part of the annotation as an additional input and propose a new representation of the recipe, hopefully correcting all the mistakes that were caused by the original error corrected by the user.

By looking closely at the graph of Fig. 4, it can be seen that the "molten" action is missing. This causes the vertices for the ingredient "butter" and the clause c 3 (visible in Fig. 6) to be isolated. This further causes butter to be missing from the input of the action "cream", and from the remainder of the recipe. If the user simply adds the missing action and connecting it to the correct input and clause, as shown in Fig. 7, the system is able to correctly repair the rest of the graph.

Conclusion

In order to adapt recipes properly, Taaable uses a highly structured case representation for cooking recipes. This representation is extracted automatically from recipe texts, using an application based on natural language processing methods, but the results are not perfect. This paper shows how a graph editing application for semantic wikis, KCatoS, can be integrated in WikiTaaable, the semantic wiki in which the knowledge required by Taaable is stored, in order to provide wiki users with a way to correct the case representation of the cooking process, thus improving the quality of the knowledge about cooking procedures stored in WikiTaaable. Making those tools available to the public will make it possible for us to gain feedback about the usability of our proposal and its ability to generate correct annotations.

Figure 1 :

 1 Figure 1: The TAAABLE interface. Queried for a dessert dish, with rice and fig, TAAABLE proposes to replace mango by fig in the "Glutinous rice with mangoes" recipe. After viewing the adapted recipe, the user can give feedback about the substitution ("OK" or "not OK").

Figure 3 :

 3 Figure 3: Simple tree representation of the recipe shown in Fig. 2.

Figure 4 :

 4 Figure 4: Recipe graph corresponding to the preparation of the recipe given Fig. 2 in KCatoS (split horizontally for readability).

Figure 2 :

 2 Figure 2: Example of a WIKITAAABLE recipe.

Figure 6 :

 6 Figure 6: "Intelligent zoom" on vertex Action:cream 4.

Figure 5 :

 5 Figure 5: Adaptation by action sequence replacement: a branch related to mango preparation is pruned (dashed arc) and a branch related to fig preparation is grafted (bold arc).

Figure 7 :

 7 Figure 7: Adding a missing action in KCatoS.

http://computercookingcontest.net/

http://wikitaaable.loria.fr

http://www.oncologik.fr

http://code.google.com/intl/en/webtoolkit/

Acknowledgements

KCatoS has been conceived and developed in the context of Kasimir project, gathering the company A2ZI, the health network Oncolor, and Université de Lorraine.


João V F Lima
email: joao.lima@imag.fr

Thierry Gautier
email: thierry.gautier@inrialpes.fr

Nicolas Maillard
email: nicolas@inf.ufrgs.br

Vincent Danjean
email: vincent.danjean@imag.fr

Exploiting Concurrent GPU Operations for Efficient Work Stealing on Multi-GPUs

The race for Exascale computing has naturally led the current technologies to converge to multi-CPU/multi-GPU computers, based on thousands of CPUs and GPUs interconnected by PCI-Express buses or interconnection networks. To exploit this high computing power, programmers have to solve the issue of scheduling parallel programs on hybrid architectures. And, since the performance of a GPU increases at a much faster rate than the throughput of a PCI bus, data transfers must be managed efficiently by the scheduler.

This paper targets multi-GPU compute nodes, where several GPUs are connected to the same machine. To overcome the data transfer limitations on such platforms, the available softwares compute, usually before the execution, a mapping of the tasks that respects their dependencies and minimizes the global data transfers. Such an approach is too rigid and it cannot adapt the execution to possible variations of the system or to the application's load.

We propose a solution that is orthogonal to the above mentioned: extensions of the XKaapi software stack that enable to exploit full performance of a multi-GPUs system through asynchronous GPU tasks. XKaapi schedules tasks by using a standard Work Stealing algorithm and the runtime efficiently exploits concurrent GPU operations. The runtime extensions make it possible to overlap the data transfers and the task executions on current generation of GPUs. We demonstrate that the overlapping capability is at least as important as computing a scheduling decision to reduce completion time of a parallel program.

Our experiments on two dense linear algebra problems (Matrix Product and Cholesky factorization) show that our solution is highly competitive with other softwares based on static scheduling. Moreover, we are able to sustain the peak performance (≈ 310 GFlop/s) on DGEMM, even for matrices that cannot be stored entirely in one GPU memory. With eight GPUs, we archive a speed-up of 6.74 with respect to single-GPU. The performance of our Cholesky factorization, with more complex dependencies between tasks, outperforms the state of the art single-GPU MAGMA code.

I. INTRODUCTION

Current processors are homogeneous chips containing many cores, whose number will increase again in the near future. The architectural trend is the emergence of hybrid systems with many tightly coupled processing units (PU) such as GPUs or accelerators (old Cell BE/Intel MIC) for high performance computing. During the last year, the number of GPU-based computers in the Top500 list has increased from 3.8% to 11% (from June, 2011 to June, 2012).

Such architectures have heterogeneous PUs in terms of computing power and programming model. The programmer is in charge of execution flow and memory consistency. Hence, in many cases, familiar algorithms need to be redesigned.

Algorithms in dense linear algebra, such as those found in the LAPACK library, and especially matrix factorizations, have already been redesigned to exploit multicore machines. The FLAME [START_REF] Quintana-Ortí | Solving dense linear systems on platforms with multiple hardware accelerators[END_REF] and PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] projects have demonstrated how much more interesting it is to exploit parallelism among the BLAS operations, than inside a given BLAS operation itself. These new algorithms are built on a software stack that allows to describe tasks with dependencies and to schedule them at runtime on multicore.

With hybrid architectures, this software has been extended to develop hybrid algorithms with multiple task implementations optimized for each kind of PU. MAGMA [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF] allows to exploit one GPU; MAGMA/StarPU [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF] reports experiments up to four GPUs; FLAME [START_REF] Quintana-Ortí | Solving dense linear systems on platforms with multiple hardware accelerators[END_REF] also shows experiments on four GPUs. The first main difference with previous homogeneous architecture is that heterogeneous machine has non-uniform processing power elements. The scheduling theory community considers a heterogeneous machine as an unrelated machine where scheduling a task graph to minimize the makespan is a well-known NP hard problem. The second difference is that a heterogeneous architecture introduces a new memory level (on the accelerator), which is non cache-coherent with the main memory. Moreover, the bandwidth of PCI-Express interconnect between the host and the accelerator memory remains low, about 8 GB/s on PCIe x16, vs. the 32 GB/s between a CPU core and the main memory.

To reach high performance, it is essential to reduce the time required to transfer data. A good schedule may take the right decision to map a task onto a GPU resource which already stores data to be reused. A good runtime may also try to exploit the high capability to overlap communication with computation on the modern GPUs. There are many publications that report experiments on multi-GPU systems on some linear algebra factorization algorithms [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF], [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF], [START_REF] Horton | A class of hybrid lapack algorithms for multicore and gpu architectures[END_REF]. The authors report higher performance by computing at runtime a static schedule of the task graph coupled with overlapping strategy during execution. Such performances come at the expense of being able to compute cost models of the task graph and of data transfers.

Nevertheless, if the communication cost can be entirely overlapped, then classical dynamic work stealing [START_REF] Frigo | The implementation of the cilk-5 multithreaded language[END_REF] with the-oretical performance guarantee should be almost as efficient, without requiring cost models anymore, and it can react well to inactivity due to system or application load variations.

This paper presents extensions of the XKaapi runtime that enable to exploit the full performance of multi-GPU systems by using a standard work stealing algorithm. We compare our design decisions with StarPU [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], that computes a static schedule of the task graph using an ingenious automatic tool to build cost models.

Our main contribution is to demonstrate that, on current GPUs, the overlapping capability is at least as important as computing a good scheduling decision to reduce the completion time of a parallel program. Moreover, because XKaapi uses an online scheduling algorithm, it can be performant for a wide range of applications, including those irregular, when the cost of the tasks depends of their input.

First, we evaluate raw performances on matrix product (DGEMM) and Cholesky factorization for single-GPU systems. Second, we study the scalability of our implementation on up to 8 GPUs, and compare it with state-of-the-art libraries. Finally, we complete our initial experiments with various grain and problem sizes, and with the volume of data transfers, in order to estimate the limitation of pure work stealing for multi-GPUs. Our experiments show that we are able to sustain the peak performance on DGEMM, even for matrices bigger than the GPU memory, with a 6.74 speed-up on 8 GPUs (relatively to 1 GPU), which is better than the other tools. For the Cholesky factorization, our XKaapi version outperforms the state of the art MAGMA single-GPU code; it suffers from the same scalability limitations as StarPU on multi-GPUs.

The remainder of this paper is organized as follows: in Sections II and III we describe the XKaapi programming model and the designed runtime extensions to support multi-GPUs. We discuss our performance results in Section IV. In Section V we present related works on runtime tools for GPUs and linear algebra (LA) libraries. We conclude and present future works in Section VI.

II. DATA FLOW TASK PROGRAMMING WITH XKAAPI

The XKaapi 1 task model [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multiprocessors[END_REF], as in Cilk [START_REF] Frigo | The implementation of the cilk-5 multithreaded language[END_REF], Intel TBB [START_REF] Reinders | Intel threading building blocks[END_REF], OpenMP-3.0 [START_REF]OpenMP Architecture Review Board[END_REF] or StarSs [START_REF] Bueno | Productive programming of gpu clusters with ompss[END_REF], [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF], enables non-blocking task creation: the caller creates the task and continues the program execution. The semantic remains sequential like in XKaapi's predecessor Athapascan [START_REF] Galilée | Athapascan-1: On-line building data flow graph in a parallel language[END_REF], but the runtime has been redesigned [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multiprocessors[END_REF] and then specialized for multi-CPUs/multi-GPUs iterative applications [START_REF] Hermann | Multigpu and multi-cpu parallelization for interactive physics simulations[END_REF]. Here we present the extension of our previous work on multi-CPUs to a general scheduling algorithm on multi-GPUs.

A. Design Choices

More than a runtime, XKaapi is a fully featured software stack to program hybrid parallel architectures. The core stack is written in C and was designed using a bottom up approach: each layer is kept as specialized as possible to fit a specific 1 http://kaapi.gforge.inria.fr need. Currently, the stack includes: a runtime supporting multicores and multiprocessors; a set of ABIs (QUARK [START_REF] Yarkhan | Quark users' guide: Queueing and runtime for kernels[END_REF], OpenMP runtime libGOMP [START_REF] Broquedis | libKOMP, an Efficient OpenMP Runtime System for Both Fork-Join and Data Flow Paradigms[END_REF]); a set of high level APIs (C, Fortran and C++); and a source to source compiler [START_REF] Lementec | The X-Kaapi's Application Programming Interface. Part I: Data Flow Programming[END_REF] based on the ROSE compiler framework. Only the C++ API currently supports multi-CPUs/multi-GPUs applications, and we will use it in the code fragments of this paper.

B. Data Flow Task Model

A XKaapi program is composed of sequential code and some annotations or runtime calls to create tasks. The parallelism in XKaapi is explicit, while the detection of synchronizations is implicit [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multiprocessors[END_REF]: the dependencies between tasks and the memory transfers are automatically managed by the runtime.

A task is a function call that returns no value except through the shared memory and the list of its effective parameters. Depending of the APIs, tasks are created using code annotation (#pragma kaapi task directive) if the XKaapi compiler [START_REF] Lementec | The X-Kaapi's Application Programming Interface. Part I: Data Flow Programming[END_REF] is used, or by library function (kaapic_spawn call using XKaapi's C API, or by calling the template function ka::Spawn), or by low level runtime function calls.

Tasks share data if they have access to the same memory region. A memory region is defined as a set of addresses in the process virtual address space. This set has the shape of a multidimensional array. The user is responsible for indicating the mode each task uses to access memory: the main access modes are read, write, reduction or exclusive [START_REF] Galilée | Athapascan-1: On-line building data flow graph in a parallel language[END_REF], [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multiprocessors[END_REF], [START_REF] Lementec | The X-Kaapi's Application Programming Interface. Part I: Data Flow Programming[END_REF]. When required ( [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multiprocessors[END_REF]), the runtime computes true dependencies (Read after Write dependencies) between tasks thanks to the access modes. At the expense of memory copy, the scheduler may solve false dependencies through variable renaming.

A thread creates tasks and pushes them on its own work queue. The work queue is represented as a stack. The enqueue operation is very fast, typically about ten cycles on the last x86/64 processors. As for Cilk, a running XKaapi task can create child tasks, which is not the case for the other data flow programming softwares previously mentioned [START_REF] Yarkhan | Quark users' guide: Queueing and runtime for kernels[END_REF], [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF], [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], except the recent StarSs extension OmpSs [START_REF] Bueno | Productive programming of gpu clusters with ompss[END_REF]. Once a task ends, the runtime executes its children following a Firstin First-out (FIFO) order. During task execution, if a thread finds a stolen task, it suspends its execution and switches to the work stealing scheduler that waits for dependencies to be met before resuming the task. Otherwise, and because sequential execution is a valid order of execution [START_REF] Galilée | Athapascan-1: On-line building data flow graph in a parallel language[END_REF], [START_REF] Gautier | KAAPI: A thread scheduling runtime system for data flow computations on cluster of multiprocessors[END_REF], tasks are performed in FIFO order without computation of data flow dependencies.

C. Blocked Linear Algebra Algorithms with XKaapi

Previous works have shown that fine granularity and asynchronism are keys to unfold parallelism on hybrid architectures [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF], [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF]. The "hybridization" strategy for linear algebra algorithms decomposes factorizations in a collection of BLASbased fine-grained tasks with dependencies among them. The tasks can then be properly scheduled on the available resources.

for ( k=0; k < N; k+= blocsize ){ ka :: Spawn<TaskPOTRF>()( A(rk,rk) ); for ( m=k+blocsize; m < N; m+= blocsize) ka :: Spawn<TaskTRSM>()( A(rk,rk), A(rm,rk) );

for ( m=k+blocsize; m < N; m+= blocsize){ ka :: Spawn<TaskSYRK>()( A(rm,rk), A(rm,rm) ); for ( n=k+blocsize; n < m; n+= blocsize ) ka :: Spawn<TaskGEMM>()( A(rm,rk), A(rn,rk), A(rm,rn) ); } } The code fragment of Figure 1 illustrates how to program a Cholesky factorization using the C++ API. The ka::Spawn<Task> creates a task of type Task. Each parameter rk,rm,rn corresponds to a range of indexes and a construction such as A(rm,rk) represents the sub-matrix of elements A(i,j) where i,j are in the range rm,rk.

Let us note that a sub-matrix does not need to be contiguous in memory as with OmpSs [START_REF] Bueno | Productive programming of gpu clusters with ompss[END_REF]. It allows faster code parallelization, even if better performance may be obtained at the cost of converting row major or column major matrix representation to block representation.

D. Versioning Task Implementations

Extensions to the C++ interface provide a high level interface for multi-versioning a task implementation [START_REF] Hermann | Multigpu and multi-cpu parallelization for interactive physics simulations[END_REF]. First, a task is associated with a signature that includes the task parameters and their access mode (read R and/or write W and/or concurrent write CW). Each CPU or GPU implementation is encapsulated in a functor object, which must respect its task signature. This concept of multi-versioning and task implementation allows a clear separation between the task definition and its implementations. StarPU uses a similar approach called codelet [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF]. Thus, architecture details are abstracted from the algorithm. Moreover, the signature allows the runtime to automatically take care of memory transfers in case of distributed memory. The application programmer does not need to code any explicit memory transfer. Figure 2 shows an example of a task with CPU (TaskBodyCPU) and GPU (TaskBodyGPU) implementations conforming to its Signature. The runtime expects at least one implementation of a task signature.

III. XKAAPI RUNTIME SUPPORT TO MULTI-GPUS

The XKaapi runtime extensions implement a programming model that offers asynchronous execution of GPU tasks and abstracts memory details. Algorithms on top of XKaapi describe the execution flow through the task dependencies and the runtime decides the target resource (CPU or GPU) and performs memory transfers as necessary. Our current version supports NVIDIA CUDA and it relies on the features of recent Fermi GPUs. This section describes the features to support multi-GPUs in XKaapi through scheduling, concurrent GPU operations and memory abstraction.

A. Work Stealing Scheduling Algorithm

The XKaapi version for multicore architectures implements a list scheduling algorithm based on Cilk's Work Stealing and has specific optimizations for fine-grained parallel algorithms, which have been sketched in section II-B. For each used GPU, the runtime launches a thread on the host machine that runs a modified work stealing algorithm.

For each task executed on a GPU, the runtime first transfers the input data and allocates memory on the GPU for the output. The runtime assumes that the GPU task implementation launches the GPU kernels asynchronously. Data transfers and kernel launch on GPU are handled asynchronously by an extension of CUDA streams described in the next section.

In comparison with original multi-CPUs work stealing, multi-GPUs work stealing adds a new state in the task state diagram that corresponds to a task for which input data are under transfer. The GPU thread polls regularly the completion of previous asynchronous GPU operations.

A task that completes its execution, when the asynchronous kernel launch has completed, activates the successor tasks (according to the data flow dependencies) that become ready. These new ready tasks are pushed on the tasks' queue attached to the current GPU and they may be stolen by one CPU or another GPU.

B. Concurrent Operations between CPU and GPU

Recent GPUs such as NVIDIA's Fermi allow new techniques to explore asynchronism in multi-GPU systems. Fermi GPUs have one execution engine and two copy engines capable of concurrent execution and transfers (two-way host-todevice and device-to-host), under the condition that no explicit nor implicit synchronization occurs. This section details how we exploit these capabilities.

XKaapi has an execution strategy for GPUs that avoids CUDA's implicit synchronizations and exploits concurrent memory transfers in the two ways along with kernel execution. It splits the execution of a GPU task in three basic operations: host-to-device input transfers (H2D), TaskBodyGPU execution (i.e. launch of CUDA kernels) (K), and device-to-host output transfers (D2H). Write-only parameters are directly allocated before launching CUDA kernel.

Since concurrency between data transfers and kernel launches must use CUDA streams, we define a new data structure, called kstream, that groups together three types of CUDA streams: a stream for host-to-device transfer, a stream for kernel execution and a stream for device-to-host transfer. The kstream structure allows to insert a request for one of the three types it handles (H2D, K, or D2H). A callback function and its argument can be specified for each request insertion. Moreover, after each request insertion, the kstream inserts a CUDA event to detect the completion of the asynchronous operation. Once the kstream detects the event completion, it calls the callback function with its argument as parameter. It is the responsibility of the client of the kstream structure to regularly poll for the completion of asynchronous requests by calling a specific function.

Figure 3 illustrates the way our kstream structure allows to pipeline concurrent operations on a Fermi GPU. This design allows concurrent execution between CUDA streams of each type. The kstream represents three flows of FIFO ordered GPU operations whose execution are independent from each other. The FIFO order is only respected among operations of the same type (H2D, K or D2H). The callback mechanism permits to compose a sequence of operations and it is typically used by the GPU work stealing algorithm, first to insert data transfers for the input of a task, and then to invoke the kernel launch when the transfer ends.

C. Memory Management

In order to enable asynchronous memory transfers with CUDA, all the user data must register to XKaapi's runtime.

XKaapi manages GPU memory through a software cache, based on the Least Recently Used (LRU) replacement policy. Each GPU thread maintains two FIFO queues in order to keep track of allocated blocks. One queue stores blocks in read-only (RO) mode and the other stores blocks in read-write (RW), write-only (WO) modes or concurrent write (CW) modes.

The first positions of the RO and RW/WO queues contain the blocks more recently accessed, and the last positions the blocks less recently accessed. When a GPU task requires to access a host memory block that is not present on the GPU, the runtime will allocate memory and insert it in one of the two queues, based on its access mode, after it has initiated the data transfer (for data in read access mode). If the GPU memory is full, it verifies first at the end of the RO queue and, then, into the RW/WO queue, respectively, if a memory block bigger or equal than the requested size is not accessed anymore. If possible, unused blocks are reused without being freed. Otherwise, it may free blocks from RO and RW/WO/CW queues as needed. This optimization avoids unnecessary CUDA calls. Furthermore, the use of two queues (RO and RW/WO/CW) ensures that data produced by one task in RW/WO/CW mode remains on the GPU if RO data can be reclaimed.

Consistency is guaranteed by a lazy strategy using a writeback policy. Data transfers to or from GPU occur only when a task accesses data and when the data is in an invalid state in the target address space. This policy avoids unnecessary transfers, unlike write-through policy. All transfer operations are asynchronous and rely on the use of our kstream data structure to signal the completion of operations. In the case of GPUto-GPU transfers, the runtime first performs a transfer deviceto-host from a GPU with a valid copy, followed by a host-todevice transfer to the GPU that owns the task. Since its version 4.1, CUDA includes transfers between two GPUs directly by DMA called Peer-to-Peer Device Access. This feature is available when the function cudaDeviceCanAccessPeer returns true. The current version of XKaapi does not make use of this feature, because of its unpredictable behaviour concerning the GPU copy engines: with it, we could not guarantee the coherence of concurrent data copies, nor their overlapping with a kernel execution.

IV. EXPERIMENTS

The experiments of this section show respectively: XKaapi's ability to overlap the communication with tasks; a comparison of XKaapi's performance with state-of-the art tools, and on Multi-GPUs platforms; an analysis of the granularity and its impact on the performance with XKaapi; and, finally, how XKaapi uses efficient data transfers.

All experiments have been conducted on an hybrid, multi-GPUs system, which is named "Idgraf". Idgraf is composed of two six-core Intel Xeon X5650 CPUs (12 CPU cores total) running at 2.66 GHz with 72 GB of memory. It is enhanced with eight NVIDIA Tesla C2050 GPUs (Fermi architecture) of 448 GPU cores running at 1.15 GHz each (2688 GPU cores total) with 3 GB GDDR5 per GPU (18 GB total). Figure 4 illustrates the hardware topology of Idgraf.

All experiments using up to 4 GPUs always run on GPUs that do not share any PEX 8647 multiplexer in order to reserve the full PCIe 16x bandwidth for each pair GPU/CPU. Experiments using more than 4 GPUs lead to having some Here the Tylesburg-36D is a QPI-PCIe bridge and the PXE8647 is a PCIe switch for two GPUs.

pairs of GPUs share the PCIe 16x between the Tylersburg chip and the PEX multiplexer.

Our experiments use the dense linear algebra problems Matrix Product (DGEMM doing C ← C + AB) and Cholesky Factorization, as PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] (the same algorithms have been implemented in XKaapi). As in [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF], we consider the peak performance from DGEMM (from CUBLAS) as an upper bound on the actual performance in double precision (DP) that one may obtain with GPUs. We used as software environment the compiler GCC 4.4 and CUDA 4.1, and the library ATLAS 3.9.39 for the CPU versions of BLAS and LAPACK. We also used MAGMA 1.1.0 for linear algebra algorithms and StarPU 1.0.1 (with its HEFT scheduling algorithm) as performance references. Each result is a mean of 30 executions, which is enough for StarPU to calibrate its internal cost model. The standard deviation is represented on the graphs.

A. Overlapping Data Transfers with GPU Kernel Executions

This section presents experiments to evaluate the capacity of our design to exploit asynchronous data transfers in concurrence with GPU kernel executions. Our experiment measures the runtime of a block matrix product algorithm. Matrices A and B are decomposed into k × k blocks of size s × s. We devised our implementation such that all the computations are performed on the GPU. Matrix computation is done with double precision, each block-matrix product launches CUDA kernels using the CUBLAS DGEMM routine.

We compare the performances of our XKaapi implementation with the performances obtained by native calls to CUBLAS DGEMM on the whole ks × ks matrices, taking into account the time to copy the matrices to/from the device or not.

Figure 5 illustrates the results of DGEMM with XKaapi; with CUBLAS when the time to copy the matrices is not accounted for (label CUBLAS (no copy)), that measures the GPU peak performance; and with CUBLAS with copy time included. Each measure of the XKaapi runtime includes all the costs of CUDA memory allocations and data transfers. XKaapi experiments always take into account all data transfers, including those required to get final results in main memory. First, the peak performance from CUBLAS nocopy is about 315 GFlop/s for square matrices of dimension 4096. Then, the performance decreases to 293 GFlop/s for bigger matrices.

Our XKaapi version, with blocks of size 1024 and 2048, reaches the GPU peak performance except for small matrices. For matrices bigger than 8192 XKaapi's implementation sustains 309 GFlop/s with block of size 2048, which is more than the GPU peak DGEMM from CUBLAS (293 GFlop/s). This algorithm with block size of about 1024 generates numerous tasks that can be exploited by our runtime to pipeline and overlap data transfers with computation. Our good performances confirm that we are able to overlap a very high percent of the data transfers with GPU kernel.

Moreover, thanks to the XKaapi software cache and to our design to exploit concurrent GPU operations, our blocked DGEMM algorithm sustains its 309 GFlop/s performance peak even after the GPU runs out of memory with matrix orders larger than 10240, which requires ≈ 2.4 GB of device memory out of 3 GB on Tesla C2050 cards.

For small matrices, because the number of tasks remains low, the data transfer is not entirely overlapped by computation. Even in this case, XKaapi presents good results. For instance, the performance of CUBLAS nocopy with matrices of size 2048 is about 312 GFlop/s. Performance drops to 152 GFlop/s if we take into account the data transfers. Our XKaapi DGEMM for this matrix dimension and with block size of 1024 generates 8 tasks for each sub-matrix product, and it reaches 240 GFlop/s, that correspond to 157% of improvement over CUBLAS when data transfers are taken into account.

B. Surpassing CUBLAS GPU Peak

A deeper look on the results presented in Figure 5 shows an interesting phenomenon. For matrix dimension 4096×4096 using block size 512, our blocked DGEMM reaches 247.5 GFlop/s. For bigger matrix dimension, using the same block size, the performance increases up to 271 GFlop/s. However, a simple analysis shows that our blocked DGEMM algorithm performance should be upper bounded by CUBLAS DGEMM, since our implementation only calls it to compute each of its blocks. Besides, the same figure shows that the performance of CUBLAS DGEMM, without taking into account data transfers, is 233 GFlop/s for a 512 × 512 matrix. Thus, the XKaapi version increases the performance by 16%.

We are not certain about the actual factor that improves the performance on the Fermi GPUs. One hypothesis is that the gain is explained by better GPU occupancy, or specific CUDA optimizations for small matrices, which includes new CUBLAS batched * GEMM routines since version 4.1.

C. Cholesky Single-CPU/Single-GPU Figure 6 reports our results using one CPU and one GPU for Cholesky factorization. We compare our work stealing based runtime XKaapi to StarPU [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] and single-GPU MAGMA [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF]. StarPU [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] schedules at runtime the entire task graph using the HEFT static scheduling algorithm. In XKaapi and StarPU, the Cholesky factorization of the diagonal block is sequential and executed on the CPU. The MAGMA version uses a more sophisticated implementation where part of the diagonal block factorization is exported on the GPU. XKaapi and StarPU, with runtime scheduling decisions, outperform MAGMA when the matrix dimension is bigger than 10240. The whole matrix size is 800 MB, and can be stored into the 3 GB of device memory. Only the last matrix of dimension 20480 can not be stored into the GPU memory. The main difference between XKaapi and StarPU vs MAGMA is that MAGMA is unable to exploit parallelism across the main iteration (see code fragment of Figure 1), as is done by XKaapi and StarPU.

For small matrix dimensions (less than 2048), the performance of MAGMA and XKaapi are similar, but StarPU seems to suffer from a higher overhead. XKaapi has a little drop and then reaches the performances of StarPU.

D. Multi-GPUs Results

For the multi-GPUs evaluation, our experiments measure the performance of DGEMM and Cholesky using from one to eight GPUs, with matrix dimension of 16384 and block size of 1024. XKaapi is compared again to StarPU. Figure 7 shows the performance results for DGEMM. XKaapi outperforms StarPU in all cases and attains 2023.14 GFlop/s (or speed-up 6.74 on 8 GPUs with respect to single-GPU).

Figure 8 reports the results for the Cholesky factorization. The algorithm corresponds to the code fragment in Figure 1. In the StarPU and XKaapi programs, all tasks, except the block factorization TaskPOTRF, are performed by a GPU. Unlike the DGEMM case, the Cholesky factorization acceleration, up to eight GPUs, is below the expected: neither XKaapi nor StarPU implementations do scale. StarPU reaches 680.82 GFlop/s (or speed-up 2.94 with respect to single-GPU). Experiments with bigger matrices (up to 20480 × 20480) show the same behaviors. This means that, when using more that 4 GPUs, communications costs can not be neglected.

E. Granularity Impact

In this set of experiments we vary the matrix dimension and the block size for our blocked DGEMM with StarPU and XKaapi. The results of Figure 9 permit to conclude that our runtime, with dynamic work stealing scheduler, has a lower overhead than StarPU. Thus, XKaapi is able to exploit more performance on DGEMM. The overhead is related to the number of tasks in the system: for matrix dimensions greater or equal than 10240 and for block size greater or equal than 1024, the difference between XKaapi and StarPU remains about the same in all the cases.

F. Memory Transfers

XKaapi and StarPU allow to monitor the execution by collecting post-mortem traces of performance counters. We have collected, for one instance of DGEMM and of the Cholesky factorization, the total number of bytes exchanged between the main memory and the GPUs.

Table I shows in GB the total memory transfers with StarPU and XKaapi on the DGEMM code. The matrix dimension is 16384. Surprisingly, StarPU generates bigger data exchanges than XKaapi up to 4 GPUs, although it uses an a priori HEFT algorithm that minimizes data transfers. For a matrix of size S bytes, if the GPU memory could store the entire data, DGEMM implies 3×S bytes from host to device transfer and S bytes of transfer to get back the result. For a DP matrix of dimension 16384 = 2 14 , the data transfer volume is 8 GB. On one GPU, XKaapi's schedules order the tasks according to the sequential order, which is memory efficient. However, StarPU's schedule seems to order the execution using a breadth-first strategy, generating too much parallelism, and thus it increases the memory consumption of the device. Then, the StarPU software cache eviction policy generates a lot of traffic that degrades the performances. On the Cholesky factorization, as illustrated in Table II, the memory transfer generated by XKaapi is larger than StarPU's. As with the DGEMM program, the scheduling strategy used by XKaapi orders the tasks close to the sequential order of their creation, which generates too much parallelism and thus too much data consumption from the device. The HEFT scheduling of StarPU minimizes the makespan and reduces the memory consumption.

DGEMM transfers (GB)

We previously showed that bad scaling of Cholesky factorization exhibited on figure 8 was due to bad overlap of communication by computation. The table II explains why XKaapi performs worst than StarPU: it has more data transfers. Moreover, when using more than 4 GPUs, the architecture share some PCIe 16x links between GPUs (Figure 4), so bottlenecks on data transfers become worst.

Using work stealing directed by data affinity would allow XKaapi to reach StarPU performances. But the amount of data transfers would still be a bottleneck. To overcome this limitation, one would need to use bigger blocks on GPU. Still, using bigger blocks means that TaskPOTRF would become a bottleneck as this task is currently run only on CPU. A parallel implementation, partially on GPU (such as MAGMA), for this task would then be required.

V. RELATED WORKS

OmpSs [START_REF] Bueno | Productive programming of gpu clusters with ompss[END_REF] is a programming tool that provides a set of OpenMP-like pragmas and a runtime system to schedule tasks while preserving dependencies. It offers different scheduling strategies and coherence protocols such as write-back and write-through. To our knowledge, OmpSs has concurrent execution and data transfers in GPUs but it shows some issues with matrix sizes that can not be entirely stored into the GPU memory. It also performs asynchronous operations in GPUs by pinned-memory buffers, which adds additional memory copies and transfer overheads. Moreover, as shown in [START_REF] Bueno | Productive programming of gpu clusters with ompss[END_REF], performance on multi-GPU systems remain difficult to compare with our results: the public downloaded version seems to have some problems on multi-GPUs and the performances reported in the paper achieve about 400 GFlop/s on single precision matrix product while we measured about 625 GFlop/s for the same problem size on similar Fermi GPU processor. The cited CUDA 3.2 version used by the authors is too old for correct comparisons.

StarPU is a runtime system for scheduling a DAG of tasks on hybrid systems optimized for numerical algorithms [START_REF] Augonnet | Starpu: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF]. In a similar way, StarPU provides a programming model for hybrid architectures and exposes an API to describe a scheduling policy which allows flexibility in work distribution. Recently StarPU has supported concurrent write (CW) access of task parameters. It has similar features to XKaapi but StarPU lacks of concurrent operations on GPUs. It uses data prefetch to predict memory transfers before task execution and provides a lazy coherence protocol. However, according to StarPU examples, each GPU task needs a synchronization at the end to ensure all kernels are finished. It forces a synchronization point in the GPU copy engines and does not allow concurrent GPU operations. This issue could not be significant in earlier GPU families, but it is crucial on Fermi GPUs. Besides, its scheduler uses the static HEFT algorithm to schedule the entire DAG thanks to cost models for data transfer and task executions. Such an approach does not allow to react to system load or task execution variations as our work stealing algorithm does.

In the case of linear algebra algorithms, PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] provides fine-grained parallel linear algebra routines with dynamic scheduling through QUARK, which was conceived specially for numerical algorithms. FLAME [START_REF] Quintana-Ortí | Solving dense linear systems on platforms with multiple hardware accelerators[END_REF] is a high-level notation to express algorithms for dense linear algebra operations on multi-CPUs/multi-GPUs. MAGMA [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF] implements static linear algebra algorithms for hybrids systems composed of GPUs. Recently it has included some methods with dynamic scheduling in multi-CPUs and multi-GPUs on top of QUARK or StarPU, in addition to static multi-GPUs version. NVIDIA CUDA Toolkit provides some BLAS routines with CUBLAS [START_REF]Nvidia developer zone[END_REF].

VI. CONCLUSION

We have presented extensions of the XKaapi runtime system to program multi-GPU architectures. The current version provides work stealing that efficiently exploit concurrent GPU operations. Our experimental results have been obtained on blocked matrix product (DGEMM) and the Cholesky factorization, with an hybrid system composed of eight Fermi-based NVIDIA GPUs.

Our main contribution is the design of an asynchronous approach of concurrent GPU operations that achieves an almost ideal overlapping of data transfer and kernel execution with DGEMM algorithm for single-GPU. Thanks to this overlapping, the use of a pure dynamic work stealing algorithm permits to reach high performance as the theory predicts for shared memory machine without communication costs. Thus, overlapping almost enables to hide the heterogeneity of the memory accesses on a multi-GPUs system.

Second, our work stealing implementation outperforms in most cases the static MAGMA approach for Cholesky factorization on single-GPU and the static HEFT scheduler from StarPU. We obtained significant performance results with both XKaapi and StarPU for the Cholesky implementation. But, our scheduling strategy makes decisions at runtime and does not rely on any cost model. Nonetheless, our scheduling strategy based on work stealing lacks of more sophisticated decisions in order to consider data locality and processing power of available PUs. In our experiments with blocked Cholesky factorization, StarPU attains better results for executions with more than 4 GPU because its scheduling strategy considers communication costs. Coupling our online scheduling with data locality or communication information might enable to get a scalable Cholesky implementation on multi-GPUs.

As future works, we plan to extend XKaapi with an OpenMP-like interface for GPUs. Also, we will design dynamic strategies to take into account the affinity between tasks and data in order to reduce the data transfers in multi-GPU executions, as well as new adaptive algorithms using XKaapi for Cholesky, QR and LU factorizations.

Fig. 1 .

 1 Fig. 1. Example of a XKaapi C++ Cholesky factorization.

/Fig. 2 .

 2 Fig. 2. Example of a XKaapi C++ program for hybrid architectures. It shows a task Signature with its parameters and access modes, as well as CPU and GPU implementations.

Fig. 3 .

 3 Fig.3. Sequential and concurrent operations in a Fermi GPU card. Fermi GPUs have two copy engines and one execution engine capable of concurrent transfer operations (host-to-device and device-to-host) and kernel execution.

Fig. 4 .

 4 Fig. 4. Idgraf hardware topology with two six-core CPUs and eight GPUs.Here the Tylesburg-36D is a QPI-PCIe bridge and the PXE8647 is a PCIe switch for two GPUs.

Fig. 5 .

 5 Fig.5.Performance results from DGEMM on Idgraf for single CPU and single GPU and different block sizes.

Fig. 6 .

 6 Fig. 6. Cholesky performance results on Idgraf for single-CPU and single-GPU with block size 1024 × 1024.

Fig. 7 .

 7 Fig. 7.DGEMM performances up to 8 GPUs. The matrix order is 16384 with block size 1024 × 1024.

Fig. 8 .

 8 Fig.8.Cholesky factorization performances up to 8 GPUs. The matrix order is 16384 with block size 1024 × 1024.

Fig. 9 .

 9 Fig. 9. Results for DGEMM using 4 GPUs on various grains.

  DGEMM IN GIGABYTES WITH MATRIX ORDER 16384 AND BLOCK SIZE 1024 × 1024. THE SUM OF INPUT AND OUTPUT DATA TRANSFERS IS 8 GB.

		1	2	4	6	8
	XKaapi	8.00 10.27 16.49 23.23 29.29
	StarPU	22.54 14.97 16.98 20.09 24.35
			TABLE I
	MEMORY TRANSFERS OF

TABLE II MEMORY

 II TRANSFERS OF CHOLESKY IN GIGABYTES WITH MATRIX ORDER 16384 AND BLOCK SIZE 1024 × 1024. THE SUM OF INPUT AND OUTPUT DATA TRANSFERS IS 4 GB.

			Cholesky transfers (GB)
	GPUs	1	2	4	6	8
	XKaapi 3.71 7.38 12.28 12.89 14.62
	StarPU	2.23 3.81	6.55	7.38	9.00

VII. ACKNOWLEDGMENTS

The authors would like to thank Fabien Lementec for providing early implementations on GPU support. This work has been partially supported by the ANR-11-BS02-013 HPAC ANR Project, and CAPES/Brazil.