# Analyzing the impact of floating-point precision adaptation in iterative programs

Guillaume Revy

► To cite this version:

Guillaume Revy. Analyzing the impact of floating-point precision adaptation in iterative programs. ARITH 2021 - IEEE 28th Symposium on Computer Arithmetic, Jun 2021, Lyngby, Denmark. pp.25-32, 10.1109/ARITH51176.2021.00016. lirmm-04939368

HAL Id: lirmm-04939368

https://hal-lirmm.ccsd.cnrs.fr/lirmm-04939368v1

Submitted on 10 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analyzing the impact of floating-point precision adaptation in iterative programs

Guillaume Revy

Univ Perpignan Via Domitia, DALI, Perpignan, France

LIRMM, Univ Montpellier, CNRS (UMR 5506), Montpellier, France

***Abstract***—The amount of floating-point computations in numerical programs is ever increasing. In order to improve their performance, the current trend is to carefully adapt the precision of some floating-point operations, to take advantage of new features of modern architectures. Such adaptation process requires tools to evaluate the impact of these transformations on the accuracy of the output. This article presents an analysis tool to help developers in evaluating the impact of adapting the precision in floating-point programs, with a particular focus on iterative programs. This tool is implemented in the LLVM compilation framework, and it includes two main modules: fp2mp that enables to instrument floating-point programs with multiple-precision computations, and a loop splitter that enables to split the iteration space of a loop into several subspaces. We illustrate the interest of this tool on various examples, and particularly how it enables to evaluate automatically the impact of adapting the precision of floating-point operations in loops by considering different precisions over the iteration space.

**Keywords:** floating-point arithmetic, multiple-precision, precision adaptation, dynamic program analysis, LLVM

## I. INTRODUCTION

Numerical applications are ever more complex, and the amount of floating-point computations is ever increasing. To design such applications, developers must choose between different formats. Indeed the IEEE 754-2008 standard defines four floating-point formats, i.e. binary16, binary32, binary64, and binary128 [1]. Various non IEEE formats also exist, like bfloat16 or Posit, for example. Yet, most of these applications use the binary64 floating-point format, that is the most accurate floating-point format directly available in hardware on most modern architectures, and the others remain underused. This choice, very often made independently of the expected result accuracy, leads to an underuse of the architecture features, and consequently to a possible loss of performance.

In order to improve application performances, the current trend is to adapt the precision of some floating-point operations. This technique enables in particular to reduce the cost of some floating-point operations, to limit the cost of memory loads, or to better take advantage of vector units of modern architectures, and consequently to improve the execution time. For example, in [2], for solving dense linear systems by iterative refinement, a speedup of up to 1.8 is observed between full binary64 precision and mixed precision programs, i.e. programs combining different arithmetic (typically binary32 and binary64), and this speedup may raise to up to 11 when the cost difference between binary32 and binary64 computations is large. In this case, the choice of using only

binary64 or combining it with binary32 is guided by the complexity of each subroutine. However in the general case, carefully and rigorously analyzing the precision requirement for each operation according to an expected output accuracy is difficult, even for numerical software experts. And this difficulty increases with the size of the program, and the combinatorics of format adaptation.

Therefore many tools have been designed to automatically adapt floating-point formats. For doing this, an alternative is to perform dynamic program analysis, where one of the main steps consists in analyzing and evaluating how some program transformations do impact the accuracy of the result, for a dataset for which a reference result is known. Hence, for this purpose, a first approach consists in comparing the result of both the transformed and the original programs, compiled and executed on the target architecture. This is implemented for example in Craft [3], Precimonious [4], or Blame [5] and HiFPTuner [6] that work on the top of Precimonious, and that target CPU architectures, while GPUMixer targets GPU [7]. Another approach is to use stochastic arithmetic to evaluate the numerical impact of the format adaptation like in Promise [8] or in Verrou [9]. Also, some tools use software libraries to emulate floating-point arithmetic to estimate this impact. This is for example the case in AMP [10] that relies on SoftFloat,<fnref n="1" /> and in fpPrecisionTuning [11] that uses the MPFR library,<fnref n="2" /> but where a limited scope of the C language is handled (e.g. mathematical function calls are not transformed).

