Pascal Havé
email: pascal.have@ifpenergiesnouvelles.fr

Roland Masson
email: roland.masson@ifpenergiesnouvelles.fr

Frédéric Nataf
email: nataf@ann.jussieu.fr

Mikolaj Szydlarski
email: mikolaj.szydlarski@gmail.com

Hua Xiang
email: hxiang@whu.edu.cn

Tao Zhao
email: zhao@ann.jussieu.fr

Algebraic Domain Decomposition Methods for Highly Heterogeneous Problems

We consider the solving of linear systems arising from porous media flow simulations with high heterogeneities. Using a Newton algorithm to handle the non-linearity leads to the solving of a sequence of linear systems with different but similar matrices and right hand sides. The parallel solver is a Schwarz domain decomposition method. The unknowns are partitioned with a criterion based on the entries of the input matrix. This leads to substantial gains compared to a partition based only on the adjacency graph of the matrix. From the information generated during the solving of the first linear system, it is possible to build a coarse space for a two-level domain decomposition algorithm that leads to an acceleration of the convergence of the subsequent linear systems. We compare two coarse spaces: a classical approach and a new one adapted to parallel implementation.

Introduction

Multiphase, compositional porous media flow models, used for example in reservoir simulations or basin modeling, lead to the solution of complex nonlinear systems of Partial Differential Equations (PDEs) accounting for the mass conservation of each component and the multiphase Darcy law, coupled with thermodynamical equilibrium and volume balance closure laws. These PDEs are typically discretized using a cell-centered finite volume scheme and a fully implicit Euler integration in time in order to allow for large time steps. After Newton type linearization, one ends up with the solution of a linear system at each Newton iteration which represents up to 90 percent of the total simulation elapsed time.

Such linear systems couple an elliptic (or parabolic) unknown, the pressure, and hyperbolic (or degenerate parabolic) unknowns, the volume or mass fractions. They are nonsymmetric, and ill-conditioned in particular due to the elliptic part of the system, and the strong heterogeneities and anisotropy of the media. Their solution using an iterative Krylov method, such as GMRES or BiCGStab, requires the construction of an efficient preconditioner which should be scalable with respect to the heterogeneities and anisotropy of the media, the mesh size, and the number of processors, and should cope with the coupling of the elliptic and hyperbolic unknowns.

Nested factorization [START_REF] Appleyard | Nested factorization[END_REF] and CPR-AMG [START_REF] Lacroix | Decoupling preconditioners in the implicit parallel accurate reservoir simulator (ipars)[END_REF], [START_REF] Scheichl | Decoupling and block preconditioning for sedimentary basin simulations[END_REF] are the main state of the art preconditioners currently used in industrial reservoir simulators. Nevertheless they still suffer from major drawbacks that should be addressed in response to the evolution of the computing architectures, and the demand for more complex physics and geology in reservoir and basin simulations. The nested factorization preconditioner is not adapted to distributed architectures and has a limited scalability with respect to the heterogeneities and anisotropy of the media. CPR-AMG is a more recent preconditioner which combines multiplicatively a parallel Algebraic MultiGrid preconditioner for a pressure block of the linear system, with a parallel incomplete factorization preconditioner (typically ILU0) for the full system. This preconditioner exhibits a very good robustness with respect to the heterogeneities and anisotropy of the media. However, its parallel scalability requires a very large number of cells per processor, say 100000, due to the coarsening step of AMG which is not strongly scalable. The robustness of the CPR-AMG also requires the definition of a pressure block which should be close to an M-matrix. This induces a sensitivity to the physics such as complex wells couplings, capillary driven flows or strong nonlinearities in pressure of the thermodynamical closure laws, and also a sensitivity to the distorsion of the mesh when multipoint flux approximations are used for the Darcy fluxes discretization.

Solving these drawbacks motivates our research for an alternative pressure block preconditioner which should remain scalable with respect to the heterogeneities and anisotropy of the media and should exhibit improved strong scalability on massively distributed architectures and improved robustness with respect to the physics and to the discretization.

The pressure block matrix is related to the discretization of a Darcy equation with high contrasts and anisotropy in the coefficients. We work in the context of overlapping Schwarz type methods on parallel computers. In order to deal with anisotropy, we force the domain decomposition to respect the strong connections between the nodes as much as possible. This is inspired by coarsening techniques in algebraic multigrid (AMG) methods. It is also well known that efficient domain decomposition methods demand a well suited coarse grid correction. For matrices arising from scalar problems with smooth coefficients, it is possible to build a priori (i.e. before any linear solve) efficient coarse spaces based on domain wise constant vectors, see [START_REF] Toselli | Domain Decomposition Methods: algorithms and theory[END_REF] and references therein. For problems with high heterogeneities, the numerical computation of these coarse spaces is often based on solving generalized eigenvalue problems in subdomains, see [START_REF] Efendiev | Multiscale finite element methods for high-contrast problems using local spectral basis functions[END_REF] and [START_REF] Nataf | A two level domain decomposition preconditioner based on local Dirichlet-to-Neumann maps[END_REF]. The corresponding local matrices are not submatrices of the global matrix A and cannot be built at the algebraic level. In this paper, we introduce a new algebraic coarse space construction based on an analysis of the Krylov space generated by the linear solve at the first iteration of a Newton-Ralphson algorithm. We take advantage of a parallel implementation to build a richer coarse space than the one proposed in[RR00, RR98].

The paper is organized as follows. In section 2, we recall the basis of the overlapping Schwarz method. In the next section, we show the benefit of an AMG style partitioning. In section 4 we consider coarse grid corrections. In § 4.1 we recall the principles of coarse grid correction in the context of Newton-Ralphson method, and in § 4.2 we introduce a new domain wise split coarse space and give numerical results on problems arising from reservoir simulations.

Parallel Schwarz method

The discretization of a linear partial differential equation on a domain Ω yields a linear system of the form

A u = f ∈ R n (1)
that we solve using a domain decomposition method. Without loss of generality, we consider here a decomposition of a domain Ω into two overlapping subdomains Ω 1 and Ω 2 . The overlap is denoted by O := Ω 1 ∩ Ω 2 . This yields a partition of the domain: Ω = Ω(1)

I ∪ Ō ∪ Ω (2)

I

where Ω (i)

I := Ω i \ Ō, i = 1, 2.
At the algebraic level this corresponds to a partition of the set of indices N into three sets: N

(1)

I , O and N (2) I .
After the discretization of the boundary value problem defined in domain Ω, we obtain a linear system of the following form

A u :=    A (1) II A (1) IO A (1) OI A OO A (2) OI A (2) IO A (2) II       u (1) I u O u (2) I    =    f (1) I f O f (2) I    .
(2)

We can also define the extended linear system by duplicating the variables

à ũ :=      A (1) II A (1) IO A (1) OI A OO A (2) OI A (1) OI A OO A (2) OI A (2) IO A (2) II           u (1) I u (1) O u (2) O u (2) I      =     f (1) I f O f O f (2) I     , (3)
where the subscript 'O' stands for 'overlap', u

O are the duplicated variables in the overlapping domain O, u

I are variables in the subdomain Ω i I . It is easy to check that if A OO is invertible, there is an equivalence between problems (2) and (3).

We introduce now three preconditioners, two for the linear system (2) and one for the extended linear system (3). A classical preconditioner to problem (2) is the additive Schwarz method. Let R j be the rectangular restriction matrix to subdomain Ω j , j = 1, 2. The additive Schwarz preconditioner is:

M -1 ASM := R T 1 A -1 1 R 1 + R T 2 A -1 2 R 2
where

A i := R i AR T i , i = 1, 2. If A is SPD, then M -1
ASM is SPD as well and the theory of this preconditioner is very well developped, see [START_REF] Toselli | Domain Decomposition Methods: algorithms and theory[END_REF] and references therein. But in the overlap, corrections are added twice. This somehow delays the convergence. In order to fix this problem, another classical preconditioner was designed: the restricted additive Schwarz (RAS) method, see [START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF]. Define Rj by setting some ones in R j to zeros, such that the operators Rj correspond to a non-overlapping decomposition,

RT 1 R 1 + RT 2 R 2 = I. (4)
Then the restricted additive Schwarz preconditioner reads

M -1 RAS := RT 1 A -1 1 R 1 + RT 2 A -1 2 R 2 . ( 5
)
Note that the RAS avoids extra correction in the overlapping zone but at the expense of the loss of the symmetry of the preconditioner. As a result, there are less theoretical results for RAS than for ASM. In practice, it was noticed that the RAS preconditioner outperforms the ASM preconditioner. The RAS preconditioner leads to an iterative method that has equivalences with the discretization of the original Schwarz method [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF], see [START_REF] Efstathiou | Why Restricted Additive Schwarz converges faster than Additive Schwarz[END_REF]. The third Schwarz type method is called the Jacobi-Schwarz method (JSM). It corresponds to the discretization of the parallel variant of the original Schwarz algorithm [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF]. It consists of solving the extended problem (3) preconditioned by a block Jacobi method:

