Christophe Bovet
email: christophe.bovet@onera.fr

Pierre Gosselet
email: gosselet@lmt.ens-cachan.fr

Nicole Spillane
email: nicole.spillane@cmap.polytechnique.fr

Multipreconditioning for nonsymmetric problems: the case of orthomin and biCG

published or not. The documents may come

Résumé

Les solveurs de Krylov préconditionnés [START_REF] Saad | Iterative methods for sparse linear systems, 2nd Edition[END_REF] sont des outils performants pour la résolution de systèmes linéaires. Il arrive cependant qu'ils convergent lentement, souvent après une phase de stagnation. Une manière naturelle de remédier à cette situation est d'agrandir l'espace dans lequel on cherche la solution à chaque itération. En suivant cette idée, nous proposons dans cette note deux algorithmes multipréconditionnés : orthomin multipréconditionné (MPorthomin) et biCG multipréconditionné (MPbiCG) avec l'objectif de résoudre des systèmes linéaires généraux en un petit nombre d'itérations. Après avoir décrit les nouveaux algorithmes, nous illustrons leur comportement sur des systèmes linéaires issus de la méthode de décomposition de domaine FETI, où au lieu d'appliquer le préconditionneur habituel, on fait agir chaque contribution locale au préconditionneur séparément.

Here is the title in French. Multipréconditionnement pour les problèmes non symétriques : le cas des solveurs orthomin et biCG.

Version française abrégée

L'objectif de cette note est de présenter brièvement deux nouveaux algorithmes pour la résolution de systèmes linéaires non symétriques. Ils appartiennent à la famille des solveurs de Krylov multipréconditionnés. Les premiers solveurs de ce type sont le gradient conjugué multipréconditionné (MPCG) [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF], FETI simultané (SFETI) [START_REF] Rixen | Substructuring and Dual Methods in Structural Analysis[END_REF] et GMRES multipréconditionné (MPGMRES) [START_REF] Greif | MPGMRES: a generalized minimum residual method with multiple preconditioners[END_REF]. Comme nous considérons les problèmes non symétriques, MPCG ne s'applique pas. MPGMRES, dans sa forme originale, conduit à une croissance exponentielle de la taille de l'espace de minimisation. Comme remarqué dans [START_REF] Greif | MPGMRES: a generalized minimum residual method with multiple preconditioners[END_REF] il est possible de se limiter à une croissance linéaire par rapport au nombre d'itérations en ne sélectionnant que certains vecteurs. Ici, nous considérons des solveurs pour lesquels cette croissance linéaire est plus naturelle. En effet, à la différence de GMRES, orthomin et biCG rendent le résidu disponible explicitement à chaque itération et c'est à celui-ci que l'on applique le (ou les) préconditionneur(s). Ceci explique la croissance linéaire de la taille de l'espace de minimisation.

Nous notons Ax * = b le système linéaire considéré, où x * ∈ R n est la solution recherchée, A est une matrice de taille n × n inversible et b est le second membre (source donnée). Une famille de N préconditionneurs (qui sont aussi des matrices inversibles de taille n × n) notée {H s } s=1,...,N est utilisée. Les algorithmes orthomin multipréconditionné (MPorthomin) et biCG multipréconditionné (MPbiCG) sont présentés dans les Algorithmes 1 et 2 respectivement.

Les propriétés qui définissent les coefficients dans MPorthomin sont : -(A A)-orthogonalité des directions de recherche (Q i Q j = 0 pour tous j = i).

x i minimise la A A-norme de l'erreur x * -y A A parmi tous les y ∈ x 0 + k<i span(P k ).

Les propriétés qui définissent les coefficients dans MPbiCG sont :

-bi-orthogonalité des directions de recherche ( P i Q j = 0 si i = j), -r i P j = 0 et r i Pj = 0 si j < i. Ceci n'est pas une propriété de minimisation. Le multipréconditionnement est particulièrement intéressant dans le cas des méthodes de décomposition de domaine où le préconditionneur usuel possède une structure additive. Il revient alors à améliorer la combinaison des termes du préconditionneur classique, ce qui peut éviter de propager des contributions locales problématiques pour la convergence [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF]. Dans le cadre symétrique défini positif, une alternative pour se prémunir contre ces mauvaises contributions est de les détecter préalablement à l'aide de problèmes aux valeurs propres généralisés et de les éliminer par un solveur de Krylov augmenté [START_REF] Spillane | Automatic spectral coarse spaces for robust FETI and BDD algorithms[END_REF]. C'est une force du multipréconditionnement de se déployer simplement sur des cas plus généraux.

Nous appliquons ces algorithmes à des systèmes linéaires issus de la méthode de décomposition de domaine FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. Le cas test utilisé est une plaque lamifiée, composée alternativement de plis à comportements élastiques et élastoviscoplastiques non associés [START_REF] Besson | [END_REF] (Figure 1). Les algorithmes multipréconditionnés sont moins sensibles à l'hétérogénéité des matériaux. Ils convergent avec un faible nombre d'itérations bien que le manque de stabilité de MPbiCG soit observé pour E 2 /E 1 10 4 , et qu'une orthogonalisation de Gram-Schmidt modifiée avec réorthogonalisation soit nécessaire pour MPorthomin (Figure 2).

Introduction

The objective of this note is to present for the first time two algorithms for solving nonsymmetric linear systems. They belong to the family of multipreconditioned Krylov subspace solvers of which MPCG (multipreconditioned CG) [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF], SFETI (Simultaneous FETI) [START_REF] Rixen | Substructuring and Dual Methods in Structural Analysis[END_REF] and MPGMRES (multipreconditioned GMRES) [START_REF] Greif | MPGMRES: a generalized minimum residual method with multiple preconditioners[END_REF] are the first. Since we consider nonsymmetric problems, MPCG is out of the question. MPGMRES, in its original form, leads to an exponential increase in the size of the minimization space.

Although it is pointed out in [START_REF] Greif | MPGMRES: a generalized minimum residual method with multiple preconditioners[END_REF] that this can be decreased to a linear growth by selecting some vectors, we propose to look at methods where this is naturally the case. Indeed, one significant difference between GMRES and both biCG and orthomin is that in the last two methods the residual is readily available at each iteration and the preconditioner (or set of preconditioners) is applied to it. This explains the linear growth of the size of the minimization space with respect to the number of iterations.

Multipreconditioning makes particular sense in domain decomposition approaches where the usual preconditioners possess an additive structure. In that case, multipreconditioning simply amounts to finding, at each iteration, an improved (or even optimal depending on the solver) combination of the terms of the classical preconditioner. It was explained in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] that such a procedure prevents the propagation of problematic local contributions which are known to penalize convergence. The alternative to prevent these bad contributions is to pre-detect them thanks to localized generalized eigenvalue computations and to remove them by augmenting the Krylov solver [START_REF] Spillane | Automatic spectral coarse spaces for robust FETI and BDD algorithms[END_REF]. Unfortunately, the theory of pre-detection is currently limited to symmetric positive semi-definite problems; it is one strength of multipreconditioning approaches to extend seamlessly to more general cases.

Multipreconditioned orthomin and multipreconditioned biCG

We consider the problem of finding x * ∈ R n such that Ax * = b where A is a non singular n × n matrix and b is a given right hand side. A family of N preconditioners (also nonsingular n × n matrices) denoted by {H s } s=1,...,N is used. The multipreconditioned orthomin (MPorthomin) and multipreconditioned biCG (MPbiCG) algorithms are presented in Algorithms 1 and 2 respectively.

The two properties that define the coefficients in MPorthomin are:

-the (A A)-conjugacy of the search directions (i.e., Q i Q j = 0 for any j = i).

-the minimization property that x i minimizes the A A-norm of the error x * -y A A over all y ∈ x 0 + k<i span(P k ) (i.e., r i Q j = 0 for any j < i).

The properties that define the coefficients in MPbiCG are:

-the biconjugacy of search directions (i.e., P i Q j = 0 if i = j), -r i P j = 0 and r i Pj = 0 if j < i. Note that this is not a minimization property: with MPbiCG, as with biCG, there is no error norm that is guaranteed to decrease at every iteration. MPbiCG is expected to inherit the lack of stability of biCG (an issue addressed in biCGstab of which we did not manage to derive a multipreconditioned version), but there exist many configurations were it remains competitive in spite of the need to store twice as many vectors at each iteration. We encountered unexpected poor conditioning of the orthogonalization step of MPorthomin; as a first answer we used a modified Gram-Schmidt method with reorthogonalization which allowed to achieve convergence, at the cost however of more operations and communications per iteration. We also note that the preconditioning step in MPorthomin could have been

Z i+1 = H 1 (Q i α i ) | . . . | H N (Q i α i )
leading to an MPorthodir algorithm.

Numerical results: application to the FETI method

The algorithms presented in the previous sections have been applied to the FETI method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. FETI is a domain decomposition-based iterative method: the domain is split into reasonably small subdomains where local solves are performed using direct solvers, and an iterative solver is used to find the Lagrange multipliers that glue all the subdomains together. FETI with multipreconditioned CG was studied in [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] Algorithm 1: Multipreconditioned orthomin (MPorthomin)

