Mathieu Clabaut

Ning Ge

Nicolas Breton

Eric Jenn

⋆⋆ Rémi Delmas
email: remi.delmas@onera.fr

Yoann Fonteneau

Industrial Grade Model Checking Use Cases, Constraints, Tools and Applications

Model checking has made a lot of progress since its infancy. For a long time, industrial applications were still limited to some very specific domains out of which the technique bumps into the state explosion wall. Nowadays things evolve and some tools are able to tackle real world use cases outside of the known domains. We give here the feedback collected when using model checking on several industrial strength use cases and give indication on how we take into account the specific domain constraints.

Model Checking for Industrial Problems

Model checking refers to the problem of exhaustively and automatically checking whether a given model of a system meets a given specification.

Model Checking is now an old technique which takes its ground in the mid 1970s as a response to concurrent problem analysis. It was until recently essentially confined to some specific areas, such as hardware analysis or protocol verification. Extension to other domains such as software verification has always been difficult due to the combinatorial explosion problem (the size of the space state grows exponentially with the size of the problem to analyze). However, recent developments in a variety of fields, ranging from symbolic model checking to SAT solver engines and including model checker parallelization lead to a broader range of application in industry including software analysis.

Industrial Use Cases

Model checking may be used for the following use cases, depending on the tool abilities: Safety Proof consists in verifying some properties on a system model. The model may be designed by hand or automatically derived from existing artefacts, such as Ada, C, Simulink or Scade code… If one of the properties is not being verified, the tool provides a counterexample explaining why the property does not hold. System Debugging. One of the great features of model checking is its ability to provide the user with counterexamples. Such a feature may be used in different ways, but is generally of a good help for debugging a system. For example, when designing a complex system which should ensure a safety property (e.g. a non collision property for an automatic train system), it is very useful to debug the root concepts which should ensure the safety upon a system model on which you can analyse counterexamples. Equivalence Checking. The ability to prove some properties upon an existing model ⋆ Seconded from Systerel, Toulouse, France ⋆⋆ Seconded from Thales Avionics, Toulouse, France may be used as a way to prove the equivalence of two models. One application of equivalence checking is to verify that a software design in C satisfies its specification, by proving the equivalence between a model of the design and a manually written model of the specification. Another application is to verify the correctness of a tool transformation -for example a translation from one formalism to another -by proving the equivalence between the before and after transformation artefacts. Constraint Solving. A somewhat more contrived use of model checking is the use of its ability to provide counterexamples as a way to solve a constraint problem. The way to do it is to ask the tool to verify that the problem has no solution. If one exists, the tool will return a counterexample which is one solution to the problem. Test Case Generation. Another use case, derived from the preceding one, is to use a test objective as a constraint. The provided counterexamples then give inputs and associated oracles which fulfill these test objectives and may be used as test scenarios.

Industrial Constraints

Putting aside the many use cases listed before, industry faces many constraints which may have prevented the use of model checking until now:

Regulatory Constraints Many companies must obtain approval from a suitable authority that the system they develop is acceptably safe to operate with regards to the applicable assurance standards (CEN-ELEC EN-50128, DO-178C,…). As such, it must be shown that the tools who have contributed to the system or to its verification have been qualified with respect to their usage and to their contribution to the global safety. Such a qualification bears a lot of constraints upon the tools and their development process

The S3 Toolbox

S3

Systerel Smart Solver4 (S3) is Systerel's response to the aforementioned use cases and industrial constraints. S3 is constructed around:

-a high level synchronous modelling language, -several frontends for C, Ada, Scade, … -a translation tool chain from the high level language towards a bit level language, -a proof engine working over the bit level language, -a proof verification engine, -several tools for animation and debugging.

The performance of the proof engine allows us and our clients to manage the proof of industrial size problems some of them we will mention in the next section. The size of those models routinely topped ten millions variables and several hundred millions of clauses.

The qualification of S3 is made possible by the use of several diversified tool-chain with some small key tools built with respect to the higher integrity constraints that may be required -namely, an equivalence model constructor, and a tool to verify the validity of the proof (see fig. 1 on the following page).

Floating-Point Arithmetic Library in S3

Critical applications used to use fixed-point arithmetic that requires less memory and less processor time then floating-point arithmetic (FPA) to perform non-integer computations on executing processors with no Floating-Point Unit (FPU), while leading to a limited-precision. Floating-point numbers support a trade-off between range and precision due to its formulaic representation which approximates a real number. Another advantage stands in the existence of a standard [1] based on solid mathematical grounds. Nowadays, FPA is more and more used in the space, aeronautics and automotive industries, due to the increasing complexity of the computations and because FPUs are becoming standard for most processors. Floatingpoint numbers are not real numbers. Floatingpoint operations behave in quite different way from the real counterparts, due, for instance, to rounding and cancellations [17]. Consequently, a software implementation of some mathematical expression usually provides results that are not strictly, mathematically, exact. As it is of-ten difficult to foresee the behavior of floatingpoint programs, formal verification of floatingpoint programs is highly desired in the industry.

The basic approaches to address formal verification of floating-point programs include abstract interpretation, formal proof and bitblasting (also called bit-flattening). Abstract interpretation partially executes program on an abstract domain. This approach performs well in program analysis with floating-point variables [START_REF] Blanchet | A static analyzer for large safety-critical software[END_REF]. Formal proof supported by proof assistants is a very powerful approach that requires to be guided by highly skilled expert to direct the reasoning towards target properties. Interactive theorem provers such as ACL2, Coq, HOL Light and PVS have been applied to floating-point verification [14]. Both of abstract interpretation and formal proof approaches lack ability to generate counterexamples when the property does not hold. Bit-blasting represents floating-point operations as circuits, which are then transformed to Boolean formula with bitwise operators to be solved by SAT solvers. Bitblasting relying on SAT solvers is a fully automatic reasoning benefiting from counterexamples for floating-point programs. It is implemented in the Satisfiability Modulo Theories (SMT) solvers such as Z3 [START_REF] De | Z3: An efficient smt solver[END_REF], MathSAT 5 [START_REF] Cimatti | The mathsat5 smt solver[END_REF], SONOLAR [15] and CBMC [START_REF] Brillout | Mixed abstractions for floating-point arithmetic[END_REF]. It is also the elementary but the most significant part of other floating-point verification strategies using SAT solvers as the back-end. Since the publication of SMT-LIB theory of binary FPA [18], solvers are starting to support it using some advanced QF_FP solving strategies, such as mixed abstraction in CBMC [START_REF] Brillout | Mixed abstractions for floating-point arithmetic[END_REF], non-conservation approximations in Z3 [13], abstraction into interval arithmetic in MathSAT [START_REF] Brain | Interpolation-based verification of floatingpoint programs with abstract cdcl[END_REF][START_REF] Brain | Deciding floating-point logic with abstract conflict driven clause learning[END_REF], translation into nonlinear reals in Realizer [16], etc.

S3 will be used to verify part of the floatingpoint software embedded in a rover platform, TwIRTee (a three-wheeled autonomous rover) used as the demonstrator for the INGEQUIP project. The INGEQUIP research project at the Institut de Recherche Technologique (IRT) Saint-Exupéry in Toulouse addresses the following equipment engineering topics: architecture modelling and evaluation, early V&V using vir-tual platforms, and formal verification. For this purpose, we implemented an FPA standalone library to enable the floating-point verification by means of bit-blasting. This optimized FPA library will establish a solid foundation and basic strategy for our future investigation on advanced FPA strategies in S3.

The

Industrial Applications

Railway Use Case -Interlocking Safety Proof

An interlocking is an arrangement of signal apparatuses that prevents conflicting movements through an arrangement of tracks such as junctions or crossings. An interlocking is designed so that it is impossible to display a signal to proceed unless the route to be used is proven safe.

The main challenge is to ensure that whatever scenario will happen, the designed interlocking will stay safe. And albeit its apparent simplicity one quickly understands the real underlying complexity -whatever rule you may imagine, it looks like one can always find a scenario breaking it (train can go backwards, can sometime appear to fly due to concurrency artefacts, …)

Instanciation Design Process In a recent work, we have considered a Computer Based Interlocking (CBI) designed through an instantiation process. In such a process, the signaling principles are captured by the engineers to produce the Generic Design in some suitable language. This design contains a set of generic descriptions of code parts that an interlocking system shall execute for some specific object of the system such as signals, switches, routes,… Each of these generic descriptions is given in a parametrized form which allows for the specificities of the object to which it applies. For example, the generic design for a route will probably be parametrized with the list of switches contained in the route. The generic design is thus specific to a set of signaling principles, but independent of a particular track layout.

