Keywords: Locality-Based Utility Computing, Peer To Peer, Self-*, Sustainability, Efficiency, Future Internet

Instead of the current trend consisting of building larger and larger data centers (DCs) in few strategic locations, the DISCOVERY initiative ‡ proposes to leverage any network point of presences (PoP, i.e., a small or medium-sized network center) available through the Internet. The key idea is to demonstrate a widely distributed Cloud platform that can better match the geographical dispersal of users. This involves radical changes in the way resources are managed, but leveraging computing resources around the endusers will enable to deliver a new generation of highly efficient and sustainable Utility Computing (UC) platforms, thus providing a strong alternative to the actual Cloud model based on mega DCs (i.e. DCs composed of tens of thousands resources). Critical to the emergence of such distributed Cloud platforms is the availability of appropriate operating mechanisms. Although, some of protagonists of Cloud federations would argue that it might be possible to federate a significant number of micro-Clouds hosted on each PoP, we emphasize that federated approaches aim at delivering a brokering service in charge of interacting with several Cloud management systems, each of them being already deployed and operated independently by at least one administrator. In other words, current federated approaches do not target to operate, remotely, a significant amount of UC resources geographically distributed but only to use them. The main objective of DISCOVERY is to design, implement, demonstrate and promote a unified system in charge of turning a complex, extremely large-scale and widely distributed infrastructure into a collection of abstracted computing resources which is efficient, reliable, secure and friendly to operate and use. After presenting the DISCOVERY vision, we explain the different choices we made, in particular the choice of revising the OpenStack solution leveraging P2P mechanisms. We believe that such a strategy is promising considering the architecture complexity of such systems and the velocity of open-source initiatives.

L'initiative DISCOVERY -les infrastructures IaaS massivement distribuées comme solution aux principales limitations des plateformes de Cloud Computing actuelles

Résumé :

La tendance actuelle pour supporter la demande croissante d'informatique utilitaire consiste à construire des centres de données de plus en plus grands, dans un nombre limité de lieux stratégiques. Cette approche permet sans aucun doute de satisfaire la demande actuelle tout en conservant une approche centralisée de la gestion de ces ressources, mais elle reste loin de pouvoir fournir des infrastructures répondant aux contraintes actuelles et futures en termes d'efficacité, de juridiction ou encore de durabilité. L'objectif de l'initiative DISCOVERY4 est de concevoir le LUC OS, un système de gestion distribuée des ressources qui permettra de tirer parti de n'importe quel noeud réseau constituant la dorsale d'Internet afin de fournir une nouvelle génération d'informatique utilitaire, plus apte à prendre en compte la dispersion géographique des utilisateurs et leur demande toujours croissante.

Après avoir rappelé les objectifs de l'initiative DISCOVERY et expliqué pourquoi les approches type fédération ne sont pas adaptées pour opérer une infrastructure d'informatique utilitaire intégrée au réseau, nous présentons les prémisses de notre système. Nous expliquerons notamment pourquoi et comment nous avons choisi de démarrer des travaux visant à revisiter la conception de la solution Openstack. De notre point de vue, choisir d'appuyer nos travaux sur cette solution est une stratégie judicieuse à la vue de la complexité des systèmes de gestion des plateformes IaaS et de la vélocité des solutions open-source.

Mots-clés : Calcul utilitaire basé sur la localité, systèmes pair-à-pair, self-*, durabilité, Internet du futur

Introduction

To satisfy the escalating demand for Cloud Computing (CC) resources while realizing economy of scale, the production of computing resources is concentrated in mega data centers (DCs) of ever-increasing size, where the number of physical resources that one DC can host is limited by the capacity of its energy supply and its cooling system. To meet these critical needs in terms of energy supply and cooling, the current trend is toward building DCs in regions with abundant and affordable electricity supplies or in regions close to the polar circle to leverage free cooling techniques [START_REF] Gary | How dirty is your data ? Greenpeace International Report[END_REF].

However, concentrating Mega-DCs in only few attractive places implies different issues. First, a disaster1 in these areas would be dramatic for IT services the DCs host as the connectivity to CC resources would not be guaranteed. Second, in addition to jurisdiction concerns, hosting computing resources in a few locations leads to useless network overheads to reach each DC. Such overheads can prevent the adoption of the UC paradigm by several kind of applications such as mobile computing or big data ones.

The concept of micro/nano DCs at the edge of the backbone [START_REF] Greenberg | The cost of a cloud: research problems in data center networks[END_REF] is a promising solution to address the aforementioned concerns. However, operating multiple small DCs breaks somehow the idea of mutualization in terms of physical resources and administration simplicity, making this approach questionable. One way to enhance mutualization is to leverage existing network centers, starting from the core nodes of the backbone to the different network access points (a.k.a.. PoPs -Points of Presence) in charge of interconnecting public and private institutions. By hosting micro/nano DCs in PoPs, it becomes possible to mutualize resources that are mandatory to operate network/data centers while delivering widely distributed CC platforms better suited to cope with disasters and to match the geographical dispersal of users and their needs (see Figure 1). A preliminary study has established the fundamentals of such an in-network distributed cloud referred by the consortium as the Locality-Based Utility Computing (LUC) concept [START_REF] Bertier | Beyond the clouds: How should next generation utility computing infrastructures be designed?[END_REF]. However, the question of how operating such an infrastructure is still under investigations. Indeed, at this level of distribution, latency and fault tolerance become primary concerns, and collaboration between servers of differents location must be organized wisely.

In this vision paper, we discuss some key-elements that motivate our choices to design and implement the LUC Operating System, a system in charge of turning a LUC infrastructure into a collection of abstracted computing facilities that are as convenient to administrate and use as available Infrastructure-as-a-Service (IaaS) managers [START_REF]Open Source Cloud Computing[END_REF][START_REF]Open Source Data Center Virtualization[END_REF][START_REF]The Open Source, Open Standards Cloud[END_REF]. We explain, in particular, why federated approaches [START_REF] Buyya | Intercloud: Utility-oriented federation of cloud computing environments for scaling of application services[END_REF] are not satisfactory and why designing a fully distributed system that operates all resources makes sense.

As the capabilities of the LUC OS are similar to existing IaaS managers and because it would be a non-sense technically speaking to develop the LUC OS from scratch we chose to revise the OpenStack solution [START_REF]The Open Source, Open Standards Cloud[END_REF], leveraging P2P mechanisms. Our current efforts focus on the validation of a distributed version of the Nova service on top of Grid'5000 [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF]. Historically, Nova relies on a MySQL centralized database, preventing it to natively scale beyond one site. To reach such a goal, we replaced the MySQL component by the REDIS backend, a distributed key/value store. Such a modification enables us to deploy several Nova controllers on distinct sites giving the illusion that there was only one global infrastructure (each controller manipulating the Nova internal states throughout the REDIS system). This first validation paves the way toward a complete LUC OS leveraging the OpenStack ecosystem.

The remaining of the article is as follows: Section 2 explains our design choices. Section 3 describes the OpenStack and gives first details of our current Proof-of-Concept. A preliminar validation of our prototype focusing on the Nova service is presented in Section 4. Finally Section 5 concludes and discusses future actions.

The LUC OS: Design Discussion

The massively distributed cloud we target is an infrastructure that is composed of up to hundreds of micro DCs, which are themselves composed of up to tens of servers. Thus the system in charge of operating such an infrastructure should be able to manage up to thousands of servers spread geographically. Delivering such a system is a tedious task where wrong design choices could prevent to achieve our goal. In this section we first discuss few conceptual considerations that led us to the LUC OS proposal and second remain the major services that the LUC OS should deliver.

From Centralized to Distributed Management