x 0 given ; // Initial Guess r 0 = b -Ax 0 ; Z 0 = H 1 r 0 | . . . | H N r i+1 ; P 0 = Z 0 ; // Initial search directions Q 0 = AP 0 ; for i = 0, 1, . . . , convergence do α i = (Q i Q i ) -1 (Q i r i ); x i+1 = x i + P i α i ; // Update approximate solution r i+1 = r i -Q i α i ; // Update residual Z i+1 = H 1 r i+1 | . . . | H N r i+1 ; // Multi Precondition β i,j = (Q j Q j ) -1 (Q j (AZ i+1 )), j = 0, . . . , i; P i+1 = Z i+1 - i j=0 P j β i,j ; // Orthogonalize Q i+1 = (AZ i+1 ) - i j=0 Q j β i,j ; end Return x i+1 ;
Algorithm 2: Multipreconditioned biconjugate gradient algorithm (MPbiCG)

x 0 given ; // Initial Guess r 0 = b -Ax 0 and r0 = r 0 ; Z 0 = H 1 r 0 | . . . | H N r i+1 ; P 0 = Z 0 ; // Initial search directions Ẑ0 = H 1 r0 | . . . | H N ri+1 ; P0 = Ẑ0 ;
// Initial conjugate search directions

for i = 0, 1, . . . , convergence do Q i = AP i and Qi = A Pi ; α i = ( P i Q i ) -1 ( P i r i
) and αi = (P i Qi ) -1 (P i ri );

x i+1 = x i + P i α i ;
// Update approximate solution r i+1 = r i -Q i α i and ri+1 = ri -Qi αi ;

// Update residual

Z i+1 = H 1 r i+1 | . . . | H N r i+1 and Ẑi+1 = H 1 ri+1 | . . . | H N ri+1 ; // Multi Precondition β i,j = ( P j Q j ) -1 ( Q j Z i+1
) and βi,j = ( P j Q j ) -T (Q j Ẑi+1 ) for j = 0, . . . , i; P i+1 = Z i+1 -i j=0 P j β i,j and Pi+1 = Ẑi+1 -i j=0 Pj βi,j ; // Orthogonalize end Return x i+1 ; where a trick was proposed to limit extra computations caused by the presence of the coarse problem which globally balances floating subdomains.

The methods have been implemented in the Finite element suite ZSet 8.6 (www.zset-software.com). The test case aims at representing a laminated material (a soft material reinforced by a stiffer one): a five-layer plate occupies the domain Ω = [0, 10] × [0, 5] × [0, 0.5]. It is clamped on the x = 0 side and a displacement is imposed on the opposite side. Stiff layers are made of a non associated elastoviscoplastic material (which causes the nonsymmetry of the tangent systems). Soft layers follows a classical linear elastic constitutive law. Small displacements are assumed. The prescribed displacement linearly grows with time, the computation ends at t = 0.05. A Newton scheme with a constant time step 0.01 is used. It is small enough for all Newton steps to converge in 1 iteration. The mechanical problem is:

   div(σ(u)) = 0 in Ω u = xt 10000 e x for x ∈ {0, 10} constitutive law σ(u) • n = 0 on ∂Ω with x / ∈ {0, 10} (1)
The soft material follows an isotropic linear elastic constitutive rule characterized by the Young's modulus E and the Poisson's coefficient ν given in Fig. 1b. The non associated elastoviscoplastic constitutive law combines a nonsymmetric yield criterion, a Norton flow rule and a Ziegler kinematic hardening [START_REF] Besson | [END_REF]. It is defined by [START_REF] Saad | Iterative methods for sparse linear systems, 2nd Edition[END_REF], where ε vp is the viscoplastic strain tensor, X is the hardening variable, ṗ is the viscoplastic multiplier. J, J 2 , I 2 are tensor invariants and σ D is the deviatoric part of σ. All material parameters are given in Fig. 1b. In order to assess the robustness of the solvers, we make the ratio of the Young's moduli vary between 1 and 10 5 . Since a prescribed displacement is applied, the increase of the heterogeneity ratio produces more (incompressible) viscoplastic flow inside the stiff layers, which contributes to the bad condition number of the tangent systems.

                 σ(u) = E 1 + ν (ε(u) -ε vp ) + νE (1 + ν)(1 -2ν) Tr(ε(u) -ε vp )I with ε(u) = 1 2 ∇u + ∇u εvp = 3 2 σ D -X J 2 (σ -X) ṗ with ṗ = f K m and f = (1 -a)J(σ -X) + aT r(σ) X = 2 3 Cα and α = ṗ σ -X I 2 (σ -X) - D C X (2)
The problem is discretized into 150 × 75 × 15 regular twenty-node brick elements (c3d20), the problem size is 2, 158, 068 degrees of freedom (dofs). The mesh is split into 64 domains using the Metis automatic graph partitioner. Subdomains are composed on average of 38, 000 dofs and 2, 600 elements. A 10 -6 relative convergence threshold is used on the Euclidean norm of the residual. All computations use a Dirichlet preconditioner with a k-scaling and symmetric coarse grid projector. Figure 2 summarizes all the results. We compare classical GMRES and orthomin with our new algorithms. MPorthomin 2 stands for MPorthomin with a modified Gram-Schmidt with reorthogonalization (see [START_REF] Saad | Iterative methods for sparse linear systems, 2nd Edition[END_REF] p. 162.). In terms of Krylov iterations, all multipreconditioned algorithms are less sensitive to the increase in the heterogeneity. The lack of stability of MPbiCG is however observed for h 10 4 . The poor conditioning of the Orthomin orthogonalization step is visible for h = 10 4 . Thanks to the reorthogonalization, MPorthomin 2 achieves convergence even for h = 10 5 . When they converge, GMRES and orthomin remain faster than the multipreconditioned algorithms. This indicates that the best algorithm should mix multipreconditioned iterations (that help with robustness) and usual preconditioned iterations (that are less expensive).

Conclusion

We have presented two new multipreconditioned algorithms for possibly non symmetric linear systems that rely on the same ideas as the existing ones, and whose complexity only grows linearly with respect to the iteration number. We assessed them in combination with the FETI domain decomposition method on non-associated viscoplasticity problems. The next step is to derive adaptive variants such as [START_REF] Spillane | An Adaptive Multipreconditioned Conjugate Gradient Algorithm[END_REF][START_REF] Spillane | Algebraic Adaptive Multipreconditioning applied to Restricted Additive Schwarz[END_REF] (developed for MPCG), in order to achieve the best compromise between robustness and efficiency.

Figure 1 .

 1 Figure 1. Heterogeneous composite

5 GMRES

 5 Total # Krylov iter.Time to solution (s)h = E 2 /E 1 10 1 10 2 10 3 10 4 10 5 10 1 10 2 10 3 10 4 10 Performance for various heterogeneity ratios. A dash means that the solver did not converge. Evolution of the residual for h = 10 4 . All Newton steps are superposed.

Figure 2 .

 2 Figure 2. Heterogeneous composite. MPomin 2 uses a modified Gram-Schmidt orthogonalization with reorthogonalization.


Iacopo Vagliano
email: iacopo.vagliano@polito.it

Cristhian Figueroa
email: cristhian.figueroa@polito.it

Oscar Rodríguez Rocha

Marco Torchiano
email: marco.torchiano@polito.it

Catherine Faron-Zucker

Maurizio Morisio
email: maurizio.morisio@polito.it

ReDyAl: A Dynamic Recommendation Algorithm based on Linked Data *

Keywords: Recommender System, Linked Data, DBpedia, Semantic Web

The Web of Data is an interconnected global dataspace in which discovering resources related to a given resource and recommend relevant ones is still an open research area. This work describes a new recommendation algorithm based on structured data published on the Web (Linked Data). The algorithm exploits existing relationships between resources by dynamically analyzing both the categories to which they belong to and their explicit references to other resources. A user study conducted to evaluate the algorithm showed that our algorithm provides more novel recommendations than other state-of-the-art algorithms and keeps a satisfying prediction accuracy. The algorithm has been applied in a mobile application to recommend movies by relying on DBpedia (the Linked Data version of Wikipedia), although it could be applied to other datasets on the Web of Data.

INTRODUCTION

The Web is evolving from an information space for sharing textual documents into a medium for publishing structured data. Linked Data 1 is a set of best practices to publish and interlink data on the Web and it is the base of the Web of Data, an interconnected global dataspace where data providers publish their content publicly.

Due to the increase in the amount of structured data pub- * The mobile application presented in Section 5 was done at the Joint Open Lab MobiLAB and was supported by a fellowship from TIM. 1 http://linkeddata.org lished on the Web through the principles of Linked Data, it is more likely to find resources that describe or represent real life concepts. The information provided by these resources may be used in many different domains. However, finding and recommending related resources is still an open research area [START_REF] Ricci | Introduction to recommender systems handbook[END_REF]. The work presented in this paper holds on the results obtained from an our previous study and is its continuation [START_REF] Figueroa | A systematic literature review of linked data-based recommender systems[END_REF]. The study stated that the problem of finding existing relationships between resources can be addressed by analyzing the categories they belong to, their explicit references to other resources and/or by combining both these approaches. The study also showed that many works aimed at resolving this problem by focusing on a specific application domain and dataset. In this paper, we address this issue and we focus on the following research questions: (i) How can we design a recommendation algorithm that exploits existing relationships between resources on the Linked Data, is independent from the application domain and may be used on different datasets on the Web of Data? (ii) How can we design a recommendation algorithm that provides novel recommendations, i.e., recommendations of resources not previously known from the user, without affecting the prediction accuracy of the algorithm?

