Ian Masliah
email: ian.masliah@lri.fr

Marc Baboulin
email: marc.baboulin@lri.fr

Joel Falcou
email: joel.falcou@lri.fr

Meta-programming and Multi-stage Programming for GPGPUs

GPGPUs and other accelerators are becoming a mainstream asset for high-performance computing. Raising the programmability of such hardware is essential to enable users to discover, master and subsequently use accelerators in day-to-day simulations. Furthermore, tools for high-level programming of parallel architectures are becoming a great way to simplify the exploitation of such systems. For this reason, we have extended NT 2 -the Numerical Template Toolbox -a C++ scientific computing library which can generate code for SIMD and multi-threading systems in a transparent way. In this paper, we study how to introduce an accelerator-based programming model into this library to allow developers to reap the benefits of such an architecture from a simple, MATLAB like code. After a brief description of the NT 2 framework, we explain how our accelerator programming model has been designed and integrated in a pure C++ library. We conclude by showing the applicability and performance of this tool on some practical applications.

I. INTRODUCTION

Developing large applications in a simple, fast and efficient way has always been an issue for software developers. As computing hardware complexity rose with the advent of SIMD, multi-processor, multi-core systems and more recently accelerators like GPUs or Intel Xeon Phi , software design methodologies did not undergo the same amount of changes. This complicates the exploitation of such hardware in mainstream applications.

Designing Domain Specific Languages (or DSL ) has been presented as a solution to these issues: DSLs allow solutions to be expressed in their programming idiom with a level of abstraction equivalent to the problem domain leading to improvements in maintainability and quality of code. One of the most popular examples is MATLAB TM which provides a large selection of toolboxes that allow a direct expression of high-level algebraic and numerical constructs in a easy-to-use imperative language. In this scope, Domain Specific Embedded Languages (or DSELs ) [START_REF] Hudak | Building domain-specific embedded languages[END_REF], [START_REF] Tratt | Model transformations and tool integration[END_REF] are languages implemented inside a general-purpose host language [START_REF] Czarnecki | DSL Implementation in MetaOCaml, Template Haskell, and C++[END_REF]. This removes the requirement of a dedicated compiler or interpreter as they are designed as a library-like component [START_REF] Czarnecki | Generative Programming and Active Libraries[END_REF], [START_REF] Veldhuizen | Active Libraries: Rethinking the roles of compilers and libraries[END_REF].

NT 2 -The Numerical Template Toolbox -is such a DSEL using C++ template meta-programming [START_REF] Abrahams | C++ Template Metaprogramming: Concepts, Tools, and Techniques from Boost and Beyond (C++ in Depth Series)[END_REF] to provide a MATLAB -inspired API while supporting a large selection of parallel architectures and keeping a high level of expressiveness [START_REF] Esterie | The numerical template toolbox: A modern c++ design for scientific computing[END_REF]. NT 2 was designed to support architectural features like SIMD extensions and multi-core programming [START_REF] Tran Tan | Automatic taskbased generation for high performance domain specific embedded language[END_REF]. However, the support for accelerators like GPGPUs was limited as GPU kernel compilers where unable to process NT 2 C++ 11 based implementation of C++ based DSEL . Furthermore, libraries such as Thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF] for CUDA can only be used in .cu files and compiled with nvcc which furthers limits the possibilities for code generation.

In this paper, we present a new extension for NT 2 that takes care of such accelerators, especially CUDA based GPUs through multi-stage programming [START_REF] Eckhardt | Implicitly heterogeneous multi-stage programming[END_REF] (or MSP ) for linear algebra and element-wise problems. MSP consists in doing multiple compilation phases allowing for type-safe program generation. Our contributions include:

• A programming model supporting both implicit or explicit data transfers to and from GPGPUs with a simple user interface • An adaptable strategy to generate CUDA kernel directly from a single C++ source file containing NT 2 statements • The integration of this kernel generator with existing CUDA libraries like cuBLAS or MAGMA. The purpose of this system is to provide the user some leeway on how to distribute the data between the host and device through a simple mechanism. As having a perfect cost model for load balancing is very complex to put in place and costly, letting the user provide some insight on data locality is beneficial.

After reviewing the concurrent state of the art software libraries (section II), we introduce NT 2 , its programming model (section III) and how we have adapted it to support GPU computation. We then describe the kernel generator process and how it integrates with the existing library (section IV). Finally, we present benchmarks assessing the generated code quality (section V) and conclude on the future work regarding NT 2 and accelerators.

II. RELATED WORKS

Software libraries for GPGPU computing try to simplify the new programming paradigm brought by many-core based systems. The general trend followed in recent years by C++ libraries is to provide a high-level interface through template meta-programming techniques. This goal is reached by providing device containers with architecture-aware generic parallel algorithms and/or a code generation based process. This approach differs from what can be seen with OpenACC [START_REF] Wienke | Openaccfirst experiences with real-world applications[END_REF] or Halide [START_REF] Ragan-Kelley | Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF] as it is a DSEL based approach and does not rely on language extensions or pragmas. We will give a detailed explanation of both type of libraries that support OpenCL/CUDA or both.

Thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF] is a header only library providing a similar interface to the C++ Standard Template Library. Its high-level interface is based on meta-programming and traits to provide efficient parallel skeletons that can run on either CPU or GPU. Container locality is expressed through an explicit container declaration limiting the abstraction but allowing for easier transfers between host and device. Locality for functions is defined by default for device vectors and can be extended with tag dispatching. Overall, it is a well rounded utility library which can be combined with CUDA Toolkit libraries such as CUBLAS, CUFFT and NPP. However, it can only be used in .cu files and compiled with nvcc which limits its usefulness.

VexCL [START_REF] Demidov | Programming CUDA and OpenCL: a case study using modern C++ libraries[END_REF] is an expression template library for OpenCL/CUDA. It provides a high-level generic interface that is suitable for both back-ends with static parameters defined within a DSEL for linear algebra. The expression template mechanism allows for code generation by lazy evaluation of vectors and elementary operations within the AST. Similarly to Thrust, it provides STL-like functions on containers that have a defined locality. It is also possible for the user to define custom functions on device that will be dynamically generated for CUDA. However, the transform used for the generation process requires a unique data locality limiting hybrid algorithms.

