# Performance on SIMD architectures of auto-tuned programs for matrix multiplication

Youssef Fakhreddine, Guillaume Revy

► To cite this version:

Youssef Fakhreddine, Guillaume Revy. Performance on SIMD architectures of auto-tuned programs for matrix multiplication. MCSoC 2024 - IEEE 17th International Symposium on Embedded Multicore/Many-core Systems-on-Chip, Dec 2024, Kuala Lumpur, Malaysia. pp.564-571, 10.1109/MCSoC64144.2024.00096. lirmm-04939346

HAL Id: lirmm-04939346

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

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
# Performance on SIMD architectures of auto-tuned programs for matrix multiplication

Youssef Fakhreddine and Guillaume Revy

Univ Perpignan Via Domitia, DALI, Perpignan, France

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

***Abstract***—The amount of numerical computations in scientific programs is ever increasing. Hence in recent years, a growing interest has emerged in dynamically adapting the precision of floating-point computations to balance performance and accuracy. We focus here on iterative routines. For that purpose, a tool was recently introduced, which enables such precision adaptation at the iteration level, leveraging multiple-precision computations and delta-debugging techniques, in order to produce a set of possible adaptations. The study presented in this article extends the exploration initiated by the aforementioned tool. For doing this, we developed a new tool to apply these adaptations to the input C program, allowing next to examine the performance characteristics of the output program. Leveraging SIMD micro-architectures, we investigate the potential for enhancing the performances delivered by this approach on the execution of precision-adapted iterative routines. By benchmarking against non-optimized versions, we assess the speedup achieved across a spectrum of the matrix multiplication program, and we illustrate how our framework of precision adaptation allows for achieving significant speedups on matrix multiplication, varying according to the accuracy threshold.

**Keywords:** floating-point arithmetic, mixed precision, precision adaptation, SIMD, performance, matrix multiplication.

## I. INTRODUCTION

Nowadays, the amount of computations in numerical programs is ever increasing. In this context many formats exist to represent real numbers on computer systems. Let us cite, for example, the standard IEEE 754 floating-point formats [1], including binary16, binary32, binary64, and binary128, but also nonstandard formats like bfloat16 developed by Google, TensorFloat-32 by NVidia, or Posit proposed by Gustafson [2]. Yet, most of the time and for accuracy reasons, in numerical computing, developers choose to use the binary64 format, which is actually the more accurate format available on most modern architectures. But this possible oversizing of computation means is also often made to the detriment of a possible degradation of the performance of the targeted programs.

In order to improve the performance of floating-point programs while not sacrificing the quality of their results, a widely used alternative consists in writing *mixed precision* programs. This consists in combining data and instructions of different types in the same programs, to reduce the load costs of data from memory, and to benefit from faster instructions in lower precision. Its interest has been shown in [3] where mixed precision programs achieve similar accuracy to the original binary64 precision program while being significantly faster and reducing memory pressure. But in this work, the appropriate

formats of the different subroutines of the programs have been decided statically based on their complexity. And this approach is not usable for any floating-point program.

Therefore several dynamic approaches and tools have been designed to help developers in implementing mixed precision programs, and to decide the data or instructions whose format can be lowered. Let us cite, for example, Precimonious [4], Blame [5], HiFPTuner [6], Promise [7], or Craft [8], [9] integrated into the FloatSmith<fnref n="1" /> project with ADAPT [11]. Most of them proceed (1) by applying a trial-and-error strategy, which consists in evaluating how some format transformations do impact the accuracy of the results, for a dataset for which a reference result is known, and (2) by reducing the combinatorics of possible transformations by using the delta-debugging algorithm [12]. However one of the limitations of existing tools is that they do not deal with the iterative structure of programs. This means that if data or instructions appear in a loop, their format is chosen and the same for all the iterations.

To overcome this issue and to handle iterative programs, a few attempts only have been proposed in [13], [14], [15]. Moreover in some previous work, we have developed an auto-tuning tool whose purpose is to choose the computation precision at the iteration level [16]. It relies on the fp2mp analysis tool [17] to propose a set of transformations that can be applied while keeping a constraint (e.g. on the accuracy of the results) satisfied. But it does not generate a new optimized program, but it just provides hints to the developers on which transformations to apply.

Therefore this article extends the exploration initiated by our auto-tuning tool for iterative routines, and it deals with the optimization of matrix multiplication through the reduction of the computation precision of instructions. The main contributions are the following:

* A new Python script that uses the hints providing by our auto-tuning tool to apply transformations on a input C program, and to generate a new optimized C program,

* And a study to measure the impact on the execution performance of this new auto-tuned program on SIMD micro-architectures, showing that a speedup of up to ≈ 1.7 can be achieved.