The first way that comes generally to the mind to pilot and use distinct clouds is to rely on classical models like federated approaches: each micro DC hosts and operates its own Cloud infrastructure and a brokering service is in charge of resources provisioning by picking on each cloud. While such approaches can be acceptable for elementary usages, advanced brokering services are mandatory to meet production environment requirements (monitoring, scheduling, automated provisioning, SLAs enforcements . . . ). In addition to dealing with scalability and single point of failure (SPOF) issues, brokering services should integrate mechanisms similar to those that are already implemented at the level of IaaS managers [START_REF] Buyya | Intercloud: Utility-oriented federation of cloud computing environments for scaling of application services[END_REF][START_REF] Houidi | Cloud service delivery across multiple cloud platforms[END_REF]. Consequently, the development of a brokering solution is as difficult as the one of a IaaS manager but with the complexity of relying only on the least common denominator APIs. While few standards such as OCCI [START_REF] Loutas | Towards a reference architecture for semantically interoperable clouds[END_REF] start to be adopted, they do not allow developers to manipulate low-level capabilities of each system, which is generally mandatory to finely administrate resources. In other words, building mechanisms on top of existing ones like in the case of federated systems prevent from going beyond the provided APIs (or require when possible, intrusive mechanisms that must be adapted to the different systems).

The other way to operate such infrastructure is to design and build a dedicated system, i.e., the LUC Operating System, in charge of operating all the geographically spread micro DCs in a distributed manner. A LUC OS will define and leverage its own software interface, thus extending capacities of traditional Clouds with its API and a set of dedicated tools. This offers a unique opportunity to go beyond classical federations of Clouds by addressing all crosscutting concerns of a software stack as complex as a IaaS manager and by revising in a fully distributed manner, mechanisms that have been traditionally implemented in a centralized one (service nodes).

The following question is now to analyze whether the collaborations between instances of the system, that is the service nodes, should be structured either in hierarchical way or in a P2P (i.e., flat) one. Few hierarchical solutions have been proposed during the last years in industry [START_REF]Cascading OpenStack[END_REF][START_REF]Scaling solutions for OpenStack[END_REF] and academia [START_REF] Farahnakian | Hierarchical vm management architecture for cloud data centers[END_REF][START_REF] Feller | Snooze: A scalable and autonomic virtual machine management framework for private clouds[END_REF]. Although they may look easier than P2P structures, hierarchical approaches require additional maintenance costs and complex operations in case of failure. Moreover, mapping and maintaining a relevant tree architecture on top of a network backbone is not meaningful (static partitioning of resources is usually performed). As a consequence, hierarchical approaches do not look to be satisfactory to operate a massively distributed IaaS infrastructure such as the one we target. On the other side, Peer-to-Peer (P2P) file sharing systems are a good example of software that works well at large scale and in a context where computing resources are geographically spread. While largely unexplored for building operating systems, peer-to-peer/decentralized mechanisms have the potential to natively handle the intrinsic distribution of LUC infrastructures as well as the scalability required to manage them. Hence, we propose to leverage advanced P2P mechanisms like overlay networks and distributed hash tables to design the LUC OS building blocks.

Cloud Capabilities

From the administrators and end-users point of views, the LUC OS should deliver a set of high level mechanisms whose assembly results in an operational IaaS system. Recent studies have showed that state of the art IaaS manager [START_REF] Peng | Comparison of several cloud computing platforms[END_REF] were constructed over the same concepts and that a reference architecture for IaaS manager can could defined [START_REF] Moreno-Vozmediano | Iaas cloud architecture: From virtualized datacenters to federated cloud infrastructures[END_REF]. This architecture covers primary services that are needed for building the LUC OS :

• The virtual machines manager is in charge of managing VMs' cycle Inria of life (configuration, scheduling, deployment, suspend/resume and shut down).

• The Image manager is in charge of VM' template files (a.k.a. VM images).

• The Network manager provides connectivity to the infrastructure: virtual networks for VMs and external access for users.

• The Storage manager provides persistent storage facilities to VMs.

• The Administrative tools provide user interfaces to operate and use the infrastructure.

• Finally, the Information manager monitors data of the infrastructure for the auditing/accounting.

The challenge is thus to propose a distributed version of the aforementioned services by relying on advanced P2P mechanisms. However, designing and developing a complete LUC OS from scratch would be an herculean work, including several non-sense actions aiming at simply providing basic mechanisms available in most IaaS solutions. Instead of reinventing the wheel, we propose to minimize both design and implementation efforts by reusing as much as possible existing piece of codes. With this in mind, we propose to investigate whether the Open-Stack solution [START_REF]The Open Source, Open Standards Cloud[END_REF] can be revised to fulfill the LUC infrastructure requirements. Concretely, we propose to determine which mechanisms can be directly used and which ones must be revisited with P2P algorithms. This strategy enables us to focus the effort on the key issues such as the distributed functioning, fault tolerance mechanisms, and the organization of efficient collaborations between service nodes of the infrastructure according to the constraints/requirements of the LUC OS.

Revising OpenStack

OpenStack is an open-source project that aims at developing a complete cloud management system. Similary to the reference architecture described in the previous Section, it is composed of several services, each one dealing with a particular aspect of a Cloud infrastructure as depicted in Figure 2. OpenStack relies on two kinds of nodes: controller and compute node. The former is in charge of managing and distributing work to the latter that provides computing/storage resources to end-users. In other words, the controllers correspond to the different services introduced in the previous section while the compute nodes host the VMs.

From the software point of view, the OpenStack architecture is based on the "shared nothing" principles: each controller (i.e., each service) is connected to the others via two different way:

• A messaging queue that enables the collaboration between sub-services of a controller. • A SQL database (DB) that stores inner states of a controller.

Finally, the controllers interact with each other through REST APIs or directly by accessing the inner-state that are stored in the differents DBs.

Considering the current structure of OpenStack, the main limitation to make it distributed is related to the SQL databases. Indeed, while OpenStack relies on the RabbitMQ messaging service, which is articulated around a centralized broker, there are few implementations of P2P messaging service such as Ac-tiveMQ [START_REF] Snyder | ActiveMQ in action[END_REF] or ZeroMQ [START_REF] Hintjens | ZeroMQ: Messaging for Many Applications[END_REF] that would be adapted to the LUC requirements. The first way to bypass the MySQL limitation is to deploy each controller DB on each location and to synchronize the different DB instances with a dedicated mechanism [START_REF] Kemme | Database replication: A tale of research across communities[END_REF]. By such a mean, when a controller processes a request and performs some actions on one site, changes in the inner-state are also propagated to all the other locations. From a certain point of view, it gives the illusion that there is only one DB for each service. Although the technique described has been used in different proof-of-concepts, current DB synchronization mechanisms are not scalable enough to cope with a LUC infrastructure deployed on large number of geographical sites.

Another approach is to replace the DBs used in OpenStack by a more suitable storage backend that would provide a better scalability. Distributed Hash Tables (DHTs) and more recently key/value systems built on top of the DHT concept such as Dynamo [START_REF] Decandia | Dynamo: amazon's highly available key-value store[END_REF] have demonstrated their efficiency in terms of scalability and fault tolerance properties. In light of this, we have revisited the Nova controller, i.e., the VM manager of OpenStack, in order to to replace the current MySQL DB system by REDIS [START_REF] Han | Survey on nosql database[END_REF], a key/value store. Technically speaking, we modified the Nova database driver. Indeed, the Nova software architecture has been organised in a way which ensures that each of its sub-services does not directly manipulate the database: they have an indirect access through a service called "nova-conductor" which in turn works with an implementation of the "nova.db.api" programming interface. Developers of Nova provide an implementation of this interface that is using SQLAlchemy to manipulate a relational database. We developed a second implementation of this interface that replaces every call to the SQLAlchemy by a call to a custom key/value store driver. This enables Nova's services to work with REDIS by only changing the database driver, limiting the level of intrusiveness in the original source code. Thanks to this modification, it is possible to instanciate a distributed cloud and operate it through a single instance of OpenStack composed of several Nova controllers deployed on distinct sites. Figure 3 depictes such a deployment. Each controller executes a REDIS instance that is configured to work in a clustering way with other instances. One or several controllers can be deployed on each site according to the expected demand in terms of end-users. Finally, a controller can be deployed either on a dedicated node or be mutalized with a compute one Inria as illustrated for Site 3. We higlight that any controller can provision VMs by orchestrating services on the whole infrastructure and not only on the site where it is deployed. Such a behavior is possible thanks to the AMQP bus and the key/value store that go through all controllers. Finally, it is noteworthy that key/value stores that focus on high-availability and partition tolerance criteria like Cassandra [START_REF] Lakshman | Cassandra: a decentralized structured storage system[END_REF] would be more appropriate than REDIS for a production deployment. We chose REDIS for its usage simplicity. Our prototype is under evaluation. However, preliminary experiments have been performed throughout 4 sites of Grid'5000 including 12 compute nodes and 4 controllers overall. While this infrastructure was rather small in comparison to our target, it aimed at validating the interconnection of several controllers WANwide and the correct behaviour of OpenStack using our noSQL backend. Our prototype suceeded to provision 500 VMs in 300 seconds (each controller creating 125 VMs in parallel). A second experiment validated the provisionning of 2000 VMs in less than 30 min. We are currently performing comparisons between OpenStack using the historical MYSQL backend v.s., using a key/value store backend. Our goal is to validate that manipulating internal states of Openstack through a noSQL deliver performances in the same order of the MySQL ones.