ViennaCL [START_REF] Rupp | ViennaCL-a high level linear algebra library for GPUs and multi-core CPUs[END_REF] is also an expression template library for OpenCL/CUDA. This library strictly follows the uBLAS programming interface and STL like algorithms for easier integration with other softwares. It has implementations for BLAS kernels and high-level solvers for sparse and dense computation that provide good performance. Through the mechanism of expression templates it can evaluate basic linear algebra operations with operator overloading. ViennaCL focuses more on OpenCL due to the necessity for separate compilation with CUDA limiting its support through the OpenCL language. It is however possible to generate CUDA code with ViennaCL at runtime. [START_REF] Lutz | Boost.Compute[END_REF] is a header only C++ library based on the OpenCL standard. Similar to other libraries, it manages device memory through a designated container. It provides an interesting concept of future for asynchronous copy on the device allowing for more versatility. Boost.Compute also supports closures and adaptable structures for device. Similarly to thrust, it is a well rounded library based on OpenCL to simplify the coding process on accelerators. It however lacks support for numerical analysis and cannot generate CUDA code.

Boost.Compute

Eigen [START_REF] Guennebaud | Eigen: A C++ Linear Algebra Library[END_REF] In Figure 1, we compare features between NT 2 and the previously described libraries. We did not include Eigen as it does not have any real support for code generation or a device API.

The purpose of the previously described libraries is usually to provide a wrapper over the C++ language for GPGPU computing. For this reason, the DSELs based on expression templates or Boost.Proto are usually lightweight and consist mostly of overloading elementary operations for containers. The code generation phase then ends up doing a dynamic compilation of the CUDA code/OpenCL kernel which adds a significant overhead on small sized or low arithmetic intensity problems. Furthermore, as mentioned in Section I, it is not possible to compile an NT 2 source code with nvcc even with the latest version (7.5). To address these issues, we have added a two step compilation in NT 2 using a serialization mechanism based on Boost.Serialization. Before describing this process, we will first detail the execution model for GPGPU computing in NT 2 . III. NT 2 EXECUTION MODEL NT 2 [START_REF] Esterie | The numerical template toolbox: A modern c++ design for scientific computing[END_REF] is a numerical computing C++ library implementing a subset of the MATLAB language as a DSEL . NT 2 simplifies the development of data-parallel applications on a large selection of architectures currently including multi-core systems [START_REF] Tran Tan | Automatic taskbased generation for high performance domain specific embedded language[END_REF] with SIMD extensions [START_REF] Estérie | Boost.SIMD: generic programming for portable SIMDization[END_REF]. Simply put, a MATLAB program can be converted to NT 2 by copying the original code into a C++ file and performing minor cosmetic changes (defining variables, calling functions in place of certain operators). NT 2 also takes great care to provide numerical precision as close to MATLAB as possible, ensuring that results between a MATLAB and an NT 2 code are sensibly equal.

Internally, NT 2 is designed to leverage the well known Expression Templates C++ idiom to build at compile time a flexible representation of the abstract syntax tree of any C++ expression containing at least one NT 2 component. This compile-time tree is then transformed in actual code to be executed on a parallel system. Contrary to other libraries based on the same technique, NT 2 relies on Boost.Proto , an external Expression Templates system to handle the creation and transformation of ASTs [START_REF] Niebler | Proto : A compiler construction toolkit for DSELs[END_REF]. Boost.Proto allows us to replace the direct walk-through of the compile-time AST done in most C++ DSELs by the execution of a mixed compiletime/runtime algorithm over the predefined AST structure generated by Boost.Proto (fig. 2). Finally, the other main difference between NT 2 and similar tools is the fact that architectural specificities of the generated code are handled by an execution model based on Algorithmic Skeletons [START_REF] Cole | Bringing skeletons out of the closet: a pragmatic manifesto for skeletal parallel programming[END_REF]. Those skeletons simplify the extension of NT 2 for various hardware by separating the concerns of optimizing the AST code for different type of architecture.

A. Support for CPU/GPU execution

If the single system computation model of NT 2 is rather classic, we needed a way to handle accelerators in a generic and extensible way. The first step is to provide a simple way to locate tables on either the host system or on a device. This is simply done by using a couple of settings in the definition of NT 2 table instances. Listing 1 show cases the nt2::host_ and device_,pinned_ settings that specifies if a table contains data stored on the host memory or device memory. Note that a special function on_device can be used to specify on which device the memory must be allocated in the case where multiple devices are available.

Semantic of operations between host and device tables is quite straightforward as they will be carried on the proper memory segment of each table. When mixing tables of different location, memory transfers are implicitly performed. This means that assigning a host/pinned table to a device table is equivalent to performing a CUDA memory transfer. This can be used for example to simplify interaction with existing GPU kernels as shown in listing 2. As streams are no assigned to tables, this transfer will be synchronous.

/ / G e n e r a t e s a h o s t t a b l e by d e f a u l t t a b l e <d o u b l e> A( o f s i z e ( 1 e3 , 1 e3 ) ) ; / / G e n e r a t e s a h o s t t a b l e w i t h c u d a h o s t p i n n e d memory

A copy function is also available to perform asynchronous memory transfers when a non-default stream is given. This semantic of transfer by assignment is a classical way of performing such operation transparently. It as been used by tools like Thrust or VexCL and have been proved to be easy enough for the user while allowing for fine grain performance tuning.

B. MAGMA Integration

Our extension to NT 2 also provides a direct high-level integration of MAGMA [START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF] for handling most linear algebra tasks as shown on listing 3 for the linsolve. Calls to linsolve can mix device and host tables. The implementation of such kernels will take care of transferring the strictly required data over the device and to perform transfer back to the host only if the result is assigned to a host table. Most MAGMA kernels are mapped onto their MATLAB equivalent API and handle the same sets of optional behaviors.

C. Kernel optimizations

As complex programs often feature multiple statements involving NT 2 tables, it may become hard in some cases to maintain a high level of performance as the cost of transfers between statements may become significant. To locally solve this issue, NT 2 provides a function called tie that allows for library-based loop fusion as depicted in Listing 4. In this scenario, the single assignment statement between call to tie will generate exactly one call to a single CUDA kernel composed of the fusion of the two loop nests over the variable A and X. Every table will be transferred in a single sequence of streaming operations and a single kernel will be executed on the device.

IV. DEVICE CODE GENERATION

