Jean-Philippe Karr

Laurent Hilico

Why three-body physics does not solve the proton radius puzzle

Keywords: numbers: 36, 10, Ee 31, 15, ac 31, 15, xt

, which could resolve the current discrepancy between determinations of the proton radius, is investigated. Using variational calculations with complex coordinate rotation, it is shown that the pµe ion, which was recently proposed as a possible candidate [2], has no resonant states in the energy region of interest. QED level shifts are included phenomenologically by including a Yukawa potential in the three-body Coulomb Hamiltonian before diagonalization. It is also shown that the ppµ molecular ion cannot play any role in the observed line.

Introduction. The recent measurement of the Lamb shift in muonic hydrogen [1], which resulted in a new determination of the proton rms charge radius, has given rise to an abundant literature attempting to resolve the 4percent discrepancy with its current CODATA value [3], mainly obtained from electron-proton (e-p) scattering and hydrogen (H) spectroscopy. The discrepancy with e-p scattering measurements [4] may be less severe than the initially announced 5 standard deviations due to a possible underestimation of error bars in e-p data analysis [5]. However, no agreement with H spectroscopic data [6] could be obtained unless some unknown systematic shift was discovered in all these experiments. An attempt to solve the discrepancy by re-evaluating the third Zemach moment of the proton [7] was shown to go against e-p scattering data [8]. Reconsideration of the assumptions on form factors used to extract the proton radius only led to a slight increase of the uncertainty [9]. Much effort has been focused on checking and improving the Lamb shift calculations used to extract the proton radius (see e.g. [2,10] and references therein, as well as [11]). In particular, the polarizability of the proton is a subject of hot debates [12], but no consensus has been reached yet regarding claims that this contribution may be able to solve the proton size puzzle. The muonic hydrogen result also triggered a search for new physics, which is however highly constrained by many low-energy data [13].

So far, no possible experimental artifacts have been proposed, with the notable exception of [2]. In that reference, Jentschura suggested that µp atoms formed in the 2S state might bind an electron present in the gas target to form a (pµe) -ion. Since the muon orbit is close to the proton, the proton charge is shielded from the outer electron, but attractive dipole interactions with the µp(2S) core could result in a quasibound (or resonant) state slightly below the 2S threshold. Jentschura gave an order-of-magnitude estimate of the associated * Electronic address: karr@spectro.jussieu.fr frequency shift, which was indeed compatible with the observed discrepancy ∆E exp -∆E th = +0.31 meV between experiment and theory. Among the many proposed solutions to the proton radius puzzle, this appears to be one of the last unexplored paths, and definitely deserves further study. The purpose of Section I is to give a conclusive -eventually negative-answer on the viability of the pµe hypothesis.

Other candidates to form such quasibound states with a µp atom can also be considered. While the possibility of forming (pµµ) -ions can be safely ignored -if only because of the low number of muons in the experiment chamber, it is known that (ppµ) + molecules can be formed in collisions of excited µp(2S) atoms with hydrogen molecules (see e.g. [14] and references therein). In this process, the favored molecular states lie ∼ 10 -20 eV below the n = 2 dissociation threshold, and their Coulomb lifetime is much too short to account for the observed transitions [15,16]. However, it is worth reconsidering the possible role of more weakly bound states which can have long lifetimes. In Section II, we discuss the relevant portion of the spectrum of ppµ and show that none of the existing states has the required characteristics to match the experimental data.

I. THE pµe HYPOTHESIS A. Coulomb Schrödinger Hamiltonian

A first insight into the physics of the pµe system is provided by the adiabatic approximation: the "slow" outer electron experiences the attractive dipole potential created by the pµ(n = 2) core. Neglecting (in a first step) the 2S-2P Lamb shift, the n = 2 degenerate level has a nonzero static dipole A = 3m e /m r in atomic units, where m r is the pµ reduced mass. The long-range behavior of the potential curve is then -A/R 2 , and it has been known for a long time [START_REF] Morse | Methods of Theoretical Physics[END_REF] that in this kind of potential, a particle of mass m has an infinite series of bound states with exponentially decreasing binding energies, provided A > A c = m e /8m. This translates into the condition m > m r /24 ≃ 7.74 m e . Since the pµe system is below that value, only a finite number of bound states can exist. This is a first indication that due to the weakness of the attractive dipole potential, the existence of resonant states is not assured.

For a more detailed investigation, we consider the fullthree body dynamics and use a numerical diagonalization of the complex rotated [START_REF] Reinhardt | [END_REF] Schrödinger Hamiltonian in a variational basis set to obtain the energies and widths of resonances. Having in mind that our goal is to determine whether resonant states do exist, we start with the pµµ negative ion, where resonances below the n = 2 threshold are known to exist, and then follow the position of the lowest resonant state (with ro-vibrational quantum numbers v = L = 0) while decreasing the mass m of the third particle -which we denote x in the following-in small steps, from the muon down to the electron mass. We used the perimetric coordinates (denoted x, y, z) and Sturmian wave functions

χ α,β nx,ny,nz (x, y, z) = (-1) nx+ny+nz αβ 2 ×L nx (αx) L ny (βy) L nz (βz) e -(αx+βy+βz)/2 (1)
where n x , n y , n z are non-negative integers and L n the Laguerre polynomials. The main idea of this method [19] is to take advantage of the dynamical symmetries of the three body Coulomb problem; in this basis set the Hamiltonian has strict coupling rules, leading to a sparse-band matrix which can be diagonalized in fewer operations with respect to a "full" matrix, maintaining better numerical stability [15,20]. This choice is well suited for the present study, in particular since the small number of variational parameters (the two length scales α -1 and β -1 ) makes it easy to readjust them for each value of m.

Our results are summarized in Fig. 1. Numerical accuracy is better than 10 -9 eV, so that even the smallest width, lying in the 10 -8 eV range, is obtained with two significant digits. The point at the extreme right corresponds to the pµµ negative ion; its binding energy E b = 120.4716427 eV and width Γ = 7.146063657 eV (obtained with the 2010 CODATA values of the proton and muon masses [3]), improve on the values reported in [21], E b = 120.39 eV and Γ = 7.12 eV. As expected, when the mass of the third particle decreases, its orbit around the pµ(n = 2) core becomes larger and its binding energy decreases.

Coulomb widths are shown mainly for illustration since the focus is on the existence of quasibound states. They decrease together with E b with decreasing m, because the system comes closer and closer to the adiabatic limit where the resonant state tends to a true bound state. The sharp drop between m = 48 m e and m = 47 m e is due to a crossing of the pµ(n = 2) and px(n = 1) thresholds which occurs at m = m p m µ /(4m p + 3m µ ) ≃ 47.67 m e . Below that value, the px(n = 1) threshold is above pµ(n = 2), so that one of the channels for Coulomb deexcitation becomes forbidden, and the width decreases.

No calculations were performed below m = 32 m e , where larger computer resources would be needed to con- verge the resonance width with sufficient accuracy. However, the results of Fig. 1 are amply sufficient to demonstrate that for m = m e the binding energy of a resonant state would be too small by many orders of magnitude to resolve the 0.31 meV discrepancy in the muonic hydrogen experiment.

B. Inclusion of QED effects

An important ingredient is still missing in our treatment, that is the Lamb shift which lifts the degeneracy of the n = 2 level, leading to a zero static dipole. The dipole potential in -A/R 2 is thus replaced by an induced dipole potential decaying like -α 0 /2R 4 where α 0 is the polarizability of the 2S (or 2P) state. As a consequence, when resonances come closer and closer to the threshold, they experience a more and more weakened attractive potential. The infinite dipole series obtained for 1/R 2 potentials becomes truncated since a 1/R 4 potential only admits a finite number of bound states [22,23]. It is therefore clear that the QED shift of energy levels does not put into question the above conclusion, since it weakens the attractive potential and can only favor the disappearance of weakly bound resonant states.

In order to show this more quantitatively, we represent QED corrections by an effective potential added to the three-body Coulomb Hamiltonian before diagonalization (the perturbative approach is not valid here due to the mixing of 2s and 2p orbitals). The most natural choice would be the Uehling vacuum polarization potential which accounts for more than 99% of the µp Lamb shift. This approach was used to obtain the energy shifts of dtµ resonances using the stabilization technique with a real scaling parameter [22]. However, that method would not yield precise enough results here because the states under consideration have a large width that is not negligible with respect to their binding energy. In principle, the complex rotation method allows for higher precision, but numerical experimentation suggests that the Uehling potential is not dilation analytic (see e.g. [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF] for a definition of this term) and thus cannot be easily combined with the complex rotation method. For this reason, we "mimic" the proton-muon Uehling potential by a Yukawa potential V ef f (r) = -λe -αr /r. First, we determined the exact shifts of the 2S and 2P levels in µp induced by the Uehling potential by diagonalization of the two-body Coulomb+Uehling Hamiltonian: ∆E(2P ) = -14.5789 meV and ∆E(2S) = -219.7374 meV [START_REF]The 2P level shift is not negligible since the radius of the muon's orbit (2([END_REF]. The deduced contribution to the Lamb shift, ∆E(2P ) -∆E(2S) = 205.1584 meV is in exact agreement with the result of [2]. Then, the parameters α = 492.86 and λ = 2.050982 10 -3 (in atomic units) were adjusted so that the diagonalization of the two-body Coulomb+Yukawa Hamiltonian yield the same exact shifts.