This article is organized as follows: After a presentation of our auto-tuning tool in Section II and the statement of the problem in Section III, Section IV deals with our new strategy

<footnote n="1">¹See https://github.com/craftthpc/floatsmith/ and [10].</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Workflow of our auto-tuning tool, presented in [16].</figcaption></figure>

to generate optimized C programs. Then Section V shows some experimental results, before a conclusion in Section VI.

## II. BACKGROUND: AN AUTO-TUNING TOOL FOR ITERATIVE ROUTINES

Our tool, called `dd-fp2mp`, is a dynamic auto-tuning tool designed to address the iterative nature of scientific computing routines [16]. This tool is implemented as a transformation pass in the LLVM pass pipeline, and its workflow is depicted in Figure 1. Unlike existing approaches that aim to minimize the number of possible transformations to consequently reduce the related combinatorics and to speed up the tuning process, the originality of `dd-fp2mp` is that it increases this number judiciously, prioritizing the quality of resulting program over speed in the optimization process.

Particularly, the tool takes as input a C program with hints indicating the loop to be tuned and a constraint to be satisfied. In this article we consider the constraint to be an output accuracy constraint, but it could also be a constraint on the behavior of the program, like the total number of iterations. Then leveraging compilation techniques, including loop splitting and unrolling, it aims to increase the combinatorics of possible transformations in the program to be tuned without altering its semantics, thereby increasing the range of possible optimizations. For doing this, it uses 3 basic blocks.

1) **Loop transformations:** Our tool relies on LLVM's loop unrolling capability and the loop splitting technique developed by Revy [17]. Even if in this article we consider loop splitting, only, both of these make instructions appear in the program without altering its semantics, consequently increasing the number of possible transformations. The way used to unroll or split the loop is passed to the LLVM core by embedding appropriate metadata into the LLVM intermediate representation (IR) through new and appropriate pragmas we have introduced in the front-end. This mainly corresponds to a splitting or unrolling factor. For example, splitting a loop by a factor of 2 gives two subloops that iterate over reduced subspaces, whose union still covers the whole iteration space of the initial loop.

2) **Accuracy analysis with fp2mp:** Then it iteratively applies a subset of possible transformations on the original program. Each configuration's impact on output

constraint is assessed using multiple-precision analysis through the `fp2mp` tool. This tool is designed to experiment with floating-point precision in computer programs. It works by duplicating in the LLVM IR each floating-point instruction used in a program with its multiple precision equivalent through the MPFR library. This allows the user to define the desired precision, and to evaluate the impact of these precision changes on the behavior of the program. But the advantage of `fp2mp` in this approach is that this module offers the possibility of adapting the precision of certain computations only, and in our particular case, of certain subloops, only.

3) **Delta-debugging algorithm:** To identify the most relevant transformations, which maintains acceptable output, our tool relies on the delta-debugging algorithm. This heuristic algorithm is used by the majority of dynamic auto-tuning tools. Initially designed for debugging as a systematic approach to isolate the minimal set of changes in a program or a program's input that cause a specific failure or bug, it has next been adapted to isolate the largest set of changes in a program or program's input, which ensures that the program still successfully works. Hence in the context of precision tuning, the delta-debugging is used to identify the largest set of transformations that still yield an acceptable output, meaning that the constraint remains satisfied.

Therefore through a trial-and-error process, our tool iteratively refines configurations to maintain acceptable output accuracy. However, rather than generating a new optimized program, the tool provides the user with a list of recommended transformations to apply. The approach implemented in this article to apply these transformations is detailed in Section IV.

## III. AUTO-TUNING OF MATRIX MULTIPLICATION PROGRAM BY LOOP SPLITTING

Let $A, B$, and $C$ be three size-$n$ matrices with coefficients in the binary64 format, and where $C = A \times B$. In this article, we consider the naive but vectorizable implementation of matrix multiplication, shown in Listing 1. Using this, the computation of the result matrix $C$ requires the evaluation of $n^2$ dot-products of size $n$. These $n^2$ dot-products are all evaluated using the same code requiring $n$ multiplications and $n$ additions, and consequently with the same computation precision. And assuming matrix coefficients in binary64, all these operations are carried out in binary64 arithmetic, as well.

In order to speed up the execution of this program, while not sacrificing the output accuracy, our alternative is to reduce the computation precision of certain well-chosen additions and multiplications. Indeed, in Listing 1, the inner loop can be vectorized, as shown in Listing 2 with a vectorization width of 4, and where the iteration step is 4, as well. In this case, 4 additions and 4 multiplications are performed at each iteration. Computing all these additions and multiplications in reduced precision, i.e. in binary32 arithmetic, with the same register size theoretically increases the vectorization width by a factor of 2. This means that 8 additions and 8 multiplications are
<caption>Listing 1: Naive but vectorizable matrix multiplication program.</caption>