The generic design is then instantiated upon a particular track layout to produce the running software.

Verification Process We graft our verification process onto the instantiation scheme by designing a Generic Safety Specification made of:

-some high level proof obligations (3 to 4 properties about absence of derailments and collisions), -some intermediate predicates modeling the domain (such as topological predicates which encode the track layout, integration predicates which associate input and output variables with their object instance, helper predicates conveying high level relations between objects such as reachability,…) -an environment model upon which the proof obligations are expressed (trains behavior,…).

and then instantiate those properties and models with the specific track layout data. The overall uncertified process is given in figure 2. As an example of intermediate predicates used in the Generic Safety Specification, the fact that "every block belonging to a given route" is free can be expressed as:

ALL b:BLOCKS (contains(route,b) -> free(b))
where ALL is the universal quantifier, b:BLOCKS indicates that b shall range on every instances of the BLOCKS, and -> is the implication operator.

Results

The proof process was applied to a CBI of more than 200 routes and permit us to pinpoint 3 safety counterexamples which where then corrected whether in the design or in the topology data. The overall proof took less than 25 mn in the worst case, which would go as high as 1 hour for a certified process implying the use of a second tool chain, an equivalence verification and a prooflog check.

The verification was successful and proved to be usable enough to assist the design team in debugging the CBI. It also paved the way for a certified verification that would be used by the safety team.

In the competitive interlocking market, this would clearly give an edge to the CBI manufacturer when comparing to other solutions integrating or not formal verification.

Aerospace and Automotive Use Case

The floating-point verification by S3 has been applied to two case studies: the avionic triplex sensor voter and the automatic rover protection system. It is also used to automatically generate the test cases.

Triplex Sensor Voter Case Study

The triplex sensor voter5 is used in a common form of redundant aircraft system Triplex Modular Redundancy (TMR), which relies on three identical sensors to compute an output value from the three input values by the voter. It is implemented using linear arithmetic operations as well as conditional expressions (such as saturation). Its formal analysis covers functional and non-functional properties including stability, absence of runtime errors, and also to parameterize certain parts of the model to help the formal analysis. The formal analysis of triplex sensor voter was first studied by Dajani-Brown et al. in [START_REF] Dajani-Brown | Formal modeling and analysis of an avionics triplex sensor voter[END_REF], where real values were abstracted by integer values and integrators were not used.

In [START_REF] Dierkes | Formal analysis of a triplex sensor voter in an industrial context[END_REF], Dierkes analyzed the Simulink model with real numbers by both simulation and formal verification, then estimated the impact of rounding errors caused by the floating-point implementation using SMT solvers and abstract interpretation. In [START_REF] Champion | Generating property-directed potential invariants by quantifier elimination in a kinduction-based framework[END_REF], Champion et al. strengthened the stability property by generating lemmas using a property-directed approach.

In our work, we start from a SCADE model of the voter and translate it to HLL model using the SCADE-translator. Thanks to our FPA library, the HLL model is then verified using the S3 solver. In parallel, we use SMT-solvers Z3 v4.4, MathSAT 5 and SONOLAR to verify the SMT model. Experiments are carried out using both simple and double precision floating-point numbers with or without subnormal numbers and with different rounding modes.

The results show that neither Z3 v4.4 (bit-blasting strategy, floating-point strategy) nor MathSAT 5 (bit-blasting strategy, abstract CDCL algorithm) or SONOLAR are able to handle the step instance, be it in simple or double precision. We managed to prove the inductive instance using a combination of SONOLAR bit-blasting to a CNF and a pure SAT solver (glucose 4.0 multithread with 8 threads and aggressive restart strategies, satellite preprocessing) in 10min of computation for the simple precision instance, and 4h15min of computation for the double precision instance (wall clock time). S3 proved the step instance in 6min using glucose 4.0 and 5mins using S3's own solver for the simple precision instance, and in 9h32min using glucose 4.0 for the double precision instance.

Test Generation

In order to conform to the domain standards, companies sometimes have to setup processes where a structural test coverage must be achieved by mean of manual testing activities. Those activities consist mainly in finding the inputs and oracles such that the structural coverage criterion is respected.

We successfully used model checking with S3 in order to setup an automatic process of finding the set of tests that will achieve the structural node coverage on a Scade model, where the criterion was that at most one entry of a given node should change at once. Automation of such a task leads to great cost reductions.

Formal Verification for ARP Case Study

The Automatic Rover Protection (ARP) system, a simple collision avoidance function, is developed for the twIRTee platform in the INGEQUIP project. It performs a predefined trajectory (a "mission") on a predefined track and avoid collisions with other rovers. The set of tracks are statically defined and embedded in the rover. The ARP is based on the following principles:

1. Missions (trajectories) are precomputed using an adapted shortest path algorithm[12]. 2. A rover shall only move on a reserved path if there is a free reserved space ahead. 3. A reserved path is a stack of reserved nodes. 4. A rover resends reservation request to get the exclusive access to the tracks located ahead of it, if there is not enough reserved space (D_REQ) ahead. 5. A rover stops if there is not enough reserved space (D_STOP) ahead. 6. A rover reserves enough space (D_RSV) ahead for each request. 7. In order to ensure a consistent management of node reservations, a distribution strategy using id-priority is implemented.

In this case, S3 is exploited to verify several properties applicable to all independent rovers, such as P 1 : rovers are never granted simultaneous access to the same area (safety), and P 2 : all rovers eventually reach their destinations (bounded liveness), etc. We present experimental results from a ARP case study with two rovers. The trajectories of rovers within an enclosed space of 10m×10m are predefined (see fig. 3). The rover R 1 proceeds with a speed v 1 = 0.4m/s from node A, through node E, F, G and H, then stops at node B, while the rover R 2 proceeds with a speed v 2 = 0.3m/s from node C, through nodes G, F, E and H, then stops at node D. The length of each trajectory is about 20m. As the first step of our study, we assume that the space occupied by rovers are ignored, and that the clocks of rovers are synchronized. Consequently, rovers are considered to perform actions synchronously. We proved the property P 1 in almost no time using induction and the property P 2 in several seconds using bounded model checking in S3.

Other Use Case -Constraint

Solving for Matrix Wiring S3 was used as a constraint solver to automate finding of solution to a wiring problem over a 100 × 100 matrix with several layers of wires and several limitations over the deformations that wires may undergo.

Once the counterexample obtained, a small and easy tool allowed to translate it into a netlist.

Conclusion and Future Work

Our experience with model checking has shown that S3 is a great tool to manage different types of industrial size problems with respect to their specific regulatory constraints. However, the tool is still at pains for some peculiar problems and has some intrinsic limitations that prevent its application on a wider range of problems, one of it being its current inability to handle designs which make use of floating-point arithmetic.

In order to tackle the floating-point limitation, we have designed an FPA library and integrated it in S3. Our next goal is to investigate how well it works on a wider range of industrial designs from automotive, aeronautic, industry, energy… and to establish benchmark (or reuse existing SMT benchmarks) to evaluate the performance of floating-point verification by S3.

The proof engine at the core of S3 is also undergoing heavy work in order to improve its efficiency for big size models. This significant improvement, associated to the heavy work ongoing on the proof engine at the core of S3, shows great promises for dealing with future industrial challenges.

Fig. 1 .

 1 Fig. 1. Example of application for proving equivalence between C sources and a HL model specification. The two tool-chains are diversified and the software in the circles are developed with the higher level of integrity (HL = High Level, BL = Bit Level).

Fig. 2 .

 2 Fig. 2.Overview of an uncertified verification process for an instantiated design. The specific topology translator tool translates topological data in an High Level Language model. The specific design translator tool translates the instantiated design in a High Level Language model. The concatenation of all models is then feed to S3 proof engine. Manual lemmas may be added into the Generic Safety Specification to help the proof.

