Clément Rossetti
email: clement.rossetti@inria.fr

Philippe Clauss
email: philippe.clauss@inria.fr

Algebraic Tiling

Keywords: Loop tiling, load balancing, parallel loop, data locality, optimizing compilers

In this paper, we present an ongoing work whose aim is to propose a new loop tiling technique where tiles are characterized by their volumes -the number of embedded iterations -instead of their sizes -the lengths of their edges. Tiles of quasi-equal volumes are dynamically generated while the tiled loops are running, whatever are the original loop bounds, which may be constant or depending linearly of surrounding loop iterators. The adopted strategy is to successively and hierarchically slice the iteration domain in parts of quasi-equal volumes, from the outermost to the innermost loop dimensions. Since the number of such slices can be exactly chosen, quasi-perfect load balancing is reached by choosing, for each parallel loop, the number of slices as being equal to the number of parallel threads, or to a multiple of this number. Moreover, the approach avoids partial tiles by construction, thus yielding a perfect covering of the iteration domain minimizing the loop control cost. Finally, algebraic tiling makes dynamic scheduling of the parallel threads fairly purposeless for the handled parallel tiled loops.

Introduction

Loop tiling [START_REF] Dongarra | Automatic Blocking of Nested Loops[END_REF][START_REF] Irigoin | Supernode Partitioning[END_REF][START_REF] Lam | A Data Locality Optimizing Algorithm[END_REF][START_REF] Ramanujam | Tiling Multidimensional Iteration Spaces for Nonshared Memory Machines[END_REF][START_REF] Wolfe | Iteration Space Tiling for Memory Hierarchies[END_REF] (variously called loop blocking or partitioning) is a well-known loop optimizing transformation often providing significant speed-ups. It explicitly addresses data locality by contracting the accessed data space of inner loops in order to promote good cache reuse. When handling parallel loops, it also provides a way to adjust the grain of parallelism or to exploit vectorization while avoiding register pressure on a single core. Loop tiling is an essential strategy used by compilers and automatic parallelizers.

However, as it was pointed out by Sato et al. in [START_REF] Sato | An Autotuning Framework for Scalable Execution of Tiled Code via Iterative Polyhedral Compilation[END_REF], traditional loop tiling techniques do not address load balancing explicitly when handling parallel loops, thus yielding bad scalability for parallelism. Moreover, rectangular tiles of constant size may generally not cover perfectly the iteration domain, since tiles overlapping borders of the domain result in partial tiles containing less iterations than the full tiles. Iooss et al. in [START_REF] Iooss | Monoparametric Tiling of Polyhedral Programs[END_REF] and Renganarayanan et al. in [START_REF] Renganarayanan | Parameterized Loop Tiling[END_REF] address the limitations of fixed size tiling. They point out that the tile IMPACT, Jan. 16, 2023, Toulouse, FR 2023. sizes have a huge influence on performance, and that parametric tiling in its full generality is known to be non-linear. Allowing tile sizes to be symbolic parameters at compile time has many benefits, including efficient auto-tuning, and run-time adaptability to system variations [START_REF] Iooss | Monoparametric Tiling of Polyhedral Programs[END_REF][START_REF] Renganarayanan | Parameterized Loop Tiling[END_REF].

In the related literature, load balancing is sometimes addressed, but from a different perspective and for the particular class of loops that implement stencil computations. While stencil loops mostly require a skewing transformation to be amenable to parallelization [START_REF] Andonov | Optimal Semioblique Tiling[END_REF][START_REF] Lam | A Data Locality Optimizing Algorithm[END_REF][START_REF] Wonnacott | Achieving Scalable Locality with Time Skewing[END_REF], dedicated techniques for data locality optimization and concurrent start for tiles are proposed. Some noticeable proposals consist in generating tiles with a specific shape (diamond tiling [START_REF] Bondhugula | Diamond Tiling: Tiling Techniques to Maximize Parallelism for Stencil Computations[END_REF], hexagonal tiling [START_REF] Grosser | Hybrid Hexagonal/Classical Tiling for GPUs[END_REF], or split and overlapped tiling [START_REF] Henretty | A Stencil Compiler for Shortvector SIMD Architectures[END_REF][START_REF] Krishnamoorthy | Effective Automatic Parallelization of Stencil Computations[END_REF]).

We propose a new tiling approach based on the volumes of the tiles, i.e., the number of iterations delimited by the tiles, instead of the sizes of standard (hyper-)rectangular tiles, i.e., the sizes of the edges of the tiles. In the proposed approach, tiles are dynamically generated and have almost equal volumes, even if their shape and edge sizes may differ. The iteration domain is well covered by a minimum number of tiles that are all almost full. Since the bounds of the generated tiles are not linear and defined by algebraic mathematical expressions, we call this loop tiling technique algebraic tiling.

Algebraic tiles are built by successive hierarchical slicing of the initial iteration domain, from the outermost to the innermost depth dimensions of the target loop nest, in a way ensuring that slices have all quasi-equal volumes. The bounds of the loop nests that are handled must be constants, or linear functions of the surrounding loop iterators and of unknown parameters -which are typically related to the data input size. Such loops are also called polyhedral loops since they may be handled using the polyhedral model [START_REF] Feautrier | Polyhedron Model[END_REF]. Quasi-perfect load balancing is achieved when each parallel loop is sliced using as many slices of quasi-equal volumes as parallel threads, and when most of the iterations have close execution times. Thus, such dynamic slicing strategy makes the resulting parallel loop scalable regarding the number of threads. Good data locality is reached by slicing profitably the non-parallelized loops, and by slicing the parallel loops in a number of parts equal to a multiple of the number of parallel threads. The number of generated slices for each dimension may stay as a parameter at compile-time, making algebraic tiling a parameterized loop tiling technique, allowing the produced code to adapt to the number of parallel threads and data layout.

We are currently implementing algebraic tiling as an extension of the automatic loop optimizer Pluto 1 [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]. Our first experiments show that algebraic tiling outperforms significantly (hyper-)rectangular tiling when parallelizing loops with OpenMP using static scheduling, and mostly provides similar or lower execution times when compared to traditionally tiled loops parallelized using dynamic scheduling of OpenMP. Thus, algebraic tiling makes dynamic scheduling fairly purposeless for the handled loop nests.

The paper is organized as follows. In the next section, we first describe the optimizing strategy of algebraic tiling and how it may be applied to any compliant loop nest. Then in Section 3, we introduce the mathematical background required for computing the bounds of algebraic tiles. Several issues regarding the applicability of algebraic tiling are addressed in Section 4. Our software that automatically generates the functions required for using algebraic tiles is described in Section 5. Experiments showing significant speedups against the (hyper-)rectangular tiling approach are exhibited in Section 6. Related work is addressed in Section 7 and conclusions are given in Section 8.

Algebraic tiling in practice

In this section, we describe algebraic tiling regarding its related programming aspects. Its mathematical foundations are presented in the next section.

Overview

Consider the loop nest shown in Figure 1, which is extracted from program syr2k of the polybench benchmark suite 2 v4.2. This loop kernel defines a triangular iteration domain. When tiling this loop nest using the standard approach, the iteration domain is partitioned into fixed-sized rectangular tiles. When also parallelizing the outermost loop among a given number of threads, each thread is traditionally assigned a slice of the iteration domain, all slices being of the same width, as illustrated in Figure 2. Two main issues can be highlighted with this example: (1) many partial tiles occur at the borders of the iteration domain and (2) slices assigned to the threads are significantly unbalanced regarding their iteration counts. In contrast, algebraic tiling first slices the outermost loop into parts of quasi-equal volumes 𝑉 𝑖 , as shown in Figure 3a, where

𝑉 0 ≃ 𝑉 1 ≃ 𝑉 2 ≃ 𝑉 3 ≃ 𝑉 4 .
Then, these slices are sliced in turn along the inner loop direction in parts of quasi-equal volumes 𝑉 𝑖 𝑗 , as shown in Figure 3b. Hence quasi-perfect load balancing is reached when parallelizing the outermost loop, while good data locality is achieved thanks to inner slicing and perfect domain covering.

A sample of the algebraic tile bounds that are dynamically generated is shown in Table 1, when slicing the outermost loop in 24 parts, and the innermost in 64 parts, and with 1 http://pluto-compiler.sourceforge.net 2 https://sourceforge.net/projects/polybench The first column of each of the three parts of the table shows the 24 outermost slices and their respective lower and upper bounds, 𝑙𝑏 𝑖 and 𝑢𝑏 𝑖 , and the resulting slice volumes (iteration count). For each outermost slice, the second column shows a sample of the related 64 inner slices and their respective lower and upper bounds, 𝑙𝑏 𝑗 and 𝑢𝑏 𝑗 , and the resulting tile volumes. Each resulting tile contains a number of iterations which is as close as possible to 𝐿𝑇𝐶/24/64 = 469140, where 𝐿𝑇𝐶 denotes the total loop trip count of the syr2k loop kernel. Note that all tiles have different sizes and may be non-rectangular along the borders of the iteration domain.