Going back to the three-body problem, the perimetric coordinates method loses its point here, because the added potential, having no selection rules, is represented by a "full" matrix. We thus used another variational expansion with Hylleraas coordinates r 1 , r 2 , r 12 and exponential basis functions

χ 1(2) n (r 1 , r 2 , r 12 ) = Re(Im) e -αnr12-βnr1-γnr2 (2)
with pseudorandom complex exponents α n , β n , γ n [START_REF] Korobov | [END_REF].

In order to test the validity of the Yukawa approximation, we computed the energy of the lowest resonances of dtµ below the n = 2 threshold. Agreement with calculations using the exact Uehling potential [22] is obtained at the meV level (see Table I). The approximation is expected to work even better for more weakly bound resonances, since the outer particle's orbit lies farther from the µp(n = 2) core and will be less sensitive to the shortrange details of the proton-muon interaction potential. The main approximation in our approach is the neglect of fine and hyperfine structure (e.g. the 2S 1/2 hyperfine splitting is 22.8 meV). This is not expected to strongly modify our results, since the polarizabilities of individual hyperfine levels differ from the structureless case by less than 2 percent. We estimate that QED-induced changes in binding energies and widths are obtained with a relative accuracy of a few percent.

dtµ v E b (eV) [15] E Y b (eV ) E U b (eV) [
Results are shown in Table I and plotted in Fig. 1 (open circles). The drop of binding energy ranges from about 62 meV for m = 80 m e down to a few meV for m = 50 m e . This is consistent with the picture of a resonant state as a µp atom in a 2S orbital interacting with the field of the distant x particle. Indeed, the 2S-2P mixing decreases as x becomes more loosely bound so that the shift of the resonance tends to that of the 2S level (220 meV). Already for m = 48 m e , the binding energy is in the 0.1 meV range. These results confirm that the pµe system can have no resonance that would be sufficiently shifted with respect to the 2S atomic level to explain even a small fraction of the 0.3 meV discrepancy.

II. MOLECULAR STATES

As mentioned in the Introduction, metastable states of the ppµ molecule are known to play a role in the muonic hydrogen deexcitation cascade [14,27]. However, their possible involvement in the observed transition signal is highly constrained by experimental data [1]. Firstly, the observed 18 GHz-wide transition cannot be linked to a photo-association (or photo-dissociation) process, which would lead to a much broader spectrum. Secondly, two hyperfine components were observed, and the second one [28] was found at the predicted frequency with the proton radius value deduced from the first one i.e. the measured line has a hyperfine structure that closely matches that of the 2S 1/2 → 2P 3/2 atomic transition. Thus the only possible scenario involving molecular states is that of a transition between resonant states, a first one lying close to the µp(2S 1/2 ) threshold and a second one close to the µp(2P 3/2 ) threshold. A prerequisite to a resolution of the proton radius puzzle would be to find such a pair of states, connected by a dipole-allowed transition, with such binding energies that the transition would be blue-shifted by a fraction of meV with respect to the 2S 1/2 -2P 3/2 atomic transition. One strong point of this hypothesis is that the lower state would have long enough Coulomb and radiative lifetimes (since its wave function mainly involves the 2S atomic orbital with extremely small admixture with 2P) to survive the ∼ 1µs delay before the spectroscopy laser pulse; similarly, the upper state wave function being mainly of 2P character, the transition width could be close to the observed one. There are three binding Born-Oppenheimer (BO) electronic curves below the n = 2 threshold (see Fig 2a). The 3dσ g and 4f σ u curves have a 1/R 2 long-range behavior (dipole potential) and support infinite series of resonant states of symmetries 1 S e , 3 P o , 1 D e ... and 3 S e , 1 P o , 3 D e ... respectively [15,16,29]. The 2pπ u curve decays like 1/R 4 (induced dipole potential) and supports a finite number of states with symmetries 1 P o , 3 P e , 1 D o , 3 D e ... [30][31][32]. The Lamb shift splits the n = 2 threshold, and the key point is that both σ curves converge to the µp(2S) threshold, while only the 2pπ u curve converges to the 2P threshold. This is clearly apparent in Fig 2b, which shows how the BO energies are modified by the Lamb shift at large internuclear distances. We obtained these curves by calculating the energy of an electron in the field of the nuclei for a grid of values of R, with a Yukawa potential added to the Coulomb potential in order to represent the Uehling potential, as described in Sect. I B. The electronic wave function was expanded on a variational set of exponential basis functions with pseudorandom exponents [33].

These considerations show that the upper level of the proposed molecular transition could only be one of the levels supported by the 2pπ u curve. These states have already been investigated in detail [30][31][32]: one 1 P o resonant state (with a binding energy E b = 13.57 eV) and one "anomalous parity" 3 P e bound state (E b = 13.54 eV) have been evidenced, and it was shown convincingly in [31] that no other state can exist [START_REF]Shape" resonances lying above the atomic thresholds are known to exist in muonic molecules[END_REF]. We can thus conclude that the hypothesis of a molecular transition fails due to the lack of a suitable upper state.

Conclusion.

We have refuted all reasonable hypotheses aiming to resolve the "proton radius puzzle" with the help of three-body physics. On the one hand, it was shown that the pµe atom has no resonant states below the µp(n = 2) threshold; on the other hand, that the spectrum of the ppµ molecular ion leaves no room for a transition lying close to the 2S-2P atomic line, since there are no long-lived resonant states in the vicinity of the 2P threshold. However, the search for possible experimental artifacts in the muonic hydrogen experiment might not be over yet; as noted in [2], formation and deexcitation of muonic hydrogen are very complicated processes [START_REF] Jensen | [END_REF] in which other many-body bound or quasibound states could play a role. Another result of this work is the demonstration of a convenient method to study the effect of QED level shifts on the position and width of few-body resonances in the framework of complex-rotated variational methods, using an adjusted Yukawa potential.

FIG. 1 :

 1 FIG.1: Binding energy (top) E b = E(pµ)n=2 -E and Coulomb width (bottom) of the lowest 1 S e resonance of the p + µ -x -system, as a function of the mass m of the third particle x. Squares: pure Coulomb Hamiltonian, circles: Coulomb Hamiltonian with an added Yukawa potential adjusted to reproduce the 2S-2P Lamb shift. For the latter, the binding energy is with respect to the shifted 2S threshold. On the top graph, the horizontal line recalls the order of magnitude of the discrepancy between theory and experiment in the proton radius experiment, while the vertical line at m/me = 7.74 signals the minimal value of m to retain an infinite series of resonant states in the adiabatic approximation. On the bottom graph, the horizontal line shows the value of the width corresponding to a Coulomb lifetime τc = 1 µs, and the vertical one at m/me = 47.67 shows the crossing point between the pµ(n = 2) and px(n = 1) thresholds.

FIG. 2 :

 2 FIG. 2: (a) Binding BO electronic energy curves below the n = 2 threshold. Distances are in units of the muonic Bohr radius aµ ≃ a0/207, and energies are with respect to the threshold. (b) Enlarged view of the region of large internuclear distances. The thin black lines are the continuation of the curves in (a), and converge to the unshifted n = 2 threshold. The vertical scale on the left (right) correspond to the σ (π) curves respectively. The σ curves are indistinguishable at such large distances; their convergence is slower because they behave like 1/R 2 instead of 1/R 4 . The thick red lines are obtained by adding a Yukawa potential Y representing the proton-muon Uehling potential. The 2pπu curve converges to the shifted 2P threshold shown on the right vertical axis. The 3dσg and 4f σu curves converge to the shifted 2S threshold shown on the left axis. Since BO curves are drawn for infinite nuclear masses, for consistency the vacuum polarization shifts and associated parameters of Y were calculated for a µp atom with infinite proton mass: ∆E(2S) = -281.5896 meV and ∆E(2P ) = -27.5666 meV yielding an exponent α = 443.7624 and coefficient λ = 1.919722 10 -3 for the Yukawa potential (see Sect. I B).

TABLE I :

 I Binding energies of resonant states of dtµ (top) and pµx (bottom) below the µp(2S) threshold. E b is obtained with the pure Coulomb Hamiltonian, E Y b with an added Yukawa potential (see text), and E U b (for dtµ only) with the exact Uehling potential. Resonance positions are given with respect to the shifted 2S atomic threshold. For pµx, the corresponding widths Γ, Γ Y are also given.

				22]
	0	217.889 8 217.828 9	217.829
			pµx
	m/me	E b (eV )	E Y b (eV)	Γ (eV)	Γ Y (eV)
	80	1.673 0	1.611 0	0.907 7	0.907 1
	70	0.623 5	0.567 1	0.423 9	0.422 6
	60	0.165 2	0.118 2	0.133 5	0.129 1
	55	0.070 2	0.032 2	0.057 7	0.050 1
	50	0.025 3	0.001 6	0.017 1	0.008 1
	48	0.016 8	7(6).10 -5	0.006 4	0.001 3