Fig. 3 .

 3 Fig. 3. Rover Trajectories for ARP

  implementation of FPA library in S3 is based on the IEEE FPA standard 754-2008 [1]. The FPA library in S3 includes:

	-Binary interchange format encoding that
	allows user-defined ranges of exponent and
	mantissa, including single and double pre-
	cisions
	-Normal, subnormal, infinity, NaN numbers
	-5 rounding directions: roundTiesToEven,
	roundTiesToAway, roundTowardPositive,
	roundTowardNegative, roundTowardZero
	-Comparison operations: Equal, NotEqual,
	Greater, GreaterEqual, Less and
	LessEqual
	-Arithmetic operations: Addition, Sub-
	traction, Multiplication, Division and
	SquareRoot
	-Conversion operations: convertIntToFloat
	and convertFloatToInt
	-Trigonometric operations: Sin, Cos, and
	Tan
	-Default exception handling: invalid opera-
	tions, division by zero, overflow and under-
	flow
	The SquareRoot and trigonometric opera-
	tions are implemented by means of both in-
	terpolation table and the function proposed by
	Cody and Waite [8].

  Methods for Industrial Critical Systems, pages 102-116. Springer, 2011. 12. Edsger W Dijkstra. A note on two problems in connexion with graphs. Numerische mathematik, 1(1):269-271, 1959. 13. Andreas Fröhlich, Armin Biere, Christoph M Wintersteiger, and Youssef Hamadi. Stochastic local search for satisfiability modulo theories. In Twenty-Ninth AAAI Conference on Artificial Intelligence, 2015. 14. John Harrison. Floating-point verification using theorem proving. In Formal Methods for Hardware Verification, pages 211-242. Springer, 2006. 15. F Lapschies, J Peleska, E Gorbachuk, and T Mangels. sonolar smt-solver. Satisfiability modulo theories competition; system description, 2012. 16. Miriam Leeser, Saoni Mukherjee, Jaideep Ramachandran, and Thomas Wahl. Make it real: Effective floating-point reasoning via exact arithmetic. In Design, Automation and Test in Europe Conference and Exhibition (DATE), 2014, pages 1-4. IEEE, 2014. 17. David Monniaux. The pitfalls of verifying floating-point computations. ACM Transactions on Programming Languages and Systems (TOPLAS), 30(3):12, 2008. 18. Philipp Rümmer and Thomas Wahl. An smt-lib theory of binary floating-point arithmetic. In International Workshop on Satisfiability Modulo Theories (SMT), 2010.

S3 is maintained, developed and distributed by Systerel.

Triplex sensor voter case study is provided by Rockwell Collins to make it publicly available to the research community.


domain: math.math-pr info math

text:
L A Abbas-Turki

S Graillat

Resolution of a large number of small random symmetric linear systems in single precision arithmetic on GPUs

Keywords: GPU, LDLt, Householder reduction, parallel cyclic reduction, divide and conquer for tridiagonal eigenproblems

This paper focuses on the resolution of a large number of small symmetric linear systems and its parallel implementation on single precision on GPUs. The computations involved by each linear system are independent from the others and the number of unknowns does not exceed 64. For this purpose, we present the adaptation to our context of largely used methods that include: LDLt, Householder reduction to a tridiagonal matrix, parallel cyclic reduction that is not a power of two and the divide and conquer algorithm for tridiagonal eigenproblems. We not only detail the implementation and optimization of each method but we also compare the sustainability of each solution and its performance which include both parallel complexity and cache memory occupation. In the context of solving a large number of small random linear systems on GPU with no information about their conditioning, we show that the best strategy seems to be the use of Householder tridiagonalization + PCR followed if necessary by a divide & conquer diagonalization.

Introdution

A quite few number of problems in physics can be divided into subproblems, solved locally and process communication steps to recover the global solution. This is also the case for simulations in mathematical finance, in particular for the challenging problem of Credit Valuation Adjustment (CVA). When American contracts are involved, the CVA can be simulated thanks to a nested Monte Carlo (NMC) that performs a Dynamic Programming Algorithm (DPA) on the inner trajectories. This procedure constitutes a straight extension of the square Monte Carlo presented in [START_REF] Abbas-Turki | Toward a coherent Monte Carlo simulation of CVA[END_REF] as a benchmark method for CVA on European contracts. Basically, the main ingredient of this extension relies on the resolution of a large number of small random symmetric linear systems.

The CVA simulation is rather the origin and not the purpose of this paper. In this contribution, we are really focused on the computational effort of adapting to our context some well known algorithms. Indeed, we are not aware of any work that deals with the resolution on GPUs of a large number of symmetric small linear systems whose size n = (number of knowns) does not exceed 64. The lack of research done in this direction could be due to the fact that large linear systems are more difficult to parallelize than smaller ones. Nevertheless, a good parallelization for small linear systems cannot be considered as straight simplifications of the work done for large systems.

Actually, when targeting performance, we have to be aware that the paradigm of the parallel implementation changes according to the size. When n ≥ 8, we should be also aware that a SIMD parallelization in our context is far from optimality. Indeed, because of the size of the cached memory available per block of threads, associating one thread per linear system reduces significantly the number of threads that can be launched in parallel. Consequently, the parallelization that we propose is neither SIMD nor straight simplification of known libraries developed for large linear systems like MAGMA [29].

The solution that we propose is really adapted for the specific problem of solving large number of small symmetric linear systems. It is based on the dependence classification of threads: The class of threads that are independent because involved in different linear systems, noted T I , and the class of those that communicate because involved in the same linear system, noted T C . The biggest part of this paper is dedicated to the organization of communicating threads T C and their use of CUDA shared memory. Regarding the independent groups of threads T I , their number will be chosen in order to saturate the use of the shared memory size available per block or at least to have a sufficient work per a streaming multiprocessors (SMs).

Our contribution can be summarized in the following points:

• We give the CUDA source code [30] associated to the adaptation of to each algorithm: LDLt, Householder reduction, parallel cyclic reduction that is not necessary a power of two and divide and conquer for eigenproblem.

As one could expect, the adaptation of the divide and conquer was the trickiest since it requires some technicality due to: Choosing the right indices for the division part, implementing the right algorithm for the solution of the secular equation, tuning the deflation parameters to get sufficiently accurate results.

• We provide an in-depth description of each implementation.

• We compare the execution time of the different methods mentioned above.

• We propose an original method to further optimize the adaptation of LDLt to our context.

• We provide an original parallel cyclic reduction that can be used for any vector size and not only a power of two that requires a zero padding.

Although we are interested by small linear systems with n ∼ 32, it is not reasonable to use LDLt for all situations. In fact, in Section 2.2, we show that even when n = 30 some random linear systems turn out to be ill-conditioned for single floating-point precision as the one used on our Geforce GPU. For this reason, we found ourselves obliged to develop both the Householder tridiagonal reduction as well as the divide & conquer diagonalization of tridiagonal matrices. Subsequently, one could ask the following legitimate question:

Must we systematically use Householder tridiagonalization with divide & conquer when we suspect the random linear systems to be ill-conditioned?

Our answer is: Perform Householder tridiagonalization and solve the linear systems cheaply using parallel cyclic reduction then take a decision according to the value of the residue error. If the residue error is small then we already have good solutions. Otherwise, we must perform divide & conquer diagonalizations and discard the smallest eigenvalues. The next time we solve this same kind of linear systems: If they used to be well-conditioned then we just process LDLt, otherwise we execute directly the combination of Householder tridiagonalization and divide & conquer diagonalization.

The answer above justifies the work detailed in this paper that is arranged as follows. In Section 2, we give a brief description of NMC for CVA and we show a realistic example where the linear systems are ill-conditioned. Afterwards, the presentation of each resolution algorithm is explained in a separate section: Section 3 for LDLt, Section 4 for Householder and parallel cyclic reductions and Section 5 for divide & conquer diagonalization. All sections 2, 3, 4 and 5 start with a subsection that describes the headlines of each problem and some references related to it. Section 6 concludes this paper with global remarks and the future work that is in preparation.

Brief description of NMC for CVA

Presentation of CVA and the references

The 2007 economic crisis raises the fear of systemic stability when the default of one financial institution could be the origin of a cascade of other defaults. To reduce this risk, several measures were established that include the calculation of the CVA as an important part of the Basel III prudential rules. Since the paper [START_REF] Brigo | Counterparty Risk and Contingent CDS under correlation between interest-rates and default[END_REF], the CVA can be viewed as an insurance contract that compensates for the no-recovered sum by the counterparty when it defaults.

For a comprehensive financial presentation of CVA, we refer to [START_REF] Brigo | Counterparty Credit Risk, Collateral and Funding: With Pricing Cases For All Asset Classes[END_REF]. Regarding the mathematical aspects, the reference [START_REF] Crépey | Counterparty Risk and Funding. A Tale of Two Puzzles[END_REF] can be viewed as the most recent and complete summary on the subject. However, little research has been dedicated to the development of numerical procedures that can be used to perform the computations. Book [START_REF] Cesari | Modelling, Pricing and Hedging Counterparty Credit Exposure[END_REF] is one of the first references that presents the industry practices in computing CVA. Among research papers, maybe the most devoted to computing CVA are [START_REF] Abbas-Turki | Toward a coherent Monte Carlo simulation of CVA[END_REF], [START_REF] Fujii | Perturbative expansion technique for nonlinear FBSDEs with interacting particle method[END_REF] and [START_REF] Henry-Labordère | Cutting CVA's complexity[END_REF].