We propose a new algorithm based on Linked Data which exploits existing relationships between resources in order to recommend related resources. It dynamically analyzes the categories they belong to and their explicit references to other resources, then combines the results. The algorithm has been applied to DBpedia2 , but it could as well be applied to other datasets on the Web of Data and it is not bound to any specific application domain.

We conducted a user study to comparatively evaluate its accuracy and novelty against three state-of-the-art algorithms, which showed that our algorithm provides a higher number of novel recommendations, while keeping a satisfying prediction accuracy. An implementation of our recommendation algorithm has been integrated into a mobile application suggesting movies based on DBpedia, which was developed in collaboration with Telecom Italia, the major network operator in Italy.

The paper is organized as follows: Section 2 reviews related works; Section 3 presents our algorithm; Section 4 describes the evaluation method and provides the results; Section 5 shows the application of our algorithm for recommending movies; Section 6 provides conclusions.

RELATED WORK

This work began by conducting a systematic literature review [START_REF] Figueroa | A systematic literature review of linked data-based recommender systems[END_REF] on the subject, which allowed us to lay the groundwork for this research. Such review listed the different approaches to exploit Linked Data in order to recommend resources. Some studies found, infer relationships between resources by taking into account the existing links between them in a dataset, and use these relationships to measure the semantic similarity of the resources. Such relationships can be direct links, paths, or shared topics between sets of items. The most important related works are summarized in the following. Damljanovic et al. [START_REF] Damljanovic | Linked data-based concept recommendation: Comparison of different methods in open innovation scenario[END_REF] recommended experts in an open innovation scenario. Their approach, named HyProximity, takes as input a description of a problem in natural language and extracts a set of relevant words that are linked with resources of DBpedia. Then it generates recommendations by combining two techniques. The first one consists in discovering resources related through hierarchical relationships, while the second one is based on traversal relationships, which connect resources without establishing a classification or hierarchy. By exploiting these two kinds of relationships, the approach identifies a set of direct or indirect topics related with potential experts to solve an innovation problem.

Passant [START_REF] Passant | dbrec -Music Recommendations Using DBpedia[END_REF] described dbrec, a recommender targeted for the music domain, which mainly relies on a distance measure named Linked Data Semantic Distance (LDSD). It takes into account the number of direct or indirect links between resources (related with the music domain) represented in DBpedia. Unlike HyProximity it does not distinguish between traversal and hierarchical links. Both Damljanovic et al. and Passant had to reduce the set of resources and links of the dataset to those belonging to a specific domain (innovation problems and music respectively), which involves a huge effort to manually define which resources or links should be considered.

Other works combine Linked Data based algorithms with other techniques of recommendation in order to improve the results. These techniques include collaborative filtering [START_REF] Kabutoya | A topic model for recommending movies via linked open data[END_REF][START_REF] Mannens | Automatic news recommendations via aggregated profiling[END_REF][START_REF] Ostuni | Top-N recommendations from implicit feedback leveraging linked open data[END_REF][START_REF] Peska | Enhancing Recommender System with Linked Open Data[END_REF], information aggregation [START_REF] Cantador | Categorising social tags to improve folksonomy-based recommendations[END_REF][START_REF] Hamdan | Experiments with DBpedia , WordNet and SentiWordNet as re-sources for sentiment analysis in micro-blogging[END_REF][START_REF] Kitaya | Music Curator Recommendations Using Linked Data[END_REF] and statistical methods like Random Indexing (RI) [START_REF] Stankovic | Discovering Relevant Topics Using DBPedia: Providing Non-obvious Recommendations[END_REF], Vector Space Model (VSM) [START_REF] Baumann | Towards a storytelling approach for novel artist recommendations[END_REF][START_REF] Ostuni | Top-N recommendations from implicit feedback leveraging linked open data[END_REF], Latent Dirichlet Allocation (LDA) [START_REF] Khrouf | Hybrid event recommendation using linked data and user diversity[END_REF], implicit feedback [START_REF] Ostuni | Top-N recommendations from implicit feedback leveraging linked open data[END_REF] and structure-based statistical semantics [START_REF] Cheng | An Empirical Study of Vocabulary Relatedness and Its Application to Recommender Systems[END_REF]. De Graaff et al. [START_REF] De Graaff | Generic knowledge-based analysis of social media for recommendations[END_REF] proposed a knowledgebased recommender system that derives the user interests from the users social media profile, which is enriched with information from DBpedia. Musto et al. [START_REF] Musto | Automatic selection of linked open data features in graph-based recommender systems[END_REF] compared several techniques to automatically feed a graph-based recommender system with features extracted from Linked Data. However, these techniques usually require additional information from the user in order to produce accurate recommendations.

We propose a new recommendation algorithm, which is cross-domain and cross-dataset. It relies only on Linked Data and does not require to reduce the set of resources and links of the dataset to those belonging to a specific domain.

ReDyAl

ReDyAl is a recommendation algorithm which takes into account the different types of relationships between the data published according to the Linked Data principles. It aims at discovering related resources from datasets that may contain either well-linked resources as well as poorly-linked resources. A resource is said to be well-linked if it has a number of links higher than the average number of links in the dataset; otherwise it is poorly-linked. The algorithm is able to dynamically adapt its behavior in order to find a set of candidate resources to be recommended, relying on the implicit knowledge contained in the Linked Data relationships.

Principles

Any dataset on the Web of Data may be seen as a tuple (R, T, L) composed by resources (R), categories (T ), and relationships (L). Categories denote types, concepts or classes. Resources are instances of concepts; they are Web resources or real world resources identified by a URI. Relationships are also known as links or properties; they are the links connecting resources or categories along the whole dataset graph. Categories often are hierarchically organized. For example, DBpedia provides information about hierarchical relationships in three different classification schemata: Wikipedia Categories, YAGO3 [START_REF] Suchanek | Yago: A core of semantic knowledge[END_REF] classes, and WordNet Synsets4 . Relationships can be of three types:

Resource-Resource (R-R) These are the traversal relationships between resources, i.e. the links between resources that do not refer to hierarchical classifications.

Resource-Category (R-T ) These are relationships between a resource and a category. They can be represented by the RDF5 property rdf:type or the dcterms:subject property from the Dublin Core vocabulary 6 .

Category-Category (T-T ) These are hierarchical relationships between categories within a hyponymy structure (a category tree). They can be represented by using the RDFS7 property rdfs:subClassOf or the SKOS8 properties skos:broader (isSubCategoryOf) and skos:narrower (isSuperCategoryOf).

Considering this model of a dataset, ReDyAl consists in three stages:

1. The first stage discovers resources by analyzing the links between the given initial resource and other resources. Only R-R relationships are considered at this stage, although they can be indirect, i.e. they can connect two resources through a third one.

2. The second stage analyzes the categorization of the given initial resource and discovers similar resources located in the same categories. It finds indirect relationships between resources through direct R-T and T-T relationships. It is possible to specify to the algorithm which specific R-T and T-T relationships to consider in this step: the choice for R-T relationships is between dcterms:subject or rdf:type, while skos:broader and skos:narrower or rdfs:subClassOf are acceptable T-T relationships.

3. The last stage intersects the results of both the previous stages and ranks them by giving priority to those found in the first stage. The algorithm computes the similarity of the initial resource with respect to any of the discovered resources, based on a similarity function which combines the Linked Data Semantic Distance (LDSD) [START_REF] Passant | dbrec -Music Recommendations Using DBpedia[END_REF] and HyProximity distance [START_REF] Damljanovic | Linked data-based concept recommendation: Comparison of different methods in open innovation scenario[END_REF], opportunely adapted and generalized.

The algorithm can be applied to any dataset on the Web of Data. In the first step, it relies only on R-R relationships: any relationship of this kind may be used, independently of the data stored on the dataset. In the second step, the algorithm can be configured to use the dcterms:subject or rdf:type properties, which are R-T relationships. DBpedia uses both to enable different categorizations; for example to rely on the Wikipedia categories, it is necessary to set dcterms:subject as R-T relationship and skos:broader and skos:narrower as T-T relationships. Any other dataset uses at least rdf:type to indicate the class which a resource is instance of. Thus, rdf:type can be used to find resources in the same class and then rdfs:subClassOf can be used to retrieve more general classes (or skos:broader and skos:narrower, if the categories are organized through SKOS properties).

The algorithm is independent on the application domain because it relies only on R-R, R-T or T-T links. If in the dataset on which the algorithm is applied there are relationships among resources in different domains the algorithm may generate cross-domain recommendations. For example, DBpedia is a general dataset which represents resources of different kind and there may be a relationship between a song and a city because the song was recorded in that city, or because is about the city. Alternatively, there may be a link between a song and a movie because the song was part of the soundtrack of the movie. Thus, a city or a movie may be recommended starting from a song. Also R-T links may generates cross-domain recommendations if resources which belong to different domains are included into the same category.