M JSM ( A) :=      A (1) II A (1) IO A (1) OI A OO A OO A (2) OI A (2) IO A (2) II      . ( 6
)
and one can easily notice that M -1 JSM ( A) can be computed in parallel. Actually, it is sufficient to factorize the diagonal blocks in parallel. When used in a Richardson algorithm, it was proved in [START_REF] Efstathiou | Why Restricted Additive Schwarz converges faster than Additive Schwarz[END_REF] that M JSM ( A) applied to (3) and M RAS applied to (2) lead to equivalent algorithms. Since RAS is easier to implement than JSM, it gives a clear advantage to RAS. But when considering two level methods, JSM has the advantage that no partition of unity is needed. It brings some benefit both in terms of implementation and efficiency as it was noticed in [START_REF] Nataf | A coarse space construction based on local Dirichlet to Neumann maps[END_REF]. Thus, in the sequel, we will only consider the JSM preconditioner M -1 JSM ( A) applied to (3).

In the general case with many subdomains, the set of indices N is decomposed into N overlapping subsets

N = ∪ N i=1 N i .
From this decomposition, we define the extended system whose right handside belongs to R N E where N E := N i=1 #N i > n. For i = 1, . . . , N we denote by R E,i the boolean matrix corresponding to the restriction operator from

R N E → R Ni : R E,i =   1 0 0 0 0 • • • 0 0 1 0 0 0 • • • 0 0 0 0 1 0 • • • 0   ( 7
)
where there is 1 on a column if the corresponding node belongs to i-th subdomain. The transpose operator R T E,i is the extension by zero from the i-th subdomain to global set of unknowns R N E . In practice, we first perform a partition of the unknowns using a graph partitioner working on the adjacency graph of the matrix. Then each partition is extended with the adjacent nodes.

Partitioning

Graph partitioning softwares such as METIS [START_REF] Karypis | A fast and highly quality multilevel scheme for partitioning irregular graphs[END_REF] or SCOTCH [START_REF] Chevalier | PT-SCOTCH: a tool for efficient parallel graph ordering[END_REF] are commonly used with the goal to compute a balanced partitioning that minimizes a cost function which for sake of efficiency has to take integer values only. The default cost function is the number of edge cuts. In order to accommodate the heterogeneous cases, it is possible to weight the costs of the edgecuts. Here, we want to get a partition that follows the anisotropy of the underlying partial differential equation (PDE). It is known to have a strong effect on the efficiency of preconditioners. Since we work at the algebraic level, this partitioning of the unknowns cannot be made using the underlying mesh and values of the parameters in the PDE. The only information on the problem comes from the entries of the matrix of the linear system. Inspired by algebraic multigrid methods, the weight of the edge (i, j) which has to be an integer is defined as follows:

weight(i, j) := γ |a i,j | |a i,i | + |a j,j | (8)
where γ is a large number. When minimizing the edge-cuts with weights, nodes that are strongly connected are kept in the same subdomain as much as possible.

In order to choose the value for γ, we have to balance two requirements. On one hand, γ must be small enough so that the value of the cost function which is the sum of the weighted edge cuts is smaller than the maximum value of an integer allowed by the partitioning software. On the other hand, γ must be large enough so that the weights take as many different values as possible in order to accurately distinguish between the different values of

|a i,j |/(|a i,i | + |a j,j |).
Typically in our tests γ = 80000.

As an example, we consider the following anisotropic problem:

-div(κ∇u) = f , discretized by a finite element method (FreeFem++, see [START_REF] Hecht | FreeFem++. Numerical Mathematics and Scientific Computation[END_REF]) on 2D unit square of size N x × N y , where N x = N y = 128 so that the number of non zero entries is 98000. The diffusion tensor κ has a strong anisotropy:

κ = κ xx 0 0 κ yy = 1 × 10 -6 0 0 1 . (9)
In Figure 2, we show on the left picture the partition obtained using the weight function of formula (8). We see that strongly connected nodes are kept as much as possible in the same domain. On the right picture, we show the unweighted partition. The problem with this partition is that strongly connected nodes belong to more than one subdomain. As a result, the Jacobi-Schwarz method converges in only 2 iterations as compared to 106 iterations for the unweighted partition. The effect of the weighted partition is now tested against an un-weighted partition on several test cases coming from "real life" reservoir simulations with problems of various sizes, see Figure 3. The first figure corresponds to 50.000 nodes whereas the last figure corresponds to the SPE10 test case with one million unknowns for the pressure block. On Figure 3, we plot the convergence histories for a weighted partition and an un-weighted one. As we see, the weighted partitioning yields a better convergence of the Schwarz method. When the iteration counts are already small (smaller than 20) for the un-weighted partition, there is no much improvement. But, when iterations counts are large, the weighted partition brings a substantial benefit. We also noticed that the cost of one iteration and the time required for partitioning are about the same in both cases.

4 Two-level Schwarz method

Background

In the previous section, we have seen the benefit of using weighted partitions that take into account the anisotropy of the problem. When the number of subdomains becomes large, this is not sufficient to prevent stagnation in the convergence of Schwarz type algorithms. Indeed, any of the three preconditioners M ASM , M RAS or M JSM removes the very large eigenvalues of the coefficient matrix, which correspond to high frequency modes. But the small eigenvalues still exist and hamper the convergence. These small eigenvalues correspond to low frequency modes and represent certain global information. We need a suitable coarse space to efficiently deal with them. This problem is closely related to deflation techniques, see for instance [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF], [PdSM + 06] or [SYEG00] and references therein. These methods are based on a knowledge of an approximation to the eigenvectors corresponding to the "bad" (small in our case) eigenvalues. Let Z denote the rectangular matrix that stores these vectors columnwise.

The number of columns of Z is the size of the coarse space. It is then classical (see for instance [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF]) to define the following matrices:

P := I -AQ, Q := ZE -1 Z T , E := Z T AZ .
Notice that if A is symmetric, we have QAZ = Z, P T Z = 0 and P T Q = 0. From these matrices, it is possible to introduce new preconditioners which will not suffer from plateaus in the convergence. Let M denote a Schwarz method, the balancing Neumann-Neumann preconditioner

P BN N := P T M -1 P + Q was introduced in[Man93]
. In [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF], a related form is introduced:

P A-DEF 2 := P T M -1 + Q . ( 10
)
These preconditioners used in any Krylov method have the interesting property that at any step n the residual r n remains orthogonal to the vector space spanned by the columns of Z: Z T r n = 0 .

Compared to the original preconditioner M , the extra cost of the new preconditioner lies in the solving of a linear system with the small matrix E. In domain decomposition methods, the resulting method is called a two-level method.

The ideal coarse space is the invariant subspace corresponding to the low part of the spectrum of the preconditioned operator. But, computing the small eigenvalues of the preconditioned operator M -1 A is of course too expensive. It is thus necessary to find a way to "guess" some good approximations to them: a priori i.e. before any solve or a posteriori i.e. after a first linear solve with the matrix A or with a matrix close to A. The a priori construction demands some analytic knowledge on the problem to be solved. For instance for Poisson or elasticity type problems, coarse spaces must contain the kernel of the operators: constant functions or rigid body motions, see [START_REF] Toselli | Domain Decomposition Methods: algorithms and theory[END_REF] and references therein. For domain decomposition methods for problems with high heterogeneities, the numerical computation of these coarse spaces is often based on solving generalized eigenvalue problems in subdomains, see [START_REF] Efendiev | Multiscale finite element methods for high-contrast problems using local spectral basis functions[END_REF] and [START_REF] Nataf | A two level domain decomposition preconditioner based on local Dirichlet-to-Neumann maps[END_REF]. The corresponding local matrices are not submatrices of the global matrix A. They cannot be built at the algebraic level. Thus, we focus on a posteriori constructions of the coarse space that can be done at the algebraic level. In order to build suitable coarse spaces, we will reuse information coming from previous solves in the context of nonlinear problems. Let F be a nonlinear mapping and suppose we solve

F (U ) = G by a Newton algorithm For k = 1, . . . , K F (u k ) • δu k+1 = G -F (u k ) u k+1 = u k + δu k+1 ( 11
)
where K is the number of iterations to reach convergence. Let us denote by A k the matrix F (u k ). The sequence of linear systems (11) are solved by the parallel Schwarz method. The idea is to use spectral information from this first solve with the matrix A 1 to build a coarse space Z and thus a coarse correction that will accelerate the solve of the subsequent linear systems with matrices A k , k ≥ 2. Although there are many variants (see for instance [RR98, RR00, GR03], [PdSM + 06] or [START_REF] Saad | h. A deflated version of the conjugate gradient algorithm[END_REF] and references therein), we can say that the principle is to proceed in the following way. From the solve of the first linear system A 1 u 1 = f 1 , a spectral analysis of the Krylov subspace enables the selection of a suitable coarse space denoted Z. More precisely, Ritz eigenvalues associated with the Krylov subspace generated by the solving of the first linear system A 1 δu 1 = . . . by the Krylov method preconditioned by M JSM are computed.

If it exists, the long plateau in the convergence is due to the presence of a few small positive eigenvalues close to zero. By computing the Ritz eigenvectors, we get an approximation to the corresponding eigenvectors. We select the nsmeig small Ritz eigenvectors, denoted (z 1 , . . . , z nsmeig ), whose real parts of the eigenvalues are smaller than a given threshold eig . A rectangular matrix

Z := |z 1 |z 2 | . . . |z nsmeig | is formed out of them.
Then, it is used to build a two-level preconditioner for solving the subsequent linear systems

A k u k = f k , k ≥ 2.

Retrieving approximate eigenvector from GMRES solver

We detail the Ritz eigencomputations after the solve of the first linear system with matrix A 1 and preconditioner M 1 := M -1 JSM ( A 1 ) into m iterations with a GMRES algorithm. The final Krylov subspace associated with the preconditioned operators reads