Due to both mathematical and computational complexity, none of the previous references provide a benchmark procedure to deal with CVA when American options are involved. This is despite the fact that American contracts are widely exchanged, especially in some markets. As it will be shown, simulating CVA on American options can be overcome thanks to an efficient way of parallelizing the resolution of a large number of small symmetric systems on GPUs. This latter point is really the heart of this work. Before detailing the background NMC algorithm in the next subsection, we introduce below the simplest formulation of the CVA

CVA 0,T = (1 -R)E P + τ 1 t<τ ≤T , (1)
where R is the recovery made by the counterparty when it defaults, E denotes the expectation operator, P t is the process of the value exposure to the counterparty, τ is the random default time of the counterparty, T is the protection time horizon and the positive part function is denoted by + . As already explained in [START_REF] Abbas-Turki | Toward a coherent Monte Carlo simulation of CVA[END_REF], one of the most important challenges of the CVA comes from the fact that the exposure P t is generally the price of a basket of different contracts that are written with our counterparty. If these contracts can be priced by closed expressions, the CVA can be calculated thanks to a onestage simulation using either the discretization of a partial differential equation or Monte Carlo method as in [START_REF] Crépey | A Lévy HJM multiple-curve model with application to CVA computation[END_REF]. However, when the underlying contracts must be simulated, it is natural from expression (1) to perform a two-stage simulation: The inner stage to compute P t and the outer stage for the CVA. This two-stage simulation leads to the square Monte Carlo simulation used in [START_REF] Abbas-Turki | Toward a coherent Monte Carlo simulation of CVA[END_REF] as a benchmark algorithm.

The square Monte Carlo proposed in [START_REF] Abbas-Turki | Toward a coherent Monte Carlo simulation of CVA[END_REF] is developed for an exposure P t of European contracts which is a category of derivatives that can be exercised only at maturity (< T ). Unlike the European case, the American derivatives allow an early exercise (< T ) of the contract which can be solved using DPA. The most used implementation of this DPA is the Longstaff-Schwartz algorithm proposed in [START_REF] Longstaff | Valuing American options by simulation: A simple least-squares approach[END_REF] and theoretically studied in [START_REF] Clément | An analysis of a least squares regression algorithm for American option pricing[END_REF]. For CVA simulation, this algorithm is implemented thanks to regressions performed on the inner trajectories providing a new NMC algorithm that extends the square Monte Carlo of [START_REF] Abbas-Turki | Toward a coherent Monte Carlo simulation of CVA[END_REF].

Description of the matrices involved by our new NMC algorithm

As usual, expression (1) is discretized using a fixed number of time steps N that introduces the sequel 0 = t 0 < t 1 < ... < t N = T and the estimation