```c
void
matrix_multiplier(double **A, double **B,
                 double **C, int n) {
    for (int k = 0; k <= n-1; k += 1) {
        for (int i = 0; i <= n-1; i += 1) {
            for (int j = 0; j <= n-1; j += 1)
                C[i][j] += A[i][k] * B[k][j];
        }
    }
}
```

<figure><img src="image_3.png" /><figcaption>Figure 2: Splitting strategies for $n = 6$.</figcaption></figure>

performed at each iteration, leading to a possible speedup of 2. However, reducing the computation precision of all these $n^2$ operations may also lead to a degradation of the accuracy of the output result $C$.

Therefore, in this work, we aim at lowering the computation precision of as many instructions as possible, by keeping coefficients of $A$, $B$, and $C$ in binary64 format, to hope for maximum gain, while guaranteeing that the constraint on the output remains satisfied. For doing this, we use the capability of `dd-fp2mp` to identify this largest set of possible transformations, by splitting or unrolling a loop in the program. The splitting strategy may lead to an improvement of the performance of the program: indeed, it creates several

<caption>Listing 2: Pseudo-code of vectorized version of Listing 1.</caption>

```c
void
matrix_multiplier(double **A, double **B,
                 double **C, int n) {
    for (int k = 0; k <= n-1; k += 1) {
        for (int i = 0; i <= n-1; i += 1) {
            for (int j = 0; j <= n-1; j += 1)
                C[i][j...j+3] += A[i][k] * B[k][j...j+3];
        }
    }
}
```

subloops operating on reduced iteration spaces, and reducing the computation precision in a subloop may allow for increasing the vectorization width and for better benefiting from SIMD micro-architectures. Conversely, the unrolling strategy does not create new subloops, thus all the instructions appear in the same loop body. And since the vectorization width is the same for all instructions in the loop body, either no instruction is transformed and no gain is observed, or the reduction involves all the instructions, but to the detriment of a degradation of the accuracy of the result. Consequently, in the sequel of the article, we consider only the splitting strategy. And to better understand the impact on the result of this loop splitting, let us assume a splitting factor of 2, and the reduction of the computation precision of the second subloop. Now there are three loops that can be split:

* Splitting the loop over $i$ means considering the first $n/2$ rows of the matrix $A$ in binary64 and the $n/2$ last ones in binary32. Consequently, it results in computing the first $n/2$ rows of $C$ using a binary64 dot-product routine, while the $n/2$ last ones are computed in binary32.

* In a similar way, splitting the loop over $j$ results in computing the first $n/2$ columns of $C$ using a binary64 dot-product routine, while the $n/2$ last ones are computed using a binary32 routine,

* And splitting loop over $k$ results in computing all the elements of $C$ using a mixed precision routine, where $n/2$ iterations are in binary64 arithmetic, and the last $n/2$ ones in binary32.

This is illustrated in Figure 2. The case presented in Figure 2c fits best what we expect of mixed precision computation since all results elements are computed using both binary64 and binary32 arithmetic. In addition, we will see in Section IV-C that, from a performance point of view, this is the one for which we can expect the best gain. For these reasons, only the third splitting will be considered in this article.

Furthermore, to decide if a set of transformations is relevant or not, we need to define an accuracy criterion. For doing this, let $\tilde{C}$ be the result matrix computed by the optimized program. Our accuracy criterion is thus defined as follows:

$$ \left\| \frac{C - \tilde{C}}{C} \right\|_{\infty} < \epsilon. \qquad (1) $$

where $\epsilon$ is the accuracy threshold to be satisfied. The objective is now to study the impact of these transformations on the performance of the optimized program.
<figure><img src="image_4.png" /><figcaption>Figure 3: New splitting tool working at the C level.</figcaption></figure>

## IV. NEW LOOP SPLITTING TOOL

This section presents the workflow we have implemented to produce auto-tuned programs for matrix multiplication. After a description of this flow, it presents the way used to produce auto-tuned C programs.

### A. Description of the flow

The flow is given in Figure 3. It takes as input a C program, as the one in Listing 3, and where the loop to be optimized is tagged using a specific pragma, and a constraint through a threshold to be satisfied. In addition, it takes a splitting factor indicating the number of subloops to be created. Using this, `dd-fp2mp` returns a list of transformations to be applied. This list contains the subloops whose computation precision is to be transformed. In this sense, our new splitting approach uses this list to transform the format of corresponding subloops.

In our case, we do not transform the format of data, but the computation precision of instructions. In C language, the computation precision of an instruction directly depends on the format of its operands. Hence applying these transformations requires first converting operands to lower format. For example, let us consider the instructions of the loop body of Listing 1. Their transformation is illustrated on the small piece of code below, where their format is lowered to binary32.

