Farah Benmouhoub

Pierre-Loic Garoche
email: pierre-loic.garoche@enac.fr

Matthieu Martel

An Efficient Summation Algorithm for the Accuracy, Convergence and Reproducibility of Parallel Numerical Methods

Keywords: floating-point arithmetic, accurate summation, numerical accuracy, numerical methods, convergence, reproducibility

Nowadays, parallel computing is ubiquitous in several application fields, both in engineering and science. The computations rely on the floating-point arithmetic specified by the IEEE754 Standard. In this context, an elementary brick of computation, used everywhere, is the sum of a sequence of numbers. This sum is subject to many numerical errors in floating-point arithmetic. To alleviate this issue, we have introduced a new parallel algorithm for summing a sequence of floating-point numbers. This algorithm which scales up easily with the number of processors, adds numbers of the same exponents first. In this article, our main contribution is an extensive analysis of its efficiency with respect to several properties: accuracy, convergence and reproducibility. In order to show the usefulness of our algorithm, we have chosen a set of representative numerical methods which are Simpson, Jacobi, LU factorization and the Iterated power method.

Introduction

Scientific computing relies heavily on floating-point arithmetic as defined by the IEEE754 Standard [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF][START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF][START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. It is therefore sensitive to round-off errors, and this problem tends to increase with parallelism. In floating-point computations, in addition to rounding errors, the order of the computations affects the accuracy of the results. For example, let us calculate in IEEE754 single precision (Binary32) the sum of three values x, y and z, where x = 10 9 , y = -10 9 and z = 10 -9 . We obtain ((x + y) + z) = ((10 9 -10 9 ) + 10 -9 ) = 10 -9 ,

(x + (y + z)) = (10 9 + (-10 9 + 10 -9 )) = 0.

(2)

Equations ( 1) and [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF] show that for the same mathematical expression, a sum of three operands, different orderings of the computations yield different results.

In the floating-point arithmetic, we note that, for the same values of x, y and z, and for the same arithmetic operation, we obtain two different results because of parsing the three values differently. In fact, many summation algorithms exist in the literature. Some of them are based on compensated summation methods [START_REF] Kahan | A survey of error analysis[END_REF][START_REF] Malcolm | On accurate floating-point summation[END_REF][START_REF] Ogita | Accurate sum and dot product[END_REF][START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Fast high precision summation[END_REF][START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF] with or without the use of the error-free transformations to compute the error introduced by each accumulation. Others are based on manipulating the exponent and/or the mantissa of the floating-point numbers in order to split data before starting computations [START_REF] Demmel | Accurate floating point summation[END_REF][START_REF] Demmel | Accurate and efficient floating point summation[END_REF].

In a similar approach, we proposed in [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF] a new algorithm for accurately summing n floating-point numbers. This algorithm performs computations only within working precision, requiring only an access to the exponents of the values. The idea is to compute the summands according to their exponents without increasing the complexity. More precisely, the complexity of the algorithm is linear in the number of elements, just like the naive summation algorithms. The main contribution of the present article is to show that this algorithm improves simultaneously the parallel execution time, the reproducibility and convergence of computations through the increase of their numerical accuracy as follows:

1. numerical accuracy improvement is illustrated on computations of MPI implementations of Simpson's rule and the LU factorization method. 2. convergence acceleration is showcased on both Jacobi's method and the iterated power method compared to versions of these methods which use a simple summation algorithm. Past results [START_REF] Damouche | Impact of accuracy optimization on the convergence of numerical iterative methods[END_REF] show that improving the accuracy of computation also leads to accelerate the convergence of iterative sequential algorithms. Our motivation is therefore to parallelize these two methods focusing first on accuracy and obtaining, as a side effect, a better convergence. 3. Last, reproducibility of numerical computation in the context of parallel summation is supported by the reduction on numerical errors. Indeed, the combination of the non-associativity of floating-point operations like addition and computations done in parallel may affect reproducibility. The intuitive solution used to ensure reproducibility is to determine a deterministic order of computation. Another method is based on reduction or elimination of round-off errors, i.e. by improving the numerical accuracy of computations that we will further see in this article. This is illustrated on Simpson's rule and a simple matrix multiplication.

A last contribution is an experimental comparison of the execution time of our algorithm with respect to the similar approach proposed by Demmel and Hida [START_REF] Demmel | Accurate floating point summation[END_REF]. Indeed, Demmel and Hida algorithm [START_REF] Demmel | Accurate floating point summation[END_REF] has a time complexity of O(n log n) because of an additional sorting step, compared to our summation algorithm which involve no explicit sorting and has a complexity in O(n). Our solution has a greater space complexity but this can be addressed with sparse datastructures. Execution time of the two approaches are compared in Sections 3.1 and 3.2.

The rest of the paper is organized as follows. Section 2 recalls elements of floating-point arithmetic and the related work on some existing summation methods proposed to improve the numerical accuracy of computations. We also present our parallel summation algorithm. Section 3 focuses on the improvement of numerical accuracy, based on two experiments, Simpson's rule and a LU factorization. Section 4 focuses on the impact of accuracy on the convergence speed. It is based on experiments on Jacobi's method and the Iterated Power Method. Last, we focus on reproducibility in Section 5, based on two experiments: the Simpson's rule and a matrix multiplication. We conclude in Section 6.

Background

This section introduces some useful notions used in the remainder of this article. Section 2.1 provides some background of the floating-point arithmetic as defined by the IEEE754 standard [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF][START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF][START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. Section 2.2 discusses related work. Our algorithm introduced in [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF], is presented in Section 2.3.

Floating-Point Arithmetic

Following the IEEE754 Standard, a floating-point number x in base β is defined by

x = s • m • β exp-f +1 (3)
where

-s ∈ {-1, 1} is the sign, -m = d 0 d 1 ....d f -1 is the mantissa with digits 0 ≤ d i < β, 0 ≤ i ≤ f -1, -f is the precision, -exp is the exponent with exp min ≤ exp ≤ exp max .
The IEEE754 Standard defines binary formats with some particular values for f , exp min and exp max which are summarized in Table 1. Moreover, the IEEE754 Standard defines four rounding modes for elementary operations over floatingpoint numbers. These modes are towards +∞, towards -∞, towards zero and towards nearest, denoted by • +∞ , • -∞ , • 0 and • ∼ , respectively. The behavior of the elementary operations ∈ {+, -, ×, ÷} between floatingpoint numbers is given by

v 1 • v 2 = •(v 1 v 2 ) ( 4
)
where • denotes the rounding mode such as

• ∈ {• +∞ , • -∞ , • 0 , • ∼ }. By Equa- tion (4)
, we illustrate that, in floating-point computations, performing an elementary operation • with rounding mode • returns the same result as the one obtained by an exact operation , then rounding the result using •. The IEEE754 Standard also specifies how the square root function must be rounded in a similar way to Equation (4) but does not specify the round-off of other functions like sin, log, etc. In this article, without loss of generality, we consider that β = 2. We assume the rounding mode to the nearest. In floating-point computations, absorption and cancellation may affect the numerical accuracy of computations. An absorption occurs when adding two floating-point numbers with different orders of magnitude. The small value is absorbed by the large one. A cancellation occurs when two nearly equal numbers are subtracted and the most significant digits cancel each other.

Related Work

Summation of floating-point numbers is one of the most basic tasks in numerical analysis. Research work has focused on improving the numerical accuracy [START_REF] Demmel | Accurate floating point summation[END_REF][START_REF] Demmel | Accurate and efficient floating point summation[END_REF][START_REF] Kahan | A survey of error analysis[END_REF][START_REF] Malcolm | On accurate floating-point summation[END_REF][START_REF] Ogita | Accurate sum and dot product[END_REF][START_REF] Rump | Ultimately fast accurate summation[END_REF] or reproducibility [START_REF] Demmel | Parallel reproducible summation[END_REF] of the computations involving summations. There are many sequential and parallel algorithms for this task. Surveys of them being presented in [START_REF] Higham | The accuracy of floating point summation[END_REF][START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. Floating-point summation is often improved by compensated summation methods [START_REF] Kahan | A survey of error analysis[END_REF][START_REF] Malcolm | On accurate floating-point summation[END_REF][START_REF] Ogita | Accurate sum and dot product[END_REF][START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Fast high precision summation[END_REF][START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF] with or without the use of error-free transformations to compute the error introduced by each accumulation. We detail some of the compensated summation algorithms further in this section. The accuracy of summation algorithms can also be improved by manipulating the exponent and the mantissa of the floating-point numbers in order to split data before starting computations [START_REF] Demmel | Accurate floating point summation[END_REF][START_REF] Demmel | Accurate and efficient floating point summation[END_REF]. This approach is the one employed by our algorithm and it is explained in details in Section 2.3.

Compensated summation methods:

The idea is to compute the exact rounding error after each addition during computations [START_REF] Kahan | A survey of error analysis[END_REF]. Compensated summation algorithms accumulate these errors and add the result to the result of the summation. The compensation process can be applied recursively yielding cascaded compensated algorithm. Malcolm [START_REF] Malcolm | On accurate floating-point summation[END_REF] describes cascading methods based on the limited exponent range of floating-point numbers. He defines an extended precision array e i where each component corresponds to an exponent. To extract and scale the exponent, Malcolm uses an integer division, without requiring the division to be a power of 2. If the extended precision has 53 + k bits in the mantissa, then, obviously, no error occurs for up to 2 k summands and

n i=1 p i = n i=1 e i .
The summands p i are added with the respect to decreasing order into the array element corresponding to their exponent. Note that such an algorithm requires twice as much running time compared to our algorithm.

Rump et al. [START_REF] Ogita | Accurate sum and dot product[END_REF][START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Fast high precision summation[END_REF] proposed several algorithms for summation and dot product of floating point numbers. These algorithms are based on iterative application of compensations. An extension of the compensation of two floating-point numbers to vectors of arbitrary length is also given and used to compute a result as if computed with twice the working precision. Various applications of compensated summation method have been proposed [START_REF] Graillat | Algorithms for Accurate, Validated and Fast Polynomial Evaluation[END_REF][START_REF] Graillat | Compensated Horner scheme in complex floating point arithmetic[END_REF]. Thévenoux et al. [START_REF] Thévenoux | Automatic sourceto-source error compensation of floating-point programs: code synthesis to optimize accuracy and time[END_REF] implement an automatic code transformation to derive a compensated programs.

Also, we mention the accurate floating-point summation algorithms introduced by Demmel and Hida [START_REF] Demmel | Accurate floating point summation[END_REF][START_REF] Demmel | Accurate and efficient floating point summation[END_REF]. Given two precision f and F with F > f , Demmel and Hida's algorithms use a fixed array accumulators A 0 , ....A N of precision F for summing n floating-point numbers of precision f such that S = n i=1 s i .

These algorithms require accessing the exponent field of each s i to decide to which accumulator A j to add it. More precisely, each A j accumulating the sum of the s i where e leading bits are j. Then, these A j are sorted in decreasing order to be summed. Consequently, complexity of these algorithms is equal to O(n log n), because of the sorting step.

Parallel approaches: In addition to the existing sequential algorithms, many other parallel algorithms have been proposed. Leuprecht and Oberaigner [START_REF] Leuprecht | Parallel algorithms for the rounding exact summation of floating point numbers[END_REF] describe parallel algorithms for summing floating-point numbers. The authors propose a pipeline version of sequential algorithms [START_REF] Bohlender | Floating-point computation of functions with maximum accuracy[END_REF][START_REF] Pichat | Correction d'une somme en arithmetique a virgule flottante[END_REF] dedicated to the computation of exact rounding summation of floating-point numbers. In order to ensure the reproducibility, Demmel and Nguyen [START_REF] Demmel | Parallel reproducible summation[END_REF] introduce a parallel floatingpoint summation method based on a technique called pre-rounding to obtain deterministic rounding errors. The main idea is to pre-round the floating-point input values to a common base, according to some boundary, so that their sum can later be computed without any rounding error. The error depends only on both input values and the boundary, contrary to the intermediate results which depend on the order of computations.

Accurate Summation Algorithm

In this section, we describe our summation algorithm introduced in [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF] for accurately summing n floating-point numbers. Our algorithm enjoys the following set of properties. First, it improves the numerical accuracy of computations without increasing the cost of complexity compared to the naive algorithm. Second, it performs all the computations in the original working precision without using accumulators of higher precision. Last, using this new algorithm, we increase the numerical accuracy and, as a side effect and shown in the next sections, we also improve the execution time and reproducibility of summation.

For the algorithm detailed hereafter, we assume that we have P processors and n summands (with n P ). We assign n/P summands to each processor. For computing the sum S = n i=1 s i , ∀0 ≤ i < P processor i computes (i+1)×n/P -1 j=i×n/P s j . Then a reduction -a last sum -is done to compute the final result. First of all, Algorithm 1 allocates an array called sum by exp which is created and initialized at 0 for all its elements before starting the summations. The array sum by exp has exp max -exp min + 1 elements. Let us assume that the exponents of summands range from exp min to exp max. The main idea is to sum all the summands whose exponent is 2 i in the cell sum by exp[i + bias] such as bias = -exp min, this avoids most absorptions while avoiding to sort explicitly the array. Let exp (s i ) denote the exponent of s i in base 2. For computing the sum S = n i=1 s i , each value s i is added to the appropriate cell sum by exp[exp s i + bias] according to its exponent. For parallelism, each processor has an sum by exp array to sum locally. In order to obtain its local final result, we add these values in increasing order. Once the final local sums are computed, a reduction is done and the processor receiving the result of the reduction gets the total sum. To emphasize the property regarding the cost of Algorithm 1 Accurate summation with local sum at each processor complexity mentioned above, we note that the cost of Algorithm 1 is O(n), no sorting being performed and the access to the data to be summed being done only once.

Let us mention that another refined implementation of Algorithm 1 has been proposed in [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF]. This second implementation shares a common idea with Algorithm 1 which is related to the way that the summands are added and differs from Algorithm 1 in that the final sum is not carried out in the same way. The advantage of this other implementation is that the summation results are more accurate than those of Algorithm 1. Concerning its drawback, the cost of complexity will be higher in the constant of the O(n) while being linear. In the following sections, we evaluate the simpler and most cost-effective version of Algorithm 1 regarding accuracy, convergence speed for iterative schema and reproducibility. We note that our implementations are done in the C programming language with MPI, compiled with MPICC 3.2, and made them run on an Intel i5 with 7.7 GB memory. Let us also note that for our experiments, we report numerical values relying on thousands separators, typesetting 1234567.89 as 1, 234, 567.89.

Numerical Accuracy

In this section, we first evaluate the numerical accuracy of our algorithm introduced in Section 2.3. Secondly, we address the issue of the compromise between numerical accuracy and running time of the studied algorithms. We take into consideration two examples, namely Simpson's rule and the LU factorization method. For each example, we have implemented two parallel versions, using MPI [START_REF] Pacheco | An Introduction to Parallel Programming[END_REF]. The first one, called original program uses simple sums: to sum n values x 1 , ...., x n it computes (((x 1 + x 2 ) + x 3 ) + .... + x n ). The second program uses Algorithm 1 and is called, the accurate program. The experiments are carried out for several numbers of processors.

Simpson's Rule

Our first example computes the integral C × b 0 f (x)dx of mathematical functions f using Simpson's rule. The Simpson's rule is a numerical method that approximates the value of a definite integral of a function f using quadratic polynomials. We measure the efficiency of our algorithm on this example by computing the absolute errors between both the results of the original and accurate programs with respect to the analytical solution of the integral, as shown in Figure 1. We integrated the following functions C × cos(x), C × (1/x 2 + 1) and C × tanh(x) with C = 10 6 , and b ranging in [2; 5]. The number of processors P ranges in [2; 8]. Each processor computes a part of the integral. The second experiment measures the execution time (in seconds) of the original program, accurate program and another program based on sorting. The choice of this last program is motivated by the main idea of Demmel and Hida algoritms [START_REF] Demmel | Accurate floating point summation[END_REF]. Let us consider the third function C × tanh(x) for P = 8, Figure 2 displays the running time in seconds taken by each program (original program, accurate program and summation by sorting) to compute the integral of this function. The results show that the summation program based on sorting like Demmel and Hida algorithms [START_REF] Demmel | Accurate floating point summation[END_REF] need more time to compute the integral of the function C × tanh(x). Contrarily to the summation by sorting program, our algorithm called accurate program requires much less time and a little more than the original program. For example, to compute the integral of the function C ×tanh(x) for b = 2 it takes 24s using the summation by sorting program, while this computation takes only 0.37s and 1.08s using the original and the accurate programs, respectively. It is well known that the summations based on sorting performed on the summands before computations are more accurate but these computations take more time (49s for b = 5) to performed than our accurate algorithm (1.15s for the same value of b) where no sorting is performed.

LU Factorization

Our second example is the parallel LU factorization method. This method consists in rewriting a matrix A as the product of a lower triangular matrix L and an upper triangular matrix U such that A = LU . The LU factorization method is a very common algorithm which can be used e.g. to solve linear systems or to compute the determinant of a matrix. In the parallel case, the matrix A is divided into blocks of rows and each processor performs its computations on a given block. For our experiments we generated square matrices of various dimensions n ∈ [200, 800] with increment of 100. These matrices contain values chosen to introduce ill-conditioned sums [START_REF] Langlois | Accuracy versus time: A case study with summation algorithms[END_REF]. In our case, we consider 30% of large values among small and medium. By small, medium and large values we mean respectively, of the order of 10 -7 , 10 0 and 10 7 . This is motivated by the IEEE754 single precision arithmetic. Also, we take vectors x with the same proportions of large values among small and medium as for the matrices. The first experimentation consists of comparing the numerical accuracy of the LU factorization carried out using the original and accurate programs. Let us consider a matrix A and vector x. We start by computing the solution vector given by Ax = b. This vector is considered as the exact solution. Next, we apply the original LU factorization program to the matrix A with P = 16 processors in order to obtain L orig and U orig . In the same way, we factorize the matrix A using the accurate LU factorization into L acc and U acc . We compare the new vector solutions b orig = L orig × U orig × x and b acc = L acc × U acc × x with the exact solution b. Figure 3 represents the absolute errors between the computed solutions after factorization and the exact solution. These experiments show significant improvements: while the difference between the absolute errors of the original and the accurate program are already up to an order of 10 5 for our smallest matrices (n = 200), it reaches up to an order of 10 7 for large ones (n = 600). More precisely, for n = 200 the absolute errors are 452, 1984 and 445, 6448 for the original and the accurate LU factorization, respectively. We obtain 996, 1472 and 655, 3600 for n = 600. Thus, we conclude from this experimentation that the accurate program shows its efficiency in terms of numerical accuracy when we handle large matrices, i.e. when various types of absorptions and cancellation have been introduced.

By the second experimentation, we want to show the running time taken by each LU factorization program for a set of matrices of size varying from 200 to 800 with P = 16. Figure 4 We notice in Figure 4 that the summation algorithm based on sorting requires much more time to compute the LU factorization. Besides, accurate program needs much less time than summation algorithm by sorting and a little more than the original algorithm. To better illustrate, let us take a matrix of size n = 300. We remark that the execution takes only 66s and 90s with original program and accurate program, respectively, while the summation by sorting program requires 224s for the same computation. Let us also remark that the running time obtained for the summation by sorting program for the large matrices are much larger than those obtained for the smallest one. In fact, the execution time increase from 101s to approximately 2h for matrices sizes 200 and 800, respectively, using the summation algorithm by sorting.

Convergence of Iterative Methods

In this section, we focus on the impact of accuracy on the number of iterations required by numerical iterative methods to converge. For our experiments, we consider two iterative methods: Jacobi's method and Iterated Power method. As in the previous section, we implemented two versions of the same algorithm, the original one and our accurate version. We observed the impact on the convergence, comparing their respective number of iterations.

Jacobi's Method

The Jacobi's method is a well known numerical method used to solve linear systems of the form Ax = b. In this method, an initial guess, an approximate solution x 0 , is selected and is iteratively updated until finding the solution x k of the linear system. More precisely, this method iterates until |x

(k+1) i -x k i | < ε.
In our case, the parallelization of the Jacobi's method is done according to the row-wise distribution. Jacobi's method is stable whenever the matrix A is strictly diagonally dominant, i.e. its satisfies the property of Equation [START_REF] Demmel | Accurate floating point summation[END_REF].

∀i ∈ 1, . . . , n,

|a ii | > j =i |a ij |. (5)
We examine the impact of accuracy on the convergence speed for the systems of sizes 10 and 100. While the chosen systems were stable with respect to the sufficient condition of the stability given by Equation ( 5), they are close to unstability with ∀i ∈ 1, . . . , n, |a ii | ≈ j =i |a ij |. Figures 5 represents the difference between the number of iterations of the original and the accurate programs.

Let us take the first system of size n = 10. We notice that for various values of ε varying from 10 -2 to 10 -5 , the convergence speed in terms of number of iterations increases from 59 to 2, 029. For the second system when n = 100, we remark that the number of iterations reduced are larger than those computed for the system n = 10 for the same values of ε. For instance the number of iterations reduced for ε = 10 -2 and ε = 10 -5 are 861 and 10, 946, respectively. From these two examples, we conclude that the smallest values of ε are those which have a large difference between the original and the accurate programs in terms of the number of iterations. Also, for a given value of ε, the accurate program shows its efficiency on the convergence speed on large matrices compared to small ones. Fig. 5: Difference between number of iterations necessary for the original and the accurate programs to achieve the convergence of the Jacobi method.

Iterated Power Method

The Iterative Power method is particularly useful for estimating numerically the largest eigenvalue and its corresponding eigenvector. The idea is to fix an arbitrary initial vector x (0) which contains a single non-zero elements. Then, we build an intermediary vector y (1) such that Ax (0) = y (1) . In order to obtain the vector x (1) , we renormalize y (1) so that the selected component is again equal to 1. For the next iteration, we use x (1) as a selected vector. The iterative process is repeated until convergence. We assume that, the parallelization of the Iterated Power method is done according to the row-wise distribution. Let us take a square matrix A of the form:

A =      d a 12 • • • a 1j a 21 d • • • a 2j . . . . . . a i1 a i2 • • • d     
We assume that a ij = 0.01 and d ∈ [300.0, 500.0] following the methodology introduced in [START_REF] Damouche | Impact of accuracy optimization on the convergence of numerical iterative methods[END_REF]. Figure 6 summarizes the difference between the number of iterations of the original and the accurate Iterated Power method. As it is observed in Figure 6, the accurate program accelerates the convergence speed of the Iterative Power method by reducing the number of iterations needed to converge. Indeed, for the matrix size n = 100 with various values of the diagonal and using P = 4 we show that the number of iterations reduced increases from 205 to 340.

Reproducibility

In this section, we aim at evaluating the efficiency of Algorithm 1 on the improvement of reproducibility. Figures 7 and8 give results of reproducibility for Simpson's rule and Matrix Multiplication, respectively. During our experiments, we consider the original and the accurate programs of each method. We compare the results of each of them on several processors and their respective results with only one processor.

Simpson's Rule

Let's take again the example of Simpson's rule already introduced in Section 3.

In practice, improving the numerical accuracy often improves reproducibility. We measure the efficiency of our algorithm on this example by computing the absolute errors between both the results of the original and the accurate programs by varying the number of processors from 2 to 8 and their respective program with only one processor, as shown in Figure 7. Let us consider several mathematical functions C × cos(x), C × (1/x 2 + 1) and C × tanh(x) with C = 10

(x)=tanh(x)*C f(x)=(1/x 2 +1)*C f(x)=cos(x)*C
Fig. 7: The reproducibility of the integral computations using Simpson method of the original program and the accurate one depending on the number of processors.

Matrix Multiplication

The computation of the matrix-matrix multiplication based on floating-point addition and multiplication which are non-associative operations is prone to accuracy problems. Moreover, the out of order execution of arithmetic operations on different or even similar parallel architectures are different, reproducibility of the results is not guaranteed. In this context, we address the problem of reproducibility in the case of matrix multiplication. To parallelize this method, we assume that each matrix is divided into sub-matrices of the size n/P . For our experiments, we consider square matrices of various dimensions n ∈ [200, 800] with increment of 200. These matrices contain a variety of floating-point values chosen with difference in magnitude. More precisely, they are made of 50% of large values (of the order of 10 7 ) among small (of the order of 10 0 ) and medium (of the order of 10 -7 ). Figure 8 represents the percentage of accuracy computed between the original and the accurate programs carried out using P = 8 and their respective result using only one processor. The results show that for different matrix sizes, the percentage of accuracy of the original program ranges from 3% to 13%. Contrarily to the original program, the percentage linked to the accurate program described in this article is equal to 100% for each matrix which confirms the usefulness of Algorithm 1.

Conclusion and Future Work

In this article, we have focused on the impact of accuracy on both the reproducibility and convergence speed of numerical algorithms. The originality of this article is to study the impact that a new accurate summation algorithm has on the accuracy of numerical methods involving sums. In this work, we have experimented our algorithm on several representative numerical methods by comparing the original and the accurate programs of each of them. The experiments show the usefulness of our algorithm on the improvement of reproducibility. Also, the results obtained show the efficiency of our algorithm to reduce the number of iterations required by numerical iterative methods to converge. More precisely, the accurate program converge more quickly than the original one without loss of accuracy. In a future work, we would like to refine our algorithm by adding a test phase after line 4 of Algorithm 1, i.e. before adding the value s i to the appropriate cell of the array sum by exp. Indeed, if we have a large set of values to be summed with the same exponent, the result produced can have a larger exponent than the initial one. Therefore, a loss of accuracy can be caused during the computations of the local sums. An interesting perspective consists in feeding our algorithms [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF] with a static analysis. Our idea is to rely on static analysis to detect the exponents range of a given set of values that will be summed. Once this range is accurately computed, one can use our summation algorithm and obtain more accurate floating-point results. Further, it will be interesting to explore the impact of our algorithms in the context of neural networks. Our goal is to improve the numerical accuracy of computations by using the accurate summation algorithms [START_REF] Benmouhoub | Parallel Accurate and Reproducible Summation[END_REF] as a replacement for their summation algorithms.

Fig. 1 :

 1 Fig. 1: The absolute errors between the original program and the accurate one for the integral computation of three different functions (C × cos(x), C × (1/x 2 + 1) and C × tanh(x)) with the corresponding analytical result by varying the upper bound of the integral b = 2, 3, 4, 5.

Fig. 2 :

 2 Fig. 2: Execution time of the original program, accurate program and summation algorithm by sorting for the integral computation of the function C × tanh(x).

Fig. 3 :

 3 Fig. 3: The absolute errors of the LU factorization by the original and accurate algorithm for matrix of different sizes.

Fig. 4 :

 4 Fig. 4: Execution time of the original program, accurate program and summation algorithm by sorting of the LU factorization.

  reduced Number of iterations reduced Size of matrix n=100 Size of matrix n=10

Fig. 6 :

 6 Fig. 6: Difference between number of iterations of original and accurate Iterated Power method (n = 100, d ∈ [300, 500] with increment of 20).

Fig. 8 :

 8 Fig. 8: The reproducibility of the matrix-matrix multiplication for different size of matrices.

Table 1 :

 1 Binary formats of the IEEE754 Standard.

	Format	#total bits f bits exp bits expmin expmax
	Half precision	16 bits	11	5	-14	+15
	Single precision	32 bits	24	8	-126	+127
	Double precision	64 bits	53	11	-1122	+1223
	Quadruple precision	128 bits	113	15	-16382 +16383

  summarizes the execution time taken by each algorithm (original program, accurate program and summation by sorting) to compute the LU factorization.

	12000	Original program,proc=16
		Accurate program,proc=16
		Summation by sorting,proc=16
	10000
	8000
	6000 Times in seconds
	4000
	2000
	0 200	300	400	500 Matrix sizes	600	700	800

Original program,proc=2 Accurate program,proc=2 Original program,proc=4 Accurate program,proc=4 Original program,proc=8 Accurate program,proc=8 Original program,proc=2 Accurate program,proc=2 Original program,proc=4 Accurate program,proc=4 Original program,proc=8 Accurate program,proc=8 Original program,proc=2 Accurate program,proc=2 Original program,proc=4 Accurate program,proc=4 Original program,proc=8 Accurate program,proc=8 f(x)=tanh(x)*C f(x)=(1/x 2 +1)*C f(x)=cos(x)*C

Acknowledgments This work was supported by a regional funding (Region Occitanie) and partially by project ANR-17-CE25-0018 FEANICSES.