f o r ( i = 0 ; i < N ; i ++ ) f o r ( j = 0 ; j <= i ; j ++ ) f o r ( k = 0 ; k < M; k ++ ) C[ i ] [ j ] += A[ j ] [ k ] * a l p h a * B [ i ] [ k ]+ B [ j ] [ k ] * a l p h a * A[ i ] [ k ] ;
i j V 0 V 1 V 2 V 3 V 4 (

Volumes of algebraic tiles

The target volumes of algebraic tiles are approached using dividers, 𝐷𝐼𝑉 1 , 𝐷𝐼𝑉 2 , etc., that are the target number of slices associated to every loop nest depth. When parallelizing the outermost loop of a loop nest, the outermost divider 𝐷𝐼𝑉 1 is typically the number of parallel threads, and the volumes of the outermost slices must be as close as possible to 𝐿𝑇𝐶/𝐷𝐼𝑉 1 . However, for data locality optimization reasons, it may be beneficial in some cases to set 𝐷𝐼𝑉 1 as equal to a multiple of the number of threads. If so, the OpenMP static schedule still yields a balanced distribution of the so-generated smaller slices of quasi-equal volumes. The next dividers for inner slices are selected to reach good cache performance. Note also that nested parallelism may also be managed efficiently by selecting next dividers according to the number of processing units of each hierarchical parallel level.

Code transformation

On Figure 4, it is shown how the original syr2k loop kernel of Figure 1 The implementation of algebraic tiling is mainly based on functions that provide the values of loop iterators where a given number of iterations has been completed (trahrhe_i and trahrhe_j in the current example). Such functions evaluate algebraic expressions using floating-point arithmetic. These algebraic expressions are symbolic roots of multivariate polynomial equations. The automatic generation of such functions is detailed in the next section.

Trahrhe expressions

Name «trahrhe» comes from «Ehrhart» with letters in reversed order. Trahrhe expressions are roots of Ehrhart polynomials, and they are also the inverses of particular Ehrhart polynomials which are called ranking polynomials. Note that these new mathematical objects were already used to collapse non-rectangular loops in [START_REF] Clauss | Automatic Collapsing of Non-Rectangular Loops[END_REF].

Ranking polynomials

Ehrhart polynomials were proposed and extended to program analysis by Clauss in [START_REF] Clauss | Counting Solutions to Linear and Nonlinear Constraints Through Ehrhart Polynomials: Applications to Analyze and Transform Scientific Programs[END_REF]. These integer-valued polynomials express the exact number of integer points contained in a finite multi-dimensional convex polyhedron which depends linearly on integer parameters. They have many applications for the quantitative analysis of loop nests whose loop bounds are linear functions of the surrounding loop indices and integer parameters, and whose statements are referencing multi-dimensional array elements through linear functions of the loop indices and parameters. Such a counting of integer points may translate to the exact number of iterations of a parameterized loop nest, the exact number of memory locations touched by a loop nest, the maximum number of parallel iterations, etc. When considering a 𝑑-dimensional polyhedron -as for example the iteration domain of a 𝑑-depth loop nest -depending linearly on integer parameters 𝑝 1 , 𝑝 2 , . . . , 𝑝 𝑚 , its Ehrhart polynomial is a polynomial of degree 𝑑 whose variables are 𝑝 1 , 𝑝 2 , . . . , 𝑝 𝑚 . Ehrhart polynomials can be automatically computed using existing algorithm implementations as the one of the barvinok library [START_REF] Verdoolaege | Counting Integer Points in Parametric Polytopes Using Barvinok's Rational Functions[END_REF].

Among their applications, Ehrhart polynomials are used by Clauss and Meister in [START_REF] Clauss | Automatic memory layout transformations to optimize spatial locality in parameterized loop nests[END_REF] to reorganize the memory layout of array elements accessed by a loop nest, in order to improve their spatial data locality: array elements are relocated in memory in the same order as they are accessed. In this approach, the new location of an array element is given by the order, or rank, of the iteration referencing it.

Such a rank of iterations is given by a polynomial, called a ranking polynomial, whose variables are the loop iterators, and whose evaluation results in the number of iterations preceding a given iteration. More formally, the ranking polynomial of a 𝑑-depth loop nest whose loop iterators, from the outermost to the innermost, are (𝑖 1 , 𝑖 2 , . . . , 𝑖 𝑑 ), is denoted 𝑟 (𝑖 1 , 𝑖 2 , . . . , 𝑖 𝑑 ). Without loss of generality, if (0, 0, . . . , 0) defines the first iteration of the loop nest, then 𝑟 (0, 0, . . . , 0) = 1, 𝑟 (0, 0, . . . , 1) = 2, and so on. If (𝑁 1 , 𝑁 2 , . . . , 𝑁 𝑑 ) are the indices of the very last iteration, then 𝑟 (𝑁 1 , 𝑁 2 , . . . , 𝑁 𝑑 ) is the total number of iterations of the loop nest (𝐿𝑇𝐶).

The computation of the ranking polynomial of a loop nest is detailed in [START_REF] Clauss | Automatic memory layout transformations to optimize spatial locality in parameterized loop nests[END_REF]. We recall this technique by applying it to the syr2k loop kernel of Figure 1, to compute the associated ranking polynomial 𝑟 (𝑖, 𝑗, 𝑘).

The rank of a given iteration (𝑖 0 , 𝑗 0 , 𝑘 0 ) is equal to the number of iterations that are executed before (𝑖 0 , 𝑗 0 , 𝑘 0 ) (included), i.e., the number of triplets (𝑖, 𝑗, 𝑘) inside the iteration domain which are lexicographically less than or equal to (𝑖 0 , 𝑗 0 , 𝑘 0 ):

∀(𝑖 0 , 𝑗 0 , 𝑘 0 ) s.t. 0 ≤ 𝑖 0 < 𝑁 and 0 ≤ 𝑗 0 ≤ 𝑖 0 and 0 ≤ 𝑘 0 < 𝑀, 𝑟 (𝑖 0 , 𝑗 0 , 𝑘 0 ) = #{(𝑖, 𝑗, 𝑘) | (𝑖, 𝑗, 𝑘) ⊴ (𝑖 0 , 𝑗 0 , 𝑘 0 ), 0 ≤ 𝑖 < 𝑁 , 0 ≤ 𝑗 ≤ 𝑖, 0 ≤ 𝑘 < 𝑀 }
where ⊴ denotes the lexicographic order. Since lexicographic inequalities are not linear, the problem is split as the conjunction of three equivalent sets of linear inequalities, according to the definition of the lexicographic order:

(𝑖, 𝑗, 𝑘) ⊴ (𝑖 0 , 𝑗 0 , 𝑘 0 ) ⇔(𝑖 < 𝑖 0 ) or (𝑖 = 𝑖 0 and 𝑗 < 𝑗 0 ) or (𝑖 = 𝑖 0 and 𝑗 = 𝑗 0 and 𝑘 ≤ 𝑘 0 )

Therefore, the sets whose integer points must be counted can be defined as the union of three disjoint convex polyhedra, and 𝑟 (𝑖 0 , 𝑗 0 , 𝑘 0 ) as the sum of three Ehrhart polynomials:

𝑟 (𝑖 0 , 𝑗 0 , 𝑘 0 ) = #{(𝑖, 𝑗, 𝑘) | 0 ≤ 𝑖 < 𝑖 0 , 0 ≤ 𝑗 ≤ 𝑖, 0 ≤ 𝑘 < 𝑀 } + #{(𝑖, 𝑗, 𝑘) | 𝑖 = 𝑖 0 , 0 ≤ 𝑗 < 𝑗 0 , 0 ≤ 𝑘 < 𝑀 } + #{(𝑖, 𝑗, 𝑘) | 𝑖 = 𝑖 0 , 𝑗 = 𝑗 0 , 0 ≤ 𝑘 ≤ 𝑘 0 } = 𝑀 𝑖 0 (𝑖 0 + 1) 2 + 𝑀 𝑗 0 + 𝑘 0 + 1 = 2 𝑘 0 + 2 𝑀 𝑗 0 + 𝑀 𝑖 2 0 + 𝑀 𝑖 0 + 2 2
One can verify that the rank of the first iteration (0, 0, 0), 𝑟 (0, 0, 0), is equal to 1, the rank of the second iteration 𝑟 (0, 0, 1) = 2, the rank of the third iteration 𝑟 (0, 0, 2) = 3 and so on. The rank of the last 𝑘-iteration when 𝑖 = 0 and 𝑗 = 0, 𝑟 (0, 0, 𝑀 -1) = 𝑀, and the rank of the first iteration when

= 0, 𝑟 (1, 0, 0) = 𝑀 + 1. The total number of iterations is 𝐿𝑇𝐶 = 𝑟 (𝑁 -1, 𝑁 -1, 𝑀 -1) = 𝑀 𝑁 (𝑁 +1) 2 .
Such a ranking polynomial associates, to each iteration index tuple, a unique integer of a continuous interval of integers starting at 1. This continuous interval is the range of integers between one and the total number of iterations. Conversely, each integer value in the interval is associated to one unique iteration index tuple. The ranking polynomial can also be seen as the one-dimensional polynomial schedule function of the iterations, which is equivalent to the original multi-dimensional linear schedule defined by the nested loops. Another important property is that such a ranking polynomial is monotonically increasing over the integers, from 1 to the total number of iterations, relatively to the lexicographic order of the loop indices. Thus, a ranking polynomial defines a bijection between the iteration domain and the interval of successive integers. It implies that in theory, it can be inverted. Such inversions result in expressions called trahrhe expressions, which are used to compute the bounds of algebraic tiles.

Generating trahrhe expressions

For a given rank 𝑝𝑐, 𝑡𝑟𝑎ℎ𝑟ℎ𝑒 (𝑝𝑐) is the tuple (𝑡 1 , 𝑡 2 , ..., 𝑡 𝑑 ) such that 𝑟 (𝑡 1 , 𝑡 2 , ..., 𝑡 𝑑 ) = 𝑝𝑐, i.e., a solution of the latter equation. The definition of function 𝑡𝑟𝑎ℎ𝑟ℎ𝑒 (𝑝𝑐) is calculated incrementally by first determining 𝑡 1 , then propagating it to determine 𝑡 2 , and so on. At each step, a symbolic uni-variate polynomial equation is solved:

Find 𝑡 1 : Solve for 𝑠: 𝑟 (𝑠, 0, ..., 0) -𝑝𝑐 = 0. Depending on the degree 𝑑 of the ranking polynomial, this equation may have 𝑑 real or complex solutions 𝑠 1 , ..., 𝑠 𝑑 . Among these solutions, only one solution 𝑠 𝑘 is such that 𝑡 1 = ⌊𝑠 𝑘 ⌋ = 0 when 𝑝𝑐 = 1, which means that for the very first iteration of the 𝑑-depth loop nest (𝑝𝑐 = 1), the value of the outermost loop iterator 𝑡 1 is, as expected, 0. This solution is propagated in the next equation. Find 𝑡 2 : Solve for 𝑠: 𝑟 (𝑡 1 , 𝑠, 0, ..., 0) -𝑝𝑐 = 0. Once again, among the 𝑑 -1 solutions, only one solution 𝑠 𝑙 is such that 𝑡 2 = ⌊𝑠 𝑙 ⌋ = 0 when 𝑝𝑐 = 1. This solution is propagated in the next equation. Find 𝑡 3 : Solve for 𝑠: 𝑟 (𝑡 1 , 𝑡 2 , 𝑠, 0, ..., 0) -𝑝𝑐 = 0. ... Find 𝑡 𝑑 : Solve for 𝑠: 𝑟 (𝑡 1 , 𝑡 2 , 𝑡 3 , ..., 𝑠) -𝑝𝑐 = 0. This last equation is obviously linear: 𝑡 𝑑 = 𝑝𝑐 -𝑟 (𝑡 1 , 𝑡 2 , ..., 0). As an example, let us compute step by step the trahrhe expressions of the syr2k loop kernel:

1. 𝑟 (𝑖, 𝑗, 𝑘) = 2 𝑘+2 𝑀 𝑗+𝑀 𝑖 2 +𝑀 𝑖+2 2 2. Solve 𝑟 (𝑠, 0, 0) -𝑝𝑐 = 𝑀 𝑠 2 +𝑀 𝑠+2

2

-𝑝𝑐 = 0. This equation has two solutions:

𝑠 1 = - √ 8 𝑀 pc + 𝑀 2 -8 𝑀 + 𝑀 2 𝑀 , 𝑠 2 = √ 8 𝑀 pc + 𝑀 2 -8 𝑀 -𝑀 2 𝑀 When 𝑝𝑐 = 1, 𝑠 2 = 0. Thus 𝑡 1 = √ 8 𝑀 pc+𝑀 2 -8 𝑀 -𝑀 2 𝑀 3. Solve 𝑟 (𝑡 1 , 𝑠, 0) -𝑝𝑐 = 𝑀 𝑡 2 1 +𝑀 𝑡 1 +2 𝑀 𝑠+2 2
-𝑝𝑐 = 0. This equation has one solution. Thus 𝑡 2 = -

𝑀 𝑡 2 1 +𝑀 𝑡 1 -2 pc+2 2 𝑀 4. Finally, 𝑡 3 = 𝑝𝑐 -𝑟 (𝑡 1 , 𝑡 2 , 0) = - 2 𝑀 𝑡 2 +𝑀 𝑡 2 1 +𝑀 𝑡 1 -2 pc+2 2
One can verify that for rank 𝑝𝑐 = 2, the trahrhe expressions result in (0, 0, 1) as expected:

• 𝑡 1 = √ 𝑀 (𝑀+8) -𝑀 2 𝑀
= 0 for any 𝑀 > 1

• 𝑡 2 = 2 2 𝑀 = 0 for any 𝑀 > 1 • 𝑡 3 = 2 2 =
1 For rank 𝑝𝑐 = 𝑀, the trahrhe expressions result in (0, 0, 𝑀 -1) as expected:

• 𝑡 1 = √ 𝑀 (9 𝑀 -8) -𝑀 2 𝑀 = 0 for any 𝑀 > 1 since 0 < √ 𝑀 (9 𝑀 -8) -𝑀 2 𝑀 < √ 9 𝑀 2 -𝑀 2 𝑀 = 1 • 𝑡 2 = 1 -1 𝑀 = 0 for any 𝑀 ≥ 1 • 𝑡 3 = 2 𝑀 -2 2 = 𝑀 -1
Note that since the degree of ranking polynomials is less or equal to the depth of the target loop nest, their roots may include radicals (square roots, cubic roots, ...), and such radicals may result in complex numbers. However, the imaginary parts of the entire numerical evaluations of trahrhe expressions are always null. Nevertheless, it means that evaluations must be performed using complex floating-point arithmetic.

Only polynomial equations whose degree is at most equal to 4 can be solved symbolically, with exact expressions for roots. The handled uni-variate polynomial equations are built from a multi-variate ranking polynomial, where one index 𝑖 𝑘 is set as the equation unknown, indices 𝑖 1 , ..., 𝑖 𝑘 -1 are set as symbolic parameters, and indices 𝑖 𝑘+1 , ..., 𝑖 𝑑 are set to their lexicographic minimum values. Thus, to ensure that such a built equation has a degree less than 4, the ranking polynomial must be such that any index 𝑖 𝑘 , in any of its monomials, has a degree less than 4, i.e., any monomial is of the form: 𝑎 𝑖

𝑝 1 1 𝑖 𝑝 2 2 ...𝑖 𝑝 𝑑
𝑑 where 𝑎 is a rational number, and every power 𝑝 𝑘 is such that 0 ≤ 𝑝 𝑘 ≤ 4.

Loop nests yielding such ranking polynomials are such that the maximum number of nested loops, whose loop trip counts depend on a given index 𝑖 𝑘 , is less than or equal to 4. For example, both outermost loops of the syr2k loop kernel in Figure 1 depend on index 𝑖, yielding a ranking polynomial where index 𝑖 is of power 2 in some monomial.

However, note that the dependence of loop trip counts regarding indices is transitive: if a loop index 𝑗 depends on a surrounding loop index 𝑖, and an inner loop index 𝑘 depends on 𝑗, then 𝑘 depends also on 𝑖: index 𝑖 is of power 3, in some monomial of the ranking polynomial. Note also that loops may depend simultaneously on several surrounding loops' indices, as for example in for(k=0;k<i+j;k++).

Hence, the loop nests that can be handled by our method may be of any depth, but are such that the number of nested loops that all depend on a given index is less than or equal to 4. This should be quite sufficient for most cases, regarding the usual loop nest depth and complexity of user codes. This equation solving process has been implemented in our software computing trahrhe expressions. In addition, several issues regarding arithmetic precision and mathematical generalization had to be solved. They are detailed in Section 5.

Bounds of algebraic tiles

As shown in Section 2, the bounds of algebraic tiles are resulting from slicing the iteration domain into parts of quasiequal volumes, these latter being sliced in turn into parts of quasi-equal volumes, and so on, from the outermost to the innermost loop dimensions.

The slicing process is driven by a target volume for the resulting parts. Nevertheless, this volume can generally never be reached exactly, since slices are built along one unique dimension, as it is illustrated by Figure 5. Moreover, the target volume cannot be less than the largest cutting hyperplane, since smallest possible slices have their lower bounds equal to their upper bounds. In the proposed approach, the target volume results from dividing the total number of iterations -which is the evaluation of an Ehrhart polynomial -by a chosen divider, as explained in Section 2. Multi-dimensional tiling is achieved dimension per dimension, from the outermost to the innermost loop depth. When handling a given dimension, the associated iterator is considered as if it was the outermost iterator, i.e., as the high-order index regarding the lexicographic order. Thus, a different lexicographic order of the iterator is considered at each step. This is illustrated by Figure 6, where it is shown on the iteration domain of kernel syr2k, that inner slices bounds must be computed "as if" the outer slices were scanned relatively to an interchanged lexicographic order.

More formally, consider a 𝑑-depth loop nest whose loop iterators are (𝑖 1 , 𝑖 2 , . . . , 𝑖 𝑑 ), from the outermost to the innermost. Let us denote by 𝐷 its iteration domain. First, 𝐷 is sliced regarding dimension 𝑖 1 . For this dimension, the lexicographic order defined by the original loop nest is considered. The associated ranking polynomial 𝑟 (𝑖 1 , 𝑖 2 , . . . , 𝑖 𝑑 ) is computed and inverted as explained in Subsection 3.2, to get the trahrhe expression, 𝑡𝑟𝑎ℎ𝑟ℎ𝑒_𝑖 1 (𝑝𝑐), defining the value of 𝑖 1 where 𝑝𝑐 iterations have been completed. The target volume, named 𝑇 𝐼𝐿𝐸_𝑉 𝑂𝐿_𝐿1, is computed by dividing the total number of iterations by a given divider 𝐷𝐼𝑉 1 . Lower and upper bounds of slices are then roughly defined as being:

𝑙𝑏 𝑖 1 (𝑖𝑡 1 ) = 𝑡𝑟𝑎ℎ𝑟ℎ𝑒_𝑖 1 (𝑖𝑡 1 × 𝑇 𝐼𝐿𝐸_𝑉 𝑂𝐿_𝐿1) 𝑢𝑏 𝑖 1 (𝑖𝑡 1 ) = 𝑡𝑟𝑎ℎ𝑟ℎ𝑒_𝑖 1 ((𝑖𝑡 1 + 1) × 𝑇 𝐼𝐿𝐸_𝑉 𝑂𝐿_𝐿1) -1 where 𝑖𝑡 1 ∈ [0...𝐷𝐼𝑉 1 -1].
Next, each resulting slice characterized by a value 𝑖𝑡 1 is sliced in turn regarding dimension 𝑖 2 . However, since each slice 𝑖𝑡 1 has to be sliced along dimension 𝑖 2 , 𝑖 2 has now to be considered as the outermost iterator regarding the lexicographic order. Hence, for each outermost slice defined by

𝐷 ∩ {𝑙𝑏 𝑖 1 (𝑖𝑡 1 ) ≤ 𝑖 1 ≤ 𝑢𝑏 𝑖 1 (𝑖𝑡 1
)}, the associated ranking polynomial 𝑟 (𝒊 2 , 𝑖 1 , 𝑖 3 , . . . , 𝑖 𝑑 ) is computed and inverted, to get the trahrhe expression, 𝑡𝑟𝑎ℎ𝑟ℎ𝑒_𝑖 2 (𝑝𝑐, 𝑙𝑏 𝑖 1 (𝑖𝑡 1 ), 𝑢𝑏 𝑖 1 (𝑖𝑡 1 )), defining the value of 𝑖 2 where 𝑝𝑐 iterations have been completed. The target volume 𝑇 𝐼𝐿𝐸_𝑉 𝑂𝐿_𝐿2 is computed by dividing the total number of iterations of the current outer slice 𝑖𝑡 1 by a given divider 𝐷𝐼𝑉 2 . Lower and upper bounds of slices are then roughly defined as being:

𝑙𝑏 𝑖 2 (𝑖𝑡 2 , 𝑙𝑏 𝑖 1 (𝑖𝑡 1 ), 𝑢𝑏 𝑖 1 (𝑖𝑡 1 )) = 𝑡𝑟𝑎ℎ𝑟ℎ𝑒_𝑖 2 𝑖𝑡 2 × 𝑇 𝐼𝐿𝐸_𝑉 𝑂𝐿_𝐿2, 𝑙𝑏 𝑖 1 (𝑖𝑡 1 ), 𝑢𝑏 𝑖 1 (𝑖𝑡 1 ) 𝑢𝑏 𝑖 2 (𝑖𝑡 2 , 𝑙𝑏 𝑖 1 (𝑖𝑡 1 ), 𝑢𝑏 𝑖 1 (𝑖𝑡 1 )) = 𝑡𝑟𝑎ℎ𝑟ℎ𝑒_𝑖 2 ((𝑖𝑡 2 + 1) × 𝑇 𝐼𝐿𝐸_𝑉 𝑂𝐿_𝐿2, 𝑙𝑏 𝑖 1 (𝑖𝑡 1 ), 𝑢𝑏 𝑖 1 (𝑖𝑡 1 ) -1 where 𝑖𝑡 2 ∈ [0...𝐷𝐼𝑉 2 -1].
The same process is repeated until the deepest tiling dimension, by considering at each step a different lexicographic order, where the current sliced dimension index is set as the outermost iterator.

Applicability of algebraic tiling

Regarding its validity, algebraic tiling follows the same rules as standard tiling: algebraic tiling may be applied when it is valid regarding data dependences. Thus a prior affine transformation of the target loop nest may be required. In practice, we use the dependence analysis process implemented in the automatic loop optimizer Pluto. Moreover, to promote vectorization and data locality, loops may be interchanged, as it is performed by Pluto.

Regarding vectorization of the innermost loops, the nonconstant algebraic tile bounds may prevent compilers to generate the fastest vector code. In such cases, a hybrid tiling approach that combines standard and algebraic tiling seems the best strategy for getting better code : the outer dimension that is linked to the innermost vectorizable loop is actually sliced using a constant value (e.g. 32 or 64), as it is achieved with standard tiling, while the other dimensions are sliced with non-constant algebraic bounds. In the future, we expect that algebraic tiling will take advantage of hardware extensions as SVE (Scalable Vector Extension) [START_REF] Stephens | The ARM Scalable Vector Extension[END_REF] to result in more efficient vector code.

Algebraic tile bounds are obviously not linear and may stay parameterized by the dividers and the problem size. Parameterized tiling in the standard approach also yields non-linear bounds, but those are polynomials while ours are algebraic expressions. In [START_REF] Iooss | Monoparametric Tiling of Polyhedral Programs[END_REF], Iooss et al. propose parametric tiling based on one unique parameter, which enables further linear code transformations. By contrast, algebraic tiling prevents any chance for further linear transformations, but opens the field to managed non-linear transformations.

Algebraic tiling is also a parametric tiling approach whose parameters are iteration count dividers rather than tile sizes. The importance of parametric tiling is exemplified by the ATLAS system [START_REF] Whaley | Automatically Tuned Linear Algebra Software[END_REF] for empirical tuning of BLAS kernels. ATLAS uses parametrically tiled BLAS kernels that are repeatedly executed on the target architecture for different problem sizes using an empirical search strategy that varies the tile sizes. A similar search could be performed by varying the dividers in place of the tile sizes.

Another main difference compared to standard tiling is the required complex floating-point arithmetic for computing the bounds, which are algebraic expressions. Such computations may require high-order precision in some cases, whenever large integer values are involved -typically large loop bounds related to problem sizes -or very small real values. Such requirements are still essential, even if the final used result is the integer part (floor) of the computed real value. For this purpose, we implemented several related features in our software: multiple-precision complex floating-point computations, dichotomous search of the polynomial roots, adjustment of the solution using the ranking polynomials, verification code, and error handling functions computing trahrhe expressions. [START_REF] Clauss | Automatic Collapsing of Non-Rectangular Loops[END_REF] The Trahrhe software and the new -atiling flag of Pluto

The Trahrhe software has been written as a bash script with a few parsers in C, that makes intensive use of the computeralgebra system Maxima3 and of iscc4 , the interactive interface to the barvinok counting library and the isl integer set library. It is freely available 5 . It takes as input a parameterized iteration domain in the iscc syntax, and several possible flags as:

• -t𝑛: computes the trahrhe expressions required to tile the 𝑛 outermost loops; • -e: generates a C header file including all the definitions of functions required to apply algebraic tiling to a loop nest, as shown in Ehrhart and ranking polynomials are computed by invoking iscc. To solve symbolic polynomial equations, our software invokes the computer algebra system Maxima.

We are currently implementing algebraic tiling in the polyhedral compiler Pluto. By using flag -atiling, algebraic tiling is automatically applied to loops identified by the pragma #pragma scop in the input C source file. Additionally, the Trahrhe software is then invoked by Pluto to generate the required header files where the relevant trahrhe functions are generated in C.

Experiments

Some first experiments were conducted on 8 programs from the polybench benchmark suite v4.2 with data size EXTRALARGE.

All programs were compiled using gcc 10.3.0 with flags -O3 -march=native -fopenmp -lm. Standard (hyper-)rectangular tiling and OpenMP parallelization were automatically applied on the 8 programs thanks to Pluto, which performed simultaneously some loop transformations to promote data locality and vectorization. On each hardware platform, tile sizes providing the lowest execution times were selected by performing an exhaustive search among sizes that are powers of 2, from 1 to 2048. When very lower execution times were reached with algebraic tiling, standard tile sizes that are as close as possible to the size of the best performing algebraic tiles were also tested, although algebraic tile sizes are not generally constant.

Note that size 1 for a tiling depth means that the related loop was actually not tiled. The best performing tile sizes were found separately by running programs with schedule static and with schedule dynamic of OpenMP. With standard tiling, schedule dynamic may provide better performance, since it aims to improve load-balancing by dynamic distribution of iterations among the threads. It is particularly significant when the iteration domain is non-rectangular, as it is for kernel syr2k (see Figure 2). However, when the iteration domain is rectangular, schedule dynamic generally yields lower performance than schedule static, due to its underlying time overhead.

Regarding algebraic tiling, dividers providing the lowest execution times were also selected through an exhaustive search among powers of 2, from 1 to 2048, excepting for the dimension of the parallel loop, where the search was performed among dividers that are multiple of the number of threads (up to 10 times the number of threads). The same loop transformations as the ones generated by Pluto for standard tiling (skewing or interchange) were also applied with algebraic tiling. Note that programs with algebraic tiling were exclusively run with schedule static of OpenMP, since quasi-perfect load balancing is obviously expected. Hybrid tiling (standard + algebraic) was applied in some cases, in order to improve vectorization.

The programs were run on 2 × 32-core AMD Zen2 EPYC 7452 @ 2.35 GHz with Turbo-Boost and Hyperthreading deactivated, using 64 parallel threads. Runs were performed after having set the environment variable OMP_PROC_BIND=true to avoid thread migration and bind the threads to processor cores.

The best performing tile sizes and dividers providing the lowest execution times are reported in Table 2. Note that only multiples of 64 are used as the outermost slices' dividers. When tiling is hybrid, symbol «×» precedes the standard tile size that was selected, additionally to the algebraic dividers.

Speed-ups were computed by comparing execution times of programs optimized with algebraic tiling (that include the time spent in computing the algebraic tile bounds) against execution times of programs optimized with standard tiling, using the best performing dividers and tile sizes:

𝑠𝑝𝑒𝑒𝑑-𝑢𝑝 = 𝑡𝑖𝑚𝑒_𝑠𝑡𝑎𝑛𝑑𝑎𝑟𝑑_𝑡𝑖𝑙𝑖𝑛𝑔 -𝑡𝑖𝑚𝑒_𝑎𝑙𝑔𝑒𝑏𝑟𝑎𝑖𝑐_𝑡𝑖𝑙𝑖𝑛𝑔 𝑡𝑖𝑚𝑒_𝑠𝑡𝑎𝑛𝑑𝑎𝑟𝑑_𝑡𝑖𝑙𝑖𝑛𝑔
The resulting speed-ups against standard rectangular tiled programs run using schedule either static or dynamic are reported in Figures 7 and8. For Figure 7, automatic vectorization was deactivated using flag -fno-tree-vectorized of gcc, while it was activated for Figure 8. One can observe that vectorization has a significant influence in the lower speedup amount that is mostly obtained with algebraic tiling, even when hybrid tiling is applied. It shows that current microprocessor SIMD units and compilers seem not well suited for inner parallel loops of varying sizes, although algebraic tiling mostly outperforms standard tiling, particularly when iteration domains are not rectangular.

Related work

Tiling has been a very productive research topic for nearly thirty years [START_REF] Dongarra | Automatic Blocking of Nested Loops[END_REF][START_REF] Irigoin | Supernode Partitioning[END_REF][START_REF] Lam | A Data Locality Optimizing Algorithm[END_REF][START_REF] Ramanujam | Tiling Multidimensional Iteration Spaces for Nonshared Memory Machines[END_REF][START_REF] Wolfe | Iteration Space Tiling for Memory Hierarchies[END_REF]. However, tiling techniques that have been proposed are all based on partitioning the iteration domain into tiles, that may be (hyper-)rectangular or nonrectangular (hexagonal, trapezoidal, etc.), but always defined by their edge directions and sizes. Thus, their bounds are always expressed as (min or max of) affine functions, or polynomial functions when parameterized. To our knowledge, tiles delimited by algebraic bounds were never proposed.

There have been many proposals regarding stencil computations, which are an important class of programs that occurs in a variety of scientific applications, and for which tiling is particularly difficult. Stencil computations mostly require skewing parallelizing transformations [START_REF] Andonov | Optimal Semioblique Tiling[END_REF][START_REF] Lam | A Data Locality Optimizing Algorithm[END_REF][START_REF] Wonnacott | Achieving Scalable Locality with Time Skewing[END_REF] that have to be handled specifically [START_REF] Bondhugula | Diamond Tiling: Tiling Techniques to Maximize Parallelism for Stencil Computations[END_REF][START_REF] Grosser | Hybrid Hexagonal/Classical Tiling for GPUs[END_REF][START_REF] Henretty | A Stencil Compiler for Shortvector SIMD Architectures[END_REF][START_REF] Krishnamoorthy | Effective Automatic Parallelization of Stencil Computations[END_REF].

In [START_REF] Hartono | Parametric Multi-level Tiling of Imperfectly Nested Loops[END_REF], Hartono et al. propose PrimeTile, a parametric multi-level tiler for imperfect loop nests. Their tiler is "multilevel" as outer tiles may also be tiled at deeper levels. The authors use this approach to distinguish full tiles from partial tiles, and thus to avoid large partial tiles by deeper tiling of areas along borders of the iteration domain. Note that such strategy may be purposeless with algebraic tiling, since (full) algebraic tiles are built by considering the volume of the whole target iteration domain from the beginning.

Jiménez et al. [START_REF] Jiménez | Register Tiling in Nonrectangular Iteration Spaces[END_REF] propose a code generation technique for register tiling of non-rectangular iteration spaces. They generate code that traverses the bounding box of the tile iteration space to enable parameterized tile sizes. They apply index-set splitting to tiled code to traverse parts of the tile space that include only full tiles. Their approach involves less overhead in the loop nest that visits the full tiles, at the price of significant code expansion.

In their work, Sato et al. [START_REF] Sato | An Autotuning Framework for Scalable Execution of Tiled Code via Iterative Polyhedral Compilation[END_REF] address explicitly load balancing through an analytic model. They model the load balancing issue analytically and combine it with empirically autotuning the loop tile size for many-core CPUs, through an iterative compilation framework. However, they only address standard rectangular tiling and obtain speed-ups against a baseline which is all tile sizes equal to 32, while in our approach, we obtain significant speed-ups against the best performing rectangular tile sizes.

Sakellariou in [START_REF] Sakellariou | A Compile-Time Partitioning Strategy for Non-Rectangular Loop Nests[END_REF] proposes a compile-time scheme for partitioning non-rectangular loop nests, where the minimization of load imbalance is based on symbolic cost estimates. In [START_REF] Kejariwal | A Geometric Approach for Partitioning Ndimensional Non-rectangular Iteration Spaces[END_REF], Kejariwal et al. present a geometric approach for partitioning N-dimensional non-rectangular iteration spaces. They partition an iteration space along the axis corresponding to the outermost loop to achieve a near-optimal partition. Kafri and Abu Sbeih in [START_REF] Kafri | Simple Near Optimal Partitioning Approach to Perfect Triangular Iteration Space[END_REF] focus on static decomposition of perfect triangular iteration spaces to achieve load balancing, by partitioning a triangular iteration space of a loop nest along the outermost loop index.

Conclusion

When compared to standard (hyper-)rectangular tiling, algebraic tiling provides significant performance improvements, despite the relatively complex algebraic expressions that have to be evaluated at runtime to determine the tile bounds. It clearly confirms that load balancing among computing hardware units is crucial for performance. Moreover, with this technique, dynamic scheduling becomes purposeless for the handled loops. Algebraic tiling may be further extended and adapted, particularly for heterogeneous computing platforms by associating performance ratios to the target computing units. More generally, trahrhe expressions open a new field for loop optimizations based on iteration counts. Besides tiling, we plan to investigate other promising directions such as algebraic scheduling.

Figure 1 . syr2k loop kernel i j thread 0 thread 1 thread 2 thread 3 thread 4 Figure 2 .

 142 Figure 1. syr2k loop kernel

Figure 3 .

 3 Figure 3. Algebraic tiling applied on syr2k iteration domain

Figure 4 .

 4 Figure 4. Algebraic tiling on syr2k loop kernel

Figure 5 .Figure 6 .

 56 Figure 5. Why exact volumes cannot be reached

Figure 4 .

 4 It may be included by any C/C++ code for algebraic tiling optimization. Complex floating-point operations are performed with C type long double. Function definitions include verification instructions that stop the program whenever incoherent values arise in the evaluations of trahrhe expressions; arithmetic precision. This program compares every tuple of original iterators to the corresponding evaluation of the trahrhe expressions. It may be used as a code pattern to algebraic-tile-transform a given target code that includes an equivalent loop nest;

Figure 7 .

 7 Figure 7. Speed-ups resulting from comparing algebraic tiling vs standard tiling with non-vectorized codes (64 threads)

Table 1 .

 1 Generated tile bounds and resulting tile volumes for kernel syr2k with dividers 24 and 64 (sample) Target: 30,025K 𝑙𝑏 𝑖 → 𝑢𝑏 𝑖 𝑙𝑏 𝑗 → 𝑢𝑏 𝑗

			𝑙𝑏 𝑖 → 𝑢𝑏 𝑖	𝑙𝑏 𝑗 → 𝑢𝑏 𝑗	𝑙𝑏 𝑖 → 𝑢𝑏 𝑖	𝑙𝑏 𝑗 → 𝑢𝑏 𝑗
		Target: 467K		Target: 468K		Target: 466K
	0→243 (29,890K)	...	692→733 (29,967K) 0→10 (462K)	979→1008 (29,835K)	0→14 (450K)
		1→2 (484K)		11→21 (462K)		15→30 (480K)
		...		...		...
		63→65 (537K)		434→444 (462K)		621→636 (480K)
		...		...		...
		191→200 (440K)		691→702 (384K)		963→977 (450K)
		...		...		978→1008 (31K)
		Target: 470K		Target: 471K		Target: 480K
	244→345 (30,141K) 0→3 (408K)	734→773 (30,180K) 0→10 (440K)	1009→1038 (30,735K) 0→15 (480K)
		...		...		...
		277→284 (496K)		707→718 (480K)		992→1007 (480K)
		...		...		1008→1038 (31K)
		Target: 463K		Target: 471K		Target: 460K
	346→422 (29,645K) 0→5 (462K)	774→811 (30,153K) 0→11 (456K)	1039→1066 (29,498K) 0→15 (448K)
		...		...		...
		362→369 (432K)		743→755 (494K)		987→1003 (476K)
		370→379 (440K)		756→767 (456K)		1004→1019 (448K)
		380→392 (403K)		768→780 (416K)		1020→1036 (476K)
		393→422 (30K)		781→811 (31K)		1037→1066 (30K)
		Target: 470K		Target: 467K		Target: 473K
	423→488 (30,129K) 0→6 (462K)	812→847 (29,898K) 0→11 (432K)	1067→1094 (30,282K) 0→15 (448K)
		...		...		...
		436→445 (440K)		791→803 (468K)		1047→1063 (476K)
		...		...		...
		Target: 469K		Target: 473K		Target: 467K
	489→546 (30,073K) 0→7 (464K)	848→882 (30,310K) 0→12 (455K)	1095→1121 (29,943K) 0→16 (459K)
		...		...		...
		494→503 (440K)		825→837 (455K)		1039→1056 (486K)
		...		...		...
		Target: 465K		Target: 464K		Target: 461K
	547→598 (29,822K) 0→7 (416K)	883→915 (29,700K) 0→13 (462K)	1122→1147 (29,523K) 0→16 (442K)
		...		...		...
		537→545 (468K)		871→885 (465K)		1082→1099 (468K)
		...		...		...
		Target: 467K		Target: 466K		Target: 471K
	599→646 (29,928K) 0→8 (432K)	916→947 (29,840K) 0→13 (448K)	1148→1173 (30,199K) 0→17 (468K)
		...		...		...
		584→593 (480K)		830→844 (480K)		1125→1142 (468K)
		...		...		1143→1173 (31K)
		Target: 471K		Target: 466K		Target: 482K
	647→691 (30,150K) 0→9 (450K)	948→978 (29,884K) 0→14 (465K)	1174→1199 (30,875K) 0→17 (468K)
		10→19 (450K)		15→29 (465K)		18→36 (494K)
		20→30 (495K)		30→44 (465K)		37→54 (468K)
		31→40 (450K)		45→59 (465K)		55→73 (494K)
		41→51 (495K)		60→74 (465K)		74→91 (468K)
		...		...		92→110 (494K)
		638→648 (484K)		948→978 (31K)		...
		...

Table 2 .

 2 Best performing tile sizes and trahrhe dividers with size EXTRALARGE_DATASET and 64 threads × 2048 × 8 8 × 64 × 128 128, 1, 2048 16 × 128 × 16

	Program	without vectorization			with vectorization
		static	dynamic	algebraic	static	dynamic	algebraic
	correlation 4 16 × 128 × 8	128, 1, 2048
	covariance 4 × 32 × 64 8 × 2048 × 1	64, 1, 32	16 × 128 × 128 16 × 128 × 128 64, 1024 × 128
	gemm	32 × 32 × 32 32 × 16 × 256 64, 128, 1	32 × 64 × 64	32 × 64 × 32	64 × 1 × 256
	gemver	5 × 2048 × 1	16 × 1 × 1	384, 1024	4 × 2048 × 1	16 × 1024 × 1	128, 128, 64
	gesummv	22 × 2 × 1	44 × 1 × 1	512, 64, 1	5 × 8 × 1	43 × 1024 × 1	192, 1024, 1
	syr2k	16 × 1 × 1	4 × 8 × 1	128, 64, 1	32 × 64 × 1	8 × 64 × 1	128, 32, 16
	syrk	16 × 16 × 1	4 × 16 × 1	64, 4, 512	32 × 128 × 1	8 × 128 × 1	64, 128 × 16
	trmm	14 × 64 × 1 8 × 2048 × 1	64, 16, 2	41 × 32 × 16	41 × 16 × 16	128, 256 × 64

Figure 8. Speed-ups resulting from comparing algebraic tiling vs standard tiling with vectorized codes (64 threads)

http://maxima.sourceforge.net

https://repo.or.cz/barvinok.git

https://webpages.gitlabpages.inria.fr/trahrhe


Damien Jaime
email: damien.jaime@lip6.fr

Joyce El Haddad
email: joyce.elhaddad@lamsade.dauphine.fr

Pascal Poizat
email: pascal.poizat@lip6.fr

A preliminary study of rhythm and speed in the maven ecosystem

Keywords: software evolution, software ecosystem, dependency management, releases, rhythm, speed, Maven BENEVOL'22:

The analysis and the evolution of dependencies are central issues in software development. They can be studied from an ego-centric perspective, i.e., focusing on a given software artifact and its dependencies. But the scope can also be enlarged to take into account the whole context in which these artifacts and their dependencies evolve, i.e., ecosystems. Our general objective is to study how the knowledge about these ecosystems can guide developers in evolving software artifacts and their dependencies. In this paper we perform a preliminary study of the Maven ecosystem with reference to two metrics related to software artifact releases, namely rhythm and speed.

Introduction

Using dependencies to third parties is popular because it reduces development time and improves software quality [START_REF] Abdalkareem | Why do developers use trivial packages? An empirical case study on npm[END_REF]. Yet, this comes at a cost as these dependencies increase the architectural complexity and maintenance costs.

Raemaekers et al. [START_REF] Raemaekers | Measuring software library stability through historical version analysis[END_REF] studied in 2012 third-party stability through version analysis, and concluded that depending on old third-party versions is a "maintenance debt". Updating dependencies may however lead to a dilemma, since it may cause breaking changes [START_REF] Raemaekers | Semantic versioning and impact of breaking changes in the maven repository[END_REF]. A trade-off has then to be found in dependency management between keeping dependencies as is, and, e.g., missing security patches, or updating them, and possibly being hit by new bugs.

The study of dependencies may be achieved with a single system (ego-centric) perspective, or with an ecosystem one. The former may increase the quality of the system of interest, through the analysis of dependencies within the system (endo-dependencies) or with third parties (exodependencies) [START_REF] Musco | A generative model of software dependency graphs to better understand software evolution[END_REF]. The latter paves the way to the empirical study of dependencies and, possibly, allows one to develop concepts and development practices based on ecosystem-wide knowledge.

In 2015, Cox et al. proposed a metric to measure the dependency "freshness" and applied it, for a single third-party, on a whole ecosystem [START_REF] Cox | Measuring dependency freshness in software systems[END_REF]. Thereafter, González-Barahona et al. proposed the "technical lag" concept, focusing on the relation to vulnerabilities and bugs [START_REF] González-Barahona | Technical lag in software compilations: Measuring how outdated a software deployment is[END_REF][START_REF] González-Barahona | Characterizing outdateness with technical lag: An exploratory study[END_REF]. Since 2018, several authors have carried out studies concerning technical lag in ecosystems. This has led to analyses for Linux distributions [START_REF] Legay | On package freshness in linux distributions[END_REF][START_REF] Legay | A quantitative assessment of package freshness in linux distributions[END_REF], Docker Hub [START_REF] Zerouali | A multi-dimensional analysis of technical lag in debian-based docker images[END_REF], package managers [START_REF] Stringer | Technical lag of dependencies in major package managers[END_REF] and more specifically for npm [START_REF] Decan | On the evolution of technical lag in the npm package dependency network[END_REF][START_REF] Zerouali | An empirical analysis of technical lag in npm package dependencies[END_REF][START_REF] Zerouali | A formal framework for measuring technical lag in component repositories -and its application to npm[END_REF][START_REF] Cogo | An empirical study of dependency downgrades in the npm ecosystem[END_REF]. In 2021, Soto-Valero et al. [START_REF] Soto-Valero | A comprehensive study of bloated dependencies in the maven ecosystem[END_REF], presented the notion of "bloated dependencies", and analyzed their presence in the Maven ecosystem.

In this paper, alike the above-mentioned works, we take an ecosystem perspective on thirdparty and dependency study. Working in the context of industrial applications with both Java and Javascript dependencies, we decide to focus, first, on the Maven ecosystem, which has been less studied than the npm one. While a snapshot of Maven Central is available for study [START_REF] Benelallam | The maven dependency graph: A temporal graph-based representation of maven central[END_REF], dated 2018, Ochoa et al. have shown how knowledge and truths could change over a few years [START_REF] Ochoa | Breaking bad? Semantic versioning and impact of breaking changes in maven central: An external and differentiated replication study[END_REF]. As a consequence, we propose first a study of Maven in 2022, differences with 2018 being clear in Table 1.

We want to define metrics according to the evolution in time. These metrics are part of a larger study related to the quality of dependencies in software ecosystems. The goal is to make software safer and more maintainable by improving its structural quality. This paper introduces two new metrics related to artifacts and releases in ecosystems, namely rhythm and speed. We will therefore introduce these terms first and then review some possible applications of them. We believe that these metrics will help to solve the trade-off between updating or not a dependency [START_REF] Meil | Log4j and the thankless high-risk task of managing software component upgrades[END_REF] by tracking the dynamics of a package or ecosystem over time (e.g. accelerating releases may indicate a need to update the dependency). The research questions that these metrics enable us to study, and the related analyses, are presented in Sections 4 and 5. Experimental data and our Python Jupyter notebooks are available at [START_REF] Jaime | Data and notebooks for the "A Preliminary Study of Rhythm and Speed in the Maven Ecosystem[END_REF].

Metrics for the Ecosystem Analysis

A software artifact release (or release for short) in the Maven ecosystem can be identified by a triple 𝑥 : 𝑦 : 𝑧, where 𝑥 is the group id, 𝑦 is the artifact id, and 𝑧 is the version. An artifact can indeed be denoted using only the first two parts of a release, i.e., 𝑥 : 𝑦. Given a release 𝑟, 𝑑𝑎𝑡𝑒(𝑟) denotes its release date. Given two subsequent releases 𝑟 𝑖-1 = 𝑥 : 𝑦 : 𝑣 𝑖-1 and 𝑟 𝑖 = 𝑥 : 𝑦 : 𝑣 𝑖 of an artifact, 𝑡𝑖𝑚𝑒(𝑟 𝑖 ) denotes the number of days between the two releases, i.e., 𝑑𝑎𝑡𝑒(𝑟 𝑖 )-𝑑𝑎𝑡𝑒(𝑟 𝑖-1 ). When 𝑟 𝑖 is the first release of an artifact, 𝑡𝑖𝑚𝑒(𝑟 𝑖 ) is 0. We use [𝑑 1 -𝑑 2 ] to denote the time interval between 𝑑 1 and 𝑑 2 . We may use [0-𝑑 2 ] to denote time up to 𝑑 2 and [0-∞] to denote the whole time. Given an artifact 𝑎 and a time interval

[𝑑 1 -𝑑 2 ], 𝑟𝑒𝑙𝑒𝑎𝑠𝑒𝑠 [𝑑 1 -𝑑 2 ] (𝑎) denotes the time-ordered sequence ⟨𝑟 1 , . . . , 𝑟 𝑖 , . . . , 𝑟 𝑛 ⟩ of releases of 𝑎 in [𝑑 1 -𝑑 2 ]
. Finally, given an ecosystem 𝜖, 𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖) denotes the set of artifacts in 𝜖. We may now introduce our two metrics: rhythm and speed.

Definition 1 (Rhythm). The release rhythm (or rhythm for short) of an artifact 𝑎 over a time interval

[𝑑 1 -𝑑 2 ], denoted by ℛ [𝑑 1 -𝑑 2 ] (𝑎)
, is given by the formula: (𝑎), is given by the formula:

ℛ [𝑑 1 -𝑑 2 ] (𝑎) = ⟨𝑡𝑖𝑚𝑒(𝑟
𝑣 𝑑 1 -𝑑 2 (𝑎) = |𝑟𝑒𝑙𝑒𝑎𝑠𝑒𝑠 [𝑑 1 -𝑑 2 ] (𝑎)| (𝑑 2 -𝑑 1 ) (2)
where |𝑠| denotes the size of sequence 𝑠. The speed is measured in releases per day (𝑟/𝑑).

Speed can be lifted up to ecosystems.

Definition 3 (Ecosystem Speed). The average speed (or speed for short) of an ecosystem 𝜖 over a time interval

[𝑑 1 -𝑑 2 ], denoted by 𝑣 [𝑑 1 -𝑑 2 ] (𝜖)
, is given by the formula:

𝑣 [𝑑 1 -𝑑 2 ] (𝜖) = ∑ 𝑎∈𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖) |𝑟𝑒𝑙𝑒𝑎𝑠𝑒𝑠 𝑑 1 -𝑑 2 (𝑎)| (𝑑 2 -𝑑 1 ) × |𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖)| (3)
One can note (i) that the speed of an artifact 𝑎 is the speed of an ecosystem that is composed only of 𝑎 and (ii) that we also have:

𝑣 [𝑑 1 -𝑑 2 ] (𝜖) = ∑ 𝑎∈𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖) 𝑣 [𝑑 1 -𝑑 2 ] (𝑎) |𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖)| (4)

Presentation of the Datasets

In this section, we present the two datasets we use to perform the analysis of the Maven ecosystem (see Tab. 1).

Maven Dependency Graph (MDG)

This dataset has been proposed by Benelallam et al. in [START_REF] Benelallam | The maven dependency graph: A temporal graph-based representation of maven central[END_REF]. It is a graph-based representation of a snapshot of the Maven Central repository at September 6 th , 2018. Its vertices include for each release its id, the packaging type (jar, war, ear), and the release date. Two kinds of edges are used: one for the precedence relationship between releases (e.g., 𝑥 : 𝑦 : 1.0.1 is next to 𝑥 : 𝑦 : 1.0.0) and one for the dependency relationship between artifact releases given a context (e.g., 𝑥 1 : 𝑦 1 : 𝑧 1 depends on 𝑥 2 : 𝑦 2 : 𝑧 2 for the compiling context). As stated in [START_REF] Benelallam | The maven dependency graph: A temporal graph-based representation of maven central[END_REF] the data collection has been achieved by "retrieving pom.xml files (at least one per artifact) from Maven Central and parsing them to retrieve metadata".