As described in section II, code generation for accelerators in C++ is based on Meta-programming techniques. This process is however limited by the fact that C++ does not natively support language extensions for runtime code generation and program execution [START_REF] Taha | A gentle introduction to multi-stage programming[END_REF]. Creating a multi-stage paradigm is therefore necessary to enable the compilation of CUDA code within a C++ framework.

A. Multi-stage programming in C++

MSP is usually applied by adding specific language extensions to trigger a new compilation phase. an example, MetaOcaml [START_REF] Taha | Metaocaml-a compiled, type-safe, multi-stage programming language[END_REF] is a multi-stage language based on Ocaml with three basic constructs for runtime code optimization.

More recent adoptions of MSP for parallel computing include techniques like Lightweight modular staging (LMS) [START_REF] Rompf | Lightweight modular staging: a pragmatic approach to runtime code generation and compiled dsls[END_REF] in SCALA, Language Virtualization [START_REF] Chafi | Language virtualization for heterogeneous parallel computing[END_REF] (base on LMS) or Terra [START_REF] Devito | Terra: a multi-stage language for high-performance computing[END_REF] using language extensions for HPC in LUA . LMS is similar to Meta-programming techniques in C++ applied to the Domain Engineering Method for Reusable Algorithmic Libraries [START_REF] Czarnecki | Generative programming -methods, tools and applications[END_REF] methodology on which NT 2 is build. Language Virtualization is an interesting concept but it is limited by the need to make the DSEL identical to the stand-alone language. This greatly increases the complexity of the generation process and makes it hard to extend.

Our approach to MSP in C++ is based on doing a compilation phase with an incomplete link to generate an object file. This object file can then be demangled (with a tool like cppfilt or nm ) to decode the C++ ABI names. Each demangled symbol will correspond to the internal C++ representation with a complete prototype of the incomplete function. By parsing this representation we can generate the CUDA/OpenCL kernel code with a corresponding host code to complete the link phase. Figure 3 describes this multi-stage process. A specific tool called symbol/code converter (Figure 3, Part 3) coupled with a serialization process was developed as an add-on to NT 2 to solve the issues described previously. The serialization process used to solve the readability of the abstract representation is based on Boost.Serialization . This representation is similar to the expression-template AST based on Boost.Proto that is parsed by the C++ compiler. This enables us to get the semantic information of the container described in section III like its data locality, data type or matrix shape. Having run-time informations like the container size is not possible as it is dynamic.