Acknowledgments. We thank V.I. Korobov for sharing his program for variational calculations of three-body systems with exponential basis functions, and V.I. Korobov, P. Indelicato, F. Nez, D. Delande, R. Pohl and F. Kottmann for helpful discussions.


Yoann Isaac
email: yoann.isaac@supelec.r

Yolaine Bourda
email: yolaine.bourda@supelec.r

Monique Grandbastien
email: monique.grandbastien@loria.fr

SemUNIT -French UNT and Linked Data

Keywords: I.2.4 [ARTIFICIAL INTELLIGENCE]: Knowledge Representation Formalisms and Methods; K.3.1 [COMPUTERS AND EDUCATION]: Computer Uses in Education E-learning, ontology for learning resources, Semantic Web, Linked Data

Over the past 15 years, the explosion in the number of learning materials available on the Web has raised the problem of their sharing. For several years, learning resources are annotated with metadata to ease this sharing. On the other hand, Semantic Web and Linked Data approach provide tools to publish metadata in a standardized way, allowing data to be shared and reused across applications, enterprises, and community boundaries.

In this paper, we describe the SemUnit project, initiated by french higher education institutions. This project aims at taking advantages of Semantic Web and Linked Data to improve e-learning services for a wide set of french higher education institutions. We present, firstly, the ontology designed to support the project: an OWL ontology taking into account the semantics of LOM elements. Afterwards, we present our architecture and some semantic services 1 .

INTRODUCTION

Today, more and more learning materials are available on the Web. The sharing of these learning materials can greatly ease the creation of quality courses, especially when they are composed of well-defined reusable learning objects [START_REF] Mohan | Learning objects on the semantic web[END_REF]. This sharing requires the creation of tools allowing to store, search and evaluate them efficiently. The full text search approach mostly used today is not adapted. Consequently, for several years, learning materials have been annotated with metadata using mainly the Learning Object Metadata (LOM) standard, one of the LOM application profiles or the Dublin Core (DC).

The creation of high quality learning objects is time consuming and requires to be supported. In 2003, the French government created the UNTs (Universités Numériques Thématiques: Thematic Digital Universities) in order to support the creation of high quality learning material and ease their sharing between higher education institutions (Universities and "Grandes Ecoles"). Each UNT gathers learning resources belonging to a specific field. Among these UNT, UNIT2 is the one dedicated to "Engineering sciences and technologies". UNIT currently provides a free access to several thousands of high quality learning objects. They are annotated with metadata complying to the SupLOMFR schema 3 , a LOM application profile dedicated to French higher education institutions. Metadata are published thanks to a network of OAI-PMH 4 portals under the open-source software ORI-OAI 5 . This system makes metadata and resources harvestable [START_REF] Van De Sompel | Resource harvesting within the oai-pmh framework[END_REF] and then allows the creation of crossrepository services. Nonetheless, only search services based on input forms are available. Although useful, this kind of service is quite restrictive and does not allow indexing of educational material by search engines. This therefore restricts the visibility of the metadata and, of the described resources and thus reduces their use (to be used a resource first needs to be found and known). Furthermore, it is difficult to automatically link those metadata with results coming from other repositories (like DBLP for example).

On the other hand, the explosion of the Semantic Web and Linked Data, in the last few years, has allowed the production of a large number of data with a clear semantic, increasing the visibility of such data. The expression "Linked Data" corresponds to a set of best practices for publishing structured data [START_REF] Heath | Linked Data: Evolving the Web into a Global Data Space[END_REF]. In the e-learning context, the respect of linked data principles can be beneficial in many ways. Indeed, if learning resources belonging to several repositories are linked together or with others resources: cross-repository services can be created (like search engines); learning re-sources can reference RDF resources close to theirs topics to offer useful information (DBpedia [START_REF] Waitelonis | Towards exploratory video search using linked data[END_REF]); learning resources can be described by linking them to specific taxonomies; learning ontologies can be mapped to others and these mapping can be used efficiently for reasoning.

Thus, the SemUnit project has been launched in order to bring the advantages of Semantic Web and Linked Data technologies to the French UNT repositories. This relies on the creation of an ontology taking into account the semantics of LOM elements. But unlike other works, the designed ontology has a structure (classes and properties) that is not based on the LOM table and reflects the meaning of metadata elements.

The main steps of the SemUnit project are: the creation of an ontology capturing the semantics of the metadata, the automatic translation of metadata from XML to RDF triples, the study of services taking advantage of the Semantic Web metadata previously created, the definition of an architecture supporting the desired services and the creation or use of the related tools.

This paper is organised as follow: section 2 describes the SemUNT ontology, section 3 presents the architecture and the tools, section 4 introduces some services grouped in a Web application, section 5 presents some related works and finally, section 6 concludes and presents some future works

THE SEMUNT ONTOLOGY

UNT learning materials are described by metadata complying with the SupLOMFR schema (an IEEE LOM application profile) and encoded in XML. Give meanings to these annotations and transform them to RDF triples requires the creation of an ontology. OWL rather than RDF(S) has been chosen to implement this ontology because of it's greater expressiveness, for example OWL cardinality restrictions are particularly important to express cardinalities. Furthermore, unlike other works, the designed ontology has a structure (classes and properties) that is not based on the LOM table and reflects the meaning of the metadata elements.

The next two subsections describe the choices made to create this ontology. The first part deals with the general structure of the ontology and the second with the management of vocabularies.

SemUNT ontology structure

Related ontologies

Several e-learning ontologies based on the LOM or RDF(S) LOM bindings have already been created. The first was defined in RDF(S) by Nilsson [START_REF] Nilsson | The LOM RDF binding-principles and implementation[END_REF]. It proposes some interesting ideas (which have been adopted for the SemUNT ontology): use language tags for lom:langString, represent VCARD text with a specific ontology, integrate Dublin Core Properties (title, format...) and DataType (DateTime, Language...) when possible, etc. Nevertheless, it does not provide solutions concerning classifications or vocabularies management and its RDF(S) implementation does not allow cardinality restrictions.

More recently, other LOM ontologies have been created (e.g. [START_REF] Evangelia | Connecting medical educational resources to the linked data cloud: the meducator rdf schema, store and api[END_REF][START_REF] García | V Simposio Pluridisciplinar sobre Diseno y Evaluación de Contenidos Educativos Reutilizables (SPDECE)[END_REF][START_REF] Ghebghoub Ouafia | Ontology Theory, Management and Design: Advanced Tools and Models, chapter Building and Use of a LOM Ontology[END_REF][START_REF] Koutsomitropoulos | The use of metadata for educational resources in digital repositories: Practices and perspectives[END_REF]). Some of these ontologies ( [START_REF] García | V Simposio Pluridisciplinar sobre Diseno y Evaluación de Contenidos Educativos Reutilizables (SPDECE)[END_REF][START_REF] Ghebghoub Ouafia | Ontology Theory, Management and Design: Advanced Tools and Models, chapter Building and Use of a LOM Ontology[END_REF]) are very closed to the LOM's table structure and then, allow an easy transformation from a LOM XML binding to an RDF one. The semantics of the LOM elements is, on the contrary, not well modelled and some ontology's classes have no real meaning. LOM categories are for instance mapped to classes even if this mapping is artificial. Only the relation category is mapped differently (e.g. in [START_REF] Ghebghoub Ouafia | Ontology Theory, Management and Design: Advanced Tools and Models, chapter Building and Use of a LOM Ontology[END_REF]) to respect its meaning: LOM relations allow to link learning resources to others and consequently are mapped to object properties. Other ontologies [START_REF] Evangelia | Connecting medical educational resources to the linked data cloud: the meducator rdf schema, store and api[END_REF][START_REF] Koutsomitropoulos | The use of metadata for educational resources in digital repositories: Practices and perspectives[END_REF] only consider a subset of the LOM and although their choices are closer to ours we have not been able to reuse them as is because we needed the entire LOM and because we wanted to use OWL (to express constraints) and not RDF(S) and SKOS to express controlled vocabularies.