Maven Central Index extraction (MCI)

The MCI is available at Maven Central 1 . It contains all artifacts, updated weekly. We will use the version at June 10 th , 2022. From the MCI we obtain a Lucene database using scripts available at Maven Central and then we extract CSV files from it. These CSV files are available at [START_REF] Jaime | Data and notebooks for the "A Preliminary Study of Rhythm and Speed in the Maven Ecosystem[END_REF].

As for the MDG, this dataset includes release ids and release dates. However, it does not (yet) include dependency information.

It is interesting to see how the number of releases has grown in almost 4 years between the MDG snapshot and the more recent MCI one. There is an important increase, which makes it interesting to have recent data. We study the increase of ecosystem speed in Section 4.1, but it would be interesting to study this behavior in more detail, especially by comparing it with npm for example.

Dependency-independent Analysis

In this section, we study the Maven ecosystem independently of the dependencies between releases. We wish to address two research questions :

• RQ1: How do rhythm and speed evolve over time in the Maven ecosystem? • RQ2: Can important events in the life of an artifact be observed through the prism of rhythm?

Analysis of the Maven ecosystem evolution

Table 1 shows an important increase in the number of releases between 2018 and 2022. This raises the question of the evolution of rhythm between these two years.

To answer this question, we took the MCI dataset and computed the distribution of times between releases (𝑑𝑡𝑏𝑟) respectively for all years up to 2018 and for all years up to 2022. For this we used the the following function:

𝑑𝑡𝑏𝑟(𝑖) = { {𝑡𝑖𝑚𝑒(𝑟) | 𝑎 ∈ 𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖) ∧ 𝑟 ∈ 𝑟𝑒𝑙𝑒𝑎𝑠𝑒𝑠 𝑖 (𝑎) ∧ 𝑡𝑖𝑚𝑒(𝑟) ̸ = 0} } (5)
where 𝑖 is a date interval. This regroups all non-null times between releases for releases with a date in the interval. Here, we use the [0 -31/12/2018] and [0 -31/12/2022] intervals and the result is shown in Table 2.

We can consider that, regarding the date of the corresponding release, a release time around Q1 is fast, around Q3 is slow, and around the median is medium compared to the ecosystem. We can see that, within 4 years, the time between releases decreased which suggests an overall acceleration.

Let us now study the general evolution of both rhythm and speed for all years between 2005 and 2022. We followed the same approach as before for aggregating the rhythm, this time focusing only on years of interest (instead of up-to), so given some year 𝑦 we use Function ( 5) with i = [01/01/y-31/12/y] (instead of [0-31/12/y]).

The result is given in Figure 1. Each box-plot corresponds to the aggregated rhythm for the year of interest, while speed is computed using Formula 2 with yearly intervals. Note that 2022 is taken into account only until June, 10 tℎ .