The role of symbol/code converter is two-fold : to parse the demangled symbols obtained and to generate the device and host code. To demangle the symbols, we use Boost.Spirit [START_REF] De Guzman | The boost spirit parser generator framework[END_REF], a C++ library to parse expressions and generate outputs based on them. These outputs correspond to the semantic information of the containers and operators. We can then generate the device and host code (Figure 3 The code in Listing 5 corresponds to Part 1 of Figure 3. During the compilation phase, the operation in Listing 5 is replaced with a call to the CUDA transform skeleton. This skeleton will then call the device kernel that is not implemented yet resulting in the incomplete link sequence (Part 2, Figure 3). The equivalent code for transform is detailed in Listing 6. This transform will analyze the informations on the input matrices and the architecture to decide if it should proceed with the generation process. From there, once we compile our code and go through phase 2 of Figure 3 we obtain the mangled code for our transform. We can then demangle the symbols resulting in the simplified code described in Listing 7. This code corresponds to the Boost.Proto AST in NT 2 that we parse to generate the host and device code. The architectural tag of the machine is depicted in purple in Listing 7 and the fma computation node in red. The function tag in gray corresponds to each element-wise operations in NT 2 . Any function can be constructed using external kernel if necessary. As an example, algorithmic skeletons (scan, reduce, zip ...) will call external kernel when the CUDA back-end is activated. Then, we can generate a thrust code in a .cu file that will then be compiled with nvcc. This enables us to overcome a weakness of the thrust library while benefiting from its optimized routines.

The representation of the AST in Figure 4 corresponds to the operation obtained from the demangled symbol that we parsed. The generated code for the .cu file is described in Listing 8 for Kepler cards. It corresponds to the AST representation with additional semantic information specific to the CUDA language. A function wrapper that calls the CUDA kernel (fma4 wrapper) is used to separate the compilation of the .cu file with nvcc from the corresponding C++ host code. The triad kernel directly calls the fma function from CUDA as the NT 2 AST recognizes every occurrence of an fma and replaces it by its function (Listing 7, see boost::simd::tag::fma ). This optimization in itself is not essential as it can be done by nvcc but demonstrate the potential of code generation. As NT 2 already optimizes the AST by replacing patterns with corresponding functions, we can benefit from this analysis. We can then call the corresponding CUDA function if available or our own implementation for each pattern.

The host code is partially described in Listing 9 and 10. Listing 9 corresponds to the initialization of parameters and data before doing the actual computation on the device. The blockSize and stream number are determined during the generation process depending on the number of parameters and architecture through prior tests. The block-size is generated by measuring the bandwidth of data transfers from host to device in a range and choosing the most optimal one. To benefit best from our Kepler GPU bandwidth for the triad kernel we use a block size of 40000 for single precision values. Double precision computations or a high number of parameters to transfer will lead to a lower block size generated. To decide the block dimension for CUDA, we use is the high value available to the card as element-wise operations are highly parallel and have no dependencies. The allocation process includes pinned memory and device memory allocation optimizations. Containers located on the GPU (defined with nt2 :: device ) do not appear in the allocation phase. To limit redundancy in data transfers between outputs and inputs ( as we have no information on the pointer) we use an unordered set.

Listing 10 describes the computation phase. It relies on block streaming with transfers to the GPU only if NT 2 tables are on the host. This streaming process is based on the overlap data transfers concept described by NVIDIA. It consists in creating multiple streams (the number depends on the architecture and problem intensity/size) and launching for each stream a transfer host to device, the CUDA kernel and the transfers device to host for a block. As the host memory has already been allocated, we must first transfer the data to pinned memory with cudaHostAlloc to benefit from GPU optimizations. Since the difference in bandwidth between pinned and page-able memory only increases with newer architectures, this optimization can give a speedup even with a mono-stream program. If the host table is declared with nt2 :: pinned , it can greatly diminish the time needed to transfer the data. As stated in section III, computations on the device only occur if some conditions are met. As of now, these conditions are limited to the problem size and data locality but can be extended as the call to transform is automatic when NT 2 has defined that CUDA is available. Due to the hierarchical tag dispatching in NT 2 , a system with an Intel processor coupled with an NVIDIA card will have a tag similar to the following : cuda < openmp < simd extension >> . Therefore, if conditions for dispatch on the GPU are not met we will call the next level of transform ( i.e. openmp). This enables us to use both the CPU and GPU in parallel depending on the problem which is a functionality rarely implemented in libraries. We further this process with our MAGMA back-end which does hybrid computations on most of its solvers.

The code generation is hidden from the user as the generation process is done during the standard compilation phase. The compilation overhead is negligible as the analysis of the AST is linear and the code generation is usually not very long as seen above. The user interface only contains the added semantic information while all the complex allocations are hidden from the user.

V. EXPERIMENTS

In this section, we will show that our generation process produces satisfactory performances in most situations. The benchmarks are realized with the following components :

• CPU : 2 x 6 cores Intel Xeon E5-2620 15MB L3, AVX • GPU : Tesla K40m

-Pageable host to device (HTD) : 3 GB/s -Pinned host to device : 9.8 GB/s

• Memory : 65 GB with a memcpy bandwidth of 5GB/s • GCC 4.9, CUDA 7.5

A. Black & Scholes kernel

The Black & Scholes algorithm represents a mathematical model that gives a theoretical estimate of the price of European call and put options on a non-divideend-paying stock. It is a bandwidth bound algorithm for GPU if we take into account the memory transfers.

The code given in Listing 11 uses the loop-fused technique described previously with the operator tie. The nt2 :: device and nt2 :: pinned tags are specific to accelerator enabled architectures. If the tag is used while no accelerator is available, NT 2 will fall back to the default architecture which is nt2 :: host . t i e ( d , d1 , d2 , r ) = t i e ( s q r t ( Ta ) , l o g ( S / X) + ( fma ( s q r ( v ) , 0 . This additional semantic information on memory locality can help to avoid useless memory transfers while staying simple enough for the user. In our experiment, we consider the worst case scenario where all entry tables are on cpu and not allocated as pinned.

This results in the .cu file in Listing 12 generated for floating point values. Since the AST does not contain the name of the parameters, the kernel generator has to give a different name to each one. This does not lead to a change in performance for the kernel as we just pass multiple times the same pointer. Memory allocation on device and data transfers between host and device do pointer checking in the host code (see triad example) to insure no redundant work is done incurring also negligible overhead. The fnms function comes from an NT 2 AST transformation and corresponds to the fused negated multiply-subtract of three values. As seen in Figure 5, the execution time of the kernel on the GPU is negligible (3 ms) compared to the overall time of the optimized version with SIMD and OPENMP in NT 2 . Most of the time is spent transferring the data between host and device memory which can be avoided with the right semantic information available ( nt2 :: device ) . In the scenario presented ( Listing 11) , which is the worse case, we still have better performance than the Boost.Compute version on which most C++ libraries are based as we use block streaming with pinned memory reaching a near-optimal throughput (average of 9.7 GB/s) for device transfers. As the bandwidth of a memcpy on the CPU (5 GB/s) is faster than page-able transfers (3GB/s) even without any overlap between transfers and computation we still increase the performance. This optimization can be disabled depending on the CUDA architecture.

As computations on the GPU are often done in great number, if the user allocates the data on the GPU he will pay no transfer cost for the rest of the computations and in this case the CUDA kernel is up to twelve times faster than the hand-optimized version for CPU. Also, using the flag nt2 :: pinned if tables must stay on cpu allows us to gain up to a x3 in transfer time. The optimized Black & Scholes version we implemented in CUDA reached the same performance in kernel and transfers time when using the same type of memory allocation for data.

B. Linsolve kernel

Linsolve is a MATLAB routine for solving linear systems. As NT 2 has its own implementation of linsolve with a LAPACK or MAGMA back-end, we can combine it with the code generation process. It also supports mixed-precision algorithms in both CPU and GPU versions through LAPACK/MAGMA or their own implementation. In this benchmark, we consider the solution of a dense linear system using the LU factorization and apply one step of iterative refinement [START_REF] Baboulin | Accelerating scientific computations with mixed precision algorithms[END_REF] :

1) Compute r = b -Ax.

2) Solve Ad = r.

3) Update y = x + d. The equivalent NT 2 code is the following : If executed on the CPU, the code in Listing 13 will call the LAPACK routines. The semantic information upper triangular allows linsolve to call the triangular solver instead of doing the classic linear solve. If executed on the GPU, the same optimizations will be applied and the iterative refinement process will trigger calls to transform for both element-wise operations.

The performance results in Figure 6 attest that the performance obtained with our model is relevant. The GPU version with MSP calls magma kernels using the CUBLAS dgemm routine without doing any transfer and reaches near peak performance of a K40m GPU which corresponds to 1.40 Tflop/s. The version that does not use MSP is slower as transfers are done during the iterative refinement step. The CPU version quickly reaches the peak performance of both CPU which is 210 Gflop/s. As we can see, there is no performance loss while call the LAPACK/MAGMA back-ends and if device pointers are passed to our code generator, there will be no memory transfers. Similar performance would also be reached using the other factorizations available in NT 2 . The development of tools for simplifying accelerator programming has been an active topic since accelerators have become a mainstream element of high-performance computing systems. In this paper, we proposed an extension of a highlevel, data-parallel scientific computing library to specifically handle GPU accelerators. Our main objectives were to keep a similar level of expressiveness in the client code with a MATLAB -like interface while supporting different use cases of accelerator programming.

To reach this goal, we have implemented a multi-stage system in which the initial C++ code is used to automatically generate the equivalent CUDA kernel by reusing our internal representation of this code. This representation is based on C++ Expression Templates and the Algorithmic Skeleton to identify and classify expressions based on the kind of loop nest that is required. Finally, we showed on a selection of examples that the performance obtained is close to the hardware capability and exhibits benefits compared to other solutions.

Work is still on-going on this system, including the final integration into the main NT 2 release and support for more specific functions on the latest GPUs. Implementing a more thorough cost model to ensure better scheduling of computation between CPU and GPU is also being studied.

Fig. 2 .

 2 Fig. 2. Templates in NT 2

  t a b l e <d o u b l e , p i n n e d > A2 ( o f s i z e ( 1 e3 , 1 e3 ) ) ; / / G e n e r a t e s a d e v i c e t a b l e t a b l e <d o u b l e , d e v i c e > D( o f s i z e ( 1 e3 , 1 e3 ) ) ; / / G e n e r a t e s a d e v i c e t a b l e on d e v i c e #2 t a b l e <d o u b l e , d e v i c e > D2 ( o f s i z e ( 1 e3 , 1 e3 ) , o n d e v i c e ( 2 ) ) ; Listing 1. NT 2 host and device specifications

  / / X i s a 1 e3 x 1 e3 m a t r i x f u l l o f 1 . t a b l e <d o u b l e> X = o n e s ( 1 e3 , 1 e3 ) ; / / T r a n s f e r t o d e v i c e t a b l e <d o u b l e , d e v i c e > Y = X ; / / cuBLAS d i r e c t c a l l c u b l a s D s c a l ( Y . s i z e ( ) , 5 . , Y . d a t a ( ) , 1 . ) ; / / T r a n s f e r b a c k t o h o s t X = Y ; Listing 2. NT 2 interaction with cuBLAS

  t a b l e <f l o a t > A, B , X ; t a b l e <f l o a t , d e v i c e > dA , dB , dX ; X = l i n s o l v e (A, B ) ; dX = l i n s o l v e (A, dB ) ; Listing 3. NT 2 interface to MAGMA kernels

  t a b l e <T> b s ( t a b l e <T> c o n s t& B , t a b l e <T> c o n s t& C , , t a b l e <T> c o n s t& Y, t a b l e <T> c o n s t& Z ) { t a b l e <T> A( B . e x t e n t ( ) ) , X( Y . e x t e n t ( ) ) ; t i e (A, X) = t i e ( B + C * s i n ( B ) , Y / ( Z + A ) ) ; r e t u r n X; } Listing 4. NT 2 merged kernels

  , Part 4-5) by de-serializing the abstract representation obtained with Boost.Spirit . B. Integration in NT 2 We use the Triad kernel to describe the generation process. This consists in doing a fused multiply-add (or fma : a = b + c * d ). The resulting code in NT 2 is : / / D e f i n e h o s t t a b l e t a b l e <f l o a t > A, B , C , D ; / / T r i a d k e r n e l A = B + C * D ; Listing 5. NT 2 Triad kernel

  t a b l e <f l o a t > A, B , C , D ; / / T r i a d k e r n e l r e p l a c e w i t h t r a n s f o r m c a l l t r a n s f o r m (A, B + C * D) ; Listing 6. NT 2 Triad transform