```c
float tmp1 = (float)C[i][j];
float tmp2 = (float)A[i][k];
float tmp3 = (float)B[k][j];
tmp1 += tmp2 * tmp3;
C[i][j] = (double)tmp1;
```

Here the three operands `C[i][j]`, `A[i][k]` and `B[k][j]` are first converted into binary32 data. Consequently, the addition and multiplication are both performed in binary32 arithmetic, as well. Finally, the result `tmp1` is converted back to binary64 arithmetic.

However the cast operations appearing in these transformations may have a negative impact on the performance of the program, for two main reasons.

1) These cast operations are costly, and they may be executed several times and redundantly. For example,

<caption>Listing 3: Annotated input C program of our new splitting tool.</caption>

```c
void
matrix_multiplier(double *A, double *B,
                  double *C, int n) {
  #pragmaOMP __host__ __host
  // SPLIT_FOR [indvar=k, start=0, end=n-1, step=1]
  // REPLACEMENT:
  //   A > A_b32 > heap (n*n)
  //   B > B_b32 > heap (n*n)
  //   C > C_b32 > heap (n*n)
  // INITIALISATION:
  //   for(i = 0; i <= n-1; i += 1) {
  //     for(j = 0; j <= n-1; j += 1) {
  //       A_b32[i*n+j] = A[i*n+j];
  //       B_b32[i*n+j] = B[i*n+j];
  //     }
  //   }
  // PREFIX:
  //   for(i = 0; i <= n-1; i += 1) {
  //     for(j = 0; j <= n-1; j += 1) {
  //       C_b32[i*n+j] = C[i*n+j];
  //     }
  // SUFFIX:
  //   for(i = 0; i <= n-1; i += 1) {
  //     for(j = 0; j <= n-1; j += 1) {
  //       C[i*n+j] = C_b32[i*n+j];
  //     }
  //   }
  //   for (k = 0; k <= n-1; k += 1) {
  //     for (i = 0; i <= n-1; i += 1) {
  //       for (j = 0; j <= n-1; j += 1) {
  //         C[i*n+j] += A[i*n+k]*B[k*n+j];
  //       }
  //     }
  // END SPLIT_FOR
}
```

in Listing 1, the coefficient $A_{0,0}$ is used $n$ times, and if inserted in the loop body, the cast operation of $A_{0,0}$ will be performed $n$ times, as well.

2) These cast operations also prevent the full benefit of vector units, which is penalizing in our case. Indeed in order to vectorize a loop, the compiler first determines its vectorization factor. It equals to the smallest vectorization factor of all its instructions, which indicates the size of vectors they can manipulate. Hence, mixing data and instructions of different types in the same loop body may leads to a reduction of its vectorization factor. For example, mixing binary32 and binary64 data and instructions on AVX2 (i.e. with 256 bits vector registers) yields a vectorization factor of 4, while if only binary32 were used, the vectorization factor would be 8.

This negative impact of cast operations was already mentioned in [6] or [18], for example. To overcome this issue, our choice is to move these cast operation once of all outside the loop.

The idea of moving cast operations outside the loop is not new, and it has already been implemented in the `vectoptim` tool manipulating LLVM IR [19]. However, we preferred to do this at the C level, which simplifies the process, and allows for breaking the dependency on the LLVM compiler.

### B. Splitting at the C level

Previous loop splitting, presented in [17], uses a compilation technique to split loops inside LLVM IR. This means that this work is fully dependent on the compiler in use. In this article, to overcome this point, we implement a new splitting tool working at C level, and based on Python. It is presented here for matrix multiplication, but it can be used for any iterative program. It uses hints inserted in the code to split the loop
<caption>Listing 4: Auto-tuned C program after loop splitting.</caption>

```c
void
matrix_multiplication(double *A, double *B,
                       double *C, int n) {
    // ALLOCATION
    float *A_b32, *B_b32, *C_b32;
    A_b32 = (float*)malloc((n*n)*sizeof(float));
    B_b32 = (float*)malloc((n*n)*sizeof(float));
    C_b32 = (float*)malloc((n*n)*sizeof(float));
    // START AUTO-TUNED LOOP
    assert((n-1)==255 && "Matrix size is wrong");
    // ... initialisation
    for(i = 0; i <= n-1; i += 1) {
        for(j = 0; j <= n-1; j += 1) {
            A_b32[i*n+j] = A[i*n+j];
            B_b32[i*n+j] = B[i*n+j];
        }
    }
    // ... loop id = {0}
    for(i = 0; i <= n-1; i += 1) {
        for(j = 0; j <= n-1; j += 1)
            C_b32[i*n+j] = C[i*n+j];
    }
    for (k = 0; k <= 63; k += 1) {
        for (i = 0; i <= n-1; i += 1) {
            for (j = 0; j <= n-1; j += 1)
                C_b32[i*n+j] += A_b32[i*n+k]*B_b32[k*n+j];
        }
    }
    for (i = 0; i <= n-1; i += 1) {
        for (j = 0; j <= n-1; j += 1)
            C[i*n+j] = C_b32[i*n+j];
    }
    // ... loop id = {1}
    for (k = 64; k <= n-1; k += 1) {
        for (i = 0; i <= n-1; i += 1) {
            for (j = 0; j <= n-1; j += 1)
                C[i*n+j] += A[i*n+k]*B[k*n+j];
        }
    }
    // END AUTO-TUNED LOOP
    // DEALLOCATION
    free(A_b32);
    free(B_b32);
    free(C_b32);
}
```