Ontology Construction

The SemUNT ontology (see figures 1 and 2) has been created with some objectives and constraints in mind: respect the meaning of metadata elements; maximize its interoperability; respect linked data rules. This last point led us to avoid some RDF constructions like collections, containers and blank nodes.

Classes.

All learning resources are instances of the class Learn-ingResource. Classes have been defined for categories representing something real for learning resources like Classification. Classes have also been defined for some structured elements like Contribution because a contribution is a grouping of several attributes: a date, a role and a contributor.

Properties.

Simple elements are mapped to properties: datatype properties for these (title, description, format, etc) taking their values in data types and object properties for those taking their values in vocabularies (see section 2.2). Furthermore, for each of the defined classes (Classification, Contribution...) an object property has been created (hasContribution, hasClassification...) whose domain is the class LearningResource and range is the defined class. The category Relation is mapped to an object property linking two learning resources. The kind of relation is defined using the subProperty construct. This representation allows a better browsing of learning resources and stays close to the meaning of the metadata element. However, with this representation, the description of the relation is lost but this element is rarely filled in with UNT learning resources. This solution has been nevertheless chosen because the pros outweigh the cons.

Interoperability.

The SupLOMFR schema, like the LOM one, defines some correspondences with Dublin Core (DC) properties. We choose to integrate them directly in the SemUNT ontology (identifier, title , language, coverage, type, rights, relation, description, format and subject).

The ontology interoperability is also improved by using the FOAF Ontology, one of the most used for people's description [START_REF] Ding | How the semantic web is being used: An analysis of foaf documents[END_REF]. Thus the foaf:Agent class is used for VCard elements. The only mapping problem is that some VCard attributes does not correspond to any FOAF properties for now but these elements are not filled in UNT metadata.

Vocabularies management

Some SupLOMFR metadata elements take their value in controlled vocabularies. Some of the possible values are coming from the LOM Schema, some of them are new values linked or not to some LOM values by a relationship (like narrower). For example in the LOM, for the 5.6 level element we have the value "Higher Education", in the SupLOMFR we have "Enseignement Supérieur" (the same than the LOM one but in French), and more specific values like "Licence" or "Doctorat". But as we also want to be able to use metadata complying with the LOM schema (and other LOM application profiles like normetic), we decided to use SKOS (Simple Knowledge Organization System) to describe vocabularies. In this way, one can define distinct controlled vocabularies sharing concepts or having concepts linked by relationships (like broader or narrower). An example of the implementation of the vocabularies is given figure 3.

ARCHITECTURE AND TRANSFORMA-TION TOOLS

Architecture

Concerning the storage of RDF triples, the system requirements were the following ones: it should accept at least 1 millions RDF triples and be scalable; it should offer a SPARQL access with relatively fast answer; it should propose a linked data interface;

To meet these requirements, the following tools have been chosen: OWLIM-Lite6 for storing RDF triples, Sesame27 for managing this repository and Pubby8 for exposing a linked data interface. OWLIM was choosen because of it's good results at the Berlin SPARQL Benchmark9 . The architecture is presented in figure 4.

Transformation tools

We defined transformation tools allowing to automatically convert XML metadata to RDF triples (complying with the ontology previously defined) and studied RDF storage solutions (1 million RDF triples have been created from 6456 UNT metadata files).

An handmade XSLT transformation has been chosen to automatically "RDFize" the XML UNT metadata although different "RDFizers" are already freely available. Indeed, Figure 4: Architecture of the SemUNIT repository these tools are very useful in some specific cases but they can only transform standard data formats, they does not allow to use a custom ontology (for most of them) and their support and maintenance are often not insured. An handmade XSLT transformation allows a better control of the transformation and does not impose heavy constraints because the SemUNIT ontology is small enough to allow an easy maintenance of the XSL file (when the ontology is modified).

The complete transformation is performed using a python script which executes different procedures before and after the launching of the XSLT transformation. These procedures are the ones which can not be done with XPath functions in XSLT transformation or which need to be done only after all metadata files have been transformed to RDF. The most important of them are:

• transformation of strings belonging to vocabularies into RDF identifiers defined in the SKOS taxonomy described previously (before XSLT);

• transformation of VCard text to RDF FOAF triples (before XSLT);

• management of lom:relation between learning resources. Indeed, relations can either link UNT documents to other UNT documents or to any Web documents (URI HTTP) and URI of UNT documents are created during the transformation. Consequently, this procedure should be done after the whole XSLT transformation.

SOME SERVICES

We have designed some services, described below, consuming the triples and taking advantages of the semantic annotations, they are grouped into a Web application 10 .

The browsing of metadata in this Web application is insured by Pubby which is integrated via a simple frame in its interface.

In addition to the services described below, we have proposed distinct visualizations: pie charts for values used in metadata elements using controlled vocabularies, timelines... 10 http://semunt.supelec.fr

Search Engine

Engine purpose

The engine's purpose is to propose a search service taking advantages of Semantic Web annotations to retrieve learning resources more efficiently. Indeed, most of the users does not know the keywords required to find the documents they need, and even when they know them, useful documents can be lost in hundreds of other documents returned as results.

In addition, a well-designed interface could also be efficient to help users finding information easily. Consequently, a Faceted Search, one of the exploratory Search techniques [START_REF] Petratos | Informing through user-centered exploratory search and human-computer interaction strategies[END_REF], have been implemented in our engine. In faceted search, facets are conceptual categories allowing to organize data from a large database through a view of clear conceptual groups. There are two kinds of facets : flat and hierarchical [START_REF] Petratos | Informing through user-centered exploratory search and human-computer interaction strategies[END_REF]. In our search engine, facets have been created to filter results obtained from the classical keyword-based approach and then perform more accurate search. Currently, users can filter results (with facets) on metadata elements taking their values in vocabularies (like structure, format or learning type). For now, only the classification based on the Dewey (mandatory in the UNT context) allows a hierarchical classification.

Detailed description

This search engine works as follows (see figure 5 for an example of a query):

• the user fills in a small form to precise : keywords corresponding to documents they want to find (separated by blank spaces); metadata fields where these keywords are searched: title, description, subject, etc;

(optional) person linked to documents they want to find; (optional) role of this person.

• the engine creates a SPARQL query to find these documents, ranks search results and return information about them to the user;

• the user is invited to filter results by selecting elements in different metadata elements presented with faceted search interface;

• the engine returns filtered results and allows the user to change filters and keywords to meet his needs.

The ranking of results is performed using the following algorithm:

• a document gets points each time a keyword appears in its metadata elements;

• the number of points earned depends on the metadata element containing the keywords: 10 for title, 4 for description, 2 for keywords. These numbers have been determined empirically;

• bonus are given when all the keywords are present in the same element (if several keywords has been entered).

Results are sorted according to points they have earned and presented as a list, giving for each document: its title, a link allowing to browse its metadata, a link allowing to get the document and its description and subjects (if their checkboxes are checked).

Expert search service

This service aims at retrieving people having an expertise in a specific topic (in fact are able to give lessons on a specific topic). The design of this kind of service needs to answer 3 questions :

• which information sources will be used to perform the search;

• which system will be offered to users for choosing the topic(s);

• which metric(s) will be chosen to rank experts found.

Concerning information sources, we used the instances of the class Contribution (which is linked to a contributor and to a kind of contribution). Thus, we can know for each learning resource what kind of contribution a person has brought: author, contributor, publisher, subject-matter expert, etc. The coefficients associated with each kind of contribution have been determined empirically.

About the choice of topic, we have chosen to let users fill in one field form to specify it rather than proposing a list of topics or a list of classification's elements. It offers a more flexible way of specifying topics and allows the reuse of search engine procedures in the ranking system.

RELATED WORK

Despite the spread of Linked Data, very few works use the linked data principles in education, the first specific workshop (Linked Learning 2011) was organized in 2011.

The Open University(OU) linked data architectures proposed [START_REF] Zablith | The ou linked open data: production and consumption[END_REF]: a workflow to RDFize data coming from different repositories; an expert search service; a service of students social network; a service allowing students to search and browse free learning materials. This is developed as part of the LUCERO 11 (Linking University Content for Education and Research Online) project. But, the schema 12 used to describe the courses is based on the courseware ontology 13 and on the AIISO ontology 14 , two ontologies designed to The mEducator project uses linked data for the publishing of medical educational resources [START_REF] Evangelia | Connecting medical educational resources to the linked data cloud: the meducator rdf schema, store and api[END_REF]. This project has defined a metadata scheme based on Healthcare LOM (HLOM). It has been designed to be interoperable with existing schemas with the integration of Dublin Core properties as super properties and the use of FOAF properties to describe people. SKOS has been chosen for controlled vocabularies [START_REF] Giordano | Developing controlled vocabularies for educational resources sharing: a case study[END_REF]. But the schema used by the mEducator project is a subset of the LOM schema, so some metadata elements like the ones used for the classification are missing. Furthermore, their ontology is only based on RDF(S), so it doesn't take advantage of OWL constructs (like constraints) or OWL tools (like reasoners).