We can observe an accelerating trend except for the period between 2009 and 2012 when there is a deceleration. We can also distinguish different parts. There is an acceleration between 2005 and 2009, a deceleration between 2009 and 2012, a new acceleration between 2012 and 2019, and finally a more important acceleration since 2019. It would be interesting to study the possible causes of these variations (new Java or Maven release? COVID-19 emergence?). While this preliminary study does not aim at finding causality relationships, these questions could be addressed in a later study.

If we observe now the evolution of rhythm at the artifact level, we get the result presented in Figure 2a. In this figure, we have computed the average time to release the 𝑥 𝑡ℎ version of artifacts, i.e., using the following formula:

𝑓 6 (𝑥) = ∑ 𝑟∈𝑥 𝑡ℎ𝑠 (𝜖,𝑥) 𝑡𝑖𝑚𝑒(𝑟) |𝑥 𝑡ℎ𝑠 (𝜖, 𝑥)| ( 6
)
where 𝑥 is the rank of the version (e.g., 4 means it is the fourth) and

𝑥 𝑡ℎ𝑠 (𝜖, 𝑥) = ⋃ 𝑎∈𝑎𝑟𝑡𝑖𝑓 𝑎𝑐𝑡𝑠(𝜖) 𝑟𝑒𝑙𝑒𝑎𝑠𝑒𝑠 [0-∞] (𝑎)[𝑥].
From the figure we can see, e.g., that on average the second release of an artifact takes about 85 days to release, while its 50 tℎ release takes only 20 days.