according to the output of `dd-fp2mp`. The C program for matrix multiplication is given in Listing 3.

The loop to be split is surrounded by tags `SPLIT_FOR` and `END SPLIT_FOR` tags, and hints are given in comment. Particularly, the following tags are used:

* **`SPLIT_FOR`**: It indicates that the surrounded loop should be a `for` loop, and it takes as parameters the induction variable, its start and end values, and the step value that the induction variable is updated at each iteration. If the surrounded loop does not match these parameters, the splitting is not performed.

* **`REPLACEMENT`**: It indicates the replacements to be made between data in binary64 and binary32 when reducing the computation precision in loop bodies. For example, in Listing 3, in loop whose computation precision is reduced, all occurrences of variable `A` will be replaced by `A_b32`. The latter variable is an array of size $n \times n$ allocated on the heap. The calls to `malloc`, respectively to `free`, are automatically inserted at the beginning, respectively the end, of the function.

* **`INITIALISATION`**: It corresponds to a C code that must be inserted before all the loops to initialize newly allocated lower precision variables. This corresponds to the cast moving presented in Section IV-A, only for variables that are read but not written in the loop body.

* **`PREFIX` / `SUFFIX`**: They correspond to C codes that must be inserted before (`PREFIX`) and after (`SUFFIX`) subloops whose precision has been lowered, and to handle the cast moving presented in Section IV-A, for variables written body but not necessarily read in the loop body.

For example, let $n = 256$, and the following splitting configuration stored as a list in Python.

```
[[0, 63, True], [64, 255, False]]
```

This means that the tagged loop is split into two subloops, where:

* subloop 1 iterates between 0 and 63, and the computation precision of its instructions is lowered to binary32,

* subloop 2 iterates between 64 and 255.

Using our new splitting tool, the resulting C program is given in Listing 4. Note that the declaration, allocation, and deallocation have been inserted once for all around the considered loop. The piece of code to ensure cast moving is inserted around the subloop whose computation precision is lowered, i.e. in our case, the first loop (`id = 0`). Note also that if `dd-fp2mp` returns a list of transformations indicating two consecutive subloops in the same precision, these are collapse in the generated programs to form a single loop. Note finally that the auto-tuning process is performed for particular input matrices, and it may no longer be valid for other input matrices. Hence an assert call has been inserted on top of the loop to ensure the sizes of the input matrices are correct.

### C. Expected performance

Before going further, this section illustrates the impact on performance that could be expected by using our new splitting tool. For doing this, we consider the size-256 matrix multiplication. Then for each loop (over $i$, $j$, and $k$), we use our splitting tool to split it into two subloops, and we apply binary64 to binary32 transformations on the first subloop, as the example above. For the loop over $k$, the result C program looks like the one in Listing 4. And finally, we vary the end index of the first subloop (and consequently, the start index of the second one) from 1 to 255 by step of 5 (i.e. by forcing at least one loop in each arithmetic). These programs are compiled using `gcc-11` with `-O3` and auto-vectorization activated. Figure 4 shows the average speedup computed as the ratio between full binary64 and auto-tuned programs on 100 random matrices, and executed on AVX2 micro-architectures.<fnref n="2" />

We observe that splitting loops over $i$ or $j$ degrade the performance of the generated programs. This is mainly due to the fact that the cast operations are moved around the tuned subloop. When we split the loop over $i$, for example, these cast operations are inserted around subloops over $i$, which remain inside the loop over $k$. Hence these cast operations are executed for each value of $k$, meaning in our case 256 times, which impacts negatively the performance of the program. These cast operations could be moved around the outer loop,

<footnote n="2">²In this article, for performance measure, the `assert` function call is disabled using the `NDEBUG` compilation option.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4: Expected speedup on auto-tuned programs for size-256 matrix multiplication.</figcaption></figure>