This article addresses the design of a LLVM tool to analyze the impact of floating-point precision adaptation. It works by instrumenting the LLVM intermediate representation (IR) with multiple-precision computations, through calls to the MPFR library. By doing this, it enables to emulate floating-point arithmetic, eventually non-standard or not available in hardware, and thus to estimate the loss of accuracy due to some program transformations. This technique of emulating precision has already been used, not only for mixed precision purpose, but also to detect floating-point anomalies. Let us cite, for example, FPDebug [13] or Herbgrind [14].<fnref n="3" /> More recently two new tools combine LLVM and MPFR to tackle this problem: PositDebug/FPSanitizer [15] and NSan [16]. However unlike most of the existing tools using MPFR as

<footnote n="1">¹See http://www.jhauser.us/arithmetic/SoftFloat.html.</footnote>

<footnote n="2">²See http://www.mpfr.org/ and [12].</footnote>

<footnote n="3">³See https://fpbench.org/community.html for a more detailed list of tools.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Number of significant bits of $u_i$, for approximating the function $1/x$ using Newton-Raphson method and $(x, u_0) = (2, 0.05)$.</figcaption></figure>

floating-point support, a peculiarity of our tool is that it enables to output a new instrumented IR, that can be further modified, at least partially, in a precision adaptation process.

Furthermore these tools work by traversing the IR, and by applying transformations to each instruction. To our knowledge, instructions of a loop body are handled regardless of the iteration at which they will be executed. However, the precision requirement for such instructions may be different at certain iterations. Some attempts have already been carried out to explore precision trade-offs in iterative programs [17]. The originality of our work is to combine our analysis tool with a loop transformation that splits the iteration space of a loop into several subspaces, providing the possibility of applying transformations on some iteration subspaces, only.

Hence the main contributions of this article are:

* A module, called fp2mp, using LLVM capabilities to instrument floating-point programs with multiple-precision computations, in order to analyze and evaluate the impact of modifying the underlying floating-point precision,

* A loop splitter, developed as a LLVM loop transformation, that enables to split the iteration space of a loop into several reduced subspaces,

* A precision update module used to adapt the precision of instructions in some parts of the IR produced by fp2mp,

* And some examples to show the interest of our approach, and to illustrate how it can be used to analyze the impact of adapting some parts of iterative programs relying on bounded or unbounded loops.

This article is organized as follows: After a motivating example in Section II, some background on LLVM and MPFR is given in Section III. Then our analysis infrastructure is presented in Section IV. It includes a description of the three main modules. Finally some illustrative examples are given in Section V, before a conclusion in Section VI.

## II. MOTIVATING EXAMPLE

Let us start here by the example of approximating the inverse function using the Newton-Raphson method. This is a well-known numerical method to find roots of a real-valued function. Starting from a first approximation $u_0$, it converges

quadratically to the solution. For the function $x \mapsto 1/x$, the sequence of $u_i$ is defined as follows:

$$u_{i+1} = u_i \cdot (2 - x \cdot u_i), \quad u_0 \in (0, 2/x).$$

Now let us consider the input $x = 2$ and a first approximation $u_0 = 0.05$. Figure 1 shows the number of significant bits of $u_i$ according the number of iterations. (The number of bits is computed using binary64 arithmetic as $-\log_2(|(\hat{x} - x)/x|)$, and binary16 curve is determined using MPFR with 11 bits of precision.) We observe that after 9 iterations, the binary64 value converges to the solution (with 53 significant bits).

But the key observation is that $u_7$ in binary32 arithmetic has the same number of significant bits as in binary64 arithmetic. We conclude that on this simple case the first 7 iterations could have been carried out in binary32 arithmetic. We observe also that this precision could have been reduced to binary16 until iteration 6.

Here this is a simple example, that can be treated by hand. But in the general case, there is a need for automated tools to reduce or extend the precision of certain instructions, at all loop iterations or only at some of them, and to analyze the impact of these transformations on the numerical quality of the result. This is particularly true if we want to test non-standard arithmetic or not available in hardware. This is the topic addressed in this article.

## III. BACKGROUND ON LLVM AND MPFR

Our work is based on LLVM and MPFR. This section presents their main features useful for our development.

### A. The LLVM compilation infrastructure

LLVM (Low Level Virtual Machine) provides a compilation infrastructure, that is particularly suitable for program analysis, transformation, and optimization [18]. Its flow follows that of a classic compiler, and it works in three main steps: the front-end (i.e. `clang` for C and C++ languages), the middle-end, and the back-end [19]. Like most modern compilers, `LLVM` middle-end embraces various analysis and transformation steps, including for example those for loop optimization, and implemented as passes. An advantage of using `LLVM` is that it provides a simple mechanism for developing new passes, which is used in this work.<fnref n="4" />