Experimental Validation

The validation of our proof-of-concept has been done via two sets of experiments. The first one aimed at measuring the impact of the use of the REDIS NoSQL solution instead of the MySQL system in a single site deployment (i.e., servers were located on the same geographical site). The second set focused on multisite scenarios by comparing the impact of the latency on our distributed Nova service with the "State of the art" approach of distributed OpenStack.

Experiments have been performed on Grid'5000 [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF], a large-scale and versatile experimental testbed that enables researchers to get an access to a large amount of computing resources with a very fine control of the experimental conditions. We deployed and configured each node involved in the experiment with a customized software stack (Ubuntu 14.04, a modified version of OpenStack "Devstack", and REDIS v3) using Python scripts and the Execo toolbox [START_REF] Imbert | Using the EXECO toolbox to perform automatic and reproducible cloud experiments[END_REF].

Impact of REDIS w.r.t MySQL

Changes made over Nova's source code to support a NoSQL database as REDIS is likely to affect its reactivity. The first reason is that a key/value store does not provide a support of operations like joining, and thus the code we developed to provide such operations, creates a computation overhead. The second reason is related to networking. Unlike a single MySQL node, data is spread over several nodes in a REDIS system. Thus, a request can lead to several network exchanges. Finally, REDIS provides a replication strategy to deal with fault tolerant aspects, leading also to possible overheads. 1 compares average response times used to satisfy API requests made during the creation of 500 VMs on an infrastructure deployed on one cluster (containing 1 controller node and 6 compute nodes), using either REDIS or the MySQL backend. The first line of Table 1 shows the overhead related to the use of the combination of a key/value store and our code instead of MySQL. Our approach leads to an average API response time that is twice higher than MySQL. While our code is probably not mature enough, we need to conduct additional experiments to see the distribution of the completion time. Indeed, we expect that some requests (in particular the ones that should deal with advanced SQL features such as join or transaction) take a significant time to be achieved in the NoSQL world. The second line of Table 1 illustrates that distributing the key/value store on several nodes does not increase the overhead related to networking (in our case, data is distributed over four nodes). The third line represents measurements with data replication enabled. The overhead is 10% compared to a clustered configuration without replication. This small overhead is achievable thanks to the asynchronous strategy used in REDIS.

While the REDIS overhead may look critical at the first sight, it must be mitigated with Table 2 which shows that higher API response time does not necessarily mean lower overall reactivity. Indeed, with REDIS, OpenStack delivers performances that are in the same order of magnitude as with MySQL. Each result of Table 1 corresponds to the average response time for all API functions, which is not necessary correlated to the VM creation time as some API functions have a more significant impact than others on the VM creation time. As explained above, such results motivate us to conduct deeper analysis to understand the time distribution of API requests.

Inria

Multi-site Scenarios

The second experiment we made consisted in evaluating a single OpenStack deployment over several locations. Our goal was to compare the behaviour of a single MySQL OpenStack with the advised Galera solution and our REDIS proposal. While the deployment of a single MySQL node is a non sense in a production infrastructure as we mentioned before, evaluating such a scenario enabled us to to get an indication regarding the maximum performance we can expect. Indeed in such a scenario, the DB is deployed on a single server located in one of the locations, without any synchronization mechanism and consequently no overhead related to communications with remote DB nodes on the contrary to a clustered Redis or an infrastructure composed of several MySQL DBs that are synchronized with the Galera mechanisms.

Regarding the experimental methodology, all executions have been conducted on servers of the same site (Rennes) in order to ensure reproducibility : distinct locations (i.e., clusters) have been emulated by adding latency between group of servers thanks to the TC unix tool. Each cluster was containing 1 controller node, 6 compute nodes, and one DB node when needed. Scenarios including 2, 4, 6, 8 clusters have been evaluated, leading to infrastructures composed of up to 8 controllers and 48 compute nodes overall. The latency between each cluster has been set to 10 ms and then 50 ms. Finally, in order to evaluate the capability of such infrastructures to distribute the workload on several controllers, and to detect concurrency problems inherent in using a non relational database backend, the creation of the 500 VMs has been fairly distributed among the available controllers in parallel. 4 present the time to create the 500 VMs. As expected, increasing the number of clusters leads to a decrease of the completion time. This is explained by the fact that a larger number of clusters means a larger number of controllers and compute nodes to handle the workload.

The results measured for a 10ms latency, show that our approach takes a rather constant time to create 500 VMs, which stabilizes around 220 seconds. While a single MySQL node has better results until 6 clusters, one can see the limitations of a single server with 8 clusters. In such a case, the single MySQL performs peforms 89% slower than our approach, while the advised Galera solution is 891% slower than our approach.

With a 50ms inter-cluster latency, the difference between REDIS and MySQL is accentuated in the 8 clusters configuration, as MySQL is 151% slower than our REDIS approach. Regarding Galera, it is noteworthy that important issues related to concurrent modifications of the databases appear with a 50 ms latency, preventing many of the 500VMs to be created (i.e., several bugs occur leading Nova to consider a lot of VMs as crashed). Such pathological behaviours are due to both the important latency between clusters and the burst mode we used to create the 500 VMs (for information, we succeeded to create 500 VMs but in a sequential manner for 2 and 4 clusters).

Conclusion and Future Work

Distributing the way Cloud are managed is one solution to favor the adoption of the distributed cloud model. In this document, we have presented our view of how such distribution can be achieved. We highlighted that it has however a design cost and it should be developed over mature and efficient solutions. With this objective in mind, we chose to design our system, the LUC Operating System, over OpenStack. This choice presents two advantages: minimizing the development efforts and maximizing the chance of being reused by a large community. As as first step, we modified the Nova SQL backend by a distributed key/value system. Although a more advanced validation of this change is required and the question of which metrics to use remains, this first prototype paves the way toward the distribution of additional OpenStack services.

Among the remaining services, the next candidate is the image service Glance. Indeed, as its images are already stored in fully distributed cloud storage software (SWIFT or CEPH), the next step to reach a fully distributed functionning with Glance is to apply the same strategy that we did with Nova. On the other hand, the situation may be different with some other services: Neutron works with drivers that may not be intented to work in a distributed way. In such situation alternatives will have to be found.

Finally, having a wan-wide infrastructure can be source of networking overheads: some objects manipulated by OpenStack are subject to be manipulated by any service of the deployed controllers, and by extension should be visible to any of the controllers. On the other hand, some objects may benefit from a restrained visibility: if a user has build an OpenStack project (tenant) that is based on few sites, appart from data-replication there is no need for storing objects related to this project on external sites. Restraining the storage of such objects according to visibility rules would enable to save network bandwidth and to settle policies for applications such as privacy and efficient data-replication.

We believe, however, that adressing all these challenges are key elements to promote a new generation of cloud computing more sustainable and efficient. Indeed, by revising OpenStack in order to make it natively cooperative, it would enable Internet Service Providers and other institutions in charge of operating a network backbone to build an extreme-scale LUC infrastructure with a limited additional cost. The interest of important actors such as Orange Labs that has officially announced its support to the initiative is an excellent sign of the importance of our action.

Inria

Figure 1 :

 1 Figure 1: The European G ÉANT backbone G ÉANT is the federation of all European Research and Educational Networks. Each black square corresponds to one network point of presence (a.k.a.. a PoP) that can host a nano/micro DC.

Figure 2 :

 2 Figure 2: Services composing OpenStack.