but this work is not straightforward, and it is not implemented in this article. The same holds for the loop over j. Conversely, the loop over k is already the outer loop. Therefore, moving the cast operations around this one does not impact too much the performance of the program, and we observe a speedup up to ≈ 1.7 compared with the untuned program. This justifies our interest only in the third splitting strategy.

## V. EXPERIMENTAL RESULTS

After some details on how to generate matrices with a pre-defined condition number, this section presents some experimental results on the performance of auto-tuned programs.

### A. Input matrix generation

Let A and B be two matrices of size n, and let κ be a predefined condition number. Recall that the matrix multiplication relies on $n^2$ dot-products $A_i \cdot B_j$, for $i, j \in \{0, \dots, n-1\}$, and where $A_i$, respectively $B_j$, is the $i+1$th row of A, respectively the $j+1$th column of B. But instead of considering the condition number of the input matrices A and B, we consider here the condition number of the dot-product $A_i \cdot B_j$. This approach was already used in [20] for matrix-vector multiplication and it is more relevant in this work, since we aim to split and optimize the computation of each dot-product according to an accuracy criterion.

Therefore we generate two matrices A and B, such that:

$$ \mathrm{cond}(A_i \cdot B_j) \le \kappa, \text{ for } i, j \in \{0, \dots, n-1\}. \quad (2) $$

For doing this, we use the GenDot routine defined in [21]. This routine can be easily extended to generate a matrix A and a vector $B_0$ for which the condition in (2) is satisfied. This is actually done in [20]. In our case, we generate such pair $(A, B_0)$, and we simply consider that all columns of B are identical, i.e. $B_i = B_0$ for all $i \in \{1, \dots, n-1\}$.

### B. Performance of auto-tuned programs

This first experiment deals with the performance of auto-tuned programs on SIMD micro-architectures. We consider the size-$n$ matrix multiplication with $n \in \{128, 256\}$, with a matrix of condition number $\kappa = 100$. The input program is implemented in binary64 floating-point arithmetic, and we consider reducing the precision of instructions to the binary32. In [16], we have shown the impact of the splitting factor. Hence we consider here a splitting factor in $\{32, 64, 128\}$. Recall that a splitting factor of 64, for example, means the

loop is split into 64 subloops, and the computation precision of each of them is decided through the auto-tuning process. Finally, as an accuracy constraint, defined in (1), we target a threshold $\epsilon$ ranging from $10^{-6}$ to $3 \cdot 10^{-4}$.

In each of these cases, we run our auto-tuning tool, and we generate optimized C programs. Then we evaluate their performance against the original program, using gcc-11 with -O3 and auto-vectorization activated. Figure 5 shows the average speedup we obtain on AVX2 and AVX512 micro-architectures on 300 performance measures, and the results of the auto-tuning, where the iterations in binary32 are in blue, and the ones in binary64 are in red.

Obviously, regarding the number of subloops transformed, as observed in [16], the tighter the precision constraint, the smaller the number of transformed subloops. For example, for $n = 128$ and a splitting factor of 32, on Figure 5a, until a threshold of $10^{-6}$, no subloop is transformed. Then from $2 \cdot 10^{-6}$, the number of transformed subloops increases, until a threshold of $8 \cdot 10^{-5}$, where all loops are transformed.

However, the important point of this experiment is that the speedup observed in practice is directly linked to the number of subloops transformed. Indeed, we can observe that, when the constraint gets less tight, the speedup increases, as well as the number of transformed loops. For example, still on Figure 5a, for $n = 128$ and a splitting factor of 32, when no subloop is transformed (threshold $\epsilon \le 10^{-6}$), a speedup of $\approx 1$ is observed. Then from $2 \cdot 10^{-6}$, the speedup increases until reaching $\approx 2$, when all subloops are transformed. Note also that, on these experiments, AVX2 micro-architectures provide slightly better performance gains than AVX512.

Further, we observe that the number of precision changes seems not impacting the performance of the generated programs. For example, for $n = 128$ and a threshold $\epsilon = 7 \cdot 10^{-5}$, a scaling factor of 32 leads to 6 changes of precision, while a scaling factor of 64 leads to 12 changes, i.e. twice as much. However, in both cases, we observe a speedup of $\approx 1.4$. This point is discussed in more details below in Section V-C.

### C. Impact of the number of precision changes

At each computation precision changes, that is, each time the program switches from a loop in binary64 to another one in binary32, or vice versa, data are duplicated and codes are executed on top and bottom of the loop in binary32, as explained in Section IV-A. This section illustrates the impact of the number changes on the performance of the programs.