Hence our passes are integrated into the middle-end, and we use the following three main features of `LLVM` IR [20]:

* It relies on an unlimited register machine instruction set in the SSA form (Static Single Assignment). Each register is thus assigned to only once, and a new register is automatically created for each instruction result.

* Most of its instructions, including arithmetic ones, are in three-address form, ensuring that they take at most two operands and produce a single result.

* It is strongly typed, meaning each register and instruction has its own type (integer, floating-point, ...).

<footnote n="4">⁴See http://llvm.org/docs/WritingAnLLVMPass.html for details.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Floating-point precision adaptation workflow.</figcaption></figure>

This work has been implemented in LLVM 10.0.0, that gives access to loop information that helped in developing our loop splitting module. This is discussed in Section IV-B.

## B. The MPFR library

In the IEEE 754-2008 standard, a floating-point number $x$ is defined by a triplet $(s, m, e)$, where $s$ is the sign, $e$ is the exponent, and $m$ is the significand in precision $p$. It defines also four formats, each of them having its own precision $p$ and exponent range, to which $e$ belongs. To emulate other arithmetic, an alternative is to use the MPFR library. It provides the capability of computing in multiple-precision, i.e. with floating-point numbers with higher (or lower) precision and with larger (or smaller) exponent range.

In MPFR, the exponent range is the same for all the MPFR variables of the program, but a precision is attached to each number. A MPFR variable may be manipulated through an object of type `mpfr_ptr`. This is a pointer to a MPFR structure (`_mpfr_struct`), that embraces especially its precision, its sign, and its exponent. This precision is set at initialization-time, as shown in Section IV-A1, while the sign and the exponent are dependent to the value of the number.

Also most MPFR functions use a prototype of the form `mpfr_op(dst, src1, src2, rnd)`, where `dst` is the output, `src1` and `src2` two inputs, and `rnd` the rounding mode. This format is particularly suitable for our work since it fits the three-address form used in the LLVM IR to manipulate floating-point instructions, which eases the process.

# IV. ANALYSIS INFRASTRUCTURE FOR ITERATIVE PROGRAMS

This section presents the infrastructure we have implemented to analyze the impact of floating-point precision adaptation, with a particular focus on iterative programs. Its flow is given in Figure 2. It has been developed into the LLVM compilation framework, and it embraces three main modules: the `fp2mp` module, a loop splitter, and a precision update module.

## A. The fp2mp instrumentation module

The `fp2mp` module is intended to help developers in executing their programs in different floating-point precisions, either in higher precision or lower precision, and eventually using non-standard floating-point arithmetic or non available in hardware. It works by instrumenting a LLVM IR with MPFR function calls. As mentioned above, such technique has already been used to detect numerical anomalies. But in our case, our aim

is to observe the impact on the accuracy of a computation of adapting floating-point precisions in a program.

1) Instrumentation process: This module works at compile-time by program instrumentation. Particularly, it traverses the LLVM IR of a program compiled with the lowest optimization level. And for each floating-point instruction, it inserts automatically the equivalent instruction in multiple-precision through a MPFR function call.

To illustrate how `fp2mp` proceeds, let us consider the toy problem of evaluating $\sqrt{2 \cdot x}$ in binary64 arithmetic. Listing 1 gives a C program, while a snippet of the IR is given in Listing 2. We can see that for each local variable and function argument, an `alloca` instruction is used to allocate memory on the stack frame of the function, that will be automatically deallocated as soon as the function returns to its caller. Also a typed register (`%0`, ...) is created for each instruction, including `load` instructions that read data from memory, and floating-point arithmetic instructions (e.g. `fadd`, `fmul`, ...). Note finally that most mathematical function calls are replaced by LLVM intrinsics, like `sqrt` by `llvm.squareRoot.f64` at Line 9.

Now in order to instrument the IR with multiple-precision computations, the `fp2mp` module works on the entire program and it proceeds in 4 main steps. The result is illustrated in the snippet of IR in Listing 3.

1) It duplicates each function (except the `main`), and it inserts in the duplicated function a multiple-precision argument for each floating-point argument, as well as one for the return value if the original function returns a floating-point value. This is shown in Listing 3: the function `eval.mpfr` is the duplication of `eval` in Listing 2, where `x.arg.mpfr` and `eval.mpfr.return` correspond to the argument `x.arg` and the return of the function, respectively. The instrumentation process now continues on this new function, so as to keep the original function unchanged.