Figure 3 :

 3 Figure 3: Nova controllers are connected through a shared key/value backend and the AMQP bus.

Table 1 :

 1 Average response time to API requests for a mono-site deployment (in ms).

	Backend configuration REDIS MySQL
	1 node	83 37
	4 nodes	82 -
	4 nodes + repl	91 -

Table 2 :

 2 Time used to create 500 VMs on a single cluster configuration (in

	sec.)
	Backend configuration REDIS MySQL
	1 node	322 298
	4 nodes	327 -
	4 nodes + repl	413 -

Table

Table 3 :

 3 Time used to create 500 VMs with a 10ms inter-site latency (in sec.).

	Nb of locations REDIS MySQL Galera
	2 clusters	271	209	2199
	4 clusters	263	139	2011
	6 clusters	229	123	1811
	8 clusters	223	422	1988

Table 4 :

 4 Time used to create 500 VMs with a 50ms inter-site latency (in sec.).

	Nb of locations REDIS MySQL Galera
	2 clusters	723	268 *
	4 clusters	427	203 *
	6 clusters	341	184 -
	8 clusters	302	759 -

Table 3 and

 3 Table

http://beyondtheclouds.github.io

On March

2014, a large crack has been found in the Wanapum Dame leading to emmergency procedures. This hydrolic plan supports the utility power supply to major data centers in central Washington.Inria

RR n°8779

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

Since July 2015, the Discovery initiative is mainly supported through the Inria Project Labs program and the I/O labs, a joint lab between Inria and Orange Labs. Further information at http://www.inria.fr/en/research/researchfields/inria-project-labs


Jianxiao Yang
email: jianxiao.yang@ensta-paristech.fr

Benoit Geller
email: benoit.geller@ensta-paristech.fr

Abdel Charbel
email: charbel.abdelnour@telecom-bretagne.eu

Catherine Nour

Douillard

Meng Li
email: meng.li@imec.be

Min Li
email: min.li@nxp.com

Charbel Abdel Nour

Catherine Douillard
email: catherine.douillard@telecom-bretagne.eu

Max-Log Demapper Architecture Design for DVB-T2 Rotated QAM Constellations

Keywords: Constellations, Log-Likelihood Ratio (LLR), Max-Log Demapper

Rotated and cyclic-Q delayed (RCQD) quadrature amplitude modulation (QAM) improve DVB-T2 system performance over highly time-frequency selective channels. However, when compared with conventional QAM demapper, the RCQD demapper requires a higher computational complexity. In this paper, a complexity-reduced max-log demapper is derived and implemented over a FPGA platform. The proposed demapper allows to find the maximum likelihood (ML) point with a search spanning only M signal constellation points and guarantees to obtain the same log-likelihood ratio (LLR) metrics as the optimum max-log soft decision demapper while spanning at most 2 M signal constellation points. The optimized hardware implementation introduces only a slight performance loss compared to the floating-point full complexity max-log performance.

Index Terms -DVB-T2, Rotated and Cyclic Q Delayed (RCQD)

as 64-QAM or 256-QAM, the computational complexity of a 2D-DEM has a non-negligible impact on receiver design. There have been many studies tackling this complexity problem. The decorrelation based method such as zero forcing (ZF) or minimum mean square error (MMSE) demapper [START_REF] Kim | One-Dimensional Soft-Demapping using Decorrelation with Interference Cancellation for Rotated QAM Constellations[END_REF] can achieve a low complexity but introduces a huge performance loss especially over channels with severe conditions such as deep fades or erasures. Various simplified 2D demappers [START_REF] Li | Design of rotated QAM mapper/demapper for the DVB-T2 standard[END_REF]- [START_REF] Bae | Low complexity two-stage soft demapper for rotated constellation in DVB-T2[END_REF] were also proposed to balance a trade-off between complexity and performance. A simplified method was first presented in [START_REF] Li | Design of rotated QAM mapper/demapper for the DVB-T2 standard[END_REF]- [START_REF] Li | Efficient iterative receiver for bit-Interleaved Coded Modulation according to the DVB-T2 standard[END_REF] based on the decomposition of signal space into sub-regions, each approximately corresponding to one quadrant of the QAM constellation; the demapping operations performed over one sub-region reduces by 61% and by 69% the computational complexity of a 64-QAM and of a 256-QAM respectively. A similar approach was presented in [12] by adaptively adjusting the sub-region with the signal-to-noise ratio (SNR). A so-called sub-region per dimension demapping (PER-DEM) method was provided in [START_REF] Butussi | Low Complexity Demapping of Rotated and Cyclic Q Delayed Constellations for DVB-T2[END_REF] with an exploration space down to 2 M for a RCQD M-QAM. This method starts by dividing the constellation space into regions limited by parallel lines along the imaginary axis; then it computes the distance to the received observation only with respect to points within the two regions closest to the received observation. Another method with sub-region decomposition was proposed in [START_REF] Bae | Low complexity two-stage soft demapper for rotated constellation in DVB-T2[END_REF] based on the observation of 2D LLR contours as a function of I/Q Rayleigh fading channel attenuations. Although this method explores only M points, no decoded performance was given so that the associated performance penalty was not reported. For all these 2D sub-region based methods, there is always a possibility of missing the closest constellation point to the channel observation due to the unbalanced (or different) I and Q channel attenuations. Therefore, all these methods are clearly not optimum. This paper proposes an exact max-log soft demapper with reduced complexity thanks to the proposed simplified detection algorithm and its corresponding hardware. The proposed exact max-log method is able to find the closest point to the received observation i.e., the optimum hard decision by exploring, at most, a space of M points. It also guarantees to perform exact max-log LLR computations i.e., soft decision for all the bits of a symbol by exploring at most a space of 21 M  points. Finally, this paper also compares objectively the C-floating point algorithm with a VHDL implementation which validates the proposed approach.

The remainder of the paper is organized as follows; system model and the conventional soft demapping process are introduced in section II. The derivation and the architecture design of the proposed complexity-reduced exact max-log demapper are detailed in section III. Logic synthesis results and the evaluated performance are given in section IV. Finally, Section V concludes the paper.

II. SYSTEM MODEL FOR RCQD CONSTELLATIONS

A. Rotated and Cyclic-Q Delayed Constellations

A conventional square M -QAM constellation can be regarded as a signal with two independent M -PAM components and takes values from the following c S set:

  cc , I Q I Q s s js s s     SA , (1)
where c A is defined as:

  c 21 1, 2 x x x x c s s s M p p                AI , (2)
where

or x I Q  ,   0,1, , 1 c M  I is an integer set and
I Q I I I Q Q Q I I Q I Q I Q Q y jy h z n j h z n h s s n j h s s n                  (3)
where IQ n n jn  represents a zero-mean circularly symmetric complex Gaussian noise term with variance 2 n  . Thus, equation (3) can be rewritten as:

  1 2 1 2 cos sin , sin cos IQ I I I I Q Q Q Q ss h h s n h h s n                           y hs n h h s n h h n (4)
where

I Q y y     y ,   12  h h h , 1,1 1 2,1 cos sin I Q h h h h          h , 1,2 2 2,2 sin cos I Q h h h h           h , I Q s s     s and I Q n n     n .

B. Soft-Demapping Process

The RCQD constellation breaks the independency between the I and Q components of the signals in the signal space plane. Indeed, both I and Q components contribute to the estimation of the log likelihood   A soft demapping solution with a negligible loss [7] can be obtained by applying the max-log approximation over (5):

                cc 22 2 01 min min 1 LLR max max 0 1 , ii i bb n ii b bb                     s S s S
y hs y hs (6) where

      c 2 min 2 1 min i i bb n bb       sS y hs
represents the log-likelihood metric of i bb  . Although (6) simplifies the computational complexity, it still implies to explore the M complex-valued constellation points. However, it should be noted from ( 6) that the LLR computation becomes an evaluation of the Euclidean distance between the observation vector y and the two closest points with bits i b taking values 0 and 1. Then if there exists such an algorithm that these two closest points can be found for each bit i b without exploring the whole constellation plane, the max-log demapper can be achieved with reduced complexity.