K m (B 1 , r 0 ) = SPAN r 0 , B 1 r 0 , B 2 1 r 0 , . . . , B m-1 1 r 0 (12)
with B 1 := M -1 1 A 1 . The computational kernel of GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] is the Arnoldi process which computes the orthonormal basis W m = |w 1 |w 2 | . . . | w m | for the Krylov space K m (B 1 , r 0 ). In the orthogonalisation process the scalars h ij are computed so that the square upper Hessenberg matrix H m ∈ R m×m satisfies the fundamental relation

B 1 W m = W m H m + h m+1,m w m+1 e H m = W m+1 H m . (13)
The rectangular upper Hessenberg matrix H m ∈ R (m+1)×m is the square upper Hessenberg matrix H m supplemented with an extra row (0, . . . , 0, h m+1,m ).

From (13) we can derive the following expression for H m :

H m = W H m B 1 W m . ( 14
)
The eigenvalues of H m are called Ritz values and they approximate the eigenvalues of B 1 . We approximate eigenvectors of B 1 by first computing eigenpairs (t i , λ i ) of matrix H m and then compute

z i := W m t i (15)
which is close to an eigenvector of B 1 , for the same eigenvalue λ i . Indeed, we left multiply

W H m B 1 W m t i = λ i t i by W m and get W m W H m B 1 W m t i B 1 W m t i = λ i W m t i .
4.2 A sparse coarse space and the two-level preconditioner

Two coarse spaces

In this paper, we shall focus on the JSM method since there is no need to build a discrete partition of unity as in (4). We use spectral information from the first solve in order to build a coarse space for the subsequent solves. More precisely, from the first solve we select nsmeig Ritz eigenvectors corresponding to the smallest Ritz eigenvalues. The rectangular matrix Z := |z 1 |z 2 | . . . |z nsmeig | stores these Ritz eigenvectors columnwise. Next, we introduce a larger coarse space Z s whose columns are defined by:

Z s i+N (j-1) := R T E,i R E,i z j , for 1 ≤ i ≤ N, 1 ≤ j ≤ nsmeig . ( 16
)
where R E,i is the restriction operator to the i-th subdomain, see (7) and z j is the j-th column of Z. The coarse space Z s is N times as large as Z, dim(Z s ) = N × nsmeig. It has a very sparse structure. For instance, for a three subdomain decomposition and nsmeig = 2, the transformation of Z into Z s has the following form:

Z =   z 11 z 12 z 21 z 22 z 31 z 32   → Z s =   z 11 0 0 z 12 0 0 0 z 21 0 0 z 22 0 0 0 z 31 0 0 z 32   . ( 17
)
where z ij := R E,i z j . If the Ritz eigenvectors are orthogonal, this is not necessarily the case for the columns of Z s . In order to improve stability of the method, we orthogonalize Z s and denote Z s ⊥ the orthogonalized basis of the coarse space. Due to the sparse structure of Z s ⊥ , this process can be performed in parallel, see § 4.2.2 for details.

We now define the coarse corrections that will be used to solve the subsequent linear systems with matrices A k , k ≥ 2. They are inspired by P A-DEF 2 (10) but different. In order to ease notations in the definition of the two-level method, we note

M -1 k := M -1 JSM ( A k ) .
First of all, we build a coarse correction for the preconditioned system

B k := M -1 k A k . ( 18
) Recall that if A k is symmetric, M -1
k is symmetric as well but the extended system A k is not symmetric. As a result B k is not symmetric and we first modify formula (10):

P k := P k + Q k , (19)
where

P k := I -Q k B k , Q k := Z s ⊥ E -1 k Z s ⊥ T , E k := Z s ⊥ T B k Z s ⊥ . (20)
It can easily be checked that we have:

(a) P k = P 2 k ; (b) P k Z s ⊥ = 0, P k Q k = 0; (c) Q k B k = I -P k , Q k B k Z s ⊥ = Z s ⊥ , Q k B k Q k = Q k . We have
Lemma 4.1 The coarse correction P k defined by (19) is invertible and has the following left-filtering property:

Z s ⊥ T B k = Z s ⊥ T P -1 k
Proof. We first prove that P k is one to one and thus invertible. Let u such that P k u = P k u + Q k u = 0. Then, we left multiply by P k and use P k Q k = 0 and P 2 k = P k to obtain P k u = 0 and thus Q k u = 0 as well. From P k u = 0, we have

u = Z s ⊥ E -1 k Z s ⊥ T B k u. Let w = E -1 k Z s ⊥ T B k u, we have u = Z s ⊥ w. Then, from Q k u = 0, we get Z s ⊥ E -1 k Z s ⊥ T Z s ⊥ w = 0. We left multiply by Z s ⊥ T B k to get Z s ⊥ T Z s ⊥ w = 0.
We take the scalar product with w and get Z s ⊥ w = 0 and so u = 0. Let us prove now the filtering property. It is equivalent to

(P T k + Q T k )B T k Z s ⊥ = Z s ⊥ . Consider first the term P T k B T k Z s
⊥ and let's prove it is null:

P T k B T k Z s ⊥ = B T k Z s ⊥ -B T k Q T k B T k Z s ⊥ = B T k Z s ⊥ -B T k Z s ⊥ (Z s ⊥ T B T k Z s ⊥ ) -1 Z s ⊥ T B T k Z s ⊥ = 0 It remains to prove that Q T k B T k Z s ⊥ = Z s ⊥ : Q T k B T k Z s ⊥ = Z s ⊥ (Z s ⊥ T B T k Z s ⊥ ) -1 Z s ⊥ T B T k Z s ⊥ = Z s ⊥
A first method would consist in using P k as a preconditioner to matrix B k or equivalently P k M -1 k as a preconditioner to the extended system A k in a Krylov method. Remark that from the left filtering property, it is easy to check that for any Krylov method with right preconditioning and at any step m the residual r m remains orthogonal to the vector space spanned by the columns of Z s ⊥ . More precisely, we have Lemma 4.2 Let x 0 be an initial guess and r 0 := b -B k x 0 be the initial residual such that Z s

⊥ T r 0 = 0. Let K m (B k P k , r 0 ) denote the Krylov space of dimension m, K m (B k P k , r 0 ) := Span{r 0 , . . . , (B k P k ) m-1 r 0 }. Then, for any x m ∈ {x 0 } ⊕ P k K m (B k P k , r 0 ) we have Z s ⊥ T (b -B k x m ) = 0 .
Proof. First remark that solving B k x = b with P k as a right preconditioner in a Krylov method amounts to solve B k P k y = b with a Krylov method and then apply P k to the resulting approximation to y. Thus, at step m of the algorithm, the approximate solution x m ∈ {x 0 } ⊕ P k K m (B k P k , r 0 ). That is,

x m = x 0 + P k m-1 i=0 α m i (B k P k ) i r 0 for some (α m i ) 0≤i≤m-1 ∈ R m . The corresponding residual reads r m = b -B k x 0 - m-1 i=0 α m i (B k P k ) i+1 r 0 = r 0 - m-1 i=0 α m i (B k P k ) i+1 r 0 .
It suffices to prove by induction that Z s ⊥ T (B k P k ) i r 0 for any 0 ≤ i ≤ m. By assumption, it is true for i = 0. Assume it holds for some i then by using Lemma 4.1, we have

Z s ⊥ T (B k P k ) i+1 r 0 = Z s ⊥ T B k P k (B k P k ) i r 0 = Z s ⊥ T (B k P k ) i r 0 = 0 .
This ends the proof of the lemma. Notice as well that the assumption Z s ⊥ T r 0 is not difficult to satisfy. Indeed, let x 0 be any initial guess, then x 0 := P k (b + P T k x 0 ) satisfies the assumption of the lemma:

Z s ⊥ T r 0 = Z s ⊥ T (b -B k P k (b + P T k x 0 )) = Z s ⊥ T b -Z s ⊥ T (b + P T k x 0 ) = 0 ,
where we have used once again Lemma 4.1 and also that P k Z s ⊥ = 0. Formula (19) demands the computation of matrix E k via formula (20) before solving each linear system with matrix A k , k ≥ 2. With a parallel implementation as explained in § 4.2.2 the cost is almost the same as nsmeig iterations of the one-level Schwarz method. When the size of the coarse space is large, the factorization of the coarse operator E k can be costly. In order to avoid this factorization cost, we reuse matrix E 2 by simplifying formula (20) :

P 2,k := I -Q 2 B k , Q 2 := Z s ⊥ E -1 2 Z s ⊥ T , E 2 := Z s ⊥ T B 2 Z s ⊥ . ( 21
)
The new coarse correction that replaces (19) reads

P 2,k := P 2,k + Q 2 . ( 22
)
Operator P 2,k is a preconditioner to matrix B k . In practice, we will use

C k := P 2,k M -1 k = [(I d + Z s ⊥ E -1 2 Z s ⊥ T ) -Z s ⊥ E -1 2 Z s ⊥ T M -1 k A k ] M -1
k as a preconditioner to the extended system A k . A naive implementation would, at each iteration, demand that one-level Schwarz preconditioner M -1 k be applied twice so that the cost of the two-level preconditioner would be at least double that of the one-level method. In order to avoid this, an equivalent definition is the following:

C k = [(I d + Z s ⊥ E -1 2 Z s ⊥ T ) -Z s ⊥ E -1 2 S T k A k ] M -1 k ( 23
)
where we precompute