2) For each floating-point instruction (`alloca`, `load`, ...) and register of the new function, it creates a MPFR pointer (`mpfr_ptr`) on the stack frame. For example, at Line 6, `x.ambr.mpfr` is the MPFR pointer corresponding to the register `x.ambr` holding the input argument `x.arg`. This pointer is then allocated, and the MPFR number is initialized, and its precision is set to 512, through a call to `__fp2mp_init`, that internally calls `mpfr_init2` (Line 7). Note that this precision value can be modified using a command-line argument.

3) For each floating-point instruction, it inserts the equivalent MPFR function call. This includes arithmetic operations (e.g. `fadd`, `fmul`, ...), and mathematical functions from standard `libm` or `LLVM` intrinsics that have a mapping in the MPFR library. For example, Line 22 of Listing 3 is the call to `mpfrMul` function that corresponds to multiple-precision version of the `fmul` instruction at Line 19, where the last argument corresponds to the rounding mode, that is by default to-nearest-even (parameter equals to `MPFR_RNDN = 0`).
<figure><img src="image_4.png" /><figcaption>Figure 3: Loop splitting description.</figcaption></figure>

latter takes this argument together with its multiple-precision counterpart, and it returns the relative error on the binary64 value assuming the MPFR value is the exact one.

However, a disadvantage of such side-by-side execution comes from deviation of branches involving floating-point values. Indeed, floating-point arithmetic is highly non-intuitive, and performing a comparison between two floating-point values may give two different results according to the underlying precision used. And if used in a conditional, it may lead to branch deviation, meaning the branch instruction may thus take a different control path in higher or lower precision than in original precision. In our approach, in this case, we choose to follow the original computation flow. As a consequence, this may lead to result inconsistency, since the results in modified precision may not correspond anymore to the results of the program obtained as if it would have been executed in this particular precision. This point is not discussed here, and we offer to the user the possibility of either inserting assert instructions in the instrumented program to be warned of such deviations, or following MPFR flow as in Section V-B.

## B. Loop splitting strategy

As mentioned in Section I, when dealing with the precision adaptation in floating-point program with iterative subroutines, the precision requirement for a variable may vary according to the iterations. To address this topic, a first approach would consist in unrolling the loop, at least partially, and in adapting the precision of each instruction. In a certain sense, this is done in [22], where the Horner rule is fully unrolled and the precision of each operation is chosen so as to achieve a given accuracy on the result. This approach could make it possible to adapt precision of iterative programs. But in the general case, it becomes unreachable as soon as the iteration space gets large or unbounded.

Therefore our alternative consists in splitting a loop into several subloops. Particularly, the approach developed here splits the iteration space of the loop into several subspaces, creating thus subloops that iterate over reduced iteration spaces, as in [23]. The underlying idea is to be able to execute these subloops with different levels of precision.

1) **Implementation in LLVM:** To illustrate this, let us consider the piece of C program below, implementing a for loop.

```c
#pragmaclang loop split_ratio(RATIO)
for (i = min; i < max; i += step) {
  /* ... */
}
```

The idea is to provide the capability to the user to split the original loop according to a given ratio. This ratio is defined here through a new pragma we have introduced in the front-end. It corresponds to the percentage of iterations executed by the first subloop (e.g., a ratio equals to 25 means the first loop iterates over 25% of the full iteration space). After splitting, the code would be the one below.

```c
bnd = min + step * floor((max-min) * RATIO / (step*100));
for (i = min; i <= bnd; i += step) {
  /* ... */
}
for (j; i < max; i += step) {
  /* ... */
}
```

This technique has been implemented in the LLVM 10.0.0 framework, that to our knowledge does not offer such a feature so far. It relies on the capabilities offered to determine the loop induction variable (*i*), the bounds of the iteration space (*min* and *max*), and the increment value (*step*). We will see in Section IV-B2 how to proceed if this information cannot be determined.

Moreover, to be transformed, this loop must be in a canonical (simplified) and LCSSA (Loop-Closed SSA) forms. Both these makes loop transformations and optimizations easier.<fnref n="6" /> Particularly, as illustrated in Figure 3a, in the canonical form, a loop is defined by the following basic blocks:

* a preheader, that dominates the loop but is not part of it,