Reducing the search space

Additionally, the algorithm may be configured with a set of forbidden links in order to restrict the kind of links the algorithm should consider. This is useful to prevent the algorithm to obtain resources over links pointing to empty nodes (i.e. resources without a URI), literals that are used to identify values such as numbers and dates, and other nodes that are not desired for the recommendation. In other words, it is a way to limit the results of the algorithm. For example the DBpedia resource dbr:Turin contains the link dbpprop:populationTotal that points to the integer value 911823: we can confgure this link as forbidden link since it does not point to a resource which can be recommended. This is also useful to increase the performance of the algorithm because limiting the number of results decreases the ranking time. All the links which are not explicitly specified as forbidden are allowed links and define a domain of interest. This may be useful when the algorithm is applied to a generic dataset as DBpedia. This dataset contains millions of links between resources, and if a developer is creating an application in the music domain then he/she may be interested only in resources of that domain, so he/she may want to consider only links pointing to those resources i.e., a set of allowed links. In fact the algorithm is cross-domain, thus it may recommend a city or a movie starting from a song, as we have already explained. While this may be an advantage in some applications, it may be confusing in others, especially if not properly explained to the user. To limit the recommendations to specific categories of resources (for example to consider only tracks and artists) it is sufficient to "allow" only the relationships which point to these kinds of resources, i.e. which have such desired category as range.

Parameter Settings

ReDyAl receives as input an initial resource by specifying its corresponding URI (inU RI), and three values (minT , minC, maxDistance) for configuring its execution. The selection of minT and minC is arbitrary and depends on the dataset and the convenience of the user who is setting up the algorithm. minT is the minimum number of links (input and output links involving the initial resource) necessary to consider a resource as well-linked. The proper value of minT depends on the dataset: if it contains resources with a high number of links between them it is expected to be higher, while if the resources have only few links it should be set to a lower value. However, this parameter impacts on the algorithm: if the initial resource is well-linked, traversal interlinking has a higher priority in the generation of candidate resources, otherwise the algorithm gives priority to the hierarchical relationships. For example, a user may consider the use of the hierarchical algorithms only if the resources are connected with less than 10 links by setting minT to 10. In a similar way, the user may arbitrary fix the value of minC, which is the minimum number of candidate resources that the algorithm is expected to generate, i.e. the number of candidate resources the user is expecting.

The value of maxDistance limits the distance (i.e. the number of hierarchical levels) that the algorithm considers in a category tree. maxDistance may be defined manually; this is particularly useful when there are not enough candidate resources from the categories found at a certain distance (i.e. the number of candidate resources retrieved is lower than minC). In this case, the algorithm increases the distances in order to find more resources and if the maxDistance value is reached with less than minC candidate resources, the algorithm ranks only the candidate resources found until that moment. Additionally, the algorithm may receive a list of forbidden links (FL) to avoid searching for candidate resources over a predefined list of undesired links. end if 25: end if 26: return CR ReDyAl (Algorithm 1) starts by retrieving a list of allowed links from the initial resource. Allowed links are those that are not specified as forbidden (F L) or that are explicitly defined in the initial resource. If there is a considerable number of allowed links (more than minT , i.e., the initial resource is well-linked) the algorithm obtains a set of candidate resources located through direct (DRl k ) or indirect traversal links (IRl k ), starting from the links explicitly defined in the initial resource (Lines 1-8). A resource is indirectly linked to the initial resource if it is linked through another resource. A resource directly linked is located at traversal distance 1 from the initial resource, while a resource indirectly linked is located at traversal distance 2 from the initial resource. With regards to the traversal links, a maximum distance of 2 is considered because for distances higher than 2 (i.e 1 direct heap plus 1 indirect heap) the number of retrieved resources is dramatically increased, therefore increasing also the number of resources that are not relevant or related with the initial resource.

Algorithm

Next, if the current number of candidate resources generated (CRtr) is greater than or equal to minC, the algorithm terminates returning the results (Lines 9-10). Otherwise, the algorithm generates a category graph (Gc) with categories of the first distance and applies iterative updates over the category graph over n distances from the initial resource, obtaining broader categories (i.e. more generic categories that are located in a higher level in a classification) until at least one of two following conditions is fulfilled: the number of candidate resources is sufficient (|CR| > minC), or the maximum distance is reached (currentDistance > maxDistance). At each iteration, candidate resources (CR hi ) are extracted from the broader categories of maximum distance (Lines 14-23). In any case, the algorithm combines these results with the results obtained in Lines 3-8 (adding CRtr and CR hi to CR). Finally, the set of candidate results is returned (Line 23).

Ranking of the recommended resources

The final operation is ranking the sets of candidate resources. The ranking process receives as input the candidate resources retrieved by the ReDyAl algorithm and ranks them according to their degree of similarity with the initial resource. This similarity is computed based on a combination of two distance measures: LDSD and HyProximity.

The LDSD distance, initially proposed by Passant [START_REF] Passant | dbrec -Music Recommendations Using DBpedia[END_REF], is based on the number of indirect and direct links between two concepts. In this measure, the similarity of two resources (r1, r2) is measured by combining four properties: the input/output direct links or the input/output indirect links between them. Equation 1 presents the basic form of the LDSD distance. Cdout is the number of direct output links (from c1 to c2), Cdin is the number of direct input links, Ciin is the number of indirect input links, and Ciout is the number of indirect output links. The implementation developed by Passant is limited to links from a specific domain, while the LDSD function implemented in ReDyAl takes into account all the concepts of the dataset unless forbidden links are specified.

LDSD(c1, c2) = 1 1 + Cdout + Cdin + Ciout + Ciin (1)
HyProximity is a similarity measure defined by Stankovic et al. [START_REF] Damljanovic | Linked data-based concept recommendation: Comparison of different methods in open innovation scenario[END_REF], which can be used to calculate both traversal and hierarchical similarities. The measure in its general form is shown in Equation 2 as the inverted distance between two concepts, balanced with a pondering function. In this equation d(r1, r2) is the distance function between the resources r1 and r2, while p(r1, r2) is the pondering function, which is used to weight different distances. Based on the structural relationships (hierarchical and traversal), different distance and pondering functions may be used to calculate the HyProximity similarity. ReDyAl reuses the HyProximity hierarchical measure, which is the quotient of a pondering function (p) and a distance (d). The distance was calculated using maxDistance such that: d(ir, ri) = maxDistance, where ir is the initial resource and ri is a candidate resources generated by the recommendation algorithm. The pondering function was calculated with an adaptation of the informational content function (Equation 3) defined by Seco et al. [START_REF] Seco | An Intrinsic Information Content Metric for Semantic Similarity in WordNet[END_REF]. In this equation hypo(C) is the number of descendants of the category C and |C| is the total number of categories in the category graph. This function was selected because it minimizes the complexity of calculation of the informational content, compared to other functions that employ an external corpus [START_REF] Hadj Taieb | New information content metric and nominalization relation for a new wordnet-based method to measure the semantic relatedness[END_REF]. Nonetheless, in ReDyAl, this measure is not limited to a specific property, and optionally can be configured to support a set of forbidden links.

hyP (r1, r2) = p(r1, r2) d(r1, r2) (2)
p(C) = 1 - log(hypo(C) + 1) log(|C|) (3)
Hybridsim = (1 -LDSD)α + (hyP (r1, r2))β (4)
Finally, the measure that combines LDSD and HyProximity used by ReDyAl is defined in Equation 4, where α and β may be set according to the convenience of the user: α is the weight for the traversal algorithm and β is the weight for the hierarchical algorithm. In this way, resources are ranked in descending order, arranged from the largest to the smallest value of Hybridsim.

USER EVALUATION

We comparatively evaluated the prediction accuracy and the novelty of the resources recommended with ReDyAl with respect to three state-of-the-art recommendation algorithms relying exclusively on Linked Data to produce recommendations: dbrec [START_REF] Passant | dbrec -Music Recommendations Using DBpedia[END_REF], HyProximity traversal and HyProximity hierarchical [START_REF] Damljanovic | Linked data-based concept recommendation: Comparison of different methods in open innovation scenario[END_REF]. This evaluation aimed to answer the following questions: (RQ1) Which of the considered algorithms is more accurate? (RQ2) Which of the considered algorithms provides the highest number of novel recommendations?

We decided to rely on a user study because we were interested in evaluating the novelty of proposed recommendations over the accuracy. Since we cannot expect that users rated all the items they already know, a user study can measure novelty more precisely than an offline study. On the other side, user studies are more expensive to conduct than an offline studies, for this reason we focus on recommendation algorithms based only on Linked Data and we did not consider algorithm exploiting traditional techniques, or combining Linked Data with traditional techniques. We plan to conduct other experiments to compare our method with other techniques and investigate on the effectiveness of our approach combined with traditional techniques.

Although our algorithm is not bound to any particular dataset, we applied it to DBpedia because it is a general dataset that offers the possibility to evaluate the results in a number of scenarios. DBpedia is one of the biggest datasets in the Web of Data and the most interlinked [START_REF] Schmachtenberg | Adoption of the linked data best practices in different topical domains[END_REF]. Furthermore, it is frequently updated and continuously grows.