S T k := Z s ⊥ T M -1 k . ( 24
)
As explained in § 4.2.2, the cost of precomputing S T k is about nsmeig applications of the one-level Schwarz method.

The corresponding algorithm for solving a series of related linear systems is summed up as follows:

Algorithm 1 For solving a series of linear systems

A k u k = f k , 1 ≤ k ≤ K: 1.
Solve system A 1 ũ1 = f1 with GMRES method preconditioned by M 1 .

Compute the Ritz eigenpairs

(z i , λ i ) for |Real(λ i )| < eig , see § 4.1.1.

2.

Orthogonalize Z s given by formula (16

) into Z s ⊥ . Compute E 2 via formula (21) Gauss factorization of E 2 3. For 2 ≤ k ≤ K Compute S T
k given by formula (24) Solve system A k ũk = fk with GMRES method preconditioned by C k , formula (23).

Numerical results in § 4.2.3 will assert the efficiency of this approach.

Parallel data structure

In this section, we explain that due to a parallel implementation, provided the number of cores is equal to the number of MPI processes, the construction of the coarse space has the same cost as in [GR03]. The only difference is in the size of the matrix E, see § 4.2.1.

The method makes use of three distributed data structures:

• vector

• sparse matrix

• coarse space

Without loss of generality, we take examples for three MPI processes and a three subdomain decomposition. A vector x ∈ R N E is distributed among the three processes, process i (1 ≤ i ≤ 3) will own R E,i x.

According to this vector distribution, the sparse matrix A is stored blockwise, A = ( A ij ) 1≤i,j≤N , see [BFF + 09] for a related data structure. For all 1 ≤ i ≤ N submatrices ( A ij ) 1≤j≤N are owned by the i-th process. Note that since A is sparse, the diagonal blocks A ii , 1 ≤ i ≤ N are sparse as well and the offdiagonal blocks are hypersparse since they correspond to interactions between neighboring subdomains. Thus, the diagonal blocks are stored in a classical CSR format. The off-diagonal blocks are stored in a compressed sparse row and columns format which means only non zero rows and columns are stored in a CSR format, see [START_REF] Buluc | On the representation and multiplication of hypersparse matrices[END_REF] for a related data structure. The matrix-vector product is then naturally parallel. The implementation of the coarse correction demands a distributed storage for sparse rectangular matrices such as Z s ⊥ and S T k (see equation ( 24)). Without entering into details, our storage is similar to that of our sparse matrix. We give some details on the parallel implementation of the coarse space operations. The first one is to split the Ritz eigenvectors following equation ( 16). Actually, since the vectors are already distributed, it is just a matter of pointers management. The next operation is to orthogonalize the split coarse space Z s into Z s ⊥ . Recall that even if the columns of Z are orthogonal that does not necessarily holds for Z s . It is easy to confirm (see formula ( 17)) that it is sufficient to orthogonalize in parallel for 1 ≤ i ≤ N the sub-vectors z ij , 1 ≤ j ≤ nsmeig and that there is no need for communication. Thus the cost of the parallel orthogonalizing of the split coarse space Z s equals the sequential cost of orthogonalizing nsmeig subvectors.

The next step in Algorithm 1 is to compute the matrix

E 2 := Z s ⊥ T M -1 2 A 2 Z s ⊥ = (M -1 2 Z s ⊥ ) T A 2 Z s ⊥ .
We shall see that the cost to compute F 2 := (M -1 2 Z) T A 2 Z for the classical coarse space is comparable to that of computing E 2 , see § 4.2.3 for wall-clock measurements as well. For simplicity, we consider the example given by formula (17): two Ritz eigenvectors and three subdomains. Remember that M is a block diagonal preconditioner so that we have:

M -1 2 Z =    M -1 11 z 11 M -1 11 z 12 M -1 22 z 21 M -1 22 z 22 M -1 33 z 31 M -1 33 z 32    ← Proc 1 ← Proc 2 ← Proc 3 and M -1 2 Z s ⊥ is given by:    M -1 11 z 11 ⊥ M -1 11 z 12 ⊥ M -1 22 z 21 ⊥ M -1 22 z 22 ⊥ M -1 33 z 31 ⊥ M -1 33 z 32 ⊥    ← Proc ← Proc ← Proc
Clearly computational costs per process are the same, they consist of factorizing once the local contributions M ii to subdomain 1 ≤ i ≤ N and then performing several backward-forward substitions. For this we use multithreaded direct solvers inside an MPI process,. As for the computation of AZ and AZ s ⊥ , we take the following example

A Z =     A 11 A 13 A 21 A 22 A 31 A 33        z 11 z 12 z 21 z 22 z 31 z 32    =     A 11 z 11 + A 13 z 31 A 11 z 12 + A 13 z 32 A 21 z 11 + A 22 z 21 A 21 z 12 + A 22 z 22 A 31 z 11 + A 33 z 31 A 31 z 12 + A 33 z 32     and A Z s ⊥ =     A 11 A 13 A 21 A 22 A 31 A 33        z 11 ⊥ z 12 ⊥ z 21 ⊥ z 22 ⊥ z 31 ⊥ z 32 ⊥    =     A 11 z 11 ⊥ A 13 z 31 ⊥ A 11 z 12 ⊥ A 13 z 32 ⊥ A 21 z 11 ⊥ A 22 z 21 ⊥ A 21 z 12 ⊥ A 22 z 22 ⊥ A 31 z 11 ⊥ A 33 z 31 ⊥ A 31 z 12 ⊥ A 33 z 32 ⊥     .
We see that the local matrix vector products computations are the same. The difference is that for AZ we sum local contributions whereas for AZ S ⊥ we store them. The extra memory requirement comes from the contribution of neighbor subdomains through the non zero off-diagonal blocks A ij , i = j. This extra storage cost is proportional to the number of neighbors of a subdomain. It is bounded as the number of subdomains increases. Similar considerations hold for the other operations like for instance computing S T k , for k ≤ 2 from formula (24). Indeed, M k is a block diagonal matrix per subdomain so that computing S T k is an embarrassingly parallel task that can be performed in nsmeig steps. Thus, the cost is smaller than nsmeig iterations of the Schwarz preconditioner. The only major difference between the classical and split coarse spaces lies in the size of matrices E 2 and F 2 . Due to the sparsity of Z S ⊥ , matrix E 2 is a square block sparse matrix of size (N × n nsmeig ) 2 with typically D × N × n 2 nsmeig non zero elements where D is the number of neighbours of a subdomain. Whereas matrix F 2 is full and of size n 2 nsmeig . The application of a direct solver to E 2 will be more costly than for matrix F 2 , see § 4.2.3 for quantitative data. In practice, each process stores one copy of matrix E 2 or F 2 even when we have for more than one subdomain per process. Usually, each process consists of several cores (four in our experiments) so that we can use multithreaded direct solvers.

Numerical performance

All runs were made on a cluster of nodes interconnected by an Infiniband network. Each node is composed of two AMD "MagnyCours" twelve-core (2.1GHz) processors. We use the optimal mapping (so called 'sweat spot') for hybrid (MPI/OpenMP) codes i.e., 4 MPI processes per node. Each MPI process has access to 6 cores, therefore we have used 6 OpenMP threads per MPI process. In practice we map either one or two subdomains per MPI process.

We consider matrices coming from oil reservoir simulations. We compare the one level JSM preconditioner (6) with two two-level methods: the classical one and the new one where the coarse space is split subdomain wise. For all experiments, the overlapping regions are made of two layers of nodes. The first matrix is the pressure block extracted from the well known SPE10 benchmark [START_REF] Christie | Tenth spe comparative solution project: A comparison of upscaling techniques[END_REF], see the convergence curves in Figure 4. This test case is characterised by large permeability variations, 8-12 orders of magnitude. The discretization is made on a regular Cartesian grid with 60 × 220 × 85 cells. The two other series of matrices come from a nonlinear black oil (BO) simulation, see the convergence curves in Figures 5 and6. In all convergence curves, No Coarse Space refers to the one-level JSM preconditioner, (Z) to the classical coarse space and (ZS) to the subdomain wise split coarse space. Notice that in the black oil simulations, even with the two level preconditioners we have some plateaus at the first iterations. This is due to the fact that the two level preconditioners are built using spectral information coming from the first linear system which is different from the subsequent ones. This is not the case for the SPE10 experiment where only the right hand side is changed from one solve to the next ones. As for numerical efficiency they are detailed in table 1 that we comment now giving more information on the test cases.

The first experiment is on the pressure block of the SPE10 (Society of Petroleum Engineers) benchmark problem [START_REF] Christie | Tenth spe comparative solution project: A comparison of upscaling techniques[END_REF]. The matrix comes from three-dimensional reservoir simulation in a highly heterogeneous and anisotropic medium with one million unknowns. In contrast with the next experiments, we have only one matrix and two solves with different right hand sides. The first solve is performed with the one level Jacobi Schwarz preconditioner, see (6) with a domain decomposition into 256 subdomains. The first solve takes 29.6s using 128 MPI processes, two subdomains per MPI process. From this first solve, we build coarse spaces for the second solve. In Table 1, we compare the classical and split coarse spaces. In terms of iterations, the coarse spaces bring huge improvements, especially for the split coarse space. This can be seen as well on the convergence curves, see Figure 4. The matrix is very ill-conditioned so that the coarse space built according to the tolerance eig = 0.1 is quite large for the split coarse space: 52 × 256 = 13, 312 degrees of freedom. The first row of Table 1 shows that it affects the cost of building the coarse space as well the overall speed-up which is close to one in this case. It motivates us to reduce the number of smallest Ritz eigenvalues (nsmeig) to ten. The results are reported in the second row of Table 1. Although the iteration counts are not as good as in the first row, we now have a substantial gain (×5.5) in terms of wall-clock time.