* a header, that contains the body of the loop,

* a latch/exiting, that has single backedge to the header, and an edge to each exit block,

* and exit blocks, that are outside the loop but dominated by the exiting block.

Using this description, Figure 3b shows the way we have implemented to split the loop. It consists in 5 main steps:

1) It inserts the computation of the bound *bnd* in the predecessor of the preheader of the loop.

2) It clones the blocks contained into the loop (header and exiting/latch) together with the preheader.

3) It updates the condition in the branch instruction of the exiting/latch block of the first loop (*i* <= *bnd*), so as to reduce its iteration space,

4) It removes the link between exiting and exit blocks, and it links carefully the exiting block of the first loop to the preheader of the second one, by taking care of propagating all the relevant data, before linking the exiting block of the second loop to the exit block.

Here the description is given to split the original loop into two subloops, but we can also split loop in more than two subloops. For this, a new pragma has also been defined to indicate the number of loops created.

<footnote n="6">⁶See https://llvm.org/docs/LoopTerminology.html for details.</footnote>
This pass produces a new LLVM IR. And in the same way as the IR produced by the fp2mp module, everything is done to enable the modification of this one later. In particular, we include all the information useful to further modify the splitting ratio. This is done thanks to the LLVM metadata mechanism, that enables data to be associated with an instruction. These metadata can then be detected when traversing the IR. This is illustrated by the IR snippet below: register %28 stores the result of `(max - min) * RATIO`, while %29 is the one of step * 100. Modifying the ratio simply consists in modifying both these instructions.

```
%28 = mul i32 %27, 25, !llvm.loop.split.ratio.ratio !281
%29 = mul i32 1, 100, !llvm.loop.split.ratio.total !282
; ....
!281 = !{!"25"} ; original value of RATIO
!282 = !{!"100"} ; original value
```

This will then avoid us having to restart this splitting process for each new ratio.

2) When loop bounds cannot be determined: As mentioned above, we need to determine the loop bounds, i.e. the start and end values, together with the increment value. These are useful to compute the split bound `bnd`. However in the general case, we have no guarantee that `LLVM` is able to compute these, and even if they are computable. This is for example the case of a `while` loop for which the number of iterations is unbounded.

In this case, instead of defining a ratio, we give the possibility of defining a number of iterations, that corresponds to the number of iterations for the first loop, as illustrated in the snippet of C program below. This is done with a new pragma.

```c
#pragmaclang loop split_iterate(ITER)
while (err < epsilon) {
  /* ... */
}
```

The idea is to split the loop, such that the first loop iterates at most `ITER` times. For doing this, the strategy follows the one presented in Section IV-B1, with the following exceptions.

1) It inserts a new counter in the predecessor of the pre-header, initialized to 0, to count the iterations.

2) It updates the condition in the branch instruction of the exiting/latch block of the first loop to be `counter < ITER && err < epsilon`.

3) And it inserts instructions to increment the counter in the exiting/latch block of the first loop.

After this, the code would be the one below.

```c
counter = 0
while (counter++ < ITER && err < epsilon) {
  /* ... */
}
while (err < epsilon) {
  /* ... */
}
```

### C. Update precision module

As mentioned above, the originality of our approach is that it produces an IR that can be further manipulated. Especially the instrumentation process is done once, and after the precision of all MPFR objects or the rounding mode of all MPFR functions can be modified directly on this output IR.

<table><caption>Table I: Functions implemented using Horner rule.</caption><thead><tr><th>Function</th><th>Degree</th><th>Interval</th></tr></thead><tbody><tr><td>log<sub>2</sub>(1 + x)</td><td>31</td><td>[-2<sup>-2</sup>; 2<sup>-2</sup>]</td></tr><tr><td>exp(x)</td><td>26</td><td>[-2<sup>-1</sup>; 2<sup>-1</sup>]</td></tr><tr><td>sin(x)</td><td>28</td><td>[-π/4; π/4]</td></tr><tr><td>sinh(x)</td><td>30</td><td>[-1; 1]</td></tr><tr><td>erf(x) - 1/2</td><td>32</td><td>[-1/4; 1/4]</td></tr></tbody></table>

For example, at a C program level, reducing the precision of an operation (here, an addition) first consists in reducing the one of its operands, as shown below.

```c
res = (float)a + (float)b;
```