Experiment

A user study was conducted involving 109 participants. The participants were mainly students of Politecnico di Torino (Italy) and University of Cauca (Colombia) enrolled in IT courses. The average age of the participants was 24 years old and they were 91 males, 14 females, and 4 of them did not provide any information about their sex. Although the proposed algorithm is not bound to any particular domain, this evaluation focused on movies because we aimed at applying our algorithm in the mobile application presented in Section 5 (which suggest movies) and in this domain a quite large amount of data is available on DBpedia. Additionally, it was easier to find participants, since no specific skills are required to express an opinion about movies. The algorithms were compared within subjects [START_REF] Shani | Evaluating recommendation systems[END_REF] since each participant evaluated recommendations from different algorithms, as it is explained in the following.

The evaluation was conducted as follows. A list of 20 recommendations generated from a given initial movie was presented to the participants. For each recommendation two questions were asked: (Q1) Did you already know this recommendation? Possible answers were: yes, yes but I haven't seen it (if it is a movie) and no. (Q2) Is it related to the movie you have chosen? Possible answers were: I strongly agree, I agree, I don't know, I disagree, I strongly disagree.

Each answer was assigned respectively a score from 5 to 1.

We developed a website9 to collect the answers from the participants. The participants were able to choose an initial movie from a list of 45 movies selected from the IMDB top 250 list 10 . The first 50 movies were considered and 5 movies were excluded because they were not available in DBpedia. Choosing these movies ensured participants to know them, but was also a limitation: the corresponding DBpedia resources are very well-linked, thus we could not properly evaluate the algorithm on poorly linked initial movies. The movies were presented to the user in a random order to avoid having most of the participants evaluating recommendations for the same initial movies (e.g. the first in the lists). When a participant selected an initial movie the tool provided the corresponding list of recommendations with the questions mentioned above. The recommendations were presented in a randomized order. Each participant was able to evaluate recommendations from as many initial movies as he wanted, but he had to answer the questions for all the recommendations, i.e. was not possible to answer only to part of the questions for the initial movie chosen. As a result, the recommendations of the lists for 40 out of 45 initial movies were evaluated by at least one participant and each movie was evaluated by an average of 6.18 participants. The dataset with the initial movies and the lists of recommendations is available online 11 .

Each list of 20 recommendations was pre-computed. In particular, recommendations were generated for each of the 45 initial movies with each of the four different algorithms.

Then, the recommendations generated by each algorithm were merged in a list of 20 recommendations to be shown to the participants. To do this, we generated a list of 40 recommendations by selecting the first 10 pre-computed recommendations for each algorithm and we ordered them by the similarity computed by each algorithms, since each algorithm ranks its recommendations by using its semantic similarity function with values between 0 and 1. Then we eliminated eventual duplicates, since the same recommendation could be provided by more than one algorithm. The final list was obtained considering the first 20 recommendations of the merged list.

With regard to the questions stated at the beginning of this section, to answer RQ1, the Root Mean Squared Error (RMSE) [START_REF] Shani | Evaluating recommendation systems[END_REF] was computed, and to answer RQ2 the ratio between the number of evaluations was computed in which the recommended item was not known by the participants and the total number of evaluations. For the RMSE measure, scores given by the participants when answering to Q2 were considered as reference and were normalized in the interval [0, 1], and these scores were compared with the similarities computed by each algorithms, since each algorithm ranks its recommendations by using its semantic similarity function.

Results

The results of the evaluation are summarized in Figure 1, which compares the algorithms with respect to their RMSE and novelty. The "sweet spot" area represents the conditions in which an algorithm has a good trade-off between novelty and prediction accuracy. In effect, presenting a high number of recommendations not known to the user is not necessarily good because it may prevent him to assess the quality of the recommendations: for example having in the provided recommendation a movie which he has seen and which he liked may increase the trust of the user in the RS.

Regarding RQ1, HyProximity accounts for the lowest RMSE measures (with 25% and about 36% for the hierarchical and traversal versions respectively), but these results are less significant due to the low number of answers to Q2 for these algorithms (this means that the RMSE was computed over a low number of recommendations). For both ReDyAl and dbrec the RMSE is roughly 45%. Concerning RQ2, the two versions of HyProximity account for the highest values (hierarchical roughly 99%, while traversal about 97%). However, such a high rate of novel recommendations may confuse the user and prevent him to judge recommendations, as we have already explained. ReDyAl has a larger rate of novel recommendations than dbrec. These two algorithms account respectively for about 60% and 45%.

The recommendations generated by HyProximity in both traversal and hierarchical version collected a low number of answers to Q2 because most of the recommendations generated by these algorithms were unknown as illustrated in Table 1. Consequently the RMSE was computed over a low number of recommendations. Thus, the results of these two algorithms related to RQ1 are less definitive than for the others, since for measuring the prediction accuracy only the evaluations for which the answer to Q1 was either "yes" or "yes but I haven't seen it (if it is a movie)" were considered.

We computed the Fleiss' kappa [START_REF] Fleiss | Measuring nominal scale agreement among many raters[END_REF] measure for assessing the agreement of the participants in answering Q2. We considered the recommendations and in particular we considered as different the same recommendation when related to a different initial movie (i.e. when appearing in different lists of recommendations). We excluded recommendations not evaluated or evaluated by only one participant. The Fleiss' kappa is 0.79; according to Landis and Koch [START_REF] Landis | The measurement of observer agreement for categorical data[END_REF], this corresponds to a substantial agreement.

In conclusion, Figure 1 illustrates that ReDyAl and dbrec provides a good trade-off between prediction accuracy and novelty (sweet spot area), although ReDyAl performs better in novelty. HyProximity hierarchical and HyProximity traversal seem to be excellent performers since the RMSE is low and the novelty is high, but the RMSE was computed on few evaluations. An additional analysis of these two algorithms is needed to verify if the user can benefit from such a high novelty and if novel recommendations are relevant.

In addition, further investigation is needed on poorly-linked resources, since the choice of the initial movies focused on selecting well known movies to make easier the evaluation from participants, but the related resources were well-linked.

On poorly-linked resources we expect ReDyAl and Hyproximity hierarchical keeping good recommendations since they can rely on categories, while dbrec and HyProximity traversal are likely to provide much less recommendations since they rely on direct links between resources.

MOBILE MOVIE RECOMMENDATIONS

An implementation of ReDyAl has been integrated into a mobile application developed in collaboration with Telecom Italia (the major network operator in Italy). This application recommends movies based on DBpedia: when the user enters the title of a movie, the application provides the Wikipedia categories to which the initial movie is related to. In this way, the user may focus on a specific scope and can receive recommendations of related resources for any category. In addition, it is possible to view any recommendation to obtain additional information.

Our algorithm can provide cross-domain recommendations because it is independent on the domain and is applied on DBpedia, which is a general dataset. Thus, the recommended resources can be movies but also other relevant entities such as actors, directors, places of recording, books on which the movie is inspired, etc. Other advantages of using DBpedia as dataset are the high number of resources that it represents, the variety of domains addressed and the continuous update and growth, since it is extracted from Wikipedia.

For example, given The Matrix as initial movie the categories which it belongs to are presented. The user may be more interested in martial arts, post-apocalyptic movies or he may prefer to consider all the movies from American directors, thus he can choose a category accordingly. a detailed view, which contains three tabs: the first contains a brief textual description, the second presents a graph view of the resource in order to show the main properties and the third summarizes the main information in a tabular form. The graph view is illustrated in Figure 2. The graph is paginated and few properties per page are presented in order to avoid information overload, since the resource can have a very high number of properties. This view can be useful also to explain the recommendation: for instance the user can understand that the recommended resource has the same director or the same main actor as the initial movie.

The graph view is based on DBpedia Mobile Explorer [START_REF] Vagliano | DBpedia Mobile Explorer[END_REF], a Linked Data visualization framework for the mobile environment, which enables the application to hide the underlying complexity of the Linked Data to the users by processing the resources to be presented received from DBpedia. The application is based on a client-server architecture and the main modules are DBpedia, a RESTful recommender service 12 which exposes our algorithm, and the mobile user interface. The main flow of interactions is represented in Figure 3. The mobile application asks for recommendations specifying an initial resource and optionally a scope such as a Wikipedia category [START_REF] Baumann | Towards a storytelling approach for novel artist recommendations[END_REF]. The recommender service answers with a list of scopes if no scope was provided or with a list of recommendations in the scope specified, otherwise [START_REF] Cantador | Categorising social tags to improve folksonomy-based recommendations[END_REF]. The recommender service relies on DBpedia to provide recommendations [START_REF] Cheng | An Empirical Study of Vocabulary Relatedness and Its Application to Recommender Systems[END_REF][START_REF] Damljanovic | Linked data-based concept recommendation: Comparison of different methods in open innovation scenario[END_REF] and the mobile application retrieves the resources to be visualized from the dataset [START_REF] De Graaff | Generic knowledge-based analysis of social media for recommendations[END_REF][START_REF] Figueroa | A systematic literature review of linked data-based recommender systems[END_REF]. The recommender service is developed in Java, while the client is an Android mobile application. The two modules use JSON as data-interchange format, while the mobile application retrieves resources from DBpedia serialized in JSON-LD 13 . The mobile application is going to be published on Google Play, but the Android Package (APK) of the first version is already available on the Web 14 .