In addition, different services depending on linked data have been proposed : see for instance [START_REF] Stankovic | Linked data metrics for flexible expert search on the open web[END_REF] which proposes an expert search based on linked data metrics or [START_REF] Kaldoudi | Educational resources as social objects in semantic social networks[END_REF] which develops the idea of integrating learning resources in a linked data social network as social object.

CONCLUSION AND FUTURE WORKS

This project has shown how the use of Semantic Web technologies coupled with linked data practices could greatly improve e-learning services by giving a formal semantics to metadata elements.

Firstly, we have created the SemUNT ontology which has given a semantic representation to all SupLOMFR metadata elements. This ontology has been designed to be compliant with the meaning of metadata elements and allow an efficient management of vocabularies thanks to the use of SKOS. Among other possibilities, this implementation could ease the management of application profiles vocabularies and the integration of vocabularies in a linked data context. In addition, the ontology has been designed to insured its interoperability with existing ones. Thus, Dublin Core properties have been integrated and FOAF has been chosen to describe persons and organizations.

Then, we have developed tools enabling to automatically transform XML metadata to RDF triples. These tools are based on handmade XSLT transformations and have pro-duced about 1 million RDF triples from 6456 metadata documents coming from several UNTs: Unit, Unisciel, Uel, Uoh, Uved.

Finally, a SPARQL endpoint 15 has ben proposed and two services 16 have been designed to consume these triples: a search engine using semantic annotations to improve its search efficiency and offer faceted search possibilities and an expert search service allowing to find people skilled in a specific area. These services have shown how a semantic binding of metadata can be used efficiently to improve classical services.

The most original part of this project is the creation of an ontology for the management of educational resources, taking into account the ontology semantics of LOM elements but having a structure (classes and properties) that is not based on that of LOM and reflects the meaning of metadata elements. this ontology is intended to be reused by teams engaged in similar work.

These results and tools will be widely spread in the UNTs community. They have already been presented to representative of different UNTs which have shown their interest. Several scenarios have been imagined as well as creation of multiples services consuming these data and taking advantage of the semantic binding. This project which for now is only a feasibility model, aims to become an institutional project involving all French UNT.

Of course many improvements can be made to our services: search engine and expert search. For example services for documents or expert. Among them: changing the flat facets by hierarchical facets when possible, improving the ranking of the search results, using reasoners reasoners to obtain new information on our data, linking our data to other repositories like DBLP, Semantic Chrunch-Base, Twarql, RDFohloh1, etc.

Figure 1 :

 1 Figure 1: Classes of the SemUNT ontology

Figure 2 :

 2 Figure 2: Properties of the SemUNT ontology

Figure 3 :

 3 Figure 3: Vocabularies in SKOS

Figure 5 :

 5 Figure 5: Search of string "Automatique"

Figure 6 :

 6 Figure 6: Expert search engine

http://www.unit.eu/fr

http://www.sup.lomfr.fr/

http://www.openarchives.org/OAI/openarchivesprotocol.html

http://www.ori-oai.org/display/ORIOAI/ORI-OAI.ORG

http://www.ontotext.com/owlim

http://www.openrdf.org/

http://www4.wiwiss.fu-berlin.de/pubby/

http://www4.wiwiss.fu-berlin.de/bizer/berlinsparqlbenchmark/

ACKNOWLEDGMENTS

This project has been funded with support from the UNT UNIT.

We want to thank particularly Marie-Hélène Comte and Julia Soyez, librarians in the SemUNIT project.


Xavier Lacoste
email: xavier.lacoste@inria.fr

Pierre Ramet
email: pierre.ramet@inria.fr

Mathieu Faverge
email: mfaverge@eecs.utk.edu

Ichitaro Yamazaki
email: iyamazak@eecs.utk.edu

Jack Dongarra
email: dongarra@eecs.utk.edu

Sparse direct solvers with accelerators over DAG runtimes

Keywords: Algèbre linéaire creuse, Calcul haute performance DAG, GPU Sparse linear algebra, High performance computing, DAG, GPU

Les architectures de calcul intègrent de plus en plus de coeurs de calcul partageant une même mémoire nécessairement hiérarchique. Les algorithmes, en particulier ceux relatifs à l'algèbre linéaire, nécessitent d'être adaptés à ces nouvelles architectures pour être efficaces. PASTIX * est un solveur direct parallèle pour matrices creuses qui intègre un ordonnanceur dynamique pour des architectures hiérarchiques de grande taille. Dans ce papier, nous étudions la possibilité de remplacer cette stratégie interne d'ordonnancement par deux supports d'exécution génériques : DAGUE † et STARPU ‡ . Ces supports d'exécution offrent la possibilité de dérouler le graphe de tâches de la factorisation numérique sur des noeuds de calcul disposant d'accélérateurs. Comme pour de précédents travaux réalisés pour les noyaux denses de l'algèbre linéaire, nous présentons nos noyaux, pour des structures creuses, issus de la librairie MAGMA et des algorithmes DAG utilisés par ces deux supports d'exécution. Nous présentons une étude comparative des performances de notre solveur supernodal avec ces trois ordonnanceurs sur des architectures multicoeurs et multigpus, et en particulier les gains obtenus sur les accélérateurs avec le support d'exécution STARPU. Ces résultats montrent qu'une approche basée sur un DAG offre une interface de programmation uniforme pour réaliser du calcul haute performance sur des problèmes irréguliers comme ceux de l'algèbre linéaire creuse.

I. INTRODUCTION

Solving a large sparse general or symmetric positive definite linear system of equations, Ax = b, is a crucial and timeconsuming part in many scientific and engineering applications. Due to their robustness, direct solvers are often used in industrial codes, despite their large memory consumption. In addition, the factorization used in the recent direct solvers can take advantage of the superscalar capabilities of the processors using blockwise algorithms and BLAS primitives. Consequently, many parallel techniques for sparse matrix factorization have been studied and implemented. For a complete survey on direct methods, we refer to [START_REF] Duff | Direct methods for sparse matrices[END_REF]- [START_REF] George | Computer solution of large sparse positive definite systems[END_REF]. The goal of this paper is to design algorithms that can plainly take advantage of the vast computing power found in modern heterogeneous computer architectures. In our current work, we focused on matrices with symmetric sparsity patterns (the sparsity patterns of A + A T is used for unsymmetric cases), and focused on the factorization with static pivoting. In this context, the block structure of the factors and the numerical operations are known in advance, and consequently, we can use a static (i.e. before the actual numerical factorization) algorithm for scheduling the communication and the computational tasks.

In previous works [START_REF] Hénon | PaStiX : A Parallel Sparse Direct Solver Based on a Static Scheduling for Mixed 1D/2D Block Distributions[END_REF], [START_REF]PaStiX : A High-Performance Parallel Direct Solver for Sparse Symmetric Definite Systems[END_REF], we have proposed a static mapping and scheduling algorithm based on a combination of 1D and 2D block distributions for sparse supernodal factorization with static pivoting. This algorithm achieves very good performance by taking into account the communication and computation requirements of each factorization as well as the communication and computation capabilities of the parallel architecture. In addition, we have developed a strategy to control the memory overhead and to reduce the communication volume needed for the message buffering. In PASTIX, this buffering corresponds to the local aggregation approach in which all local contributions for the same non-local block are summed into a temporary block buffer before being sent. The new strategy improves the mechanism of this local aggregation that often lead to great reduction in the memory consumption, especially for 3D problems.

Emerging supercomputers consist of many microprocessors, each of which may have many computational cores. Hence, these emerging architectures exhibit strongly hierarchical topologies, both in terms of memories and processor interconnect. Achieving good performance requires a mapping of the algorithms on the computational resources, and scheduling algorithms specifically designed for NUMA architectures.

In PASTIX the internal data structures of the solver, as well as the communication patterns, have been modified for the dynamic scheduling on these architectures [START_REF] Faverge | Dynamic scheduling for sparse direct solver on NUMA architectures[END_REF]. A dynamic scheduler based on a work-stealing algorithm has been also developed to fill in communication idle times while preserving a good locality for data mapping [START_REF]Fine grain scheduling for sparse solver on manycore architectures[END_REF]. Furthermore, we have integrated the dynamic adaptation of the computational task grain to efficiently use multi-core architectures and shared memory. Experiments on several numerical test cases have been performed to prove the efficiency of these approach on different architectures.