III. COMPLEXITY-REDUCED MAX-LOG DEMAPPER DESIGN

A. The Derivation of the Proposed Algorithm

For a given constellation component

c Q s  A (resp. c I s  A ) in (4), the equivalent observation I r (resp. Q r ) of the other component I s (resp. Q s ) becomes: 21 12
, .

I Q I Q I Q ss ss            r y h h n r y h h n (7)
Then the Euclidean distance between the equivalent observation and the component 2 1

II s

 rh can be expressed as:

         22 2 1 ,1 1,1 ,2 2,1 2 2 1,1 ,2 2,1 ,1 22 1,1 2,1 22
1,1 2,1

I I I I I I II I s r h s r h s h r h r h h s v hh            rh (8)
where In order to further reduce the computational complexity, the following transformations over (4) are proposed.

1,1 ,1 2,1 ,2 22 1,1 2,1 II h r h r v hh    , so
First, both sides of (4) are divided by the factor n  and are extended as below: Second, if the last constant term in ( 9) is moved to the left-hand side, an equivalent observation model y can be derived as: [START_REF] Li | Design and FPGA prototyping of a bit-interleaved coded modulation receiver for the DV B-T2 standard[END_REF] where the equivalent observation terms I y , Q y , the equivalent channel response terms

    cos sin
I I I I Q Q Q Q n n n I I I I s Q Q Q nn Q s II QQ ns h h s n h h s n Mp h h n h h n Mp h h p hh                                                                     y       1 cos sin 1 1, sin cos II QQ n I Q ns n pn h M h                         (
              1,1 1,2 2,
I I Q Q n n s I Q n n s I I I I Q Q Q Q n s n I Q y h M y h h M h h h p n h h p n p hh p hh                                                          yy hs n , I Q n n              hp n
Q Q Q ns y y M h         (12) 1,1 1 2cos , I ns hh    (13) 2,1 1 2sin , Q ns hh    (14)
,

I Q n n    n (20)
,.

I Q p p    p (21)
Finally, based on the results of ( 8), the local optimum solution ,min I p (resp. Q p ) with the given value Qc p  I (resp.

Ic p  I ) can be derived as well:

      1,1 1,2 2,1 2,2 ,min 22 1,1 2,1 1,1 2,1 1,1 1,2 2,1 2,2 2 2 2 2 1,1 2,1 1,1 2,1
arg min round round ,

Ic I Q Q Q II p IQ Q I Q Q h

y h p h y h p pp hh h y h y h h h h p h

h h h r c p                I (22)       1,2 1,1 2,2 2,1 ,min 22 1,2 2,2 1,2 2,2 1,2 1,1 2,2 2,1 2 2 2 2 1,2 2,2 1,2 2,2
arg min round round ,

Qc I I Q I QQ p IQ I Q I I h

y h p h y h p pp hh h y h y h h h h p h

h h h r c p                I ( 23
)
where I r , Q c , Q r and I c are:

1,1 2,1 22 1,1 2,1
,

IQ I h y h y r hh    (24) 1,1 1,2 2,1 2,2 22 1,1 2,1 , Q h h h h c hh    (25) 1,2 2,2 22 1,2 2,2 , IQ Q h y h y r hh    (26) 1,2 1,1 2,2 2,1 22 1,2 2,2 . I h h h h c hh    (27)
Therefore, the corresponding local minimum distance with the corresponding value of Qc p  I (or Ic p  I ) can be computed as:

      22 ,min 1,1 ,min 1,2 2,1 ,min 2,2
,,

I Q I I Q Q I Q p p y h p h p y h p h p        (28)       22 ,min 1,1 1,2 ,min 2,1 2,2 ,min
,.

I Q I I Q Q I Q p p y h p h p y h p h p       
(29) Thanks to these previous derivations, the optimization of the proposed complexity-reduced exact max-log demapper is achieved at the algorithm level; we summarize it as follows:

The complexity-reduced max-log (CRML) demapping algorithm: 0.

Initialize the metric terms h for the equivalent observation model in [START_REF] Li | Design and FPGA prototyping of a bit-interleaved coded modulation receiver for the DV B-T2 standard[END_REF] by using ( 11)-( 16); b. Compute the terms I r , Q c , Q r and I c for the equivalent Euclidean distance terms of ( 22) and ( 23) by using ( 24)-( 27); 2.

  min i bb     for 2 0,1,
For each 0,1, , 1

 

,min ,

IQ pp : { If       min ,min ,min
,,

i i I Q I Q b b p p p p   , then       min ,min ,min
,,

i i I Q I Q b b p p p p   }; 4.
For each bit i b with 2 0,1, ,log 1 iM  , perform: As for the soft bit information i.e., the LLR value, the proposed algorithm allows that these values are obtained, when step 4 is finished after exploring at most (step 2 and step 3 may have overlapped points) 2 M integer points. In this way, (6) can be performed without exploring the whole signal space.

      min min LLR 1 0 i i i b b b      (30)

B. Architecture Design for Hardware Implementation

The CRML demapper presented in the previous section has not yet been optimized for hardware implementation and needs some further modifications to minimize the amount of required hardware complexity, and to maximize the active time of each hardware component, while respecting the constrains in terms of system throughput.

Since the largest constellation supported by DVB-T2 is 256-QAM with 8 bits per symbol decoded by a LDPC decoder, it is reasonable to assume there are 8 clock cycles between two consecutive QAM symbols for the demapper. Among the different elementary steps, the computations in ( 11)-( 16), ( 24)-( 27) are highly dependent and have a low parallelism degree of 2 processes corresponding to the two components I and Q. Differently from step 1, steps 2 and 3 perform computations for different candidate constellation points and have a high parallelism of M . Finally, step 4 has a 2 log M level of parallelism since it deals with pure bit-level computations. Therefore the proposed CRML demapper is divided into 3 parts: preprocessing (step 1), computing Euclidean distance metrics (step 2 and step 3), and bit LLRs generation (step 4).

In the preprocessing step, the inversions in ( 24)-( 27) can be implemented by using Newton's method. For any positive number c , its inverse 1 c  can be calculated iteratively as follows [START_REF] Karp | High Precision Division and Square Root[END_REF]:

        11 2 n n n x x c x     , 1, , nN  , ( 31
)
where n is the iteration index and  . Moreover, in order to reduce the applied number of iterations, the initial value   0 x is set as [START_REF] Wagner | Newton's Method to Perform Division[END_REF]:

  0 00 1.8823 2.8235 xc     . ( 32
)
In this way, the approximate reciprocal   0 n

x of 0 c can be obtained by substituting 0 c and   0 0

x into (31), i.e.,

        11 0 0 0 0 2 n n n x x c x     , 1, , nN  , ( 33
)
where   0 N x is sufficiently accurate after 2 or 3 iterations.

Finally, the approximate reciprocal   N x can be achieved by compensating the factor 2 S :

      00 00 2 , if 2 , 2 , if 2 . N SS N N SS x c c x x c c          (34)
During the second step dedicated to compute Euclidean distance metrics for all the candidates, the operation 22 IQ dd  in (28) and (29) needs 2 multiplications and 1 addition. A possible low-complexity approximation [START_REF] Crawford | Approximation for RMS[END_REF] can be provided by:

    2 22 1 1 1 max , min , . 2 16 32 I Q I Q I Q d d d d d d              (35)
In this way, the 2 multiplications and 1 addition are replaced by 1 multiplication and 3 additions. Since the parallelism in step 2 and step 3 is M , this approximation can save at least M multipliers.

Therefore, the step for computing Euclidean distance metrics for all the candidates includes the following sub-steps where each sub-step represents 1 system clock cycle:

i.

The candidate pairs    

,min , | 0,1, , 1

I Q I p p p M  and     ,min , | 0,1, , 1
I Q Q p p p M
 are computed by using [START_REF] Yang | A Low-Complexity 2D Signal Space Diversity Solution for Future Broadcasting Systems[END_REF] and [START_REF] Brossier | Self adaptative PLL for general QAM constellations[END_REF]. The weighted channel attenuation coefficients   , | , 1,2 and 0,1, ,

1 ij h p i j p M    
can be obtained by only shifting and adding operations in first clock period; ii.

According to (35), the terms

C. Computational Complexity Analysis of the Proposed Demapper

In order to simplify the complexity comparison and to get a direct knowledge of the proposed CRML algorithm, a computational complexity analysis is performed in terms of real multiplications (RMs) and real additions (RAs) which are in general the most important metrics for the complexity evaluation. It should be mentioned that a real subtraction is considered as a RA in this paper. RAs for one RCQD constellation symbol.

IV. LOGIC SYNTHESIS RESULTS AND PERFORMANCE EVALUATION

A. Logic Synthesis Results

The general QAM demapper for rotated constellation was synthesized and implemented by using Xilinx ISE. Computational resources of the demapper take up 7637 slice Flip-Flops and 32764 slice LUTs. The corresponding occupation rates are about 3% and (resp. 15%) of a Xilinx FPGA 5VLX330FF1760-2 for slice registers (resp. slice LUTs). In addition, multiplication resources for the demapper module take up 16 DSP blocks. It represents 8% of the total DSP blocks available in the chosen device. The maximum system clock frequency reaches 96 MHz (the clock period is 10.4 ns) and there are 8 system clocks between two consecutive symbols. So an output LLR rate of 96 MLLR/s for 256-QAM can be achieved at the input I/Q symbol rate of 12 Msymbol/s. Most importantly, only 8 multipliers and 213 adders/subtractors are used.

B. Performance Evaluation and Numerical Computational Complexity Comparisons

The proposed demapping algorithm is compared with other methods in terms of BER and computational complexity, such as the max-log method (see (6)), the MMSE method [START_REF] Kim | One-Dimensional Soft-Demapping using Decorrelation with Interference Cancellation for Rotated QAM Constellations[END_REF], the sub-region method [START_REF] Li | Design of rotated QAM mapper/demapper for the DVB-T2 standard[END_REF], and the PD-DEM method [START_REF] Butussi | Low Complexity Demapping of Rotated and Cyclic Q Delayed Constellations for DVB-T2[END_REF]. In these simulations, perfect synchronization and channel estimation are assumed, which is different from the practical case [START_REF] Bay | On the Hybrid Cramé r-Rao bound and its Application to Dynamical Phase Estimation[END_REF]- [START_REF] Brossier | Self adaptative PLL for general QAM constellations[END_REF]. It should be mentioned that the low RM number has to be attributed to the equivalent observation model in [START_REF] Li | Design and FPGA prototyping of a bit-interleaved coded modulation receiver for the DV B-T2 standard[END_REF] which greatly simplifies the computation of Euclidean distance metrics. Moreover, although the MMSE demapper seems also attractive over fading channel, it has a very high error floor over fading erasure channel [START_REF] Yang | A Low-Complexity 2D Signal Space Diversity Solution for Future Broadcasting Systems[END_REF], [START_REF] Yang | Uniformly Projected RCQD QAM: A Low-Complexity Signal Space Diversity Solution for Broadcasting Systems and Even More[END_REF]. . Moreover, a hardware architecture design is detailed with additional simplifications over the original computations. Since the non-rotated constellations are special cases of the rotated ones with rotation angle 0   , this proposed demapper can also be applied to the demapping of conventional QAM constellations. Some future work concerns the interface of the proposed scheme with other iterative channel coding [START_REF] Vanstraceele | Block Turbo Codes for Multicarrier Local Loop Transmission[END_REF]- [START_REF] Geller | Block Turbo Codes: From Architecture to Application[END_REF].

s.I h  and 0 Qh

 0 is a QAM normalization factor (e.g., for 256-QAM 170 s  ). In order to obtain a RCQD constellation, the conventional square symbol is first rotated by an angle  to obtain a rotated symbol Then the imaginary part of z is delayed by one symbol period to build the symbol x such that the real and imaginary parts of symbol z are transmitted over two different independent and identically distributed fading events. Equivalently, for orthogonal frequency-division multiplexing (OFDM)-based systems such as DVB-T2, these parts are transmitted over two sufficiently spaced subcarriers of an OFDM symbol, therefore doubling the diversity order. Let 0  denote Rayleigh distributed fading coefficients that affect the two subcarriers where symbol z has been transmitted. The observed symbol IQ y y jy  received by the demapper can be expressed as:

1 iM

 1  . Therefore the accurate LLR computation requires an exploration of a signal space containing all the possible M complex-valued constellation points:

1 h

 1 observation I r (resp. Q r ) can be computed without being compared with other Euclidean distance terms.Therefore, the essential idea of this 2-D global minimization problem can be simplified as follows: for each given equivalent observation solution I r (resp. Q r ) in (7) with the known equivalent channel response (resp. 2 h ) by rounding v . The global optimum solution must then belong to one of the M local optimum solutions many parallel dividers, multipliers and comparators.

  9) where I p (resp. Q p ) is defined in (2) and represents the equivalent constellation component of I s (resp. Q s )

  end of step 2 after exploring M integer (constellation) points in signal space.

  needs only two multiplications and one subtraction without any division operation.The accuracy of the inverse value is highly dependent on the iteration number, while the iteration number relies on the initial value   0 x . In order to increase computation stability, the fixed point value c is normalized by continuously left-shifting S ( 2 S  ) or right-shifting S ( 2 S  ) so that the normalized value

  the step of bit LLR generation, since the maximum number of bits per symbol is 8 (256-QAM), only one LLR   LLR i b is generated within each system clock cycle for each bit i b by using (30).