CVA 0,T = N -1 k=0 E P + t k+1 1 τ ∈(t k ,t k+1 ] . (2)
In order to approximate the expectation E in (2), we simulate M 0 outer stage trajectories of the underlying asset S = (S established. In order to compute the exposition P in (2) at each time t k of the outer trajectories, we simulate M k inner stage trajectories of the same underlying asset S. In these outer and inner simulations, the vector S is a given Markov process whose realisations could be drawn thanks to a given random number generator like those presented in [START_REF] Abbas-Turki | Pricing derivatives on graphics processing units using Monte Carlo simulation[END_REF]. An illustration of this NMC algorithm is given in Figure 1.

When the exposure P includes American contracts, we perform the Longstaff-Schwartz algorithm on the inner trajectories. This requires N -k -1 regressions at each time step k ∈ {1, ..., N -1} and for each outer trajectory l ∈ {1, ..., M 0 }. For a fixed couple (k, l), the regression is performed using a projection on the space generated by ψ l (S t k ) = (ψ l 1 (S t k ), ..., ψ l n (S t k )). The choice of the latter family should obviously depend on the considered problem, but generally practitioners use some family of polynomials.

The regression matrix is the Monte Carlo approximation

A k,l of A k,l = E(ψ l (S t k )ψ l (S t k ) t ) given by A k,l = 1 M k M k j=1 ψ l (S (j) t k )ψ l (S (j) t k ) t (3)
where (j) is the inner trajectory index and t is the transpose operator. By definition, the correlation matrices A k,l are symmetric. Moreover, the family ψ(S) is always chosen to make all {A k,l } 1≤k≤N -1, 1≤l≤M0 positive definite i.e.

∀ A ∈ {A k,l , 1 ≤ k ≤ N -1, 1 ≤ l ≤ M 0 } X t AX > 0, for every X ∈ R n -{(0, ..., 0)}. (4)
However, the Monte Carlo approximations A k,l , that are symmetric, do not necessary fulfill condition (4) since they depend on the convergence parameter M k . Indeed, although the values taken by M k can be sufficient to have a good overall convergence of NMC, some of the small values produce either numerical indefiniteness or even negativity. As already introduced in [START_REF] Gordy | Nested Simulation in Portfolio Risk Measurement[END_REF] and adapted to CVA in [START_REF] Abbas-Turki | TVA on American Derivatives[END_REF], M k should be of the order of √ M 0 . In Figure 2, we give some condition values for the benchmark model of Black & Scholes (d = 29) with independent coordinates. The parameters of this model are its volatility σ, interest rate µ and spot value S 0 . The regression performed in this figure are linear and includes the constant ψ(S) = (1, S 1 , ..., S 29 ) which makes n = 30.

Although the value n = 30 could be considered high by some practitioners, it is possible to use it especially for sufficiently large expositions to the counterparty, for instance the exposition of a bank to another bank. Figure 2 shows then an example of matrices that can corrupt the DPA when implemented in a single precision. In fact, the number of trusted decimals is not sufficient to make a decision on the early exercise strategy computed by the DPA. When this kind of situation is confronted, one has to discard some smallest eigenvalues before resolving any linear system.

LDLt decomposition

Presentation of the algorithm and the references

The resolution of the linear system AX = Y with A symmetric is divided into two steps: The factorization of the matrix A that leads to A = LDL t and the resolution of LZ = Y as well as DL t X = Z where L t is the transpose of L. The matrix D is diagonal and the matrix L is lower triangular with 1 on its diagonal. The factorization is performed thanks to the following expressions

A = LDL t , D j,j = A j,j - j-1 k=1 L 2 j,k D k,k , L i,j = 1 D j,j A i,j - j-1 k=1 L i,k L j,k D k,k if i > j.
(

) 5
Because it prevents the computation of square roots, LDLt is generally considered as a better alternative to the Cholesky decomposition. Both methods share the same important stability for symmetric positive definite matrices A characterized by [START_REF] Ballard | Communication-Optimal Parallel and Sequential Cholesky Decomposition[END_REF]. Furthermore, when the positive definiteness is numerically questioned, one should avoid the use of either LDLt or Cholesky decomposition. Also, these two methods share the same complexity order and the same memory space occupation. Due to all these similarities and the fact that Cholesky's literature is larger than LDLt's, we will not distinguish between the references of each method.

To our knowledge, [START_REF] Volkov | LU, QR and Cholesky Factorizations using Vector Capabilities of GPUs[END_REF] is the first reference that implements Cholesky decomposition on GPUs. Paper [START_REF] Ballard | Communication-Optimal Parallel and Sequential Cholesky Decomposition[END_REF] comes after and it theorizes the minimization of the communication cost involved in Cholesky factorization. Even though both papers are very interesting, the extent of the work developed their is adapted to large matrices n ≥ 64. The same can be said on the Cholesky's code of MAGMA [29]. Indeed, MAGMA library is even dedicated to heterogeneous CPU/GPU implementations that are generally justified for sufficiently large sizes.

As said previously, the stability property of LDLt and Cholesky is quite important. In fact, when the correlation matrix is not numerically singular, these two methods are so stable that they do not need any pivoting like those performed for LU decomposition [START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF]. Escaping the pivoting phases makes a great advantage for the GPU implementation since it reduces communications between threads. Besides, LDLt and Cholesky are the most efficient methods of factorization with a complexity given by O(n 3 /6) where n × n is the size of the matrix. They are also the ones that use the least the memory space as they involve only n(n + 1)/2 values.

Once the LDLt factorization performed, the resolution of LZ = Y then of DL t X = Z are quite straightforward. These resolutions are even quadratic in complexity with respect to n.

Adaptation and optimization

We present three different versions of the LDLt factorization:

1. An SIMD version that requires only threads of T I , one for each linear system. 2. A collaborative version that involves n threads of T C for each linear system with n unknowns. 3. An optimal hybrid solution that involves n * (n * < n) threads of T C for each linear system with n unknowns. The number of the T I -independent groups of threads is chosen in order to saturate the shared memory or at least to have a sufficient work per an SMs. The SIMD version is straightforward and used only to convince the reader of its inefficiency. Regarding the collaborative and the hybrid versions, they are both based on a column after column processing. In fact, as shown on Figure 3, for a fixed value of j, the different coefficients {L i,j } j+1≤i≤n can be computed by at most n -j independent threads. Thus, {L i,1 } 2≤i≤n involves the largest number of possible independent threads equal to n -1. In the collaborative version, we use the maximum n -1 threads +1 additional thread that is involved in the copy from global to shared and in the solution of the linear system after factorization. This makes n threads for the collaborative version and one of these threads is also involved in the computation D j,j which needs a synchronization before calculating L i,j . For j > n/2 in the collaborative version, more than the half number of threads are in a wait state. This is not a problem when n is large enough, because the shared memory is sufficiently filled which limits the possibility of launching independent threads T I on other linear systems. Nevertheless, when n is small, the communicating threads T C in a wait state prevent the schedular from the execution of independent threads T I even though there is a sufficient shared memory space for other linear systems. This situation motivates the use of an hybrid solution that either employs about the half number n * (n) ≃ n/2 of maximum number of communicating threads or all of it n * (n) ≃ n. As we will see in Figure 4 for some n, this hybrid solution is significantly better than the simple collaborative one.

Finally, we precise that we use n(n + 1)/2 memory space for the LDLt factorization +n for the resolution. In order to keep a coalesced access to the shared memory and reduce the bank conflicts, the matrices are also stored column after column and not row after row as it is generally done.

Comparison of the different versions

Ξ ρ =          1 ρ • • • ρ ρ ρ 1 ρ . . . ρ . . . . . . . . . ρ . . . ρ 1 ρ ρ ρ • • • ρ 1          with 0 < ρ < 1 (6)
are positive definite. Because these matrices are strongly structured, we prefer to use a randomized version given by A = Q Ξρ RQ ′ Ξρ , where R is a diagonally dominant tridiagonal symmetric random matrix and Q Ξρ is the orthogonal matrix that results from the Householder tridiagonalization of Ξ ρ . The components of R are set using uniform random variables and the multiplication of the diagonal elements by the appropriate factor to make R diagonally dominant.

We point out that all the results are obtained from an implementation on an NVIDIA Geforce 970.

From Figure 4, we notice that the hybrid solution outperforms the collaborative one when n < 40. Moreover, the SIMD version is clearly unsatisfactory for all sizes even when n = 4. The speedup obtained when using communicating threads gets relatively high according to the size n and it exceeds the speedup of 14 per linear system when n = 64.

The optimal number of communicating threads in the hybrid version depends on the GPU used. In Figure 5, we give the experimental values obtained for different sizes n when the implementation is performed on the Geforce 970. We distinguish two regimes,

n * (n) ≃ n/2 if n < 40, n otherwise. ( 7
)
Using the LDLt hybrid implementation, we establish Figure 6 that shows the number of linear systems that can be solved per second. The curve obtained is almost proportional to n -3 which coincides with the theoretical result.

Householder and parallel cyclic reductions

Presentation of the algorithms and the references

Householder tridiagonalization

Similar to the method based on LDLt, we propose here to solve the linear system AX = Y , with A symmetric, through two steps:

• The tridiagonal Householder decomposition A = QU Q t where Q is orthogonal and U is symmetric tridiagonal.

• The Parallel Cyclic Reduction (PCR) associated to the problem U Z = Q t Y that allows to recover X thanks to X = QZ.

When the linear system is symmetric, the Householder tridiagonalization is generally used as the first step of a diagonalization algorithm which could employ: QR factorization, bisection method, multiple relatively robust representations or divide & conquer. We refer to [START_REF] Demmel | Performance and Accuracy of LAPACK's Symmetric Tridiagonal Eigensolvers[END_REF] for a sequential comparison between these four algorithms according to speed and accuracy.

As advised in the introduction, we would like to use the Householder tridiagonalization with PCR and check the residue error before looking for a diagonalization of the system. This procedure is justified by the fact that PCR is less complex than any of the four algorithms cited above with a theoretical ratio equal at least to n/ log 2 (n) in favour of PCR. Moreover, as already shown in [START_REF] Zhang | Fast Tridiagonal Solvers on the GPU[END_REF], PCR is quite stable for symmetric and positive definite matrices and is suited to parallel architectures like GPUs.

Without going through details that can be found for instance in [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF][START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF], let us present the main points of the Householder tridiagonalization. The basic ingredient is the Householder matrix H whose expression, for some vector u different from the zero vector, is given by

H = I -uu t /b, b = u t u/2. ( 8
)
The idea then is to choose the right vectors u n ,..., u 3 associated to H n ,..., H 3 . The product of these matrices yields to the orthogonal matrix Q = H n H n-1 ...H 3 and their successive applications on A provide:

U = Q t AQ = H t 3 ...H t n AH n ...H 3 .
In [START_REF] Vömel | Divide & Conquer on Hybrid GPU-Accelerated Multicore Systems[END_REF], the authors use a Level 3 BLAS that involves a reduction to intermediate banded and subsequent tridiagonal form.

Parallel Cyclic Reduction

Let us also give the highlights of PCR and refer to [START_REF] Hockney | Parallel Computers: Architecture, Programming and Algorithms[END_REF][START_REF] Zhang | Fast Tridiagonal Solvers on the GPU[END_REF] for more details on the subject. At this stage, we are interested by the resolution of the linear system U Z = Q t Y where the value of V = Q t Y = (v 1 , ..., v n ) is known and U is tridiagonal and symmetric i.e.

U =             d 1 c 1 c 1 d 2 c 2 0 c 2 d 3 . . . . . . . . . . . . 0 . . . . . . c n-1 c n-1 d n             . ( 9
)
PCR comes from a simple modification of Cyclic Reduction (CR) which is schematized in Figure 7 where the linear equations e1,...,e8 constitute the equality U Z = V when n = 8. Applied to this symmetric linear system of equations (e1,...,e8), the first step of CR reduces the number of 8 equations to 4 equations defined by

             e2': d ′ 2 z 2 + c ′ 2 z 4 = v ′ 2 e4':c ′ 2 z 2 + d ′ 4 z 4 + c ′ 4 z 6 = v ′ 4 e6':c ′ 4 z 4 + d ′ 6 z 6 + c ′ 6 z 8 = v ′ 8 e8': c ′ 6 z 6 + d ′ 8 z 8 = v ′ 8 with              d ′ i = d i - c 2 i-1 d i-1 - c 2 i d i+1 c ′ i = -c i c i+1 /d i+1 v ′ i = d i - c i-1 v i-1 d i-1 - c i v i+1 d i+1 . (10)
On this new system, we perform another similar reduction that yields to a system of two equations that only involve z 4 and z 8 . The resolution of the latter system makes possible the backward stage of resolving [START_REF] Crépey | A Lévy HJM multiple-curve model with application to CVA computation[END_REF] and finally the original system of 8 equations.

CR is specified by a forward reduction phase then a backward phase to recover the solution. CR suitability to GPU and its implementation was already studied in both [START_REF] Goddeke | Cyclic Reduction Tridiagonal Solvers on GPUs Applied to Mixed Precision Multigrid[END_REF][START_REF] Zhang | Fast Tridiagonal Solvers on the GPU[END_REF]. The authors of [START_REF] Goddeke | Cyclic Reduction Tridiagonal Solvers on GPUs Applied to Mixed Precision Multigrid[END_REF] propose a method to overcome shared memory bank conflicts during CR, but it uses 50% more on-chip storage. Because of this extra-storage, this trick must not be used in our case since it reduces the number of independent groups of threads involved on different linear systems.

Unlike CR, PCR requires only forward reductions. In Figure 7 for instance, the PCR version would apply two simultaneous reductions in Step 1: A reduction to obtain a new system involving (z 2 , z 4 , z 6 , z 8 ) and another reduction that provides a new system involving (z 1 , z 3 , z 5 , z 7 ). Generally speaking, for a system of size n, PCR reduces the system of n equations to 2 systems of n/2 equations, then to 4 systems of n/4 equations and so on until reaching n/2 systems of 2 e1 e7 e5 z2 z8 z4 z6

Step 1: Forward reduction to a 4-unknown system involving z2, z4, z6 and z8

Step 2: Forward reduction to a 2-unknown system involving z4 and z8

Step 3: Solve 2-unknown system

Step 4: Backward substitution to solve the rest 2 unknowns

Step 5: Backward substitution to solve the rest 4 unknowns equations that can be simply solved. This process makes PCR more suited to parallel architecture and prevent the bank conflicts of shared memory. Nevertheless, PCR can be improved for large systems n > 64 by a combination with CR as detailed in [START_REF] Zhang | Fast Tridiagonal Solvers on the GPU[END_REF].

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A e3 z3 z1

Adaptation and optimization

Householder tridiagonalization

We present two different versions of the Householder tridiagonal factorization.

1. An SIMD version that requires only threads of T I , one for each linear system. 2. A collaborative version that involves n threads of T C for each linear system with n unknowns.

The number of the T I -independent groups of threads is taken to be the one that saturates the shared memory or that executes sufficient work per SMs.

The SIMD version is a single-threaded CUDA adaptation of the procedure proposed in [START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF]. The collaborative version is also based on [START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF], but it provides a multi-threaded implementation of independent tasks which makes it much more efficient than the SIMD version. We begin by explaining the algorithmic steps of the common procedure. The first stage is to compute the tridiagonal form U through successive zeroing of the columns of matrix A = (A i,j ) i,j=1,...,n . This stage is processed at each step i = n, ..., 3 beginning by the vector

u t i = A i,1 , ..., A i,i-1 ± √ σ, 0, ..., 0 , σ = A 2 i,1 + ... + A 2 i,i-1 , (11)
then calculating the intermediary variables

b i = u t i u i 2 , , p i = U i+1 u i b i , B i = u t i p i 2b i , q i = p i -B i u i ( 12
)
Step 1: Reduced to 2 systems of 4 unknowns

Step 2: Reduced to 4 systems of 2 unknowns

Step 3: Solve which allow us to set

A A A A A A A A A A A A A A A A A A A A A A e'7 e'5 e'3 e'1 A A A A A A A A A A A A A A A A A A A A e''
U = U 3 with U i = U i+1 -q i u t i -u i q t i and U n+1 = A. (13)
Now that we have the tridiagonal form U , the second stage is to compute the orthogonal matrix Q defined by Q = H n H n-1 ...H 3 . Besides, we remind that a Householder matrix H i is completely defined thanks to u i . Consequently, during the first stage, the nonzero components of u i are stored in the ith row of the shared memory space allocated for A and u i /b i in the ith column. Thus, the computation of Q is performed in the second stage using Q = Q n and the induction

Q i = H i Q i-1 for i = 4, ..., n with Q 3 = H 3 . ( 14
)
By definition, Q i is an identity matrix in the last i rows and columns and only its elements up to row and column i -1 need to be computed. These then overwrite u i and u i /b i stored in A in the first stage.

As far as the first stage is concerned, in addition to the n × n shared memory space allocated for A we need 2n + 1 extra shared memory space. The latter space is used to store the diagonal and the off-diagonal plus 1 value needed for the synchronization between phases where only one thread can be used and the other phases. Also, since p i is of size i its components can be stored temporarily in the place of undetermined elements of the off-diagonal. Regarding q i , it overwrites p i in the off-diagonal.

Let us now take a look at the multi-threaded parts of the collaborative version. For the second stage, we can use i -1 threads of type T C that need synchronization only when the calculation of Q i is finished. As for the first stage, the computations of p i and q i in ( 12) and the induction performed in (13) are all parallelized using i -1 threads. The other parts of this stage are executed using only one thread.

Parallel Cyclic Reduction

It is important to point out that CR and PCR implementations proposed in [START_REF] Goddeke | Cyclic Reduction Tridiagonal Solvers on GPUs Applied to Mixed Precision Multigrid[END_REF][START_REF] Zhang | Fast Tridiagonal Solvers on the GPU[END_REF] cannot be used directly for any system size. Indeed, the first reference requires a size that is equal to a power of two plus one and the versions of the second paper are done for a size equal to a power of two. The simplest way to use both implementations for any size would be to perform a zero padding. In our case, using this latter technique is not a good option since we fill a part of the shared memory with zeros instead of using it for the resolution of other systems. This fact motivates our version of the PCR presented bellow.

We propose to implement the PCR with permutations of equations in the way illustrated in Figure 8. The idea is to gather the equations involved in the same system and to separate those that are independent. Indeed, using this simple idea one can deal with any size. Let us take the example n = 7, the changes that occur on the matrix of the system are the following where (R) and (P) mean respectively Reductions and Permutations and the empty parts of the matrices represent the null part. Obviously, one should also set the same kind of reductions and permutations on the vector of unknowns Z as well as on the vector of values V . To reorder the final solution, we use an array of size n in addition to the usual 3n memory space needed for PCR.

          d 1 c 1 c 1 d 2 c 2 c 2 d 3 c 3 c 3 d 4 c 4 c 4 d 5 c 5 c 5 d 6 c 6 c 6 d 7           (R) -→           d ′ 1 0 c ′ 2 0 d ′ 2 0 c ′ 3 c ′ 2 0 d ′ 3 0 c ′ 4 c ′ 3 0 d ′ 4 0 c ′ 5 c ′ 4 0 d ′ 5 0 c ′ 6 c ′ 5 0 d ′ 6 0 c ′ 6 0 d ′ 7           (P ) -→           d ′ 1 c ′ 2 c ′ 2 d ′ 3 c ′ 4 c ′ 4 d ′ 5 c ′ 6 c ′ 6 d ′ 7 0 0 d ′ 2 c ′ 3 c ′ 3 d ′ 4 c ′ 5 c ′ 5 d ′ 6           (R) -→           d ′′ 1 0 c ′′ 2 0 d ′′ 3 0 c ′′ 4 c ′′ 2 0 d ′′ 5 0 c ′ 4 0 d ′ 7 0 0 d ′′ 2 0 c ′′ 3 0 d ′′ 4 0 c ′′ 3 0 d ′′ 6           (P ) -→           d ′′ 1 c ′′ 2 c ′′ 2 d ′′ 5 0 0 d ′′ 3 c ′′ 4 c ′′ 4 d ′ 7 0 0 d ′′ 2 c ′′ 3 c ′′ 3 d ′′ 6 0 0 d ′′ 4          

Comparison of versions and comparison with LDLt

In this section, we reuse the matrices A = Q Ξρ RQ ′ Ξρ introduced in Section 3.3. Unlike LDLt and PCR, the Householder factorization does not require A to be definite positive and thus one can even take R to be only tridiagonal symmetric random matrix and not diagonally dominant.

As stated before, we ought to tridiagonalize A: A = QU Q t then use the PCR as well as two matrix/vector multiplications U Z = Q t Y , X = QZ to recover the solution. The PCR and the multiplications are much faster than the tridiagonal factorization and, like for resolutions LZ = Y and DL t X = Z in LDLt, can be neglected when compared to the overall execution time.

To make the previous affirmation more quantitative, Figure 9 shows a huge numbers of PCRs and matrix/vector multiplications that can be computed per second. These numbers generally exceeds the number of systems solved using LDLt shown in Figure 6 of the previous section. However, they coincide for very small systems as we do much less computations on the shared memory compared to the time spent in the access to the GPU global memory.

Regarding the benefits of using the collaborative solution instead of the SIMD version of Householder tridiagonalization, the continuous line in Figure 10 shows a quite significant speedup that increases with respect to the size n. Besides, the dashed line in Figure 10 shows the execution time superiority of the LDLt hybrid solution when compared to the collaborative version of the Householder tridigonalization + PCR. Basically, the LDLt hybrid solution is

Divide and conquer algorithm for tridiagonal eigenproblems

Presentation of the algorithm and the references

We reuse here the Householder tridiagonalization of Section 4 as a first stage for the resolution of AX = Y with A symmetric. This new resolution procedure is implemented through the following steps:

• We perform the tridiagonal Householder decomposition A = QU Q t where Q is orthogonal and U is symmetric tridiagonal.

• We use the divide & conquer algorithm for symmetric tridiagonal eigenproblems to establish the eigenvalues and eigenvectors of U = ODO t where O is orthogonal and D is diagonal. Consequently, we have A = N DN t where the orthogonal matrix N = QO.

• In the way that is usually done to solve a linear system with numerical singularities [START_REF] Press | Numerical Recipes in C++: The Art of Scientific Computing[END_REF], we discard the smallest eigenvalues of D that provide a condition number larger than 10 5 .

Because the first step was already studied and the third step is quite standard, we are only interested by the divide & conquer part. This latter method goes back to the reference [START_REF] Cuppen | A divide and conquer method for the symmetric tridi-agonal eigenproblem[END_REF] and it became numerically sustainable since the work presented in [START_REF] Gu | A stable algorithm for the rank-1 modification of the symmetric eigenproblem[END_REF][START_REF] Gu | A divide-and-conquer algorithm for the symmetric tridiagonal eigenproblem[END_REF]. In [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF], one can find a quite detailed presentation of divide & conquer algorithm for eigenproblems. As far as we are concerned, we study the important points that should be explored in our adaptation, which are:

1. Divide the diagonalization problem in two diagonalization subproblems with known diagonal factorization. 2. Solve the secular equation. 3. Use Löwner's Theorem [START_REF] Löwner | Über monotone Matrixfunctionen[END_REF][START_REF] Demmel | Applied Numerical Linear Algebra[END_REF] for the stability of the overall procedure. 4. Perform a matrix multiplication to conquer the diagonalization problem from the diagonalized subproblems.

Let a matrix U given by ( 9), the first point would be to make the following division

U =                 d 1 c 1 c 1 . . . . . . . . . . . . c m-1 c m-1 d m -c m 0 0 d m+1 -c m c m+1 c m+1 . . . . . . . . . . . . c n-1 c n-1 d n                 +c m 1 m,m+1 1 t m,m+1 (15)
= U 1 0 0 U 2 + c m 1 m,m+1 1 t m,m+1
where 1 m,m+1 = (0, ..., 0, 1, 1, 0, ..., 0) with only the (m) th and (m + 1) th coordinates equal to 1, all the other coordinates are null. As assumed, U 1 and U 2 have a known diagonal factorization i.e. there exists diagonal matrices D 1 , D 2 and orthogonal matrices

O 1 , O 2 such that U 1 = O 1 D 1 O t 1 and U 2 = O 2 D 2 O t 2 . Consequently, one can rewrite U as U = O 1 0 0 O 2 D 1 0 0 D 2 + c m uu t O t 1 0 0 O t 2 where u = O t 1 0 0 O t 2 1 m,m+1 = last column of O t 1 first column of O t 2 .
Denote now Λ = {λ 1 , ..., λ n } the ordered family of eigenvalues of

D 1 0 0 D 2 .
One can show that if c m = 0 and the eigenvalue λ of U satisfies λ / ∈ Λ, then its value is obtained as a solution of the secular equation

n i=1 u 2 i λ i -λ + 1 c m = 0. ( 16
)
The reference [START_REF] Li | Solving Secular Equations Stably and Efficiently[END_REF] provides a good summary on the different methods used for the solution of [START_REF] Gordy | Nested Simulation in Portfolio Risk Measurement[END_REF]. It also proposes an hybrid procedure whose performances compete even with Gragg's scheme [START_REF] Gragg | Parallel divide and conquer algorithms for the symmetric tridiagonal eigenproblem and bidiagonal singular value problem[END_REF] that has a cubic convergence. The major advantage of the hybrid scheme comes from the fact that it prevents additional computations due to the second order differentiation of the left term of equality [START_REF] Gordy | Nested Simulation in Portfolio Risk Measurement[END_REF].

Once we fix an eigenvalue λ which is a solution of ( 16), the eigenvector V λ of D 1 0 0 D 2 + c m uu t can be computed by

V λ = D 1 0 0 D 2 -λI -1 u (17)
where the vector u is defined in [START_REF] Gu | A stable algorithm for the rank-1 modification of the symmetric eigenproblem[END_REF][START_REF] Gu | A divide-and-conquer algorithm for the symmetric tridiagonal eigenproblem[END_REF] thanks to Löwner's Theorem. Replacing u by u is quite important in order to ensure stability and sustainability of the algorithm, especially when some eigenvalues are almost equal. Let us assume now that all eigenvectors W = (V λ ) λ eigenvalue of U are known which brings us to the final point. To conquer, we need to compute the eigenvectors of U using the product

Q 1 0 0 Q 2 W
. This last step is the heaviest numerically in the whole algorithm. Finally, we voluntarily did not present the deflation that appears when u i or λ i -λ i+1 vanish numerically. This is due to the fact that deflation is already well presented in the references cited above and to the fact that deflation is not quite important when matrices are small.

Adaptation and optimization

Let us study how the steps 1.→4. should be and are implemented in our source code [30]. Because 4. is the heaviest part, we start with it then we go decreasingly until the first step.

From step 3., we have at our disposal on the shared memory: the eigenvalues of U , the transpose matrix W t as well as

Q 1 0 0 Q 2 . Consequently, the computation of Q 1 0 0 Q 2 W , at step 4., is done using W t Q t 1 0 0 Q t 2 then
processing a transpose operation. Nevertheless, we do not need to transpose

Q 1 0 0 Q 2 since W t Q t 1 0 0 Q t 2
involves the dot product of the rows of W t with the rows of Q 1 0 0 Q 2 . Associating one thread for each row of the latter matrix, the memory acccess is coalescent during the successive multiplications.

The result of the dot product performed by each thread is stored in the register memory then it is copied, after a synchronization, to the shared memory space of each row of W t . The complexity of step 4. is then O(n 3 ) and we use 2n 2 of the shared memory to perform the matrix multiplication. Because of the zeros in Q 1 0 0 Q 2 , one can decrease further the memory occupation but to the detriment of the readability of the code.

From step 2., we have the eigenvalues of U and Q 1 0 0 Q 2 at our disposal on the shared memory. We would like to use one thread for the computation of each eigenvector of D 1 0 0 D 2 +c m uu t to be stored in the column of W . This is performed thanks to expression [START_REF] Gragg | Parallel divide and conquer algorithms for the symmetric tridiagonal eigenproblem and bidiagonal singular value problem[END_REF] and each resulted eigenvector is saved in a row-form to keep the coalescent access of each thread. Subsequently, we obtain W t instead of W . The complexity of this step is O(n 2 ) and it needs 3n shared memory space: n for the eigenvalues and 2n for both u and the diagonal (D 1 , D 2 ) involved in [START_REF] Gragg | Parallel divide and conquer algorithms for the symmetric tridiagonal eigenproblem and bidiagonal singular value problem[END_REF].

From step 1., we have the diagonal (D 1 , D 2 ), u and Q 1 0 0 Q 2 at our disposal on the shared memory. Before starting the resolution of ( 16), we need to sort (D 1 , D 2 ) and build the ordered set Λ = {λ 1 , ..., λ 2 }. Although multithreaded, this sorting does not need to be optimized because its complexity is at most O(n 2 ). Moreover, the sorting result is stored in a new shared memory space of size n plus some variables stored temporarily in the memory space of W t (Not used yet). Afterwards, we solve (16) using Gragg's scheme [START_REF] Gragg | Parallel divide and conquer algorithms for the symmetric tridiagonal eigenproblem and bidiagonal singular value problem[END_REF] that has the advantage of a cubic and monotonous convergence. The fact that this scheme requires a second order differentiation, of the left term of equality [START_REF] Gordy | Nested Simulation in Portfolio Risk Measurement[END_REF], is not an important drawback when the size n is small. Also, because n is small, the complexity of the iterative Gragg's scheme is rather O(n 3 ) instead of O(n 2 ), considered for big values of n. Finally, we arrive to step 1. that represents the heart of the algorithm as it sets the division. One has then to choose a constant m and perform [START_REF] Goddeke | Cyclic Reduction Tridiagonal Solvers on GPUs Applied to Mixed Precision Multigrid[END_REF] such that U 1 and U 2 are already diagonalized. Otherwise, we have to reiterate the division (15) for U 1 and U 2 separately and so on till reaching a division that has a diagonal factorization. Assume now that for all m ∈ {2, ..., n -2}, we do not know yet the diagonal factorization of U 1 and the diagonal factorization of U 2 . What is then the best choice to do on the value of m?

There are two answers to the previous question depending on whether there exists m 0 such that c m0 = 0 or not:

• If yes, then set m = m 0 . Moreover, we re-apply directly (15) on the sub-matrices.

• If no, then set m = ⌊n/2⌋.

The yes case is obvious since the larger linear system can be decomposed, without conquering, in two independent systems. Regarding the no case, the choice is set in such a way that the computations induced by each subproblem is comparable to the other. To make this latter fact possible, the best way is to impose the same size (±1) for both problems.

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A n L n/2 L n-m m= L m/2 L L m/2 L m- L (n-m)/2 L L (n-m)/2 L n-m- 1 2 2 1 2 1
In order to simplify the comparison between the Householder tridiagonalization and the divide and conquer algorithm, we consider only matrices with c m = 0 for all m ∈ {1, ..., n -1}. Our division is then performed using the scheme given in Figure 11 till reaching matrices of dimension 1 × 1 or 2 × 2 for which the diagonal form can be obtained easily. This division scheme requires an extra shared memory storage of size 2 1+⌊log 2 (n-1)⌋ , but it provides a pure divide and conquer algorithm (Not a combination of divide and conquer with another method like QR). In particular, this pure divide and conquer prevents to have eigenvalues of multiplicity bigger than two at each conquering step.

Thanks to what is explained above, it is not difficult to conclude that the overall shared memory occupation is given by 2n(n + 2) + 2 1+⌊log 2 (n-1)⌋ . In addition, the complexity t of the proposed implementation can be computed thanks to the induction t(m) = 2t(m/2) + α(m)m 3 , m ∈ {2, ..., n}

with t(1) = 1 and α(m) is a decreasing sequence that is bigger than 2 for the sizes considered in this paper. Using [START_REF] Gu | A stable algorithm for the rank-1 modification of the symmetric eigenproblem[END_REF], we check that t(n) ≤ α(n)4 3 n 3 .

Comparison with Householder tridiagonalization

The complexity of the divide and conquer algorithm is generally considered as similar to the one of Householder tridiagonalization for large matrices (see [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF]). This is not the case in our examples because we deal with small ones (dimension of at most 64).

Moreover, the divide and conquer algorithm suffers from divergence problems when implemented on GPU. Indeed, the need for deflation in some cases can lead numerous threads to wait. The necessary use of more synchronization in this algorithm also reduces the performance. For instance, the resolution of the secular equation is iterative and so makes some threads wait for the others.

All the facts related above justify the results obtained in Figure 12. In particular, we see that for the largest matrices, 48 ≤ n ≤ 64, the divide and conquer takes at most three × the execution time of a Householder factorization.

Conclusion and future work

In this work, we presented the adaptation of well-known methods to the resolution of large number of small symmetric linear systems on single precision GPUs. We also provided original ideas for further optimization using LDLt and PCR. Our goal was to know if the use of Householder tridiagonalization with divide & conquer is the best solution when we suspect the linear systems to be ill-conditioned.

We have shown that it is better to first check the residual using Householder tridiagonalization + PCR. If the really fast PCR is not sufficient, performing a divide & conquer diagonalizations and discard the smallest eigenvalues becomes mandatory.

If we are sure that the system is well-conditioned then we just process an LDLt decomposition. If we are sure of the converse, we execute directly a combination of Householder tridiagonalization and divide & conquer diagonalization.

The conclusion is definitely related to the context of the CVA computation. As we have to solve a large number of small random linear system, it is difficult to assert that all the systems will be well-conditioned. As a consequence, it would be naive to implement directly a simple LDLt or Cholesky decomposition.

As a future work, we plan to explore the accuracy of each method by studying the rounding errors and error propagation. For that, we aim to present a sufficiently consistent study of the residue errors as well as compare the results of CADNA [31] software obtained from the various solutions.

Figure 1 :

 1 Figure 1: An example of a two-stage simulation with M 0 = 2, M 6 = 8 and M 8 = 4.

Figure 2 :

 2 Figure 2: Condition numbers for linear regression associated to Black & Scholes model: σ = 0.1, µ = 0.1, S 0 = (1, ..., 1), t k = 0.1, n = 30, M k = 300 and l ∈ {1, ..., 1000}.

Figure 3 :

 3 Figure 3: Standard LDLt parallel strategy.

First

  , let us introduce the set of matrices used for the tests. The matrices introduced in sections 4.3 and 5.3 are related to the one presented here which

Figure 4 :

 4 Figure 4: The speedup of the collaborative and the hybrid versions when compared to the SIMD implementation.

Figure 5 :

 5 Figure 5: The optimal number of communicating threads in the hybrid version.

Figure 6 :

 6 Figure 6: Number of linear systems that can be solved within a second.

Figure 7 :

 7 Figure 7: Cyclic Reduction for n = 8 unknowns: Communication pattern showing the dataflow between equations. Letters e' and e" stand for updated equation.

Figure 8 :

 8 Figure 8: Modified Parallel Cyclic Reduction for n = 8: Communication pattern showing the dataflow and permutation of equations. Letters e' and e" stand for updated equation.

Figure 9 :

 9 Figure 9: Number of PCRs + two matrix/vector multiplications that can be performed per second.

Figure 10 :

 10 Figure 10: The execution time ratio of: SIMD/Collaborative and (tridiagonal + PCR)/LDLt

Figure 11 :

 11 Figure 11: The division scheme.

Figure 12 :

 12 Figure 12: The execution time of the divide and conquer when compared to the Householder tridiagonalization.

Acknowledgements

This work was funded by project ARRAND (ANR-15-CE39-0002-01) and partially supported by the project FastRelax (ANR-14-CE25-0018-01).


Edito #7: Poïétique des ambiances architecturales et urbaines : les mises en oeuvre de l'expérience_The poietics of architectural and urban ambiences: making use of experience Par_By Grégoire Chelkoff Email: gregoire.chelkoff [at] grenoble.archi.fr, architecte, professeur à l'ENSA de Grenoble, directeur du Cresson.

. Les multiples dimensions que l'idée d'ambiance articulent peuvent dérouter celui qui en interroge la faisabilité et les différentes modalités de création. Point de recettes ni de voies toutes tracées. Partant que l'environnement est à concevoir comme façonnant et façonné, la perspective qualitative proposée dans ce court texte est d'envisager une approche poïétique dont l'objet est d'étudier l'ambiance mise en oeuvre par l'usage ordinaire des espaces. C'est en même temps la manière dont le lieu s'imagine en nous et « s'énacte » dans nos pratiques, qui forment l'ambiance, on pourrait dire : qui font oeuvre d'ambiance.

The many dimensions contained in the idea of "ambience" can be confusing for anybody considering the feasibility and various conditions of creativity. There are no recipes or predefined paths. Based on the idea that the environment should be designed as a structuring and structured element, the qualitative perspective proposed in this short text requires the use of a poietic approach to study the ambience introduced by the ordinary utilisation of space. It is both the way in which the place is imagined in our minds and the way in which it is "enacted" in our everyday life that form ambience. In fact, one might almost say that they create a work of ambience.

Penser l'ambiance en architecte, et penser l'architecture comme ambiance, c'est alors considérer les manières dont les dispositifs spatiaux et matériels initient ou prédisposent certaines formes d'expérience composant avec les flux Thinking of ambience as an architect and thinking of architecture as an ambience means considering the ways in which spatial and material devices initiate or predispose certain types of experience that combine flows of light, (lumière, chaleur, son, mouvement) et qui sont prégnantes et significatives dans l'accomplissement sensible de la vie quotidienne. Les formes d'expérience sont à envisager dans toute leur banalité mais aussi dans toute leur force de façonnage.

Monter un escalier, traverser une rue, se mettre au balcon, marcher sur un sol, s'appuyer contre une paroi, appeler de sa fenêtre, tourner un angle de rue, passer une porte, se caler dans une embrasure, chacune de ces situations peut être pensée comme forme d'expérience spécifique faisant de l'environnement une oeuvre en cours qui nous implique.

A l'idée plastique séduisante de l'oeuvre comme « jeu savant et correct des volumes sous la lumière » vient se greffer une autre dimension oeuvrante de l'architecture qui pourrait être celle des jeux savants du corps et de la lumière dans l'exécution d'un mouvement humain.

Cette

poïétique des ambiances ordinaires demande l'accumulation d'expériences directes qualitatives pour établir un socle de connaissance empirique argumenté. L'observation en est un enjeu primordial, tant par les phénomènes saisis que par le processus heat, sound and movement. These flows are relevant and significant in the sensitive accomplishment of everyday life. Forms of experience should be envisaged in all their banality and also with all their strength as structuring elements.

Climbing a flight of steps, crossing a street, going out onto a balcony, walking, leaning against a wall, calling from a window, turning the corner of a street, going through a door, standing in the doorway -all these situations can be seen as a specific form of experience, making the environment a work in progress, a work that involves us.

The attractive aesthetic idea of a work as the "skilful, correct play of volumes under light" is extended by another structuring dimension of architecture that might be called the skilful play of the body and light in the performance of a human movement.

The poietics of ordinary ambiences demand an accumulation of direct