For doing this, we consider matrices of size 256, with a splitting factor of 64. To observe the impact of the number of changes, we modify our auto-tuning tool to generate a transformation set where only at most $m \in \{1, 2\}$ changes appear. Note that $m = 1$ means that at most 1 precision change appears, that is, the generated program contains at most two subloops, one in binary64 and one in binary32. For each configuration, i.e. no maximum number of changes ($m = \infty$), and at most one and two changes, Table I shows the number of changes in the generated programs and the number of iterations in binary32 (for the loop over k), for a threshold
<figure><img src="image_6.png" /><figcaption>(a) $n = 128$, splitting factor = 32</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) $n = 128$, splitting factor = 64</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) $n = 128$, splitting factor = 128</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) $n = 256$, splitting factor = 32</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) $n = 256$, splitting factor = 64</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f) $n = 256$, splitting factor = 128</figcaption></figure>

<caption>Figure 5: Speedup and iteration patterns of programs, for $n \in \{128, 256\}$ and a splitting factor $\in \{32, 64, 128\}$.</caption>

from $10^{-6}$ to $3 \cdot 10^{-4}$ (which are thresholds for which most of changes appear). Finally Figure 6 shows the speedup obtained on these programs.

We observe that the threshold and the maximum number of changes allowed impact both the number of changes in the generated programs, and the number of iterations transformed in binary32. Moreover the number of changes appearing in the programs is linked to the number of iterations in binary32. For example, for a threshold of $3 \cdot 10^{-5}$, the auto-tuning tool finds a transformation set with 12 precision changes for $m = \infty$, but with only 1 change for $m = 1$ and 2 changes for $m = 2$. Regarding the number of iterations in binary32, we obtain 72, 40, 48 iterations, respectively, out of the 256 in total.

An important point is that, on a first hand, we observe in Figure 6 that this maximum number of changes impacts, in the same way, the speedup of the generated programs. Indeed, for example, for a threshold $\epsilon = 10^{-4}$, with $m = 1$, the speedup is slightly worst than for the other cases, i.e. $\approx 1.1$ for AVX2 and AVX512, while a speedup of $\approx 1.3$ is obtained with $m = 2$, and of $\approx 1.4$ with $m = \infty$. On the other hand, when the threshold decreases, the speedup increases more quickly for $m = \infty$ than for $m \in \{1, 2\}$. This comes from the fact that in this case, there are more opportunities for transformations, but these opportunities involve more than one or two changes.

### D. Impact of the condition number

The condition number allows for estimating the impact on the result of perturbations on the data. This section illustrates the impact of this condition number on the auto-tuning of matrix multiplication, and the performance of the result programs. Here we consider matrices of size 256, a splitting factor of 64, a maximum number of changes $m = 1$, and a condition

<figure><img src="image_12.png" /><figcaption>Figure 6: Speedup of auto-tuned programs, for $n = 256$, a splitting factor = 64, and a maximum number of changes $m \in \{1, 2\}$.</figcaption></figure>

number $\kappa \in \{100, 1000, 10000\}$. Figure 7 shows the speedup observed on the programs we obtained.

As expected, the higher the condition number, the fewer transformations are applied. For example, for $\kappa = 10000$, almost no transformation is applied, and we observe a speedup close to 1. This is due to the fact that these transformations imply data cast operations, which can be seen as data perturbations, and the higher the condition number, the more sensitive the computation of each dot-product is to these perturbations.

## VI. CONCLUSION

This article builds on previous work in the field of precision adaptation for iterative routines, particularly in the context of matrix multiplication. By extending the capabilities of an existing auto-tuning tool, we introduced a Python script that applies precision-reducing transformations to an input C program, generating a new, optimized version. Our experimental results demonstrate that this approach can achieve significant performance improvements, with speedups of up to
<table><caption>Table I: Number of changes and number of iterations in binary32, for $n = 256$, a splitting factor of 64, and $m \in \{\infty, 1, 2\}$.</caption><thead><tr><th></th><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th><th>15</th><th>16</th></tr></thead><tbody><tr><td rowspan="2">m = ∞</td><td>threshold ε</td><td>0</td><td>1</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>4</td><td>4</td><td>1</td><td>2</td><td>2</td><td>2</td><td>1</td><td>0</td></tr><tr><td># changes</td><td>0</td><td>4</td><td>8</td><td>8</td><td>12</td><td>12</td><td>12</td><td>12</td><td>16</td><td>16</td><td>32</td><td>12</td><td>6</td><td>10</td><td>14</td><td>10</td></tr><tr><td rowspan="2">m = 1</td><td># iterations in b32</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td></tr><tr><td># changes</td><td>0</td><td>4</td><td>4</td><td>4</td><td>8</td><td>8</td><td>8</td><td>8</td><td>8</td><td>8</td><td>24</td><td>40</td><td>48</td><td>52</td><td>52</td><td>256</td></tr><tr><td rowspan="2">m = 2</td><td># changes</td><td>0</td><td>1</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>32</td><td>2</td><td>2</td><td>2</td><td>2</td><td>0</td></tr><tr><td># iterations in b32</td><td>0</td><td>4</td><td>8</td><td>8</td><td>12</td><td>12</td><td>12</td><td>12</td><td>12</td><td>12</td><td>48</td><td>80</td><td>88</td><td>124</td><td>172</td><td>256</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 7: Speedup and iteration patterns of programs according to the condition number, for $n = 256$ and a splitting factor of 64.</figcaption></figure>