Modern GPUs can substantially outperform high-end multicore CPUs both in terms of data processing rate and memory bandwidth. In the past, porting general purpose codes onto GPUs required a considerable programming effort, mostly due to the lack of tools and interfaces. However, API for GPUs such as CUDA or OpenCL have been rapidly evolving in the last few years, bringing accelerator programming into the mainstream. Hence, GPUs are becoming a more and more attractive alternative to traditional CPUs, in particular for the more interesting ratio of cost-per-flop and watts-per-flop. However, efficient GPU programming remains a laborious challenge. In this paper, in order to exploit the computing power of GPUs, we have integrated a runtime system into the sparse direct solver. We divided the algorithm into computational tasks, and use a Directed Acyclic Graph (DAG) to represent the dependencies between these tasks. This DAG representation of an algorithm enables a clear separation between the flow of data between tasks and the data distribution. Then, a runtime system is used to schedule these tasks while respecting the dependencies. Since the runtime system offers a uniform programming interface for a specific subset of hardware or low-level software entities, applications can use these uniform programming interfaces in a portable manner. Furthermore, the runtime system can optimize application requests by dynamically mapping the tasks onto resources as efficiently as possible. The work presented in this paper is part of the MORSE project1 to design dense and sparse linear algebra methods that achieve the fastest possible time to an accurate solution on large-scale heterogeneous multicore systems with GPU accelerators.

Concerning accelerator-based platforms for sparse direct solvers, a lot of attention has been recently paid to design new algorithms that can exploit the huge potential of GPUs. For a multifrontal sparse direct solver, some preliminary works have been proposed in the community [START_REF] George | Multifrontal Factorization of Sparse SPD Matrices on GPUs[END_REF], [START_REF] Yu | A CPU-GPU Hybrid Approach for the Unsymmetric Multifrontal Method[END_REF], resulting in single-GPU implementations based on off-loading parts of the computations to the GPU. The main idea is to treat some parts of the task dependency graph entirely on the GPU. Therefore, the main originality of these efforts was in the methods and algorithms used to decide whether or not a task can be processed on a GPU. In most cases this was achieved through a threshold based criterion on the size of the computational tasks. From the software point of view, most of these studies have only produced software prototypes and few sparse direct solvers exploiting GPUs have been made available to the users, most of them being developed by private software companies such as MatrixPro2 , Acceleware3 and BCSLib-GPU 4 . As far as we know, there are no publications nor reports where the algorithmic choices are depicted. A recent progress towards a multifrontal sparse QR factorization on GPU have been presented in [START_REF] Davis | Multifrontal sparse qr factorization : Multicore, and gpu work in progress[END_REF].

The rest of the paper is organized as follows : In the next section, we present the libraries and software used for our current study. There, both DAGUE and STARPU runtimes are introduced, and we summarize the main features of the PAS-TIX solver. Then, in Section III, we describe the algorithms and main ideas that have been implemented to optimize the task scheduling and the granularity of the sparse factorization. In Section IV, we present experiments we performed with our sparse direct solver to validate our approach on industrial challenging matrices. Finally, in Section V, we conclude with some prospects of the current work.

II. RUNTIME AND SOLVER

DAGUE [START_REF] Bosilca | DAGuE : A generic distributed DAG engine for High Performance Computing[END_REF] is a distributed runtime system designed to achieve extreme computational scalability by exploiting an algebraic representation of Direct Acyclic Graphs that efficiently captures the totality of tasks involved in a computation and the flow of data between them. Its primary goal is to maximize parallelism while automatically orchestrating task execution so as to minimize both communication and load imbalance. Unlike other available DAG-based runtimes, the concise symbolic representation of the algorithm that DAGUE uses minimizes the memory required to express the map of tasks ; at the same time, it provides extreme flexibility during the scheduling process. This algebraic representation allows the DAG to be traversed at very high speed, while tracking any flow of data from task to task. By combining this underlying mechanism with an understanding of the specific hardware capabilities of the target architecture, DAGUE is able to schedule tasks in ways that creatively adapt alternative work-stealing strategies to the unique features of the system. These capabilities enable DAGUE to optimize data movement between available computational resources, including both different nodes of the full system and different accelerators on the same node.

The DAGUE runtime aimed first at providing a scheduler for large distributed system of multicore nodes and is able to handle heterogeneous architectures to relocate the most compute intensive kernel on GPUs. A full description of DAGUE, and the implementation of classical linear algebra factorizations in this environment, can be found in [START_REF] Bosilca | DAGuE : A generic distributed DAG engine for High Performance Computing[END_REF], [START_REF] Bosilca | Flexible development of dense linear algebra algorithms on massively parallel architectures with DPLASMA[END_REF].

STARPU [START_REF] Augonnet | StarPU : A unified platform for task scheduling on heterogeneous multicore architectures[END_REF] is a software tool aiming to allow programmers to exploit the computing power of the available CPUs and GPUs, while relieving them from the need to specially adapt their programs to the target machine and processing units. The STARPU run-time supports a task-based programming model. Applications submit computational tasks, with CPU and/or GPU implementations, and StarPU schedules these tasks and associated data transfers on available CPUs and GPUs. The data that a task manipulates is automatically transferred among accelerators and the main memory, so that programmers are freed from the scheduling issues and technical details associated with these transfers. STARPU takes particular care of scheduling tasks efficiently, using well-known algorithms from the literature. In addition, it allows scheduling experts, such as compiler or computational library developers, to implement custom scheduling policies in a portable fashion.

The main differences between DAGUE and STARPU are the tasks submission process, the centralized scheduling and the data movement strategy. DAGUE uses its own parametrized language to describe the DAG in comparison to the simple sequential submission loops used by STARPU. Therefore STARPU relies on a centralized strategy which analyzes the dependencies between tasks and schedules these tasks on the available resources, while each computational unit of DAGUE immediately release the dependencies of the completed task solely using the local knowledge of the DAG. At last, STARPU scheduling strategy exploits cost models of the computation and data movements to schedule a task to the right resource (CPU or GPU) in order to minimize overall execution time.

However it has no data movement policy on shared memory systems resulting in lower efficiencies when no GPUs are used compared to the data-reuse heuristic of DAGUE.

Hence, the research around STARPU has focused mainly on the case of an heterogeneous multicore node enhanced with multiple GPUs, while research around DAGUE has more focused on scalability issues on a large number of homogeneous nodes.

PASTIX is a scientific library that provides a high performance parallel direct solver for very large sparse linear systems. Numerical algorithms are implemented in single or double precision (real or complex) using LL T or LDL T factorizations for symmetric matrices, and LU factorizations with static pivoting for non symmetric matrices having symmetric patterns. Theres is a version of PASTIX for multicore node architectures, which uses a hybrid MPI-thread programming to fully exploit the advantage of shared memory and to reduce the memory overhead. Direct methods are numerically robust methods, but very large three dimensional problems may require a large amount of memory, even with any memory optimization. For this type of problems, PASTIX provides an adaptive blockwise incomplete factorization that is much more accurate (and numerically more robust) than the scalar incomplete factorizations which are commonly used as preconditioner for iterative solvers. Such incomplete factorization can take advantage of the latest breakthroughs in sparse blocked direct methods, and particularly should be very competitive in CPU time (effective power used from processors and good scalability), while avoiding the memory limitation encountered by direct methods.

III. ALGORITHMS

In this paper, we study and design algorithms and parallel programming models for implementing sparse supernodal direct methods on an emerging computers equipped with GPU accelerators. Our ultimate goal is to release a manycores (CPU+GPU) version of the PASTIX solver [START_REF]PaStiX : A High-Performance Parallel Direct Solver for Sparse Symmetric Definite Systems[END_REF]. In this paper, we consider the Cholesky, LDL T , and LU factorization algorithms which are already present in PASTIX.

The task dependencies of numerical factorization can be represented by a tree whose nodes represent computational tasks and whose edges represent transfer of data between the tasks. The distributed memory version of PASTIX uses a right-looking formulation which, having computed the factorization of a column-block corresponding to a node of the tree, immediately sends the data to update the column-blocks corresponding to its ancestors in the tree. In the actual implementation, we locally aggregate contributions to the same block before sending the contributions. This can significantly reduce the number of messages, and is known to limit the memory overhead induced by the direct methods.

Either a static or dynamic scheduling of block computations can be used, independently from these different approaches described above. For homogeneous parallel architectures, it is useful to have an efficient static scheduling scheme. For the PASTIX solver, we have recently developed a dynamic scheduling scheme, specifically designed for modern supercomputers, that have many microprocessors, each of which consists of one or many computational cores, and induces a strong hierarchical topology. To address the needs of dynamic scheduling, the data structures of the solver, as well as the patterns of communication, have been also modified [START_REF] Faverge | Dynamic scheduling for sparse direct solver on NUMA architectures[END_REF]. Thanks to these efforts on the multicore implementation [START_REF] Hénon | On using an hybrid MPI-Thread programming for the implementation of a parallel sparse direct solver on a network of SMP nodes[END_REF], we have a platform ready to deal with bigger problem sizes on today's supercomputers.