Fig. 1 .

 1 Fig. 1. BER evaluations of floating-point C simulations for DVB-T2 RCQD 256-QAM.

Fig. 1

 1 Fig. 1 further compares the floating-point bit error rate (BER) performance of the algorithms in section IV.B. All these methods are evaluated over a Rayleigh fading channel [7] for the RCQD 256-QAM with the 64800-bit long size and 4/5 rate low density parity check (LDPC) channel code defined in the DVB-T2 standard. Without any surprise, the CRML algorithm achieves exactly the same best floating-point performance as the full-complexity max-log

Fig. 2 .

 2 Fig. 2. BER evaluations of the proposed demappers in floating-point C, fixed-point C, and fixed-point hardware simulations for DVB-T2 RCQD 256-QAM.The C and VHDL fixed-point BER performance of the proposed CRML algorithm is also presented on Fig.2. In this figure, the C floating-point BER performance of the proposed CRML algorithm serves as reference. It can be observed that there is only a 0.4 dB loss difference between the floating-point version and the VHDL Implemented version including the loss introduced by the hardware of the parallelized LDPC decoder.

TABLE I Complexity

 I Comparison of the considered algorithms for the DVB-T2 RCQD 256-QAM constellation.

	Algorithm	CP	RM	RA
	Max-Log	256	1032	776
	Sub-region	81	332	251
	MMSE	16	64	48
	PD-DEM	80	390	279
	Proposed CRML	32	60	371

TABLE I

 I gives the complexity comparison in terms of candidate point explored (CP), RM and RA to demap the RCQD 256-QAM signal.

  16.2 16.3 16.4 16.5 16.6 16.7 16.8 16.

		10 -2
	BER	10 -4
			Proposed C Floating
			Proposed C Fixed
		10 -6	Proposed VHDL Fixed
			9	17	17.1 17.2
			SNR (dB)

ANR Greencocom project.


Keywords: BLAS library, reproducibility, accuracy, superaccumulator, floating-point expansions, error-free transformations, multi-and many-core architectures

Due to non-associativity of floating-point operations and dynamic scheduling on parallel architectures, getting a bitwise reproducible floating-point result for multiple executions of the same code on different or even similar parallel architectures is challenging. We address the problem of reproducibility in the context of fundamental linear algebra operations -like the ones included in the BLAS library -and propose algorithms that yield both reproducible and accurate results (correct rounding, except for triangular solver). We present implementations of these algorithms for the BLAS routines along with the performance results in parallel environments such as Intel desktop and server CPUs, Intel Xeon Phi, and both NVIDIA and AMD GPUs.

I. INTRODUCTION