Therefore, this module takes as input the instrumented `LLVM` IR and a list of `MPFR` objects, and it updates first the precision of these objects. This list is built using the `JSON` output file produced by `fp2mp`. Second, for the particular case of instructions, it updates also the precision of its operands, by eventually inserting temporary `MPFR` variables to handle these reduced precision arguments.

## V. EXPERIMENTAL RESULTS

This section presents two examples to show how our tool can help in evaluating the impact of modifying the precision in iterative programs based on bounded and unbounded loops.

### A. Case of bounded loops: polynomial evaluation

Let us consider the case of bounded loops, through the example of the polynomial evaluation using Horner rule. The C program is given below, where the number of iterations is the degree of the polynomial.

```c
double
evaluate(double *a, int n, double x)
{
  double res = a[n];
  #pragmaclang loop split_ratio(RATIO)
  for (int i = n-1; i >= 0; i--)
    res = res * x + a[i];
  return res;
}
```

For this experiment, we used 5 functions, whose descriptions are given in Table I. The polynomials were computed using `fpminimax` from the `Sollya` tool,<fnref n="7" /> with coefficients in the binary64 format. Then we used our tool to split this loop using the strategy presented in Section IV-B1, and to apply binary64 → binary32 transformations on the first loop. The number of iterations of this loop was determined using the `RATIO` parameter, whose value varied from 0% to 100% by step of 5% (0% and 100% are the binary64 and binary32 versions, respectively). Figure 4a shows the maximum relative error, when evaluating these polynomials on 10000 inputs equally spaced on the approximation interval. Here and hereafter, in the instrumented instructions, binary64 and binary32 are emulated with `MPFR` in precision 53 and 24, respectively, the result in binary64 being considered as the "exact" one.

As expected, the percentage of iterations executed in binary32 impacts the accuracy of the result. Especially it decreases when this percentage increases. Furthermore we can

<footnote n="7">⁷See https://www.sollya.org and [24]</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Maximum relative error according to the percentage of iterations in binary32 (left) or binary16 (right).</figcaption></figure>

observe that up to ≈ 60%, this percentage has a relatively small impact on this accuracy. Indeed for these 5 functions, until this bound, the maximum relative error is about $2^{-50}$ and it does not vary too much. Conversely, after this bound, the degradation of the result increases with this percentage.

We performed the same experiment with the binary16 format emulated in MPFR with precision 11, instead of binary32. The results are given in Figure 4b. Obviously, the impact of reducing some operations in precision 11 is more important than in the previous case. However, we still observe the same curve behavior: until a ratio of ≈ 60%, the error does not vary too much, and after the degradation increases with the ratio.

Our approach does not provide any guarantee on this splitting bound, since this experiment is performed on 10000 sample inputs, only. This can still be used as a hint in the implementation process of a function. But the key point here is twofold: (1) all this work was done in a fully automated way, and (2) it enables to evaluate the impact of using arithmetic not implemented in hardware.

## B. Case of unbounded loops: conjugate gradient

Let us now consider the case of unbounded loops, through the example of the Conjugate Gradient method. This is a well-known iterative solver for linear systems. Particularly, given a symmetric definite-positive matrix A and a vector b, it enables to solve $Ax = b$. The method is given in Algorithm 1, where we can see that the loop stops as soon $||r_k||$ is no greater than $\epsilon$

or $k$ has reached maxiter, meaning the number of iterations is not known in advance.

This is an iterative method, where the number of iterations is directly linked to the precision of the computations. For example, the impact of the computing precision on the number of iterations is illustrated in [25]. But in this work, the precision of all the program is modified. Conversely, in our case, we are interested in modifying the precision of certain iterations only, to observe the impact on the speed of convergence.

In this second experiment, we considered a binary64 implementation of the conjugate gradient, and we used our tool to split the loop using the strategy presented in Section IV-B2. Since the number of iterations in not known in advance, the ratio is replaced by the number of iterations for the first loop. Finally we applied binary64 → binary32 transformations on this first loop, and asked fp2mp to follow the MPFR flow.

This method is defined for symmetric definite-positive matrices. Hence as a matrix example, we used the 494_bus matrix from the Suite Sparse Matrix Collection.<fnref n="8" /> And we built a second member $b$, such that the solution $x$ is a vector of ones.

Now to evaluate the impact on the number of iterations of modifying the computing precision of some instructions, we executed this program with $\epsilon = 10^{-6}$ (and maxiter = 100000). Solving this problem in binary64 requires 1315 iterations, while in binary32, it requires 2494 iterations. We then varied the number of iterations in binary32 from 10 to 2490 by step of 10. Figure 5 shows how evolves the total number of iterations.