It seems that the more mature an artifact is, the more its rhythm tends to accelerate. But the reality is more complex than that. It is important to distinguish several groups in the ecosystem and to compare them independently. Indeed the artifacts with many releases have a different rhythm than those with few releases, and the proportion of artifacts with 10 releases is not at all the same as those with 400 releases.

Figure 2b, focuses only on the (931) artifacts with more than 400 releases. We can see that they have immediately a very low time between releases. It would be interesting to define groups of artifacts in the ecosystem to be able to associate a given artifact to the right group and compare it (based on rhythm and speed) more precisely.

Analysis of the evolution of an artifact: log4j-core

We will now see if rhythm allows us to observe changes in an artifact's life. Each artifact has a different rhythm, so to know if a release was faster than the others we have to perform an analysis at the artifact level.

Take the example of Apache Log4j2 (more precisely, artifact org.apache.logging.log4j:log4jcore), the most used logger library in the Java world. On December 9 tℎ , 2021, a level 10 security fault (CVE-2021-44228) was discovered which impacted a large number of artifacts depending on it. The rhythm distribution of Log4j2 before the issue discovery (last version on March 7 tℎ , 2021) was medians: 63.0, Q1: 26.5, Q3: 124.0, which position it as a slow library compared to the Maven ecosystem (see Sect. 4.1). At the time of writing this paper, the last version is dated February 23 r𝑑 , 2022, and the rhythm distribution is medians: 48.5, Q1: 19.0, Q3: 114.2. We can see that the security issue has had an impact.