The two other series of matrices, denoted BO, come from a black oil simulation. This computation implies the solving of large-scale nonlinear problems arising from the finite-volume discretization in which the non-linearity is handled by a Newton-Raphson algorithm. Solving these problems leads to a succession of linear problems the solution to which converges towards the solution to the considered problem. In our numerical experiment we consider a series of linear systems extracted at some time step in the simulation. Following the strategy described in § 4.2, we solve the first linear system with a one-level Schwarz algorithm. Then, we reuse the coarse operator built from eigenvectors approximated during the first resolution in the solutions of the linear systems for the remaining Newton-Raphson iterations. To be more precise, we consider two series of five linear systems. In the first serie we have 60 × 60 × 32 grid points and in the second one 120 × 120 × 64. At each linear solve (except for the first one), we use the solution of the previous linear system as an initial guess.

In Table 2, we report results for the black oil test case with 60 × 60 × 32 = 115, 200 grid points. Due to the nonlinearity of the problem, the number of non zero entries is not the same for the various matrices in the series. The average nnz entries is 791, 550. We compare the classical coarse space made of nsmeig Ritz eigenvectors corresponding to small eigenvalues of the preconditioned system with the coarse space introduced in § 4.2 based on a domain wise splitting of the previous coarse space. We have 64 subdomains and we use either 32 or 64 MPI processes. The average iteration counts for the one-level Schwarz method is 54 iterations. The solving time of one system is 2.07s when using 32 MPI processes (two subdomains per MPI process) and 1.23s with 64 MPI processes (one subdomain per MPI process). When we solve the next four linear systems in the Newton-Raphson algorithm by the one-level Schwarz method, the iteration counts (see Figure 5) and solving time slightly decrease due to the use of the previous solution as an initial guess. Therefore we do not report them in the table. In the third and fourth columns, we give the average iteration counts for the next four solves using the classical (denoted by Z) or split (denoted by Z S ) coarse spaces. When using the split coarse space, we gain a factor 2.7 in iteration counts. We automatically selected nsmeig = 7 Ritz eigenvectors whose eigenvalues are smaller than a given threshold eig = 0.1. In column 6, we report the time needed to build and factorize the square matrix E 2 of size 448 × 448 for the split coarse space. Recall that this operation is performed only once for the series of matrices, see Algorithm 1. In the last two columns, we report the average speedup for the next subsequent four solves taking into account all extra costs due to the coarse space: the construction of the coarse space matrix E 2 and at each iteration the cost of precomputing matrices S T k , see (24). Since the classical coarse space is much smaller than the split coarse space if the same threshold eig is used, we also tried to use larger nsmeig's. This yields a classical coarse space which is slightly larger but richer. Results are reported in Table 3. We see that enlarging nsmeig has a small impact on both iteration counts and speedup.

Table 4 is organized as Table 2 except that we consider the large black oil test case with 120 × 120 × 64 grid points. The average nnz entries is 6, 391, 950. We have 128 subdomains and we use either 64 or 128 MPI processes. The solve performed by the one-level Schwarz method require in average 102 iterations. The solving time is 21.05s when using 64 MPI processes and 11.97s with 128 MPI processes. When using the split coarse space, we gain a factor of almost four in iteration counts. We automatically selected nsmeig = 14 small Ritz eigenvectors. The overall speed-up is better than two for the split coarse space Z S . This speed-up computed from wall-clock times is in agreement with the iteration counts. Indeed, the one-level solves for the five matrices take 5 × 102 = 510 iterations. Whereas, the split coarse space solves take 104 iterations for the first (one-level) solve plus 14 iterations for building the coarse operator plus 4 × 25 iterations for the four subsequent two-level solves which means 218 iterations for the whole process. This is about half the total number of iterations for the one-level method. As in the previous test case, we also tried to use larger nsmeig's for the classical coarse space. Results are reported in Table 5 which is organized as Table 3. We see that enlarging nsmeig has a marginal impact on both the iteration counts and wall-clock times.

Conclusion

We studied the solving of linear systems arising from porous media flows simulations with high heterogeneities by Schwarz type methods. We first investigated the influence of the partition into subdomains. Using an AMG type partitioning in Metis or Scotch improves the convergence with almost no extra cost. Then, we introduced two two-level preconditioners based on coarse space corrections. They are algebraic in the sense that they only make use of information generated during the solving of the first linear system. Taking advantage of a parallel implementation, it is possible, at nearly no extra cost, to use larger but sparse coarse space (Z S ) obtained by splitting subdomain wise a classical coarse space. The iteration counts and convergence curves are then always substantially better than with the classical coarse space. As long as the coarse space is not too large, we also have a gain in wall-clock solving times. But when the coarse space is too large, the direct solver used to invert E 2 in (23) takes too much time, see Table 1. Our current solution is to reduce the size of the coarse space by limiting the number of Ritz eigenvectors that span the classical coarse space. Another way to bypass this limitation would be to design iterative methods for this specific type of problem. This requires further investigation.

BO

Figure 1 :

 1 Figure 1: Decomposition into two overlapping subdomains.

Figure 2 :Figure 3 :

 23 Figure 2: Weighted vs. un-weighted partitions for a problem with a strong anisotropy (see equation (9))

Figure 4 :

 4 Figure 4: Convergence curves for SPE10 benchmark -domain decomposition into 256 subdomains

Table 1 :

 1 Comparison of the coarse spaces for SPE10 benchmark with 256 subdomains

	SRA	Iterations		C.S.	Speedup
	Matrix nsmeig	1-lvl sol.	Z	Z S	Z S	Z	Z S
	SPE10	52 TOL 10 FV	399	165 192	28 81	11.19s 2.03 1.37 0.97s 2.63 5.54
	BO	Iterations		C.S.	Speedup
	# MPI	1-lvl sol.	Z	Z S	nsmeig	Z S	Z	Z S
	32 64	54	32	20	7	0.46s 1.09 1.40 0.25s 1.10 1.44

Table 2 :

 2 Comparison of the coarse spaces, 60 × 60 × 32 unknowns for a black oil simulation and 64 subdomains. Speedup refers to wall-clock timings of the overall process.

Table 3 :

 3 Effect of the size of the classical coarse space, 60 × 60 × 32 unknowns for a black oil simulation and 64 subdomains.

		Iterations	Speedup
	# MPI	Z(12)	Z(24)	Z(12)	Z(24)
	32 64	31	30	1.10 1.10	1.10 1.09

Table 4 :

 4 Comparison of the coarse spaces, 120 × 120 × 64 unknowns for a black oil simulation and 128 subdomains. Speedup refers to wall-clock timings of the overall process.

	BO	Iterations	Speedup
	# MPI	Z(28)	Z(56)	Z(28)	Z(56)
	64 128	80	79	1.05 1.05	1.05 1.06

Table 5 :

 5 Effect of the size of the classical coarse space and 128 subdomains, 120 × 120 × 64 unknowns for a black oil simulation

		10 0							No Coarse Space
										ZS (10v)
										ZS (tol, 52v)
		10	-1							Z (10v) Z (tol, 52v)
	|| b -Ax n || 2 / || b || 2	10 10	-3 -2
		10	-4
		10	-5
		10	-6	0	50	100	150	200	250	300	350	400
								Number of iterations

supported by the National Natural Science Foundation of China under grant 10901125,


Sylvie Boldo
email: sylvie.boldo@inria.fr

How to Compute the Area of a Triangle: a Formal Revisit

Keywords: floating-point arithmetic, formal proof, Coq, triangle, underflow

Mathematical values are usually computed using well-known mathematical formulas without thinking about their accuracy, which may turn awful with particular instances. This is the case for the computation of the area of a triangle. When the triangle is needle-like, the common formula has a very poor accuracy. Kahan proposed in 1986 an algorithm he claimed correct within a few ulps. Goldberg took over this algorithm in 1991 and gave a precise error bound. This article presents a formal proof of this algorithm, an improvement of its error bound and new investigations in case of underflow.

I. INTRODUCTION

Floating-point (FP) arithmetic is seen as intricate because too few people have sufficient knowledge to understand how it works. For people having been only trained with mathematics, facts such that (x+y)+z may be different from x+(y +z) for certain values or the fact that there exists x such that x = 0, but x 2 = 0 is beyond comprehension. This is the reason why mathematical formulas are most of the time programmed as they stand in mathematical textbooks.

We are interested here in computing the area of a triangle, given its side lengths as FP numbers. This is especially difficult (using FP computations) for needle-like triangles like the one in Figure 1. The common formula to compute the area is two millennia old and is attributed to Heron of Alexandria:

∆ = s (s -a) (s -b) (s -c)
where s = a + b + c 2 This formula is known to be inaccurate using floatingpoint arithmetic since the 80s. This has been first studied by Kahan [START_REF] Kahan | Miscalculating Area and Angles of a Needle-like Triangle[END_REF]. He gave examples of incorrect computations: either the result was very wrong or the computation was stopped due to a negative square root, created by round-off errors. Kahan also proposed an algorithm that behaves correctly using floating-point arithmetic, that is to use the following formula:

1 4 (a + (b + c)) (c -(a -b)) (c + (a -b)) (a + (b -c))
The parenthesis are not to be removed to guarantee that the square root will carry out on a non negative number. Kahan claimed the safety and the accuracy of the result within a few ulps by giving concise and precise arguments.

A little later, Goldberg presented this algorithm in his famous article "What Every Computer Scientist Should Know About Floating-Point Arithmetic" [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]. More than the algorithm, he gave a precise error bound, but gave no proof or hint of the reason why it behaves so well. He claimed that, given a machine epsilon ε such that ε ≤ .005, the rounding error was at most 11 ε.

To get a high guarantee on mathematical results or programs, formal methods have recently developed. This of course includes floating-point arithmetic that has been formalized since 1989 in order to formally prove hardware components or algorithms [START_REF] Carreño | Specification of the IEEE-854 floatingpoint standard in HOL and PVS[END_REF], [START_REF] Russinoff | A mechanically checked proof of IEEE compliance of the floating point multiplication, division and square root algorithms of the AMD-K7 processor[END_REF], [START_REF] Harrison | Formal verification of floating point trigonometric functions[END_REF]. This algorithm is a good test case for formal proof checking. The first reason is that this is common knowledge, so it is believed both that it cannot be false and that it can hardly be enhanced. A recent example of the computation of the discriminant [START_REF] Kahan | On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic[END_REF] has shown that pen-and-paper proofs may miss difficulties: here the fact that the floating-point test may be mistaken [START_REF] Boldo | Kahan's algorithm for a correct discriminant computation at last formally proven[END_REF]. The second reason is that this is an uncommon kind of proof compared to what has been proved before. It is somewhat easier: this is forward error analysis with few floating-point cunning facts. Formal proofs are usually applied to more complex and trickier algorithms. Anyway, even if this algorithm and these proofs are decade old, we were able to notably improve the error bound.

From the algorithm point of view, this work gives a high guarantee of its correctness and gives precise hypotheses on the radix and the needed precision. We are also sure to take into account second-order terms for the round-off error, which are usually dismissed. Here, they cannot be dismissed and must therefore be bounded, even coarsely. From the formal methods point of view, we will base our proof on the recent Flocq library [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF] that has not yet been thoroughly used for floating-point algorithms. Flocq is a formalization in Coq that offers a multi-radix and multi-precision formalization for various floating-and fixed-point formats (including FP with or without gradual underflow) with a comprehensive library of theorems. Its usability and practicality have to be established against test-cases.

Another point is the difficulty in handling subnormals. Most pen-and-paper proofs assume there is neither underflow, nor overflow. This is a very strong hypothesis that greatly simplifies the proofs. Unfortunately, it may be difficult to give conditions beforehand or to check afterwards that this hypothesis is fulfilled. This is the reason why we take care of subnormals and their consequences for round-off errors.

This article is organized as follows. Section II presents the proof when no underflow occur, the improvement of the error bound and the formal demonstration. Section III presents the modification when taking gradual underflow into account. Section IV presents the proved C program with annotations stating its precise specifications, including overflow considerations.

Notations:

The side lengths of the triangle will be denoted by a, b and c and are assumed to be exact FP numbers. We have ordered these lengths beforehand, so that 0 ≤ c ≤ b ≤ a. As these values represent a (possibly degenerate) triangle, basic geometry requires that a ≤ b + c.

We will denote the radix by β, an integer greater than 1. The precision of the floating-point format will be denoted by p and will be greater than 1. The minimal exponent will be denoted by E i : this means that the smallest positive FP number (the smallest subnormal) is β Ei . We will denote by • the default rounding mode (rounding to nearest, ties to even), and we will denote by ⊕ the FP addition, ⊖ the FP subtraction, and ⊗ the FP multiplication. We will denote by ε the machine epsilon, that is to say the relative error bound for normal numbers in rounding to nearest, which is ε = β 1-p 2 . We will denote by C(x) (for correct) the exact value of a floating-point x, meaning the value it would have had without any rounding. In particular, C(a) = a, and C(a ⊕ b) = a + b.

All the theorems stated in this article correspond to one or several Coq theorems. This development is available on http://www.lri.fr/ ∼ sboldo/research.html.

Algorithm

For the sake of readability, here is Kahan's algorithm with temporary variables:

t 1 = a ⊕ (b ⊕ c) t 2 = a ⊕ (b ⊖ c) t 3 = c ⊕ (a ⊖ b) t 4 = c ⊖ (a ⊖ b) M = ((t 1 ⊗ t 2 ) ⊗ t 3 ) ⊗ t 4 ∆ = • • 1 4 • ( √ M )
Note that Kahan's algorithm computes t 1 ⊗ t 2 ⊗ t 3 ⊗ t 4 without parenthesis, but we chose to parenthesize the multiplication for reasons explained later in Section III-C.

II. WHEN NO UNDERFLOW OCCURS

A. Hypotheses

We first consider a floating-point format on p bits with unbounded exponent range. It exactly corresponds to the very common "provided no underflow, or overflow occur". This greatly simplifies the proof for a beginning, and helps to get a tight error bound without having more to take subnormal into account. This corresponds to the FLX format defined in Flocq. A value is said to fit in the format when it can be represented by a floating-point number, and is therefore computed without rounding error.

Even if we consider a generic radix, we have in thought that β will be 2 or 10 (or even 4 or 16), therefore we have assumed that 1 4 fits in the format. We also proved this was correct for β = 2 and β = 10 as soon as p ≥ 2.

We also assume that the precision is not too small. More precisely, we require that the machine epsilon, that is to say ε = β 1-p 2 is smaller or equal to 1 100 . This is guaranteed for a precision greater than 6 with β = 2 and greater than 2 with β = 10. This is summarized as follows:

Set of Hypotheses 1 We assume in this section:

• the exponent range is unbounded,

• 1 4 fits in the format, • ε ≤ 1 100 , • 0 ≤ c ≤ b ≤ a ≤ b + c.

B. Non-negativity and Exact Subtraction

First, we will prove that the computation will not fail due to taking the square root of a negative number. This requires to prove that M ≥ 0. As roundings are monotone, it is sufficient to prove that all t i s are non-negative. The monotonicity of the rounding will be used thoroughly and is sufficient to prove all the required inequalities, given the assumptions on a, b, and c. For example, to prove that t 4 ≥ 0, it is sufficient to prove that c -(a ⊖ b) ≥ 0, which is equivalent to c ≥ a ⊖ b. By the monotonicity of the rounding and as c fits in the format, it is sufficient that c ≥ ab which is exactly the assumption a ≤ b + c.

One of the key point in the error bound proof is the fact that a ⊖ b = ab. This subtraction could have created a large relative error, but is in fact exact due to the assumptions on the inputs. More precisely, we use Sterbenz theorem on exact subtraction [START_REF] Sterbenz | Floating point computation[END_REF]. There is left to prove that

b 2 ≤ a ≤ 2 • b. We know that b ≤ a and that a ≤ b + c ≤ 2 • b as b ≤ c. So a ⊖ b is computed without error.

C. Error Lemmas

All the remaining operations may not be exact and we have to bound their relative errors to get a final relative error. We will denote by err(x, y, e) the mathematical inequality |x -y| ≤ e • |y|. It will mean that x is an approximation of the exact value y with a relative error e.

The first two theorems may seem silly. Their use will be explained in Section II-D.

Theorem 1 (err aux) Given x, y, e 1 , e 2 , if e 1 ≤ e 2 and err(x, y, e1), then err(x, y, e2).

Theorem 2 (err 0) Given x, we have err(x, x, 0).

We will now express floating-point operations with this error. For now, this is naive forward error analysis.

Theorem 3 (err init flx) We assume the set of hypotheses 1. Given x, we have err(•(x), x, ε).

Remember we assume an unbounded exponent range. For all real numbers, the error can be seen as a relative error bound bounded by ε. Note this was an application of a standard theorem from Flocq.

Theorem 4 (err add flx) We assume the set of hypotheses 1. Given x 1 , y 1 , e 1 , x 2 , y 2 , e 2 , if err(x 1 , y 1 , e 1 ), and err(x 2 , y 2 , e 2 ), and 0 ≤ y 1 and 0 ≤ y 2 , then

err(x 1 ⊕ x 2 , y 1 + y 2 , ε + (1 + ε) • max(e 1 ,

e 2 )).

There is nothing new here: this is a typical theorem from forward analysis. Given the fact that y 1 and y 2 are nonnegative, the error of an addition is ε plus the maximum of the errors of the inputs multiplied by 1 + ε. The proof is straightforward.

Theorem 5 (err mult flx) We assume the set of hypotheses 1. Given x 1 , y 1 , e 1 , x 2 , y 2 , e 2 , if err(x 1 , y 1 , e 1 ), and err(x 2 , y 2 , e 2 ) then

err(x 1 ⊗ x 2 , y 1 • y 2 , ε + (1 + ε) • (e 1 + e 2 + e 1 • e 2 )).
There is nothing new here either. We take care of the second order terms e 1 • e 2 even if it will probably be negligible. The proof is also straightforward.

Theorem 6 (err sqrt flx) We assume the set of hypotheses 1. Given x, y, e, if 0 ≤ y, and e ≤ 0.5 and err(x, y, e), then err(

•( √ x), √ y, ε + (1 + ε) • 5 8 • e).
This is the theorem that allows us to improve over Goldberg's bound. Probably, his point of view was to consider that | √ 1 + h -1| ≤ |h| for |h| ≤ 1 and this gives an error bound which is ε + (1 + ε) • e and gives 11 ε at the end. But in fact, | √ 1 + h -1| ≈ |h| 2 but may be greater for negative h. Therefore, we assume h is small, and prove that for all |h| ≤ 0.5, we have