As expected, starting the computation in binary32 arithmetic leads to an increase of the total number of iterations. For example, when using 10 iterations in binary32, the total number of iterations raises to 1769. But the most interesting point is that, in a fully automated way, we are able to evaluate the impact of reducing the precision at some iterations of the program, even on unbounded loops. Hence, on this example, we observe that the precision of the first iterations significantly impacts the convergence of the method. Indeed, between 10 and 100 iterations in binary32, the total number increases quickly to reach ≈ 2400. Then from 100 iterations, the curve reaches a plateau, that evolves slowly until 2494 iterations.

### Algorithm 1 Conjugate Gradient method.

1: $r_0 := p_0 := b - Ax_0$, and $k = 0$

2: **while** $||r_k|| \ge \epsilon$ and $k < \text{maxiter}$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;$\alpha_k := \frac{r_k^T r_k}{p_k^T A p_k}$

4: &nbsp;&nbsp;&nbsp;&nbsp;$x_{k+1} := x_k + \alpha_k p_k$

5: &nbsp;&nbsp;&nbsp;&nbsp;$r_{k+1} := r_k - \alpha_k A p_k$

6: &nbsp;&nbsp;&nbsp;&nbsp;$\beta_k := \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$

7: &nbsp;&nbsp;&nbsp;&nbsp;$p_{k+1} := r_{k+1} + \beta_k p_k$

8: &nbsp;&nbsp;&nbsp;&nbsp;$k = k + 1$

9: **end while**

<footnote n="8">⁸See https://sparse.tamu.edu/ and [26].</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5: Total number of iterations according to the number of iterations in binary 32, for the conjugate gradient method.</figcaption></figure>

## VI. CONCLUSION AND FUTURE WORK

This article addresses the topic of analyzing the impact of the adaptation of precision in iterative programs. We have designed an analysis infrastructure, that mainly relies on the fp2mp module, that instruments floating-point programs with multiple-precision computations to enable to execute these programs in different precisions, and a loop transformation, that splits the iteration space of a loop into several reduced subspaces. This work has been developed in the LLVM 10.0.0 compilation infrastructure, using the MPFR library to handle multiple-precision. We finally illustrated how this works on two examples, both for bounded and unbounded loops.

The research direction is threefold: First, we have illustrated the use of this tool on relatively small routines. Since all these techniques have been implemented in a compiler infrastructure, it seems to be reachable to handle any programs. Hence it could be interesting to see how it behaves on larger real life scientific applications. Second, in all our experiments, we have considered the numerical quality of the computed result. However, reducing the precision of some instructions may also lead to an improvement of the performance of the application. It could be interesting to systematically evaluate the gain in performance associated with each transformation, to be able to evaluate if these transformations are relevant. Third, we have seen in Section I that dynamic tools dedicated to the design of mixed precision programs heavily rely on such analysis technique. In this sense, this approach can be seen as a first step toward auto-adaptation process. And it could be interesting to see how it can be used for producing mixed-precision programs.

## REFERENCES

[1] "IEEE standard for floating-point arithmetic," IEEE Std. 754-2008, pages 1-58, 2008.

[2] M. Baboulin, A. Buttari, J. Dongarra, J. Kurzak, J. Langou, J. Langou, P. Luszczek, and S. Tomov, "Accelerating scientific computations with mixed precision algorithms," *Computer Physics Communications*, vol. 180, no. 12, pp. 2526-2533, 2009.

[3] M. O. Lam, J. K. Hollingsworth, B. R. de Supinski, and M. P. Legendre, "Automatically adapting programs for mixed-precision floating-point computation," in *ICS'13*, 2013.

[4] C. Rubio-González, C. Nguyen, H. D. Nguyen, J. Demmel, W. Kahan, K. Sen, D. H. Bailey, C. Iancu, and D. Hough, "Precimonious: tuning assistant for floating-point precision," in *SC'13*, 2013, pp. 1-12.

[5] C. Rubio-González, C. Nguyen, B. Mehne, K. Sen, J. Demmel, W. Kahan, C. Iancu, W. Lavrijsen, D. H. Bailey, and D. Hough, "Floating-point precision tuning using Blame analysis," in *ICSE 2016*, 2016, pp. 1074-1085.