The main contributions of our current paper can be subdivided into three subtasks described hereafter :

-Kernels for sparse factorization. Our first step of developing a sparse direct solver for clusters of multi-GPU nodes is to develop supernodal computational kernels designed specifically for this purpose. In particular, one of the most important kernels is the matrix/matrix-product involved in the updates of the trailing matrix. In the case of supernodes, data structures can be compacted in order to improve the efficiency of the computations, but we still need to improve the matrix-add operations that correspond with each block of the sparse updates. This kernel must also prefetch the data to minimize the data transfers on the GPU devices. -Scheduling of computations and data transfers. The main aim of this paper is to evaluate the potential of replacing the current scheduler of PASTIX with two generic runtime frameworks DAGuE and STARPU for executing the task graph corresponding to a sparse factorization.

In the current scheduling scheme of PASTIX, a task corresponds to a supernode (1D block distribution), see Fig. 1. To improve the efficiency of the sparse factorization on a multicore implementation, we introduced a way of controlling the granularity of the BLAS operations (referred to as an ESP option for Enhanced Sparse Parallelism). This functionality dynamically splits a single task of computing the contribution to the trailing submatrix, using the current panel into subtasks, so that the critical path of the algorithm can be reduced.

In this paper, for both DAGUE and STARPU scheduling, one computational task corresponds to the computation FIGURE 1: Task graph for supernodal factorization of the contributions to a remaining column-block using a single block in a panel. Hence, the number of tasks to schedule is equal to the number of blocks in the factored matrix. -Sparse algorithms. In order to control the granularity of the computational tasks, the automatic criteria used to set adaptive block sizes has to be extended to heterogeneous architectures. Some references for benchmarking dense linear algebra kernels are described in [START_REF] Volkov | Benchmarking GPUs to Tune Dense Linear Algebra[END_REF] and show that efficiency could be obtained on GPU devices only on relatively large blocks -that can be found on top of the elimination tree. Similarly, the amalgamation algorithm [START_REF] Hénon | On finding approximate supernodes for an efficient ILU(k) factorization[END_REF], reused from the implementation of an incomplete factorization, is a crucial step to obtaining larger supernodes and efficiency on GPU devices. The default parameter for amalgamation has been slightly increased ; we allow up to 12% more fill-in to build larger blocks while keeping a good a high level of parallelism.

A. Sparse GEMM on GPU Fig. 2a illustrates the data structure of PASTIX storing the sparse matrix. The nonzero entries are grouped into supernodes which are composed of several dense blocks. Each column block is stored as a dense block -represented by S 1 storage in the figure. When we update the column block S 2 facing the first extradiagonal block of S 1 , all blocks of S 1 are included in the blocks of S 2 . Thus, in PASTIX, we first factorize the diagonal block of S 1 , and we then update the off-diagonal blocks of S 1 . Next, for each extradiagonal block B i of S 1 , we compute the contribution B k B T i to the (k, i)-th block for k > i, and store the contribution in a temporary buffer. Finally, we update S i for each extradiagonal block B i . On a GPU, the computation is exactly the same except that our buffer is in the shared memory of the GPU.

To have an efficient usage of the GPU during sparse factorization, we wrote a new GPU kernel that takes the sparsity of our matrix into account. Our kernel is based on the framework proposed in [START_REF] Kurzak | Autotuning gemm kernels for the fermi gpu[END_REF], which also allows us to use an auto-tuning script to select near-optimal block sizes for the kernel. Specifically, this kernel computes C = αAB +βC, where C is divided into a 2D grid of tiles (see Fig. 2b), and each thread block computes a tile of the matrix-multiplication, AB, which is stored in the shared GPU memory. Each thread is in charge of computing several entries of this tile. Each of these entries is computed one by one, and to compute each entry, the corresponding row of A and the corresponding column of B are loaded into the thread registers. Once the tile is computed, it is added into C. In the supernodal solver PASTIX, C may span several discontinuous supernodal blocks in a columnblock (see Fig. 2b). Hence, during the summation of αAB into βC, we use an offset that specifies the row of C, to which the computed values of αAB should be accumulated into. For this, we added two pairs of integer arrays as input arguments to the kernel, which store the first and last row indexes of each blocks of A and C.

B. DAG scheduling

In this section, we describe the different approaches used to describe the DAGs in DAGUE and STARPU. For both codes, we define two kinds of tasks in our LL T factorization :

panel(s) : factorization of the diagonal block of the supernode s (POTRF) and update all the extradiagonal blocks of this supernode (TRSM) ; -gemm(b,s) : perform the GEMM product of all the blocks beneath B b in its supernode by B T b and subtract the result from the supernode s which has its diagonal block facing B b . For the LU factorization, we are using the same tasks with addition of the update of U in both kernels. The panel task does one extra solve on the upper part, and the gemm task does two matrix-matrix multiplications : one for the lower part and one for the upper part. Finally, the LDL t factorization follows the same scheme as Cholesky with different kernels to integrate the storage and the computation of the diagonal. The three factorizations have then the same DAG representation with different kernels, LU having more dataflow to move the upper part.

1) JDF representation of a DAG: In DAGUE, the data distribution and dependencies are specified using the Job Data Flow (JDF) format. Fig. 3 shows our JDF representation of the sparse Cholesky factorization using the tasks panel and gemm described previously. The second one is based only on the block id parameter and computes internally the supernode (fcblk) in which to apply the update. On Line 2 of panel(j), cblknbr is the number of block columns in the Cholesky factor. Once the j-th panel is factorized, the trailing submatrix can be updated using the j-th panel. This data dependency of the submatrix update on the panel factorization is specified on Line 6, where firstblock is the block index of the j-th diagonal block, and lastblock is the block index of the last block in the j-th block column. The output dependency on Line 7 indicates that the j-th panel is written to memory at the completion of the panel factorization. The input dependency of the j-th panel factorization is specified on Lines 4 and 5, where leaf is true if the j-th panel is a leaf in the elimination-tree, and lastbrow is the index of the last block updating the j-th panel. Hence, if the j-th panel is a leaf, the panel is read from memory. Otherwise, the panel is passed in as the output of the last update on the panel.

Similarly, gemm(k) updates the fcblk-th block column Hence, the data dependency of gemm(k) is resolved once the cblk-th panel is factorized, and the fcblk-th column is updated using the prev-th block. Notice that the diagonal blocks are not used to update the trailing submatrix, but it is included in the code to have a continuous space of execution for the task required by DAGUE. Finally, Lines 7 through 9 specify the output dependencies of gemm(k), where next is false if the k-th block is a diagonal block, and it is the index of the next block after the k-th block in the fcblk-th row otherwise. Hence, the completion of gemm(k) resolves the data dependency of the fcblk-th panel factorization if this is the last update on the panel, or it resolves the dependency of updating the fcblk-th block column using the next-th block otherwise.

2) STARPU tasks submission: STARPU builds its DAG following the tasks ordering provided by the user and by using data dependencies. The following pseudocode shows the STARPU tasks submission loop for the LL T decomposition. for all extra diagonal block b i of S 1 do 4:

S 2 ← supernode in f ront of (B i ) 5: submit gemm(S 1 , S 2 ) {sparse GEMM B k,k≥i × B T i substracted from S 2 } 6:
end for 7: end for By default, STARPU generates the tasks graph following sequential consistency (Fig. 4a). Since the order in which the facing supernode receives contribution is not relevant in our factorization, we disabled STARPU sequential consistency (Fig. 4b). Then, we notify the scheduler that a panel task depends on all the GEMM updates on the supernode and that this task triggers the GEMM on the trailling supernodes.

We are not using a reduction on the target panel, because for each update only a subarea is impacted. Moreover, a reduction requires a copy of the destination supernode on each computing device, that could represent a lot of memory overhead. This allows only an out of order execution of the updates while guarantying the mutual exclusion of the computations. DAGUE in its actual release doesn't allow those reduction operations and keep the order defined by the dependencies. As the order of the reduction is imposed by the data dependencies, the potential concurrency is reduced but the accuracy of the numerical results is maintained between successive runs.

C. Sparse mapping

Our mapping algorithm is based on a static scheduling based on a performance model. Thus, the partitioning and mapping step generates a fully ordered schedule that can be used in the parallel factorization. This schedule aims at statically regulating all of the issues that are classically managed at runtime. To make our scheme very reliable, we estimate the workload and message passing latency by using a BLAS and communication network time model, which is automatically calibrated on the target architecture.