| √ 1 + h -1| ≤ 5 8 • |h|.
The proof is only based on interval arithmetic using bisections. This could have been very tedious using Coq, but the proof was lightened by the interval tactic [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF] that was able to solve the cases for positive and negative h. We unfortunately had to split the cases and remove the absolute values by hand.

D. Main Proof

We now have all the necessary theorems to do forward error analysis in Coq on this algorithm. This was expected to be very tedious as formulas get rather complicated, but this was rather easy. The idea is to use the eapply Coq tactic that applies partly a theorem: this means that some variables remain uninstantiated for some time during the proof, and are noted with a ? followed by a number. For example to prove x ≤ y, we first prove that x ≤ ?1375 and then that ?1375 ≤ y. Here, we want to bound the error of the computation of M . As we do not know beforehand the error of M and we do not want to compute it exactly by hand, we use the err aux theorem to get an unknown value as the error bound of M . We then have to prove that err(M, C(M ), ?1352). Then, as M is the result of a multiplication, we apply the err mult theorem. So the unknown value is partly instantiated: the error bound of M is ε + (1 + ε) • (?1359 + ?1360 + ?1359 • ?1360)) and we have to prove that err(t4, c -(ab), ?1360) and that err((t

1 ⊗t 2 )⊗t 3 , (C(t 1 )•C(t 2 ))•C(t 3 )), ?1359).
Step by step, depending on the last FP operation, we instantiate unknown values, and sometimes create them. We use Theorems err init and err 0 to solve the simplest goals and we take advantage of the fact that a ⊖ b = ab. At the end, we get a large formula for the error of M :

ε + (1 + ε) • (ε + (1 + ε) • (ε + (1 + ε) • (2 • ε + ε • ε+ (2 • ε + ε • ε) + (2 • ε + ε • ε) • (2 • ε + ε • ε)) + ε + (ε+ (1 + ε) • (2 • ε + ε • ε + (2 • ε + ε • ε) + (2 • ε + ε • ε)• (2 • ε + ε • ε))) • ε) + ε + (ε + (1 + ε) • (ε + (1 + ε)• (2 • ε + ε • ε + (2 • ε + ε • ε) + (2 • ε + ε • ε) • (2 • ε + ε • ε)) +ε + (ε + (1 + ε) • (2 • ε + ε • ε + (2 • ε + ε • ε)+ (2 • ε + ε • ε) • (2 • ε + ε • ε))) • ε)) • ε)
but we did not have to give it directly to Coq, as was expected! The value is indeed needed to be exactly put in the right form for theorems to be applied.

We can then make the prover reorganize the formula and we get:

ε 9 +9•ε 8 +36•ε 7 +84•ε 6 +126•ε 5 +126•ε 4 +84•ε 3 +36•ε 2 +9•ε
We know that ε is small and that the terms with the highest exponents will be negligible, but proving this was cumbersome. As we required ε to be smaller than 0.01, we were able with some effort to bound this value by 9ε + 37ε 2 .

Here is an idea of the respective sizes. To prove that the very long formula was the error bound of M from the previous lemmas, the proof is 34 lines long (this includes the errors of all the t i ). To prove that this value is smaller than 9ε+37ε 2 , it takes 80 lines, more than twice the length of the "intelligent" proof.

Theorem 7 (err M flx) We assume the set of hypotheses 1. We have err(M, C(M ), 9ε + 37ε 2 ).

Note that the order of the t i does not matter here: we may have chosen M = ((t 2 ⊗ t 4 ) ⊗ t 1 ) ⊗ t 3 , it would have given the same error bound.

E. Correctness Theorem with unbounded exponent range

To end the proof, we just have to go on with the forward analysis with the square root computation and the multiplication by 1 4 . This last value is computed correctly as it is assumed to fit in the format. Thus the bound on the round-off error of ∆: Theorem 8 (err ∆ flx) We assume the set of hypotheses 1. We have

err(∆, C(∆), 61 8 ε + 36ε 2 ).
Instead of the 11 ε, we are able to formally guarantee that the relative error is 7.625 ε (plus the second-order terms).

We can still tighten this bound in radix 2: in this case, multiplying by 1 4 is exact, therefore the last computation does not create any round-off error, and we can prove: Theorem 9 (err ∆ flx radix2) We assume the set of hypotheses 1 and β = 2. We have

err(∆, C(∆), 53 8 ε + 29ε 2 ).
Instead of the 11 ε, we formally guarantee that the relative error is 6.625 ε (plus the second-order terms). This means a 40 % better bound on the relative round-off error.

III. TAKING GRADUAL UNDERFLOW INTO ACCOUNT

Unfortunately, the exponent range is limited: the IEEE-754 standard [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] precisely states what are the minimal and maximal exponents allowed in the binary32 and binary64 formats. In this Section, we will only consider gradual underflow. Overflows will be looked into in the next Section. We will have a minimal exponent and take into account the fact that subnormal results may appear and produce a huge relative error bound. This corresponds to the FLT format of Flocq.

A. Hypotheses

Set of Hypotheses 2 We assume in this section:

• gradual underflow with E i as minimal exponent,

• E i ≤ -3 -p,
• no upper bound on the exponent,

• 1 4 fits in the format, • ε ≤ 1 100 , • 0 ≤ c ≤ b ≤ a ≤ b + c.
The last four hypotheses were assumed in the previous Section. The second hypothesis is true in all reasonable formats. It only states that subnormal numbers are small and is equivalent to have the smallest positive normal number smaller or equal to β -4 .

B. Error Lemmas

First and for the same reason as before, we have proved that M ≥ 0 and that a ⊖ b = ab. These proofs are based on the monotonicity of the rounding and on Sterbenz theorem, which are still valid with gradual underflow.

We are of course trying to use as much as possible the previous formal proofs, but some theorems are not valid any more (for example err init). Here are the changes when subnormal may appear.

For the addition, it is well-known that, if the result is subnormal, it is exact. Therefore the ε bound also holds here. Thus we exactly have the same formula as in Theorem err add flx: Theorem 10 (err add flt) We assume the set of hypotheses 2. Given x 1 , y 1 , e 1 , x 2 , y 2 , e 2 , if err(x 1 , y 1 , e 1 ), and err(x 2 , y 2 , e 2 ), and 0 ≤ y 1 and 0 ≤ y 2 , then

err(x 1 ⊕ x 2 , y 1 + y 2 , ε + (1 + ε) • max(e 1 , e 2 )).
Concerning the multiplication, we cannot dodge the subnormals as with the addition. In the case where the output of the FP multiplication is subnormal, the relative error bound gets huge. To prevent this, we will have to prove that gradual underflow will not happen to guarantee the wanted error bound. We recall that β Ei is the smallest positive subnormal number.

Theorem 11 (err mult flt) We assume the set of hypotheses 2. Given x 1 , y 1 , e 1 , x 2 , y 2 , e 2 , if x 1 and x 2 fit in the format, if err(x 1 , y 1 , e 1 ), and err(x 2 , y 2 , e 2 ), and

if β Ei+p-1 < |x 1 ⊗ x 2 |, then err(x 1 ⊗ x 2 , y 1 • y 2 , ε + (1 + ε) • (e 1 + e 2 + e 1 • e 2 )).
As for the square root, there is nearly no subnormal problem. If x is in the format, its square root cannot be subnormal, except for zero. When x = 0, the next theorem is not valid, therefore we put an hypothesis to prevent that, knowing it will be easy to prove in our context. We also could have proved that x = 0 was sufficient.

Theorem 12 (err sqrt flt) We assume the set of hypotheses 2. Given x, y, e, if 0 ≤ y, and e ≤ 0.5, and β Ei+p-1 < •( √ x), and err(x, y, e), then

err(•( √ x), √ y, ε + (1 + ε) • 5 8 • e).
We have all the basic blocks to do forward analysis on a FP format with gradual underflow. But, as far as the multiplication are concerned, we have to prove that no subnormal will appear.

C. Ordering the t i s

As explained, we will need to prove that no subnormal is created. This means we need either to know beforehand or to detect afterwards if subnormals appear. In this case, the detection afterwards was a better solution. The reason is that the hypothesis needed on a, b and c to guarantee this fact would have been very strong, and much stronger than the chosen hypothesis on the result. More precisely, we give a limit value for the result of the algorithm: if it is above this limit, we are sure that no subnormal was created and that the error bound holds. If it is under the limit, a subnormal may have appeared and the relative error bound may be much greater than the expected bound.

We then need to detect any underflow. As the multiplications are the only problems, we have to look into the computation of M and we may try to organize the t i as we want. As explained before, the original algorithm by Kahan did not put any parenthesis and was stated as t 1 ⊗ t 4 ⊗ t 3 ⊗ t 2 while we choose M = ((t 1 ⊗t 2 )⊗t 3 )⊗t 4 . The reason is that we ordered the t i by magnitude. More precisely, we proved that:

Theorem 13 We assume the set of hypotheses 2. We have