[6] H. Guo and C. Rubio-González, "Exploiting Community Structure for Floating-Point Precision Tuning," in *Proceedings of the 27th ACM SIGSOFT International Symposium on Software Testing and Analysis. Association for Computing Machinery*, 2018, p. 333-343.

[7] I. Laguna, P. C. Wood, R. Singh, and S. Bagchi, "GPUMixer: Performance-Driven Floating-Point Tuning for GPU Scientific Applications," in *Proceedings of 34th International Conference, ISC High Performance 2019*, M. Weiland, G. Juckeland, C. Trinitis, and P. Sadayappan, Eds., vol. 11501. Springer, 2019, pp. 227-246.

[8] S. Graillat., F. Jézéquel, R. Picot, F. Févotte, and B. Lathuilière, "Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic," 2016.

[9] F. Févotte and B. Lathuilière, "Debugging and optimization of HPC programs with the Verrou tool," Denver, CO, USA, Nov. 2019.

[10] R. Nathan, H. Naeimi, D. J. Sorin, and X. Sun, "Profile-driven automated mixed precision," CoRR, vol. abs/1606.00251, 2016.

[11] N. Ho, E. Manogaran, W. Wong, and A. Anoosheh, "Efficient floating point precision tuning for approximate computing," in *2017 22nd Asia and South Pacific Design Automation Conference (ASP-DAC)*, 2017, pp. 63-68.

[12] L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, and P. Zimmermann, "MPFR: A Multiple-Precision Binary Floating-Point Library with Correct Rounding," *ACM Transactions on Mathematical Software*, 2007.

[13] F. Benz, A. Hildebrandt, and S. Hack, "A dynamic program analysis to find floating-point accuracy problems," in *PLDI'2012*, 2012.

[14] A. Sanchez-Stern, P. Panchekha, S. Lerner, and Z. Tatlock, "Finding Root Causes of Floating Point Error," *SIGPLAN Not.*, vol. 53, no. 4, p. 256-269, 2018.

[15] S. Chowdhary, J. P. Lim, and S. Nagarakatte, "Debugging and Detecting Numerical Errors in Computation with Posits," in *Proceedings of the 41st ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI'2020*. Association for Computing Machinery, 2020, p. 731-746.

[16] C. Courbet, "NSan: A Floating-Point Numerical Sanitizer," in *Proceedings of the 30th ACM SIGPLAN 2021 International Conference on Compiler Construction*, 2021.

[17] Y. Chatelain, E. Petit, P. De Oliveira Castro, G. Lartigue, and D. Defour, "Automatic Exploration of Reduced Floating-Point Representations in Iterative Methods," in *25th International Conference Euro-Par 2019 Parallel Processing*, 2019, pp. 481-494.

[18] C. Lattner and V. Adve, "LLVM: A Compilation Framework for Lifelong Program Analysis & Transformation," in *CGO '04*, 2004, pp. 75-86.

[19] A. V. Aho, R. Sethi, and J. D. Ullman, *Compilers: principles, techniques and tools*. Addison Wesley, 1986.

[20] Lattner, C. and Adve, V., "The LLVM Instruction Set and Compilation Strategy," CS Dept., Univ. of Illinois at Urbana-Champaign, Tech. Rep., 2002.

[21] N. Revol and F. Rouillier, "Motivations for an arbitrary precision interval arithmetic and the MPFI library," *Reliable Computing*, pp. 275-290, 2005.

[22] C. Lauter, "Arrondi correct de fonctions mathématiques - fonctions variées et bivariées, certification et automatisation," Ph.D. dissertation, Univ. de Lyon - ÉNS Lyon, 2008.

[23] L. Thévenoux, "Synthèse de code avec compromis entre performance et précision en arithmétique flottante IEEE 754," Ph.D. dissertation, Université de Perpignan Via Domitia, Jul. 2014.

[24] S. Chevillard, M. Joldeş, and C. Lauter, "Sollya: An environment for the development of numerical codes," in *Mathematical Software - ICMS 2010*, vol. 6327, 2010, pp. 28-31.

[25] A. Cohen, C. Fabre, F. Pétrot, T. T. Jost, and Y. Durand, "Seamless Compiler Integration of Variable Precision Floating Point Arithmetic," in *CGO 2021*, 2021.

[26] T. A. Davis and Y. Hu, "The University of Florida Sparse Matrix Collection," *ACM Transactions on Mathematical Software*, vol. 38, 2011.