Exascale computing (10 18 operations per second) is likely to be reached within a decade. This increased computational power of modern supercomputers enables us to solve more complex and larger problems. That, consequently, leads to the higher number of floating-point (FP) operations to be performed, each of them potentially causing a round-off error. As FP operations like addition and multiplication are nonassociative, so that, for instance, parallel implementations of the BLAS routines become non-reproducible. We define reproducibility as an ability to obtain a bit-wise identical FP result from multiple runs of the code on the same input data. However, the result often varies from one parallel machine to another or even from one run to another. These discrepancies worsen on heterogeneous architectures such as clusters with co-processors or GPUs. Such non-determinism of FP calculations in parallel programs causes validation and debugging issues and may even lead to deadlocks. The problem of both non-determinism and non-reproducibility of FP computations becomes so important that it is listed among the top ten ExaScale challenges [START_REF] Lucas | Top Ten ExaScale Research Challenges[END_REF].

Numerical reproducibility can be addressed by targeting either the order of operations or the error resulting from the finite arithmetic. One solution consists in providing the deterministic control over rounding errors by, for example, enforcing the execution order for each operation. However, this approach is not portable and/or does not scale well with the number of processing cores. The other solution aims at avoiding cancellation and rounding errors by using an exact accumulation method such as the superaccumulator proposed by Kulisch [8]. This solution increases the accuracy at the price of more operations and memory transfers. Because of that, for a long time, it was considered too costly.

To enhance reproducibility, Intel proposed a "Conditional Numerical Reproducibility" (CNR) in its MKL. However, CNR does not ensure correct rounding and it induces large performance overhead. For instance, for large arrays the MKL's summation with CNR is 85-93 % slower than both the regular MKL's and our reproducible summation. Demmel and Nguyen introduced a family of algorithms for reproducible summation in FP arithmetic [START_REF] Demmel | Fast reproducible floating-point summation[END_REF]. Their solution induces a twofold slowdown as data transfers and reductions need to be performed twice. They improved the algorithm [START_REF] Demmel | Parallel Reproducible Summation[END_REF] using a single reduction among nodes, which minimized the overhead to roughly 20 %. Demmel and Nguyen have extended their concept to the other BLAS-1 routines, distributed as ReproBLAS. Arteaga et al. [START_REF] Arteaga | Designing bitreproducible portable high-performance applications[END_REF] used this approach with improved communication for summation and obtained the same accuracy with roughly 10 % overhead. Neal [START_REF] Radford | Fast exact summation using small and large superaccumulators[END_REF] considered scalar superaccumulators of different sizes for summation in the R software. Alternatively, we introduced in [START_REF] Sylvain Collange | Numerical Reproducibility for the Parallel Reduction on Multi-and Many-Core Architectures[END_REF] a multi-level approach to compute reproducible sums. This approach is based on FP expansions (FPEs) and superaccumulators. The proposed implementations showed that the numerical reproducibility and bit-perfect accuracy can be achieved at no additional cost for large sums with dynamic ranges of up to 90 orders of magnitude.

In this work, we address the problem of reproducibility of the BLAS routines due to the cancellations and rounding errors that occur during FP operations. We begin with parallel reduction using our multi-level summation approach, which allows to perform accumulation in any order without loosing a bit of information. Then, we extend this approach to some of the BLAS routines from level-1 to level-3. The paper is organized as follows. Section II reviews aspects of computer arithmetic. Section III presents our multi-level reproducible approach. We evaluate our implementations in Section IV and discuss conclusions in Section V.

II. COMPUTER ARITHMETIC

FP arithmetic consists in approximating real numbers with a mantissa, an exponent, and a sign:

x = ± x 0 .x 1 . . . x M -1 mantissa ×b e , 0 ≤ x i ≤ b -1, x 0 = 0
where b is the basis, M is the precision, and e is the exponent (e min ≤ e ≤ e max ). The IEEE-754 standard specifies FP formats and operations as well as requires correctly rounded results for the basic FP operations (+, -, ×, /, √ ) and fma 1 .

This means the operations are performed as if the result was computed with the infinite precision and rounded to the FP format. In this paper, we extend the correct rounding criterion to the BLAS routines. We consider double precision format (binary64) and assume rounding-to-nearest.

Two approaches exist to perform one FP addition without introducing rounding error. The first solution aims at computing the error which occurred during rounding using FP expansions in conjunction with error-free transformations (EFTs). FPEs represent the result as an unevaluated sum of FP numbers, whose components are ordered by magnitude with minimal overlap to cover a wide range of exponents. FPEs of sizes 2 are described in [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF]. When working with the rounding-to-nearest mode, the rounding error of addition and multiplication can be represented as a FP number. This error can be computed in FP arithmetic using EFTs: TwoSum, Alg. 1, for addition; TwoProd, Alg. 2, for multiplication.

Algorithm 1: EFT for the sum of two FP numbers.

Function [r, s] = TwoSum(a, b) r ← a + b z ← r -a s ← (a -(r -z)) + (b -z)
Algorithm 2: EFT for the product of two FP numbers.

Function [r, s] = TwoProd(a, b) r ← a × b s ← f ma(a, b, -r)
Adding one FP number to an expansion is an iterative operation. The FP number is first added to the head of the expansion. The rounding error is recovered as a FP number, using an EFT such as TwoSum, and recursively accumulated to the remainder of the expansion. FPE computations are free of conditional branches and memory accesses, making fast pipelined and vectorized implementations practical. However, their complexity is linear with the number of terms in the FPE, so they are best suited for low dynamic range sums.

The second solution exploits the finite range of representable FP numbers by storing every bit of the sum. It can be implemented as a large fixed-point accumulator (superaccumulator). This accumulator covers the range from the minimum representable FP value to the maximum value independently of the sign. We use a superaccumulator of 2098 (emin + emax + mantissa = 1022 + 1023 + 53) bits for double precision FP accumulation. Superaccumulator implementations involve irregular memory accesses, but the amortized complexity of accumulation does not depend on the accumulator length, making it suitable for large dynamic range sums.

III. MULTI-LEVEL REPRODUCIBLE APPROACH

We split the reproducible parallel reduction algorithm into five stages [START_REF] Sylvain Collange | Numerical Reproducibility for the Parallel Reduction on Multi-and Many-Core Architectures[END_REF]:filtering, private superaccumulation, local super-

1 fma(a, b, c) computes a × b + c with a single rounding.
accumulation, parallel reduction, and rounding. This decomposition is suitable for nested parallelism on modern architectures. The idea is to accumulate numbers with commonlyoccurred exponents using a fast small FPE cache and to use slower superaccumulators only for numbers of unlikely exponents.

The first stage uses FPEs with EFTs for addition of two FP numbers, see Alg. 1. Each vector lane or GPU thread maintains its own FPE of size p. In case the accuracy provided by FPEs is not enough, the remaining rounding error is accumulated to private superaccumulators on the second stage. Private superaccumulators belong to a CPU thread or a set of GPU threads. Each private superaccumulator also receives the contents of the FPEs at the end of the summation. On the third stage, private superaccumulators are merged into a single local superaccumulator (one per a node or a GPU warp). The fourth stage performs a parallel reduction of local superaccumulators to a final superaccumulator. Finally, on the fifth stage, the result is rounded to the desired FP format.

We extended and adapted this approach to dot product, triangular solver [START_REF] Iakymchuk | Reproducible triangular solvers for high-performance computing[END_REF], and matrix-matrix multiplication [START_REF] Iakymchuk | Reproducible and Accurate Matrix Multiplication for GPU Accelerators[END_REF].

IV. PERFORMANCE EXPERIMENTS

This section details our implementations of the multi-level approach and presents their evaluation on a range of parallel platforms, which are listed in Tab. I. We verify the correctness of our implementations against the generic multiple precision library MPFR, which is sequential and runs on CPUs only. Our implementations profit from the technical features of the considered architectures: SIMD instructions, fusedmultiply-and-add (FMA) instructions, and multi-threading on CPUs and Xeon Phi; local memory and atomic instructions on GPUs. For example, on the Intel Haswell architecture, we use AVX intrinsics to benefit from the 4-way SIMD and implement TwoSum by replacing half of the 6 additions/subtractions by FMAs (multiplying by one). The latter optimization allows to use both the FP adder and FMA units that can run in parallel on Haswell. Thread-level parallelism is exposed using OpenMP to benefit from multi-core and hardware multi-threading. The final inter-node reduction within a cluster is performed using MPI. On the Intel Xeon Phi co-processor, we benefit from 8-way SIMD by using 512bit vector intrinsics. We perform explicit memory prefetching in order to maximize memory throughput. We develop handtuned OpenCL implementations for NVIDIA and AMD GPUs. They use multiple superaccumulators per work group, which are stored in local or global memory. In order to avoid bank conflicts, superaccumulators are interleaved together to spread their digits among different memory banks. Concurrency between multiple threads, which share single superaccumulator, is handled through atomics.