Unlike usual algorithms, our partitioning and distribution strategy is divided into two distinct phases. The partitioning algorithm is based on a recursive top-down strategy over the block elimination tree provided by block symbolic factorization. Pothen and Sun presented such a strategy in [START_REF] Pothen | A mapping algorithm for parallel sparse Cholesky factorization[END_REF]. The partition phase splits column-blocks associated with large supernodes, and then for each column-block, it builds a set of candidate threads for its mapping. Once the partitioning step is over, the task graph is built. In this graph, each task is associated with the set of candidate threads for its columnblock. The mapping and scheduling phase then try to optimally map each task onto one of these sets of threads. An important constraint is that once a task has been mapped to a thread, then all the data accessed by this thread are also mapped on the process associated with the thread. This means that an unmapped task that accesses a block that has already been mapped should be mapped on the same thread to preserve the data locality.

IV. RESULTS

To compare the factorization times of PASTIX using generic runtime systems DAGUE and STARPU, with that using the internal scheduler of PASTIX, we present experimental results on challenging matrices from industrial applications. Some properties of the three different test matrices used for our scaling studies are shown in Table I. For all three approaches, the same partitioning and mapping were used, while the fully ordered schedule was used only for experiments using the internal scheduler of PASTIX. 256GB of total RAM, and is equipped with two Tesla T20 GPUs.

Fig. 5 and 6 compare the performance of PASTIX using DAGUE and STARPU runtimes, with PASTIX using its original scheduler on Audi and MHD test matrices. A good scalability was obtained using all three approaches. On these two test cases, generic schedulers behaved quite well, and obtained the performances similar to the fine-tuned PASTIX scheduler. On a small number of processors, we could even obtain better results using the generic schedulers. DAGUE lost some performance when more than a socket (12 cores) was used, but recovered as soon as the computation spawned over more than 2 sockets. The DAGUE scheduler seemed to not be able to extract any performance gain between 12 and 24 cores.

With the 10 Millions test matrix (Fig. 7), the finely-tuned scheduler of PASTIX outperformed the generic runtimes. Specifically, STARPU showed its limitation as the number of threads increased, while DAGUE maintained a good scalability. This can be explained by the fact that STARPU does not take data locality into account.

B. Experiments with GPUs

Fig. 8 and9 show the results of using one and two GPUs. For these experiments, we did not use all 48 available cores for computation as one core was dedicated to each GPU. The CUDA kernel gave good acceleration on both single and double precision. The factorization time was reduced significantly using the first GPU when the number of cores was small. The speedups of up to 5 was obtained using one GPU. However, the second GPU was relevant only with a small number of cores (less than 4 threads in single precision (Fig. 8), and less than 12 threads in double precision (Fig. 9)). With one GPU, once the cores on a socket (12 core on Romulus, and 6 on Mirage) were fully utilized, the GPU had no effect on the factorization time.

We also conducted additional GPU experiments using a compute node of the Mirage machine from INRIA -Bordeaux. Mirage nodes are composed of 2 Hexa-core (Westmere Intel® Xeon® X5650), with 36GB of RAM. The results on Mirage (Fig. 10) were similar to the ones on Romulus. V. CONCLUSION

In this paper, we examined the potential benefits of using generic runtime systems, DAGUE and STARPU, in a parallel sparse direct solver PASTIX. The experimental results using up to 48 cores and two NVIDIA Tesla GPUs demonstrated the potential of this approach to design a sparse direct solver on heterogeneous manycore architectures with accelerators using a uniform interface.

Through the study presented in this paper, we have identified three potential research path, which we plan to continue to investigate in the future.

First, in order to minimize the overhead induced by the scheduler, we need to increase the granularity of the tasks at the bottom of the elimination tree. These leaves or subtrees may be merged into bigger tasks to achieve this goal.

Second, we would like to pursue a similar experiment in distributed heterogeneous environments, composed by manycores nodes with multiGPUs. On such an platform, when a supernode update another non-local supernode, the update blocks are stored in a local extra-memory space (this is called "fan-in" approach). In order to reduce communication time overhead, we delay sending these updates until the last updates to the supernodes are accumulated, trading latency for bandwidth. We will study potential approaches to implement such an challenging optimization using generic runtime systems, STARPU and DAGUE.

More specifically, in the context of STARPU, for successive solution steps, the performance models used by our direct solver for the static mapping and scheduling step, could benefit from an online model refined during the execution.

Finally, the availability of extra computational resources, highlight the potential to dynamically build or rebuild the supernodal structures according to the load on the cores and the GPUs. A first approach will be to adjust the task granularity at runtime. Simultaneously, we will work on the distributed version of our solver, and on the challenging problem of defining an initial mapping of the data compatible with heterogeneous capabilities of the distributed memory architectures.

  FIGURE 2: Data management

1: for all Supernode S 1 do 2 :

 2 submit panel(S 1 ) {update of the panel} 3:

FIGURE 4 :

 4 FIGURE 4: Task graph

FIGURE 9 :

 9 FIGURE 9: GPU results on Romulus (double precision)

  the k-th block, where fcblk is the index of the block row that the k-th block belongs to, and blocknbr on Line 2 is the number of blocks in the Cholesky factor. The input dependencies of gemm are specified on Lines 4 through 6, where the cblk-th panel A is being used to update the fcblkth column C. Specifically, on these lines, diag is true if the k-th block is a diagonal block, and it is false otherwise, and prev is false if the k-th block is the first block in the fcblkth block row, and it is the index of the block in the block row just before the k-th block otherwise. Hence, the prev-th block updated the fcblk-th column just before the k-th block does.

	1. panel(j) [high priority = on]	1. gemm(k)
	2. j = 0 .. cblknbr-1	2. k = 0 .. blocknbr-1
	3. ... set up parameters for the j-th task ...	3. ... set up parameters for the k-th task ...
	4. :A(j)		4. :A(fcblk)
	5. RW A ← ( leaf ) ? A(j) : C gemm( lastbrow )	5. READ A ← diag ? A(fcblk) : A panel(cblk)
	6.	→ A gemm(firstblock+1..lastblock)	6. RW C ← first ? A(fcblk) : C gemm( prev )
	7.	→ A(j)	7.	→ diag ? A(fcblk)
			8.	→ (( !diag) && (next == 0)) ? A panel( fcblk )
			9.	→ (( !diag) && (next != 0)) ? C gemm( next )
		(a) Panel factorization		(b) Trailing submatrix update
		FIGURE 3: JDF representation of Cholesky.

using

  The first set of experiments were conducted on the Romulus machine at the university of Tennessee. Romulus has four twelve-core AMD Opteron 6180 SE CPUs (2.5 GHz) with

		Name	N		NNZ A		Fill ratio	OPC		Type	Factorization	Source
		MHD Audi	485,597 943,695		12,359,369 39,297,771		61.20 31.28	9.84e+12 5.23e+12		Real 1,600 Real	LU LL T	University of Minnesota PASTIX PARASOL Collection
		10M		10,423737		89,072,871		75.66	1.72e+14	Complex	LDL T	PASTIX with STARPU French CEA-Cesta
											1,400	PASTIX with DAGUE
										Factorization Time (s)	600 800 1,000 1,200
											400
											200
											0
											1	2	4	6	12	24	36	48
											Number of Threads
										FIGURE 6: LU decomposition on MHD (double precision)
											40,000	PASTIX PASTIX with STARPU
											PASTIX with DAGUE
	Factorization Time (s)	400 600 800				PASTIX PASTIX with STARPU PASTIX with DAGUE	Factorization Time (s)	10,000 20,000 30,000
		200									0	1	2	4	6	12	24	36	48
											Number of Threads
										FIGURE 7: LDL T decomposition on 10M (double complex)
		0
		1	2	4	6	12	24	36	48
				Number of Threads			A. Experiments on multicores
	FIGURE 5: LL T decomposition on Audi (double precision)

TABLE I

 I

							: Matrices description
		500
						0 CUDA device 1 CUDA devices		800				0 CUDA device 1 CUDA devices
		400				2 CUDA devices		700				2 CUDA devices
	Factorization Time (s)	200 300						Factorization Time (s)	300 400 500 600
		100							200
									100
		0							0
		1	2	4	6	12	24	36	1	2	4	6	12	24	36
				Number of Threads					Number of Threads
			(a) Audi LL T decomposition				(b) MHD LU decomposition
					FIGURE 8: GPU results on Romulus (single precision)

http ://icl.eecs.utk.edu/morse

http://www.matrixprosoftware.com/

http://www.acceleware.com/matrix-solvers

http://www.boeing.com/phantom/bcslib/ RR n°7972

RR n°7972

ACKNOWLEDGMENT

The authors would like to thank the DAGUE and STARPU teams for their support and assistance with this project. Special thanks also go to Sam Crawford for the valuable comments on our submitted manuscript and to Abdou Guermouche for his advice.

R ÉF ÉRENCES