Figure 3: The interactions between the main modules of the application

CONCLUSIONS AND FUTURE WORK

We presented ReDyAl which is a hybrid algorithm that dynamically uses both the traversal and hierarchical approach for discovering resources. It is independent from the application domain and, although we applied it to DBpedia, it could be easily adapted to other dataset in the Web of Data. It relies only on Linked Data and does not require to reduce the set of resources and links of the dataset to those belonging to a specific domain.

We evaluated and compared our algorithm against three state-of-the-art algorithms by conducting a user study and we also showed a practical application of the algorithm by presenting a mobile application that provides movie recommendations relying on DBpedia. Although the algorithm could be applied to other datasets in the Web of Data, we selected DBpedia because it is a general dataset, thus crossdomain recommendations were possible. In addition, there is a high number of resources represented, a variety of domains addressed and it is continuously updated, since it is extracted from Wikipedia. The user study demonstrated that ReDyAl improves in the novelty of the results discovered, although the accuracy of the algorithm is not the highest (due to its inherent complexity). Although ReDyAl is not bound to any particular domain, the study focused on movies as for this domain there is a quite large amount of data available on DBpedia and participants were not required to have specific skills.

Future work includes studying the relevance under different domains and improving the accuracy of ReDyAl while maintaining its novelty. We plan to conduct other studies to compare it with traditional techniques and with approaches which combine Linked Data with traditional techniques. We are also working on combining ReDyAl with collaborative filtering techniques in order to take user preferences into account while providing recommendations.

14 https://www.dropbox.com/sh/0q8d2mcbko9e2oj/ AAASh-YHGz0MmG_Z8hH6mfWOa?dl=0

Figure 1 :

 1 Figure 1: Prediction accuracy and novelty of the algorithms evaluated

Figure 2 :

 2 Figure 2: The graph view of V for Vendetta

Table 1 :

 1 Percentage of answers for Q1 by algorithm

	By se-

http://dbpedia.org

http://www.mpi-inf.mpg.de/yago-naga/yago/

https://wordnet.princeton.edu

http://www.w3.org/TR/rdf11-concepts/

http://dublincore.org/documents/dcmi-terms/

http://www.w3.org/TR/rdf-schema/

https://www.w3.org/TR/skos-reference

http://natasha.polito.it/RSEvaluation/

http://www.imdb.com/chart/top

http://natasha.polito.it/RSEvaluation/faces/ resultsdownload.xhtml

http://natasha.polito.it/LDRecommenderWeb/

http://json-ld.org/


Pacôme Eberhart
email: pacome.eberhart@upmc.fr

Julien Brajard
email: julien.brajard@upmc.fr

Pierre Fortin
email: pierre.fortin@upmc.fr

Fabienne Jézéquel
email: fabienne.jezequel@upmc.fr

Estimation of Round-off Errors in OpenMP Codes

Keywords: CADNA library, Discrete Stochastic Arithmetic, floatingpoint arithmetic, numerical validation, multi-core architectures, manycore architectures, OpenMP, round-off error

come

Introduction

The power of computational resources tends to increase and so does the number of floating-point operations performed in numerical simulations. Unfortunately each floating-point operation may generate a round-off error. These errors can accumulate and significantly affect the computed results. In parallel applications, more computations can lead to more round-off errors, and the order of computations can change with respect to the sequential execution or to another parallel execution: this can cause different round-off error propagations and hence different numerical results

The CADNA library5 [START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF][START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF] that implements DSA (Discrete Stochastic Arithmetic) [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF][START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF] enables one to estimate round-off errors in numerical simulation codes. In this paper we present a CADNA version for the numerical validation of codes using the OpenMP standard which is widely used for parallel numerical simulations. We describe how this new CADNA version has been developed taking into account constraints related to both OpenMP and DSA. We show the benefit of CADNA and analyse its cost in terms of execution time for various OpenMP programs: several benchmark codes and also a real-life application.

As far as related work is concerned, DSA is sometimes compared to interval arithmetic [START_REF] Alefeld | Introduction to interval analysis[END_REF][START_REF] Kulisch | Advanced Arithmetic for the Digital Computer[END_REF] that takes into accout round-off errors and computes, instead of each floating-point result, an interval with guaranteed bounds. However, with a naive implementation of interval arithmetic in a code, these bounds may overestimate the effective error. Therefore the use of interval arithmetic in an application requires specific algorithms and methods. Interval versions of linear algebra kernels have been proposed for multicore architectures [START_REF] Hölbig | High accuracy and interval arithmetic on multicore processors[END_REF][START_REF] Milani | Solving dense interval linear systems with verified computing on multicore architectures[END_REF][START_REF] Revol | Parallel implementation of interval matrix multiplication[END_REF][START_REF] Zimmer | Using C-XSC in a multi-threaded environment[END_REF] and have been implemented using OpenMP pragmas or a standard POSIX threads [START_REF]The Open Group Base Specifications[END_REF] library.

The rest of the paper is organized as follows. Section 2 presents the CADNA library and its use for the numerical validation of sequential codes. Section 3 describes the new OpenMP-enabled CADNA version developed to control the numerical quality of OpenMP codes. Section 4 presents, for various OpenMP applications, the interest of this CADNA version, as well as some performance and numerical results. Finally concluding remarks are given in Sect. 5.

The CADNA library

Principles of DSA (Discrete Stochastic Arithmetic)

The CADNA (Control of Accuracy and Debugging for Numerical Applications) library [START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF][START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF] that implements DSA (Discrete Stochastic Arithmetic) [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF][START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF] enables one to estimate round-off errors in any scientific code written in C, C++ or Fortran. CADNA uses a random rounding mode: at each elementary operation, the result is rounded up or down with the probability 0.5. The computer's deterministic arithmetic, therefore, is replaced by a stochastic arithmetic where each arithmetic operation is performed N times before the next one is executed, thereby propagating the round-off error differently each time. For each computed result, we obtain N samples R 1 , . . . , R N . The value of the computed result R is chosen to be the mean value of {R i } and, if no overflow occurs, the number of exact significant digits in R, C R , can be estimated as

C R = log 10 √ N R στ β with R = 1 N N i=1 R i and σ 2 = 1 N -1 N i=1 R i -R 2 .
τ β is the value of Student's distribution for N -1 degrees of freedom and a probability level 1 -β. In practice β = 0.05 and N = 3. Thus we can get a 95% confidence interval on the number of exact significant digits of the computed result. DSA is based on an hypothesis on the distribution of round-off errors.

It has been shown that in case this hypothesis is not rigourously satisfied the accuracy estimation is not altered if it is considered as exact up to one digit. The complete theory can be found in [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF]. The validity of C R is compromised if the two operands in a multiplication or the divisor in a division are not significant. Therefore CADNA performs a socalled self-validation that consists in controlling all multiplications and divisions during the execution of the code. This self-validation has led to the parallel computation of the N samples R i and also to the concept of computational zero [START_REF] Vignes | Zéro mathématique et zéro informatique[END_REF] defined as either a mathematical zero or a number without any significance, i.e. numerical noise. From this concept, relational operators that take into accout round-off error propagation have been defined [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF].

Numerical validation of sequential codes using CADNA

The CADNA library allows one to use new numerical types: the stochastic types. In practice, classic floating-point variables are replaced by the corresponding stochastic variables, which are composed of three floating-point values and one integer to store the accuracy. The library contains the definition of all arithmetic operations and order relations for the stochastic types. Only the exact significant digits of a stochastic variable are printed or "@.0" if it is numerical noise. Because all operators are redefined for stochastic variables, the use of CADNA in a program requires only a few modifications: essentially changes in the declarations of variables and in input/output statements. CADNA has been successfully used for the numerical validation of academic and industrial simulation codes in various domains such as astrophysics, atomic physics, chemistry, climate science, fluid dynamics, geophysics [START_REF] Jézéquel | Estimation of numerical reproducibility on CPU and GPU[END_REF].

We rely here on the latest CADNA version [START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF], which performs better with sequential codes and also enables the numerical validation of vectorized codes. Unlike in the previous CADNA version, the rounding mode is not explicitly changed in this version. It is set to rounding towards +∞ once, in the cadna init initialisation function that must be called at the beginning of the user program. The random rounding mode is then emulated by taking advantage of arithmetic properties of rounding towards +∞.

CADNA can detect numerical instabilities which occur during the execution of the code. These instabilities are of several types.

-Self-validation: both operands in a multiplication, the divisor in a division, or the first argument of the power function are not significant. -Instability in an intrinsic or a mathematical function: non significant argument in such a function. -Branching instability: indeterminism in a branching test.

-Cancellation: sudden loss of accuracy on an addition or a subtraction.

At the end of the run, the total number of instabilities and each type of instability together with their occurrences are printed. If instabilities have occurred, their source need to be identified and, if necessary, the code changed. The user can specify the instabilities to be detected. One may choose, for instance, to activate only self-validation, to detect all types of instabilities or to deactivate the detection of instabilities.

CADNA for OpenMP codes