The blue plot in Figure 3 shows the rhythm of the project since 2020. If we focus on [07/03/2021-23/02/2022], we get [277, 4, 1, 4, 2, 1, 6, 1, 1, 56]. We observe a sharp drop on the same day as the issue was discovered. Speed had already suffered decelerations in the past, but this one is distinguished by the fact that the drop is vertical and is followed by a plateau. We can therefore say that there is a strong reactivity, a period of stabilization and then a beginning of a return to normal.

Rhythm allowed us to observe a change in the life of Log4j2, through the imprint it had. Yet, it would be interesting to analyze the whole Maven ecosystem to search for such patterns and see, e.g., if there is a correlation between CVEs and rhythm changes.

Dependency-aware analysis

In this section, we bring dependencies into the big picture. We wish to address one research question:

• RQ3: Does the rhythm of an artifact have an impact on the rhythm of artifacts that depend on it?

Since we have seen in Section 4.2 how the rhythm of Log4j2 evolved concerning CVE-2021-44228, we will address this question on this artifact.

To recover dependency information, we searched, in the Maven Repository2 , for the uses of the last version of log4j-core before the discovery of the issue, namely version 2.14.1 released on March 12 tℎ , 2021. Not having found either an API or export mechanisms for this, we created a script for browsing and parsing the Web pages to find the dependent artifacts. This gave us 500 different releases, for 149 distinct artifacts. The data is, again, available at [START_REF] Jaime | Data and notebooks for the "A Preliminary Study of Rhythm and Speed in the Maven Ecosystem[END_REF].

In Figure 3, we have seen that the blue plot is the Log4j2 rhythm from 2020 to 2022. The red plot shows the average rhythm of the Log4j2 dependents on the same period. Finally, the green plot shows the total number of releases of these dependents. All three plots are computed monthly.

We notice an increase in the number of releases of dependents in the month of the CVE-2021-44228 discovery. In December 2021, 224 artifacts have made a release, the release time median is 5 days, with a minimum value at 1 day and a maximum value at 594 days. Yet, the average rhythm of these dependents does not show a significant variation. Possibly, artifacts that had not been released for a long time and those that had just been released compensate for each other. Thus, we can conclude that the average rhythm is not a tool precise enough to study the effect of an artifact on its dependents.

Limitations & Threats to Validity

First of all, our analyses are based on extractions from Maven Central. Results may therefore not apply to all of the Maven or Java ecosystem. Broadly, we need to study if our concepts are amenable to software ecosystems in general.

We made artifact-centered experiments on Log4j2, with time intervals centered around the CVE-2021-44228 discovery. This required to use the MCI dataset and own-made dependency extraction. To generalize these experiments, we would have either to use the MDG dataset and limit to 2018, update MDG, or add dependencies in our MCI extraction.

In experiments based on the time between release values or rhythm, we do not take into account releases with a time value of 0 (any 𝑟 such that 𝑡𝑖𝑚𝑒(𝑟) = 0). This occurs for the first releases of artifacts or if an artifact has several releases on the same day. We should study the impact of this choice.

Finally, a need for categorization arises when comparing a specific artifact to the ecosystem. This is due to possible heterogeneity between artifacts (we have seen an example of this in Fig. 2a and2b). It would be interesting to use labeling information (e.g., the Maven categories of artifacts) or to use clustering to be able to relate a given artifact to close ones.

Conclusion & Perspectives

In this paper, we have studied the rhythm and speed of Maven Central until June 10 tℎ , 2022. We have explored possible applications of these metrics. It will be necessary to deepen their applications, but we believe that they can be useful in the management and analysis of dependencies. We have observed a general trend toward a decrease in time between releases, and accordingly an increase in speed, even as the number of artifacts in the ecosystem increases. We were also able to observe a correlation between the number of releases of an artifact and its rhythm. Finally, we have seen that an event in the life of an artifact could have a visible imprint on its rhythm and that we could observe a change of behavior in the artifacts depending on it, even if not on their rhythm.

This preliminary study paves the way to further analyses such as clustering based on rhythm and speed. We have first results showing a possible relation between the number of dependents of an artifact and its speed. We have not been able to talk about it in this paper, due to room running short, but our experimentation on the MDG dataset is available at [START_REF] Jaime | Data and notebooks for the "A Preliminary Study of Rhythm and Speed in the Maven Ecosystem[END_REF].

A first perspective is to study how the knowledge of the rhythm and speed of one's system dependencies, including relative to the ecosystem, may help in managing these dependencies, e.g., prioritizing upgrades or selecting alternatives. Another perspective is related to the comparison of the Maven and npm ecosystems, first from the rhythm and speed perspective. We also want to extend the MDG dataset dependency graph model (or define our own) in three directions. First, we want to support metrics related to aging dependencies, such as freshness [START_REF] Cox | Measuring dependency freshness in software systems[END_REF], and study how they compose transitively (along dependencies of dependencies). The second extension is concerned with the expressiveness of dependency requirements. Maven for example supports version ranges. Multiple dependency versions are also possible up to transitivity and Maven has its own resolution algorithm for this. A consequence is that the effective dependency used at some date 𝑑 1 may be different from the one used at 𝑑 2 , even if the POM file has not changed in between. A dependency graph should support expressive dependency constraints rather than hardcoded dependency relations as in the MDG dataset, and realization of constraints (fixing the version to be used) only when needed. The last extension is to include other parts of the ecosystem in its dependency graph. Examples include the compiler, virtual machine, and build system versions, which have an impact, e.g., on the success of a system build.

Figure 1 :

 1 Figure 1: Evolution of aggregated rhythm and speed in the Maven ecosystem between 2005 and 2022.

  (a) all artifacts (b) only artifacts with more than 400 releases

Figure 2 :

 2 Figure 2: Average time to release 𝑖 𝑡ℎ versions.

Figure 3 :

 3 Figure 3: Evolution of the Log4j2 rhythm and relation to its dependents.

Definition 2 (Speed). The

  1 ), . . . , 𝑡𝑖𝑚𝑒(𝑟 𝑖 ), . . . , 𝑡𝑖𝑚𝑒(𝑟 𝑛 )⟩ (1) with ⟨𝑟 1 , . . . , 𝑟 𝑖 , . . . , 𝑟 𝑛 ⟩ = 𝑟𝑒𝑙𝑒𝑎𝑠𝑒𝑠 [𝑑 1 -𝑑 2 ] (𝑎). average speed (or speed for short) of an artifact 𝑎 over a time interval [𝑑 1 -𝑑 2 ], denoted by 𝑣 [𝑑 1 -𝑑 2 ]

Table 1

 1 Descriptive statistics about the MDG (2018-09-06) and the MCI (2022-06-10) datasets.

	Property	MDG	MCI Increase
	Releases (unique 𝑥 : 𝑦 : 𝑧)	2,407,335 9,060,212	276%
	Artifacts (unique 𝑥 : 𝑦 in 𝑥 : 𝑦 : 𝑧)	223,478	477,707	113%
	Groups (unique 𝑥 in 𝑥 : 𝑦 : 𝑧)	35,699	66,139	85%
	Upgrades (release evolutions)	2,183,845 8,836,734	304%
	Dependency relationships	9,715,669	-	-

Table 2

 2 Comparison of time between releases in days up to 2018 and up to 2022.

		MCI up to 2018 MCI up to 2022
	Medians	14.0	9.0
	Means	43.9	35.6
	Q1	4.0	2.0
	Q3	43.0	32.0
	Minimums	1	1
	Maximums	101	77
	Outlier maximums	4341	5206

https://maven.apache.org/repository/central-index.html

https://mvnrepository.com/

Acknowledgments

This work is funded by an ANRT PhD grant, number 2021/0047.