≈ 1.7 on SIMD micro-architectures. These findings highlight the potential of precision adaptation as a viable strategy for enhancing computational efficiency while maintaining the required accuracy, offering insights for developers aiming to optimize numerical computations in scientific programs.

Future work could explore the application of this method to other vectorizable codes, such as matrix inversion. Additionally, we aim to estimate the performance gain a *primary* to provide an option for users to request a targeted speedup. Furthermore, it could be interesting to investigate how this approach can be adapted using the blocking method for matrix multiplication. Such expansions and refinements could broaden the applicability and effectiveness of precision adaptation, making it a powerful tool in the optimization of computational performance across various domains.

## ACKNOWLEDGMENT

This work was supported by the PADOC ANR project under grant n°ANR-18-CE25-0004.

## REFERENCES

[1] IEEE Computer Society, *IEEE Standard for Floating-Point Arithmetic*, 2019.

[2] J. L. Gustafson and I. T. Yonemoto, “Beating Floating Point at Its Own Game: Posit Arithmetic,” *Supercomputing Frontiers and Innovations*, vol. 4, no. 2, p. 71–86, 2017.

[3] 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.

[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 *27th ACM SIGSOFT International Symposium on Software Testing and Analysis*. Association for Computing Machinery, 2018, p. 333–343.

[7] S. Graillat, F. Jézéquel, R. Picot, F. Févotte, and B. Lathuilière, “Autotuning for floating-point precision with Discrete Stochastic Arithmetic,” *Journal of Computational Science*, vol. 36, p. 101017, 2019.

[8] 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, pp. 369–378.

[9] M. O. Lam and J. K. Hollingsworth, “Fine-grained floating-point precision analysis,” vol. 32, no. 2, 2018, pp. 231–245.

[10] M. O. Lam, T. Vanderbruggen, H. Menon, and M. Schordan, “Tool integration for source-level mixed precision,” in *2019 IEEE/ACM 3rd International Workshop on Software Correctness for HPC Applications (Correctness)*, 2019, pp. 27–35.

[11] H. Menon, M. O. Lam, D. Osei-Kuffuor, M. Schordan, S. Lloyd, K. Mohror, and J. Hittinger, “ADAPT: Algorithmic Differentiation Applied to Floating-Point Precision Tuning,” in *International Conference for High Performance Computing, Networking, Storage, and Analysis*. IEEE Press, 2018.

[12] A. Zeller and R. Hildebrandt, “Simplifying and Isolating Failure-Inducing Input,” *IEEE Transactions on Software Engineering*, pp. 183–200, 2002.

[13] 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.

[14] H. Brunie, C. Iancu, K. Z. Ibrahim, P. Brisk, and B. Cook, “Tuning Floating-Point Precision Using Dynamic Program Information and Temporal Locality,” in *SC20: International Conference for High Performance Computing, Networking, Storage and Analysis*, 2020, pp. 1–14.

[15] R. Gu, P. Beata, and M. Becchi, “A Loop-Aware Autotuner for High-Precision Floating-Point Applications,” in *International Symposium on Performance Analysis of Systems and Software (ISPASS 2020)*, 2020, pp. 285–295.

[16] Y. Fakhreddine and G. Revy, “Using loop transformations for precision tuning in iterative programs,” in *30th IEEE Symposium on Computer Arithmetic (ARITH 2023)*, 2023, pp. 159–166.

[17] G. Revy, “Analyzing the impact of floating-point precision adaptation in iterative programs,” in *2021 IEEE 28th Symposium on Computer Arithmetic (ARITH)*, 2021, pp. 25–32.

[18] ——, “Contribution to the automated generation of programs in computer arithmetic,” Habilitation à diriger des recherches, Université de Perpignan Via Domitia, 2023.

[19] H. L. Boudi, P. Langlois, and G. Revy, “Increasing the impact of vectorization by automatic transformation of programs,” Tech. Rep., 2020, unpublished (available upon request).

[20] C. Chohra, “Towards Reproducible, Accurately Rounded and Efficient BLAS,” Ph.D. dissertation, Université de Perpignan Via Domitia, 2017.

[21] T. Ogita, S. M. Rump, and S. Oishi, “Accurate Sum and Dot Product,” *SIAM Journal on Scientific Computing*, vol. 26, no. 6, pp. 1955–1988, 2005.