Fig. 4 .

 4 Fig. 4. Triad kernel transform AST

  u s i n g b o o s t : : p r o t o : : c h i l d c ; u s i n g b o o s t : : p r o t o : : v a l u e ; s i z e t s i z e = numel ( b o o s t : : p r o t o : : c h i l d c <0>(a1 ) ) ; s i z e t b l o c k S i z e = s t d : : min ( s t d : : s i z e t ( 4 0 0 0 0 ) , s i z e ) ; s i z e t n S t r e a m s = s t d : : min ( s t d : : s i z e t ( 2 ) , s i z e / b l o c k S i z e ) ; s i z e t n = s i z e / b l o c k S i z e ; s i z e t l e f t o v e r = s i z e % b l o c k S i z e ; dim3 blockDim = s t d : : min ( s t d : : s i z e t ( 1 0 2 4 ) , s i z e ) ; dim3 d i m G r i d = b l o c k S i z e / s i z e t ( 1 0 2 4 ) ; c u d a S t r e a m t s t r e a m [ n S t r e a m s ] ; / / A l l o c a t i n g memory on t h e d e v i c e v a l u e ( a0 ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , n S t r e a m s , s i z e , t r u e ) ; v a l u e ( c h i l d c <0>(a1 ) ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , n S t r e a m s , s i z e ) ; v a l u e ( c h i l d c <1>(a1 ) ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , n S t r e a m s , s i z e ) ; v a l u e ( c h i l d c <2>(a1 ) ) . s p e c i f i c s ( ) . a l l o c a t e ( b l o c k S i z e , n S t r e a m s , s i z e ) ; / / c h e c k s r e d u n d a n c y b e t w e e n i n p u t s and o u t p u t s s t d : : u n o r d e r e d s e t <c o n s t f l o a t * > a d d r ; a d d r . i n s e r t ( c h i l d c <0>(a0 ) . d a t a ( ) ) ; f o r ( s t d : : s i z e t i = 0 ; i < n S t r e a m s ; ++ i ) { c u d a S t r e a m C r e a t e (& s t r e a m [ i ] ) ; } Listing 9. NT 2 CUDA Triad Host Code 1

  f o r ( s t d : : s i z e t i = 0 ; i < n ; ++ i ) { s t d : : s i z e t j = i % n S t r e a m s ; v a l u e ( a0 ) . s p e c i f i c s ( ) . t r a n s f e r h t d ( a0 , i , s t r e a m [ j ] , j ) ; v a l u e ( c h i l d c <0>(a1 ) ) . s p e c i f i c s ( ) . t r a n s f e r h t d ( c h i l d c <0>(a1 ) , i , s t r e a m [ j ] , j , a d d r ) ; v a l u e ( c h i l d c <1>(a1 ) ) . s p e c i f i c s ( ) . t r a n s f e r h t d ( c h i l d c <1>(a1 ) , i , s t r e a m [ j ] , j , a d d r ) ; v a l u e ( c h i l d c <2>(a1 ) ) . s p e c i f i c s ( ) . t r a n s f e r h t d ( c h i l d c <2>(a1 ) , i , s t r e a m [ j ] , j , a d d r ) ; f m a 4 w r a p p e r ( v a l u e ( a0 ) . s p e c i f i c s ( ) . d a t a ( j ) , v a l u e ( c h i l d c <0>(a1 ) ) . s p e c i f i c s ( ) . d a t a ( j ) , v a l u e ( c h i l d c <1>(a1 ) ) . s p e c i f i c s ( ) . d a t a ( j ) , v a l u e ( c h i l d c <2>( a1 ) ) . s p e c i f i c s ( ) . d a t a ( j ) , dimGrid , blockDim , s t r e a m [ j ] ) ; b o o s t : : p r o t o : : v a l u e ( a0 ) . s p e c i f i c s ( ) . t r a n s f e r d t h ( a0 , i , s t r e a m [ j ] , j ) ; } i f ( l e f t o v e r ! = 0 ) { . . . } Listing 10. NT 2 CUDA Triad Host Code 2

  t a b l e <T> b l a c k s c h o l e s ( t a b l e <T> c o n s t& S , t a b l e <T> c o n s t& X , t a b l e <T> c o n s t& Ta , T c o n s t r , T c o n s t v ) { a u t o s = e x t e n t ( Ta ) ; t a b l e <T , d e v i c e > d ( s ) , d1 ( s ) , d2 ( s ) ; t a b l e <T> r ;

  5 f , r ) * Ta ) / ( v * d ) , fma(-v , d , d1 ) , S * n o r m c d f ( d1 )-X * exp(-r * Ta ) * n o r m c d f ( d2 ) ) ; r e t u r n r ; } Listing 11. NT 2 black and scholes