As baselines, we consider our vectorized and parallelized non-deterministic algorithms. In case of summation, we also use: the deterministic Intel TBB reduction (referred as "TBB deterministic"); the Fast Deterministic Parallel Sum algorithm [START_REF] Demmel | Fast reproducible floating-point summation[END_REF] ("Demmel fast"); the one reduction reproducible summation [START_REF] Demmel | Parallel Reproducible Summation[END_REF] ("ReproBLAS") from the ReproBLAS library; the single-sweep reduction [START_REF] Arteaga | Designing bitreproducible portable high-performance applications[END_REF] with two and three levels ("bitrep2" and "bitrep3", accordingly) from the bitrep library.

Figs. 1a and 1b present the timings achieved by the summation algorithms as a function of the input dataset size N of double precision FP numbers. In the legends of all figures, "Superacc" corresponds to our algorithm that is solely based on superaccumulators; "FPEp + Superacc" stands for our algorithm with FPEs of size p (p = 2 : 8) in conjunction with superaccumulators when needed; "FPEpEE + Superacc" represents our algorithm based on the expansion of size p and the early-exit optimization technique, meaning stop propagating errors when they are zeros. These plots show that on large datasets the performance of the baseline unordered sum as well as the expansion cache is constrained by the memory bandwidth on all platforms. For instance, these algorithms achieve 23.3 GB/s on platform A over the theoretical peak bandwidth of 25.6 GB/s. As expected, for small dataset that fits in cache, the non-reproducible parallel sum outperforms our implementations by a factor of 2 to 5. "TBB deterministic" and MPFR (not shown on the plot) run two and three orders of magnitude, respectively, slower than our implementations. The Fast Deterministic Parallel Sum algorithm requires two passes on the data. In memory-bounded scenarios, its execution is two times slower than the unordered reduction and the other singlepass algorithms. In addition, the overhead of the second kernel call impacts the performance of the GPU implementation. Fig. 2 compare the performance scaling of the reduction algorithms on platform B. The number of MPI process varies from 1 to 64, each of them performing the summation of 16M double-precision FP numbers. This dataset size ensure that we fall in the out-of-cache case. For each process, we measure the local summation time, which is colored red, and the MPI reduction time, which is colored green. For the whole range of processors, the execution time of each algorithm is dominated by the local summation time because of the dataset size. In addition, due to the equal distribution of computations among MPI processes, the computation time is roughly equivalent on the whole range of MPI processes, while the reduction time changes according to the number of MPI processes involved. We normalize the total runtime of each algorithm by the total execution time of the parallel FP summation. The "bitrep2", "bitrep3", "Demmel fast", and "ReproBLAS" algorithms are 2.9 %, 3.1 %, 95.5 %, and 20.6 %, accordingly, slower than the conventional parallel summation on 64 CPUs. In contrast, our "FPE2 + Superacc" and "FPE8EE + Superacc" deliver both correctly rounded and bit-wise reproducible results with the overhead of 9.2 % and 7.8 %, respectively.

We extend the multi-level summation algorithm to dot product by applying TwoProd for exact multiplication of two FP numbers. Fig. 1c shows the performance results of dot product on platform D, see Tab. I. TwoProd induces small performance overhead due to summation of two arrays: one for the result and another for the error. Therefore, the performance penalties are slightly higher than in case of ExSUM. Nevertheless, for large array sizes, ExDOT delivers numerically reproducible results with no performance loses.

On the BLAS level-2, we propose a blocked multi-level reproducible algorithm for triangular solver (ExTRSV). The idea behind is to divide matrix by blocks and apply Ex-TRSV on diagonal blocks and matrix-vector multiplication (ExGEMV) on sub/upper diagonal blocks; the later is benefi- cial for the performance. Both blocked ExTRSV and ExGEMV are built upon ExSUM and ExDOT. Fig. 1d presents the measured time achieved by the substitution algorithms as a function of the matrix size n on platform D, see Tab. I. The implementations with expansions deliver better performance than with superaccumulators only. However, due to switches to superaccumulators at the end of computing each element of the solution as well as when the accuracy provided by expansions is not sufficient, the benefit of using expansions is small. Moreover, the division by diagonal elements and the following accumulated errors have negative impact on the final accuracy of the solution, see Fig. 1e, but not on reproducibility. As an enhancement of the accuracy, we propose to use one step of iterative refinement, which is based on ExGEMV, ExTRSV, and ExAXPY. Although this solution will induce larger performance overhead, the accuracy of the results will be improved significantly.

We also tackle the core of BLAS -matrix-matrix multiplication (GEMM) kernel -by providing a blocked multilevel reproducible ExGEMM algorithm. Fig. 1f presents the performance results achieved by the GEMM implementations as a function of the matrix size n on platform E, see Tab. I. The performance of ExGEMM, even with FPEs in conjunction with superaccumulators is limited, because the non-deterministic double precision GEMM already squeezes the architecture performance and does not leave much resources for our approach. That leads to 12-16 times performance overhead on NVIDIA and AMD GPUs. We consider to ameliorate the performance of both ExTRSV and ExGEMM in order to be within the 10 times performance overhead.

V. CONCLUSIONS AND FUTURE WORK

We presented a multi-level approach that delivers both accurate and bit-wise reproducible results for fundamental linear algebra operations such as the ones included in the BLAS library. We derived an algorithm for parallel FP reduction, which is a common operation among the BLAS routines, and provided implementations on multi-and many-core architectures. For large sums with more than 10 7 elements of moderate dynamic range, which cover the most common cases in practice, the proposed implementations deliver comparable performance to the highly optimized parallel FP summations along with bit-wise reproducible results on a large range of platforms. Moreover, our summation algorithm either comparable to or outperforms the existing alternatives.

We extended the multi-level approach to dot product, triangular solver with one right-hand side, and matrix-matrix multiplication. ExDOT delivers similar performance results to ExSUM. Even though the performance of ExTRSV and ExGEMM lag behind their non-reproducible counterparts, their outputs are consistently reproducible and accurate, in terms of rounding-to-nearest in case of the later, independently of threads scheduling and data partitioning. We proposed to use the iterative refinement technique in order to improve accuracy of ExTRSV, which is saturated by the division.

Our roadmap includes deriving a complete set of reproducible, accurate, and fast BLAS routines on parallel architectures from desktop and server processors to Intel Xeon Phi co-processors and GPU accelerators. We plan to conduct a priori error analysis of the derived ExBLAS (Exact BLAS) routines. More information on the ExBLAS project as well as its sources can be found at https://exblas.lip6.fr/.

Fig. 1 :

 1 Fig. 1: Performance results on various architectures from Tab. I.

Fig. 2 :

 2 Fig. 2: Performance scaling of parallel reduction on platform B.

TABLE I :

 I Hardware platforms employed in the evaluation.

	A	Intel Core i7-4770	4 cores with Hyper-Threading	3.400 GHz
	B	Intel Xeon E5-4650L	64 nodes with 2 × 8 cores	2.600 GHz
	C	Intel Xeon Phi 5110P	60 cores × 4-way Multi-Threading	1.053 GHz
	D	NVIDIA Tesla K20c	13 SMs × 192 CUDA cores	0.705 GHz
	E	AMD Radeon HD 7970	32 CUs × 64 units	0.925 GHz

ACKNOWLEDGEMENT

This work undertaken in the framework of CALSIMLAB is supported by the public grant ANR-11-LABX-0037-01 and used the HPC resources of ICS funded by Region Île-de-France and the project Equip@Meso (ANR-10-EQPX-29-01) both overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (ANR-11-IDEX-0004-02). This work was partially supported by the project FastRelax (ANR-14-CE25-0018-01).