We detail here how we have designed an OpenMP-enabled CADNA version for C/C++ codes, although up to our knowledge such a version could also be developed for Fortran codes.

Compatibility check. Our CADNA implementation requires the rounding mode of each thread being set to rounding towards +∞. But to our knowledge, there is no guarantee in the OpenMP standard that the worker threads have the same rounding mode as their master thread (though vendor specific extension may ensure this such as the KMP_INHERIT_FP_CONTROL environment variable for the Intel compiler [START_REF] Intel | User and Reference Guide for the Intel C++ Compiler 15[END_REF]).

We thus have to check the compatibility between each OpenMP implementation used and CADNA. Within the cadna init call (considered as not performed inside an OpenMP parallel region), we hence first set the rounding mode of the master thread to rounding towards +∞, and then check in an OpenMP parallel region that all worker threads have inherited their rounding mode from the master thread. This is performed thanks to a software test: in double precision, only rounding towards +∞ ensures that 1.0 + 1.0e-20 != 1.0 This enables us to check OpenMP implementations that fork new worker threads at each OpenMP parallel region. If the parallel region within the cadna init call is the first one in the program, this test enables us to also check for OpenMP implementations that create worker threads once at the first OpenMP parallel region, and then retrieve worker threads from a thread pool.

However, when partially analyzing a large OpenMP code with CADNA, it is convenient to be able to consider that the cadna init call is performed after some parallel regions. In such case, and considering that all parallel regions (the ones within the cadna init call and the ones in the CADNA instrumented code) have the same number of threads (condition hereafter referred to as (*)), we check thread pool based OpenMP implementations with the following test within the cadna init call: we set to rounding towards +∞ the rounding mode of all threads in a second parallel region, and we check in a third parallel region that the rounding mode of each thread has been correctly saved. If this second check fails, we consider that the OpenMP execution environment is not compatible with CADNA.

In our tests, both GOMP and Intel implementations have found to be compatible with CADNA, whether cadna init is called prior to any OpenMP parallel region or not.

Random rounding mode. The next issue lies in the random number generation. Thanks to the new version of CADNA (see Sect. 2.2) used here, we now rely on a random number generator that was designed to be easily replicated in each execution flow; this was indeed required for a SIMD execution of CADNA where the random generator was replicated for each scalar lane of the SIMD unit. It is thus straightforward to have a distinct random generator in each OpenMP thread: the three required variables are defined as threadprivate and are initialized in a parallel region within the cadna init call.

However, in order to ensure that the values of these threadprivate variables will persist from one parallel region to the next one, we have to check the conditions listed in Sect. 2.14.2 of the OpenMP 4.0 reference [START_REF]OpenMP Application Program Interface[END_REF], which we hereafter denote as (**) and which we summarize here as: same number of threads and same affinity policies for all parallel regions, no nested parallelism and no dynamic adjustment of the number of threads. It can be noticed that the conditions (**) encompass the condition (*).

Instability detection. In an OpenMP application instrumented with CADNA, the counters dedicated to the numerical instabilities (see Sect. 2.2) can be concurrently incremented by multiple threads, which can lead to data race conditions. We have thus protected the updates of these counters with OpenMP atomic constructs.

OpenMP reductions. When an OpenMP application includes an OpenMP reduction on floating-point variables, instrumenting such a reduction with CADNA implies to perform the corresponding reduction operation with stochastic variables. This is possible since OpenMP 4.0 thanks to the declare reduction construct and the redefinition of all arithmetic operators in the CADNA codes. Currently the +, -and * operators are supported.

Atomic constructs on stochastic variables. If an OpenMP application contains atomic constructs, the code cannot be instrumented as is with CADNA. As specified in the OpenMP 4.0 reference [START_REF]OpenMP Application Program Interface[END_REF], an atomic statement should involve scalar types, and cannot thus be applied on the CADNA stochastic types.

A CADNA-redefined arithmetic operation implies three floating-point IEEE operations, some bit manipulations (depending on the type of the arithmetic operation and on the rounding mode), and some instability detections (depending on the user specification). During this sequence of operations, we have to ensure that the variables are accessed by only one thread. The instability detection depends indeed on the result of the three floating-point computations. We have therefore chosen to include such CADNA-overloaded operation in a critical block. In practice, each atomic construct is replaced by a critical construct in the user code, and one can use distinct critical regions (with distinct names) when accessing distinct variables or arrays.

Provided that the conditions (**) are ensured, CADNA can now be used in OpenMP codes. We emphasize that, except for the atomic constructs, the OpenMP-enabled CADNA version does not require additional modifications to the user code other than those required by the sequential CADNA library.

Performance tests and application to OpenMP codes

For our tests, we use two different HPC servers. The first one, named CPUserver, is composed of two Intel E5-2660 CPUs (each having 8 cores with 2-way SMT at 2.20 GHz) and 32 GB of memory. The second one is a 5110P Xeon Phi (Knight Corner), used in native mode as a distant server and named Xeon-Phi, with 60 cores (4-way SMT at 1.053 GHz) and 8 GB of dedicated memory.

In all our tests, the number of threads in each parallel region is determined according to the OMP NUM THREADS environment variable. Moreover, the same affinity policies are set for all parallel region, and no nested parallelism as well as no dynamic adjustment of the number of threads are used on both servers in order to meet the conditions (**) (see Sect. 3).

Unless otherwise mentioned, we use here CADNA with only self-validation activated, since this is its classic usage. It has to be noticed that all runs are scalar ones (no vectorization used): as shown in [START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF] it is possible to use CADNA for SIMD codes but this requires to access the SIMD lane identifier. To our knowledge this is currently not possible with OpenMP directives, and we therefore rely on the SPMD-on-SIMD programming paradigm (e.g. with the Intel SPMD Program Compilerispc6 ). In this paper we only focus on OpenMP programming: we thus consider scalar IEEE (i.e. not instrumented with CADNA) codes and scalar CADNA codes in order to precisely measure the performance impact of CADNA for OpenMP multi-threaded applications.

A reduction code

This first experiment shows that a sequential code and its parallel version may exhibit a very different numerical quality.

The sum S of all elements of an array A of size 2n with n = 10 6 is computed in single precision. For i = 0, ..., 2n -1, each element A[i] is initialised as A[i] = -i if i is even, and A[i] = i otherwise. Therefore the exact sum is 10 6 . The sum S is computed using an OpenMP for loop construct with a reduction clause and two possible scheduling clauses: (static) so as the loop is divided into chunks of size 2n over the number of threads, and (static,1) so as the chunk size is 1. Results obtained for 1, 32, and 240 threads, with and without CADNA, are reported in Table 1. As already mentioned in Sect. 2.2, CADNA prints only the exact significant digits of results and @.0 if they are numerical noise.

Using the (static) scheduling option, whatever the number of threads, the result S is 1.000000E+06 with or without CADNA. Therefore the result has the maximum possible accuracy in single precision, i.e. 7 exact significant digits. Indeed, in this case each thread has in charge contiguous elements of A and alternatively sums positive and negative values: all computed results are accurate.

Using the (static,1) scheduling option, with 32 or 240 threads, the result has no correct digits and is provided by CADNA as @.0. Without CADNA the result depends on the number of threads and one can notice that for 240 threads (static) (static,1) # threads without CADNA with CADNA without CADNA with CADNA 1 1.000000E+06 1.000000E+06 1.000000E+06 1.000000E+06 32 1.000000E+06 1.000000E+06 1.892352E+06 @.0 240 1.000000E+06 1.000000E+06 1.617920E+05 @.0 Table 1. Results S obtained with and without CADNA using two possible schedulings.

it has a wrong order of magnitude. Each thread sums indeed values of the same sign, the absolute value of the sum thus keeps on increasing, and its accuracy regularly decreases due to the limited floating-point precision. The reduction involves inaccurate with close absolute values and different signs and thus provides numerical noise. In this case, CADNA detects cancellations, i.e. losses of accuracy due to the subtraction of close inaccurate results.

One advantage of this simple numerical example is that its exact result can be easily compared to the computed one. It points out that the accuracy of results may change from a sequential to a parallel execution of a code and also among several parallel executions with different schedulings. This numerical example emphasizes the need for a CADNA version available for OpenMP codes.

Performance tests

We now focus on the performance impact of CADNA on OpenMP codes by studying different micro-applications or benchmarks. In order to handle NUMA effects on CPU-server, we use the following OpenMP settings for thread affinity policies:

when the number of threads is lower or equal to 8, OMP PLACES is set to sockets [START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF] and OMP PROC BIND to master in order to have all threads on the same socket as the master thread; otherwise (i.e. for 16 or 2 × 16 threads) no thread affinity policy is used (OMP PROC BIND set to false): we let the OS scheduler run the threads on all physical CPU cores (with 16 threads, i.e. without SMT) or on all logical CPU cores (with 2 × 16 threads, i.e. with 2-way SMT).