Fig. 5 .

 5 Fig. 5. Black and Scholes Performance Comparison (in ms)

t a b l

  e <T , d e v i c e > A, b , x , e ; t a b l e <T , s e t t i n g s ( d e v i c e , u p p e r t r i a n g u l a r )> u ; t a b l e <T , s e t t i n g s ( d e v i c e , l o w e r t r i a n g u l a r )> l ; t i e ( l , u ) = l u (A) x = l i n s o l v e ( l , b ) ; / / l o w e r t r i a n g u l a r s o l v e x = l i n s o l v e ( u , x ) ; / / u p p e r t r i a n g u l a r s o l v e / / One-s t e p r e f i n e m e n t e = bn t 2 : : mtimes (A, x ) ; e = n t 2 : : l i n s o l v e ( l , e ) ; e = n t 2 : : l i n s o l v e ( u , e ) ; x = x + e ; Listing 13. NT 2 LU linear solve with iterative refinement

Fig. 6 .

 6 Fig. 6. Performance comparison of NT 2 linear solve (in Gflop/s)

  is a popular library to solve linear systems using expression templates. It should be able to support accelerator code with CUDA by writing Eigen code in a .cu file. It however does not provide any high-level interface or special containers for GPU computing. Eigen does not support OpenCL.

	Feature	Thrust	VexCL	ViennaCL	Boost.C	NT 2
	MATLAB API	-	-	-	-
	AST optimization	-
	Device Arrays
	Cuda code gen			-	-
	OpenCL code gen	-
	parallel skeletons			-
	CUBLAS support		-		-
	Static code gen	-	-	-	-
	dense LA solvers	-	-		-
	sparse LA solvers	-	-		-	-
	Fig. 1. Feature set comparison between NT 2 and similar libraries

  Listing 12. NT 2 black and scholes fused cuda kernelThe Black & Scholes algorithm involves high latency and high register count operations. This results in sub-optimal performance on SIMD for the CPU due to spilled registers while a Kepler GPU has no such problem.

	g l o b a l	v o i d b s ( f l o a t * t 0 , f l o a t * t 1 , f l o a t * t 2 ,
	f l o a t * t 3 , c o n s t f l o a t * t 4 , c o n s t f l o a t * t 5 , c o n s t
	f l o a t * t 6 , c o n s t f l o a t t 7 , c o n s t f l o a t * t 8 , c o n s t
	f l o a t t 9 , c o n s t f l o a t * t 1 0 , c o n s t f l o a t t 1 1 , c o n s t
	f l o a t * t 1 2 , c o n s t f l o a t * t 1 3 , c o n s t f l o a t * t 1 4 ,
	c o n s t f l o a t t 1 5 , c o n s t f l o a t * t 1 6 , c o n s t f l o a t * t 1 7 ,
	c o n s t f l o a t * t 1 8 , c o n s t f l o a t * t 1 9 )
	{
	i n t i = b l o c k I d x . x * blockDim . x+ t h r e a d I d x . x ;
	t 0 [ i ] = s q r t f ( t 4 [ i ] ) ;
	t 1 [ i ] = p l u s ( l o g f ( d i v i d e s ( t 5 [ i ] , t 6 [ i ] ) ) , d i v i d e s ( m u l t i p l i e s
	( t 7 , t 8 [ i ] ) , m u l t i p l i e s ( t 9 , t 1 0 [ i ] ) ) ) ;
	t 2 [ i ] = fnms ( t 1 1 , t 1 2 [ i ] , t 1 3 [ i ] ) ;
	t 3 [ i ] = fnms ( m u l t i p l i e s ( t 1 4 [ i ] , e x p f ( m u l t i p l i e s ( t 1 5 , t 1 6 [ i ] )
	) ) , f a s t n o r m c d f ( t 1 7 [ i ] ) , m u l t i p l i e s ( t 1 8 [ i ] , f a s t n o r m c d f (
	t 1 9 [ i ] ) ) ) ;
	}


Arnaud Ralko

Didier Frédéric Mila

Didier Poilblanc

Phase separation and flux quantization in the doped quantum dimer model on the square and triangular lattices

Keywords: numbers: 75.10.Jm, 05.50.+q, 05.30.-d

HAL is

Understanding electron pairing in high temperature superconductors is a major challenge in strongly correlated systems. In his milestone paper, Anderson proposed a simple connection between high temperature superconductors and Mott insulators [1]. Electron pairs "hidden" in the strongly correlated insulating parent state as Valence Bond (VB) singlets lead, once fried to move at finite doping, to a superconducting behavior. A very good candidate of the insulating parent state is the resonating VB state (RVB), a state with only exponentially decaying correlations and no lattice symmetry breaking. A simple realization of RVB has been proposed by Rokhsar and Kivelson (RK) in the framework of an effective quantum dimer model (QDM) with only local processes and orthogonal dimer coverings [2]. Even though the relevance of these models for the description of SU(2) Heisenberg models is not fully established, this approach is expected to capture the physics of systems that naturally possess singlet ground states (GS). Indeed, specific quantum dimer models have recently been derived from a spin-orbital model describing LiNiO 2 [3], or from the trimerized kagome antiferromagnet [4]. In a recent work, a family of doped QDMs (at T=0) generalizing the so-called RK point of Ref. [2] has been constructed [5], taking advantage of a mapping to classical dimer models [6] that extends the mapping of the RK model onto a classical model at infinite temperature. This investigation has suggested that phase separation into two macroscopic constituents, an undoped phase and a phase with larger doping, could be quite generic in doped QDMs [5]. Despite their ad-hoc construction, the "quasi-classical" QDMs of Ref. [5] are expected to reflect many features of more realistic models such as the two-dimensional quantum hard-core dimer-gas defined by the Hamiltonian:

H = v c N c |c c| -J (c,c ′ ) |c ′ c| -t (c,c ′′ ) |c ′′ c| (1)
where the sum on (c) runs over all configurations of the Hilbert space, N c is the number of flippable plaquettes, the sum on (c ′ , c) runs over all configurations |c and |c ′ that differ by a single plaquette dimer flip, and the sum on (c ′′ , c) runs over all configurations |c and |c ′′ that differ by a single hole hopping (see Ref. [5] for more details). Troughout the energy scale is set by J = 1. A schematic phase diagram for the two lattices is depicted in Fig. 1 in the undoped case. Remarkably, these lattices lead to quite different insulating states. Indeed, an ordered plaquette phase appears on the square lattice immediately away from the special RK point, whereas a RVB liquid phase is present in the triangular lattice In this Letter, we investigate in details the properties of model (1) on the square and triangular lattices at finite doping. Building on the differences between the two lattices in the undoped case, we investigate to which extent the properties of the doped system are governed by the nature of the insulating parent state. This investigation is based on exact Diagonalisations and extensive Green's Function Monte-Carlo (GFMC) simulations [7,8], particularly well adapted for QDM for which all off-diagonal elements are non-positive. [9].

Phase separation: At small t, it is expected that holes experience an effective attractive potential. It is there-fore natural to first address the issue of phase separation (PS), i.e. the possibility for the system to spontaneously undergo a macroscopic segregation into two phases with different hole concentrations. We analyze the problem as a function of the hopping parameter t and hole concentration x = n h /N , where n h is the number of holes in the system and N the number of sites. In order to perform a Maxwell construction (see for example [10]) we define:

s(x) = e(x) -e(0) x (2)
where e(x) is the energy per site at doping x. This quantity corresponds to the slope of the line passing through e(0) and e(x). If the system exhibits PS, the energy will present a change of curvature implying s(x) to have a minimum at a critical doping x c . The fact that the local curvature of e(x) at x = 0 is negative then implies that the two seperated phases will have x = x c and x = 0 (the undoped insulator). In Fig. noticeable differences. While for the square lattice (lower panel) the critical hole concentration x c is roughly size independent, there is a strong size dependence for the triangular lattice (upper panel). This size effect can be traced back to the nature of the parent undoped state. On the square lattice, the crystalline phase (for v < 1) at zero doping is very robust and for increasing size, its local order changes only weakly. On the triangular lattice, it has been shown that size effects play an important role [9], especially in the RVB liquid phase for 0.8 v ≤ 1. Periodic boundary conditions tend to stabilize the socalled √ 12 × √ 12 phase on small clusters, and clusters with more than 192 sites are necessary to significantly reduce finite-size effects, in particular, as in Fig. 2, close to the transition point with the crystalline phase. Hence the PS observed around x = 0.075 for the 3×6×6 cluster is not representative of the thermodynamic limit.

To obtain the phase diagram in the (v, x) plane, we have performed a systematic size-scaling analysis at fixed t and for various v's depicted in Fig. 3 the previous discussion, a significant size dependence is only present for the triangular lattice, in which case PS disappears for large clusters in the RVB phase in the vicinity of the RK point [START_REF]For small values of t, the GFMC algorithm is no[END_REF]. In Fig. 4, we report the thermodynamic limit of x c for the two lattices as a function of v, and for different values of t. For the square lattice, calculations have been done from the RK point down to the expected phase transition between the plaquette phase and the columnar phase, namely v ≃ 0.6 [START_REF] Syljuasen | [END_REF]. For the triangular lattice, the range between the RK point down to the RVB-√ 12 × √ 12 transition point at v ≃ 0.8 has been covered [9]. These results clearly demonstrate the difference between the square and triangular lattices. In the first case, as soon as v = 1, PS occurs for x < x c . Moreover, upon decreasing v, crystalline order strengthens and, for fixed t, it is necessary to consider a higher concentration of holes to reach a stable conducting phase. Similarly, the bigger t, the lower x c . On the triangular lattice, a finite size-scaling analysis shows that no phase separation appears down to a critical value v ∼ 0.9, well above the critical value v ∼ 0.8 below which plaquette order sets in. Although numerical limitations prevent computations for smaller v and t, our results up to the 3 × 12 × 12-site cluster provide clear evidence for a region of PS inside the RVB region, between v ∼ 0.8 and v ∼ 0.9.

Dimer ordering on the square lattice: Next we investi- The dashed lines correspond to the approximate location of the phase transition between plaquette and columnar phases [START_REF] Syljuasen | [END_REF] for the square lattice and between plaquette and RVB phases [9,13] for the triangular lattice.

gate how dimer order, known to exist at x = 0, evolves under finite doping. Two scenarii are a priori possible: i) the dimer order vanishes in the stable conducting phase immediately at x c ; ii) dimer order survives above x c in a narrow region of the conducting phase. To solve this problem, we have calculated the squared order parameter D 2 ( k) defined by:

D 2 ( k) = 1 N Ψ 0 |d(-k)d( k)|Ψ 0 Ψ 0 |Ψ 0 (3)
along the path Γ → M → K → Γ in the first Brillouin zone of the square lattice, using standard notations, where d( k) is the Fourier transform of the dimer operator d i in the u x direction which multiplies a configuration by 1 if there is dimer between site i and its right neighbor, and 0 otherwise. Note that this calculation has not been tried for the triangular lattice since no Bragg peak is present in the RVB phase, and the algorithm is losing efficiency for v 0.8 [START_REF]For small values of t, the GFMC algorithm is no[END_REF]. In the pure plaquette phase on the square lattice, a Bragg peak develops at point k M = (π, 0), the middle of the side of the BZ. We show in the inset of Fig. 5 a typical result for the squared order parameter on the 196-site cluster for different values of x. Clearly, the Bragg peak disappears upon doping. A finite size scaling of the order parameter can be performed thanks to the linear behaviour at low concentration. Within our data,