0 ≤ t 4 ≤ t 3 ≤ t 2 ≤ t 1 .
This is useful as it means that a subnormal result will not disappear. More precisely, if the result of a multiplication is a subnormal, the results of the following FP multiplications will also be subnormal. This means that, if a subnormal appear, then M will be a subnormal, which is very easy to check. The idea is to forbid the cases where an internal multiplication creates a subnormal, but this fact cannot be detected at the end of the computation.

For that, we will use the following theorem:

Theorem 14 (subnormal aux) We assume the set of hypotheses 2. Given x and y, we assume that x fits in the format and that β Ei+p-1 < |x ⊗ y|. We also assume that,

if |x| ≤ 1, then |y| ≤ 1. Then β Ei+p-1 < |x|.
The idea of the proof is quite simple. By the absurd, we assume that |x| ≤ β Ei+p-1 . From the hypothesis on E i , it means that |x| ≤ 1 and so that |y| ≤ 1. Then |x • y| ≤ β Ei+p-1 • 1 = β Ei+p-1 . As β Ei+p-1 fits in the format and by monotonicity of the rounding, we have that |x ⊗ y| ≤ β Ei+p-1 which is absurd.

This theorem allows us to prove that, if M is normal, then (t 1 ⊗ t 2 ) ⊗ t 3 is also normal, which also implies that t 1 ⊗ t 2 is normal. then all the multiplication results are normal and the err mult flt Theorem can be applied. In this proof, there are several other goals, namely the ones corresponding to "if |x| ≤ 1, then |y| ≤ 1". They are indeed straightforward as the t i are ordered.

D. Correctness Theorem with Gradual Underflow

We can now apply the same kind of proof as in Section II-D with forward error analysis, helped by the proof assistant. We of course get the same error bound, provided M is not a subnormal:

Theorem 15 (err M flt) We assume the set of hypotheses 2 and that β Ei+p-1 < M . We have err(M, C(M ), 9ε + 37ε 2 ).

To end the proof, we just have to go on with the forward analysis with the square root computation and the multiplication by 1 4 . This last value is computed correctly as it is assumed to fit in the format. Thus the bound on the round-off error of ∆, provided ∆ is big enough (so that M is big enough):

Theorem 16 (err ∆ flt) We assume the set of hypotheses 2 and that 1 4 β

⌈ E i +p-1 2 ⌉ < ∆. We have err(∆, C(∆), 61 8 ε + 36ε 2 ).
Instead of the 11 ε, we are able to formally guarantee that the relative error is 7.625ε (plus the second-order terms) while taking into account gradual underflow.

We can still tighten this bound in radix 2: in this case, multiplying by 1 4 is exact, therefore the last computation does not create any round-off error, as we are far from the underflow threshold:

Theorem 17 (err ∆ flt radix2) We assume the set of hypotheses 2, that β = 2, and that 2

⌈ E i +p-1 2 ⌉-2 < ∆. We have err(∆, C(∆), 53 8 ε + 29ε 2 ).
Instead of the 11 ε, we formally guarantee that the relative error is 6.625 ε (plus the second-order terms) while taking into account gradual underflow.

/ * @ requires 0 <= x; @ ensures \result==\round_double(\NearestEven,\sqrt(x)); @ * / double sqrt(double x); / * @ logic real S(real a, real b, real c) = @ \let s = (a+b+c)/2; @ \sqrt(s * (s We use the Frama-C platform 1 to perform formal verification of C programs at the source-code level. Frama-C is an extensible framework that combines static analyzers for C programs, written as plug-ins, within a single tool. In this work, we use the Jessie plug-in for deductive verification. C programs are annotated with behavioral contracts written using the ANSI C Specification Language [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language, version 1.5[END_REF] that tries to be as near C statements as possible. The Jessie plugin translates them to the Why3 verification platform [START_REF] Bobot | Why3: Shepherd your herd of provers[END_REF]. Finally, the Why3 platform computes verification conditions from these programs, using traditional techniques of weakest preconditions, and emits them to a wide set of existing theorem provers, ranging from interactive proof assistants to automated theorem provers. In this work, we use the Coq proof assistant, and the automated theorem prover Gappa that uses interval arithmetic to prove properties that occur when verifying numerical applications [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF].

The full annotated program is in Figure 2. Here are some details about the annotations. We only consider the double type meaning the binary64 type of the IEEE-754. First, the square root is defined as an external function with a specification: it requires the input to be non negative and produces the rounding to nearest of the exact square root. Then S is the mathematical exact value of the area of the triangle, 1 http://www.frama-c.cea.fr/ computed with Heron's formula. We then require the inputs of the function are such that 0 ≤ c ≤ b ≤ a ≤ b + c as explained before.

We also require that a ≤ 2 255 . The reason is to prevent overflows. This is sufficient to guarantee that no operation will overflow. It may seem a strong hypothesis, but if you consider a = b = c = 2 256 , then M is the rounding of 3 • 2 1024 and thus overflows.

The last annotation is what the function ensures, meaning what it guarantees: if the result is greater than 2 -513 , then the relative error is smaller than 53 8 2 -53 + 29 • 2 -106 . Let us now detail the proofs. There are 3 kinds of proofs. The first one is the precondition of the square root function that requires the input to be non negative. This was already proved in Coq and we just had to plug the given proof. The second kind concerns the overflows. All those were automatic thanks to Gappa: the hypothesis a ≤ 2 255 and the facts that 0 ≤ c ≤ b ≤ a were sufficient for Gappa to prove no exceptional behavior (infinities here) will occur.

The last proof is the one of what the function ensures. First we wanted to compare our algorithm to Heron's formula, so we first prove that, for all real numbers a, b and c

, if s = a+b+c 2 , then s (s -a) (s -b) (s -c) = 1 4 (a + (b + c)) (c -(a -b)) (c + (a -b)) (a + (b -c))
. This do not require any hypothesis of the non negativity of the input as the square root would give the same answer if the inputs were negative. The proof is straightforward using the ring tactic.

For the round-off error, we use the err ∆ flt radix2 theorem. We have several things to prove in terms of precision, minimal exponent,. . . The hypothesis 2 -513 < ∆ is of course exactly the 2 ⌈ E i +p-1 2 ⌉-2 required by the theorem, as p = 53 and E i = -1074 in binary64.

All proof obligations were done either with Coq or with Gappa as shown in Figure 3. We proved the C program respects its specification, and that it will not fail, due to overflow or to a negative square root.

V. CONCLUSION

With this case study, we have shown several interesting facts. The first one is that Kahan's algorithm for computing the area of a triangle is correct, and that its error bound is much better than what Goldberg gave in [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]. The second one is that the formal proof of this fact was not as cumbersome as expected, as features from Coq were really helpful to get the error bound without much effort. Unfortunately, bounding the higher order terms was tedious and this should be automatized in the future. These proofs have also shown that the Flocq library was both complete (all needed theorems were there) and helpful. We did not find useless goals to prove or tedious stating of theorems, as can be found sometimes in the standard library of Coq.

Another fact is that taking underflow and overflow into account gives only a small range of inputs where this program is correct. As soon as a is greater than 2 255 or if the result is smaller than 2 -513 , we do not guarantee anything. This correct range was much smaller than what we expected. As for the proofs, overflow was very easy as it is entirely handled by Gappa. But managing possible underflows was rather complex. It lead us to reorder the computations to be sure to get a hint that a subnormal did appear. More generally, it does not seem easy to give a recipe to handle gradual underflow and each example seem to come with its own subtlety that helps to prove it, but that cannot be applied to another.

There are a few generalizations of this work. If anyone ever needs another radix than 2 or 10 or a power of 2, the proof can be re-run to get a slightly increased error bound taking into account the fact that 1 4 does not fit in the format. As for the other rounding to nearest, this could probably apply to rounding to nearest, ties away from zero. The proof was unpractical with rounding to nearest with an arbitrary tie, as it causes problems with the symmetry of the rounding: we cannot prove that •(|x|) = | • (x)| so we decided this complexity was not worth it and we chose ties to even.

Another perspective would be to consider the side lengths as real numbers. They will therefore come as FP numbers with an error. The previous algorithm is not safe in this case as it may lead to take the square root of a negative number. The reason is that the inequality a ≤ b + c will hold on the real side lengths, but not always for their roundings. A long-term perspective is to consider all (or most?) the algorithms from the floating-point literature and formally prove them. It could be done under the common assumptions, that is to say no underflow and no overflow. But it would be more interesting to handle gradual underflow by either giving constraints for subnormal not to appear or by giving correct results even in this case. Overflow has also to be considered, but our experience shows that constraints on initial values are usually enough. Going from well-known facts to formally proved facts would be decisive step towards a high guarantee of our scientific results.

Fig. 1 .

 1 Fig. 1. A Needle-Like Triangle

  This work was supported by the VERASCO project (ANR-11-INSE-003) of the French National Agency for Research (ANR).

Fig. 2 .

 2 Fig. 2. Annotated and Proved C program for the computation of the area of a triangle

Fig. 3 .

 3 Fig.3. Snapshot of Why3Ide: all goals were proved (green dots). The hidden goals concern overflows and are proved using Gappa. The shown goals are the only one proved using Coq and their compilation time. The first one is the post-condition of the function while the last one is the non-negativity of the square root (precondition of the square root function).

  Consider a = 1 + 3 • 2 -53 , b = 1 and c = 3 • 2 -53 . In binary64, the value a will be rounded in ã = 1 + 4 • 2 -53 while b and c are unchanged. Then,•(c -•(ãb)) = -2 -53 and the algorithm will fail. Another algorithm should be created to also handle these cases.