On Xeon-Phi, no thread affinity policy is used. By default, we use gcc/g++ version 4.9.2, and Intel icc/icpc version 2017 (beta). Figures 1 and2 first present performance results for a Mandelbrot set computation (4096 × 4096 pixels, with a maximum number of 4096 iterations per pixel) and for a 3D finite difference (FD) stencil (6 iterations over a 256×256×256 grid with a six-order central stencil). First, for serial runs the CADNA over IEEE overheads match the ones previously obtained [START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF]. When considering OpenMP runs with only one thread, one can see an important additional overhead with gcc (see Figs. 1(a) and 2(a)). This seems to be related to the handling of threadprivate variables by gcc. Such a problem does not exist with icc (see Figs. As far as Mandelbrot computations are concerned (see Figs. 1(a) and 1(b)), the CADNA overheads remain stable when the number of physical CPUs used increases. This shows no scaling overhead with CADNA for compute bound applications. One can however notice that, contrary to the IEEE version, the CADNA version shows no performance gain with the SMT capability.

As far as 3D FD stencil computations are concerned (see Figs. 2(a) and 2(b)), the CADNA overheads remain stable (or decrease somewhat) when the number of physical CPUs used increases. This also shows no scaling overhead for CADNA for such an application which is memory-bound. Here again the CADNA version shows no performance gain with the SMT capability, but this is also the case for the IEEE version since this application is memory-bound.

Figure 3 presents results for the + reduction. The decrease in CADNA overheads between the serial execution and the OpenMP with 1 thread is only due to a slower computation for the IEEE code with 1 OpenMP thread compared to the serial IEEE code without OpenMP. When using multiple threads, one can see that the CADNA overheads are lower. This is due to the fact that the arithmetic intensity of the benchmarks is modified by the use of CADNA. Using CADNA, we have indeed at least 3 times more computations (for each sample, and for the extra operations associated with the rounding mode emulation and the random number generation), while requiring 4 memory accesses. But these memory accesses can be performed at once with the same cache line (as the members of the stochastic types are contiguous in memory): the arithmetic intensity thus increases with CADNA. Since this reduction benchmark presents an even lower arithmetic intensity than the 3D FD stencil one, its performance is even more limited by the memory bandwidth. Using CADNA thus leads to performance less limited by the memory bandwidth and to greater speedups when increasing the number of threads and when using SMT. Similar conclusions have been obtained for the * reduction.

In order to study the impact of atomic constructs required by the instability detection, we present in Fig. 4 performance results for the Mandelbrot set computation with detection of no instability and of all instabilities. The increase of CADNA overheads between serial and 1-thread OpenMP executions is thus 4.5% for no instability (hence no atomic construct) in Fig. 4(a) and 6% for self-validation only (some atomic constructs used) in Fig. 1 Fig. 3. Performance results (using icc) for reductions (+ operator on 2 28 elements), with CADNA over IEEE overheads (Ovh.). OMP-XT denotes OpenMP runs with X threads.

construct impact is thus not significant for 1 OpenMP thread. When detecting all instabilities (all atomic constructs used) in Fig. 4(b), this increase is even lower (2.1%) since the CADNA overheads are much more important. Moreover, the atomic construct impact does not increase for multiple OpenMP threads as shown in Figs. Finally, Fig. 5 presents the CADNA over IEEE overheads of our benchmarks for serial and OpenMP executions on Xeon-Phi. The CADNA overheads are greater on this architecture for Mandelbrot and reduction computations due to the use of the strict and no-except floating-point models (required by CADNA with icc): contrary to CPU-server, such models imply an important overhead on the Knight Corner architecture (e.g. 5.1x for Mandelbrot computations). Nevertheless these tests confirm the previous conclusions (obtained on CPUserver), here for a higher number of cores and threads: no scaling overhead for CADNA (and even lower overheads for memory-bound applications with SMT), and no extra overhead due to the atomic constructs used for the instability detection.

The overhead associated with the atomic constructs on stochastic variables will be studied in the next section.

A shallow-water application

Presentation. The CADNA library is finally applied to a "real-life" application: a Shallow-Water model (denoted SW). A SW model is a numerical integration of the SW equations [START_REF] Vreugdenhil | Numerical Methods for Shallow-Water Flow[END_REF] that describe the flow in a fluid. The application presented here describes the evolution of the water height and velocities in a twodimensional oceanic basin (the water velocity is constant on the water depth). Fundamentally, a SW computation can be viewed as a combination of finite difference stencils.

In our application, the model is integrated for 500 time steps on a 256 × 256 grid. There are two integration modes: the direct model (denoted forward ) that computes the output of the model, and the adjoint model (denoted backward ) that computes the sensitivity of the output to the initial conditions. Both modes (forward and backward) involve the same set of computations, the only difference lies in the direction of the time loop: from past to future in the forward mode, from future to past in the backward mode. As a consequence, the forward mode implementation is fully parallel in space using a space domain decomposition, whereas the backward mode requires the use of numerous atomic operations. The successive computations of the forward and backward modes allow the estimation of a scalar quantity called the "residual" which is a key parameter for validating adjoint codes [START_REF] Brajard | Numerical Validation of Data Assimilation Codes Generated by the YAO Software[END_REF]. The code was compiled using icc/icpc version 2015, and run on CPU-server.

Performance and numerical results. We first present the performance impact of CADNA on this SW code. Figure 6 presents performance results of the forward mode. It shows that there is no scaling overhead in CADNA. The overhead is even lower with multiple threads (9.9x for 16 threads) than for one Fig. 5. CADNA over IEEE overheads for various benchmarks on the Xeon Phi, using icc and the same parameters than on CPU-server (except for the reduction here with 2 27 elements due to memory constraints). thread (15.0x) as it was already explained in Sect. 4.2 for such memory-bound application. Figure 7 presents performance results of the backward mode. In this case, parallel speedups are limited by the use of numerous atomic operations on floating point numbers, as already discussed in [START_REF] Nardi | YAO: a generator of parallel code for variational data assimilation applications[END_REF]. With CADNA, these operations must be converted into critical regions as shown in Sect. 3 (in this SW model, we use three distinct critical regions). As a consequence, the computation time and therefore the overhead are increasing with the number of threads, because of increased lock contention. It should thus be noticed that the numerical validation of a parallel code using (numerous) atomic directives remains possible using CADNA but with a higher cost.

Finally we focus on the numerical stability of the residual of the SW model (see Table 2). For clarity only the serial and 16-thread OpenMP results are presented. It can be noticed that there are only 6 common digits between the results computed without CADNA by the serial and parallel versions. CADNA estimates that 5 digits in serial and 4 digits in parallel are correct. These digits (except the last one in serial which may not be reliable in CADNA) are consistent with the IEEE versions. Thanks to CADNA, we can therefore state that the OpenMP version of this SW application is reliable with a 4-digit accuracy. Moreover, by comparing the serial and parallel numerical results, we can deduce that the differences in the IEEE results are likely due to round-off error propagations and not to a bug in the OpenMP parallelization.

Conclusion

In this paper, we have presented a new OpenMP-enabled CADNA version that can be used to validate parallel codes by providing the number of exact significant digits of their results, and that presents similar or lower CADNA overheads with respect to the serial ones. This library can also help to explain differences in the numerical results between serial and (multiple) parallel executions.

It is worth noting that the control of accuracy can also be performed on MPI codes [START_REF] Montan | Numerical verification of industrial numerical codes[END_REF]. Therefore the work presented on this paper could be extended to enable the numerical validation of large scale applications using both MPI and OpenMP. Besides, the CADNA instrumentation of codes could be eased thanks to OpenMP-like compiler directives that would avoid the insertion and modifications of C/C++ instructions in the user program. A preliminary study has been performed on such directives. However work has still to be carried out to propose a complete set of directives that enable all CADNA functionalities and can control the numerical quality of any C/C++ program.

  1(b) and 2(b)).

Fig. 1 .

 1 Fig. 1. Performance results for Mandelbrot computations, along with CADNA over IEEE overheads (Ovh.). OMP-XT denotes OpenMP runs with X threads.

Fig. 2 .

 2 Fig.2. Performance results for FD stencil computations, along with CADNA over IEEE overheads (Ovh.). OMP-XT denotes OpenMP runs with X threads.

  1(b), 2(b), 4(a) and 4(b).

Fig. 4 .

 4 Fig.[START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF]. Performance results for Mandelbrot computations (using icc), along with CADNA over IEEE overheads (Ovh.). OMP-XT denotes OpenMP runs with X threads.

Table 2 .

 2 Performance results (using icc) for Shallow-Water (backward mode), with CADNA over IEEE overheads (Ovh.). OMP-XT denotes OpenMP runs with X threads. Comparison between IEEE and CADNA computations of the Shallow-Water residual for serial and parallel executions.

		1000
							IEEE
							CADNA
	Computation time (s, log scale)	10 100	8.9x Ovh.	Ovh. 6.8x	Ovh. 20.5x	Ovh. 43.9x	Ovh. 57.9x
		1
			Serial	OMP-1T	OMP-8T	OMP-16T	OMP-2x16T
	Fig. 7. IEEE		CADNA
			Serial	3.446611873236805E-06	3.4461E-06
		OpenMP -16 treads	3.446619149194419E-06	3.446E-06

URL address: http://cadna.lip6.fr

http://ispc.github.io

Acknowledgments

The authors thank Pierre-Emmanuel Le Roux (LIP6) for managing the compute servers, and Philippe Thierry (Intel) for providing the Intel compiler (2017-beta). They also thank the reviewers for their helpful comments.