D 2 ( k = k M ) ≡ D 2
M (L, x) behaves like a L x + b L in the linear region. In this case, one can determine rather precisely x +∞ such that D 2 M (+∞, x +∞ ) = 0 i.e. x +∞ is the concentration in the thermodynamic limit where the Bragg peak vanishes. By definition, x +∞ = -b +∞ /a +∞ ≃ 0.05 (8). This value, and the linear behaviour in the thermodynamic limit, are displayed in Fig. 5. If we compare x +∞ to the corresponding x c from Fig. 4, we can conclude that the Bragg peak indeed vanishes at x c ≃ 0.067(5) within error bars. Note that numerical errors increase for larger clusters, and we are not able to use the same analysis for clusters larger than 256 sites (the results for the 18 × 18 cluster have not been used). Although the determination of x +∞ is delicate, the linear behavior of D 2 vs x is consistent with the physical behavior expected for the binary system of Fig. 4. No dimer order is present above x c , showing that the system is simply "conducting "in this case.

Flux quantization: Finally, let us turn to a better characterization of the "conducting phase". In the model defined by the Hamiltonian of Eq.( 1), the holes are bosonic, and we therefore expect the conducting phase to be a superfluid (through Bose condensation) [14]. However, since these bosons move in a dimer environment, pair correlations might be present. To investigate this issue, we pierce the torus by an Aharonov-Bohm flux of strength φ = 2πξ with 0 ≤ ξ ≤ 1. The flux is simply implemented by the Peierls substitution, changing the hole hopping term into t ′ = exp(±i2πξ/L x )t, where the ± depends on whether the hole hops to the right or to the left, and where L x is the linear size of the system. Obviously, the whole spectrum should be periodic in ξ with period 1. We show in Fig. 6 the spectrum of the 4 × 4 cluster on the square lattice, with 4 holes. It turns out that the spectrum is rigorously periodic with period ξ = 1/2, which means that there is flux quantization in units of half the flux quantum (red curve). This was suggested quite some time ago in the context of a more general QDM by Kivelson[14], who also predicted that, in the cylinder geometry, one can in principle tunnel between the two branches of Fig. 6 by creating finite energy vortexantivortex excitations that are finally pulled off through the edges of the cylinder, lifting the degeneracy at the level crossing. This lifting of degeneracy is not present in our case, neither in the torus geometry, nor in the cylinder geometry. However, this is just a consequence of the translational symmetry, which puts the two states that are degenerate at ξ = 1/4 into different symmetry sectors. Removing the translational symmetry by changing the amplitude of a local dimer flip indeed removes the degeneracy (upper panel of Fig. 6). In that case, the flux quantization in units of half the flux quantum could indeed be detected in an experiment in which the flux is sweeped. Thus, in our model, the ground-state energy has periodicity hc/2e, suggesting the presence of mobile pairs (or elementary particles of charge Q = 2e) in the system [10]. Note that, unlike what was recently found in a bosonic model with correlated hopping [15], these particles are probably not boson pairs (although a definite proof would require the knowledge of the two particle Greens's function not accessible with GFMC, boson pairs are unlikely, at least far enough from the phase separation): From the bosonic point of view, it is the statistical flux of the dimer background that leads to the half-flux quantization. If dimers are interpreted as SU(2) electron singlets, these singlets are the physical pairs that lead to half-flux quantization.

Summary and conclusions: The numerical investigation with Green's function Quantum Monte Carlo and exact diagonalizations of the doped two-dimensional quantum hard-core dimer model on the square and triangular lattices has led to a number of interesting conclusions regarding hole motion in a dimer background. Phase separation is often present at low doping, as suggested by earlier investigations, but our results indicate that it related to the presence of valence bond order: In the RVB phase of the triangular lattice, PS only occurs close to the plaquette phase, where short-range dimer correlations are already strong enough. Close to the RK point, doping the RVB phase leads directly to a conducting (most likely superfluid) phase. To characterize this phase, we have studied the effect of an Aharonov-Bohm flux, with the conclusion that the flux quantization is in units of half a flux quantum, consistent with the idea that the dimer background leads to effective particles of charge 2e. All these results are in qualitative agreement with the gauge theories of high T c superconductivity in strongly correlated systems [16]. Finally, let us emphasize that, thanks to the large sizes that we could simulate with GFMC, all these conclusions are essentially free of the usual finitesize limitations of the numerical investigation of strongly correlated electron models.

FIG. 1 :

 1 FIG.1:(color online) Schematic phase diagrams for the triangular and the square lattice.

10 FIG. 2 :

 102 FIG. 2: (color online) Slope of energy density (Eq.(2)) vs doping for different sizes. (a) Triangular lattice. (b) Square lattice.

15 FIG. 3 :

 153 FIG. 3: (color online) Scaling of the critical doping xc defined by Maxwell construction with the inverse total number of sites. (a) Triangular lattice. (b) Square lattice.

FIG. 4 :

 4 FIG.4:(color online) Phase separation boundaries for the square and triangular lattices in the thermodynamic limit. The dashed lines correspond to the approximate location of the phase transition between plaquette and columnar phases[START_REF] Syljuasen | [END_REF] for the square lattice and between plaquette and RVB phases[9,13] for the triangular lattice.

FIG. 6 :

 6 FIG. 6: (color online) Energy spectrum vs (reduced) Aharonov-Bohm flux ξ, for a 4 × 4 cluster with 4 holes at v = 0.70 and t = 1.00. (a) Torus geometry. The momenta of the two low-energy branches are specified and the complete spectrum is displayed up to ξ = 0.2. (b) The two low-energy branches for the case of a cylinder (dashed line) and including a small bond disorder (full red line).

  2, typical results are shown for both square and triangular lattices and for different sizes. Interestingly, PS appears in both cases, but with

	-0.02
	(x)			3x6x6 3x8x8
	s t			3x10x10
	-0.025
	0	0.05	0.1	0.15
	-0.134
	-0.136 s s (x)			12x12 16x16 14x14
	0 -0.138	0.05	0.1	0.15

  . In agreement with

	0.12
	0.09	v=0.90
	0.06	v=0.85 v=0.80
	0.03
	0	0	0.002	0.004	0.006	0.008	0.01
	0.1
	0.08				v=0.95
						v=0.90
	0.06				v=0.85
	0.04
	0.003	0.004	0.005	0.006	0.007	0.008

We acknowledge useful discussions with Federico Becca. This work was supported by the Swiss National Fund, by MaNEP, and by the Agence Nationale de la Recherche (France).