Patrick Amestoy

Je Souhaite Remercier

Alexei Mikchevitch

Roméo, Matthieu Bastien

Dimitri

Keywords: sparse matrices, direct methods for linear systems, low-rank approximations, floatingpoint arithmetic, rounding error analysis, high-performance computing

Je souhaite tout d'abord remercier mes encadrants pour leur aide durant ces trois années.

Introduction

Many industrial and scientific simulations require to solve linear systems of the form Ax = b. Being able to perform such an algorithmic step with enough accuracy is crucial for of feasibility and performance of the simulation and may take a huge part of the computation time. It is also the case for the codes developed at EDF R&D, such as code_aster, code_carmel and code_saturne.

In order to solve a linear system in a robust and accurate way, one might need to use a direct method such as performing an LU factorization A = LU , before solving triangular systems Ly = b and U x = y. The multifrontal method and the supernodal method are both part of this category of algorithms, while handling sparse matrices, whose coefficients are mostly zero.

However, direct methods can be very costly in terms of memory consumption and computation time. As an addition to the use of massive parallelism, one may try to take advantage of the properties of data sparsity of the matrices manipulated in order to further reduce the complexity. Indeed, carefully chosen sub-matrices of A, L and U can often be approximated by matrices of low ranks. By doing this we obtain a Block Low-Rank compression (BLR) that can be exploited to reduce the computational complexity of direct solvers [START_REF] Amestoy | On the complexity of the Block Low-Rank multifrontal factorization[END_REF] and therefore reduce their time and memory consumption significantly. The BLR format is notably used in the MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF](Amestoy et al., , 2019a)), PaStiX [START_REF] Hénon | PaStiX: A high-performance parallel direct solver for sparse symmetric definite systems[END_REF][START_REF] Pichon | Sparse supernodal solver using Block Low-Rank compression: Design, performance and analysis[END_REF][START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization techniques[END_REF], and STRUMPACK [START_REF] Rouet | A distributed-memory package for dense hierarchically semi-separable matrix computations using randomization[END_REF][START_REF] Ghysels | An efficient multicore implementation of a novel HSS-structured multifrontal solver using randomized sampling[END_REF] solvers.

Another perspective to improve the performance of numerical methods is the use of mixed-precision algorithms. Low precision is used to improve the performance of the computations: speed, memory, and energy consumption. However, those low-precision formats are carefully combined with higher ones so that the result will have a high accuracy. This category of algorithms has recently known a renewed interest due to the growing availability of low-precision formats on modern hardware, such as fp16 and bfloat16.

The goal of this work is to develop new techniques that aim at further improving the performance gains of using BLR compression within a direct method, in particular based on mixed precision.

In chapter 2 we investigate the possibility of combining low-rank approximations with mixed-precision arithmetic. The main idea is that not all coefficients of a low-rank approximation are equally important. We try to store different components in different precision formats, the components with the lowest norms being kept and used in low CONTENTS precision. We justify this approach with an error analysis in the case where an arbitrary number of floating-point formats are used. As a result, we obtain a criterion determining what precision format should be used for each component. Our numerical experiments suggest that the proposed mixed-precision low-rank representation presents a high potential: a large fraction of both the entries needed to represent BLR matrices and the floating-point operations needed to compute their LU factorization could be switched to lower precisions.

In chapter 3 we adapt the approach described in chapter 2 for dense matrices to the factorization of sparse matrices based on a multifrontal method. Indeed, in the multifrontal method we process a sequence of dense frontal matrices on which the algorithms presented on chapter 2 for the dense LU factorization and its BLR variant can be applied. We propose two main contributions, while being more focused on the practical aspects of the algorithms. First, we describe how mixed-precision BLR (MPBLR) can be used as a storage format only, and could benefit from using custom floating-point formats. We implemented it in the multifrontal solver MUMPS. Our experiments show that the memory peak is significantly reduced at the price of a small time overhead, without compromising the accuracy. Second, we describe algorithms that perform most computations using low precision in the phases of factorization and triangular solution. We wrote a high-performance implementation of the solution phase within the MUMPS solver. The computation time is reduced, while preserving the previous memory gains, in the case there is only one right-hand side (RHS) to the linear system.

Our experiments revealed a weakness of the BLR triangular solution phase in case we have several RHS: in MUMPS, the time gains of both the classical and the mixedprecision BLR variants diminish when increasing the number of RHS. After analyzing the algorithm in case of multiple RHS, we conclude that the data locality is rather poor, which motivates us to rethink the communication and memory access patterns. In chapter 4 we propose a new algorithm in which the order of the block operations has been modified. We carry out a communication volume analysis that shows that the data movement is reduced and therefore the data locality is improved. We implement several variants of the algorithms, including one that is better adapted to the parallelism scheme of MUMPS. The latter obtains time gains compared to the standard algorithm, in case of multiple RHS.

In order to make this thesis self-contained, we also provide a chapter of general background (chapter 1). In particular we explore the techniques used in the multifrontal method, how to handle BLR compressions, the basics of rounding error analysis, as well as the state of the art of the mixed-precision algorithms that are somewhat related to this thesis.

Gaussian elimination

Gaussian elimination (Algorithm 1.1) is one of the most basic direct methods for solving a linear system Ax = b. Its main step consists in performing a series of updates on each coefficient: overall, we will subtract k<min(i,j) l ik u kj from each coefficient a ij . The factors L and U are progressively filled with their final values, one column at a time for L, and one row at a time for U . After the last step of the algorithm (should there be no division by 0), we obtain a decomposition A = LU of the original matrix, where L is a lower triangular matrix and U an upper triangular matrix. Then, we solve the two triangular systems Ly = b and U x = y using Algorithms 1.2 and 1.3 respectively. Algorithm 1.1 Gaussian elimination (right-looking)

1: L ← 0 2: U ← 0 3: for k=1 to n do 4:

for i=k to n do 5:

l i,k ← a i,k /a k,k 6: u k,i ← a k,i 7:
end for 8:

for i=k+1 to n do 9:

for j=k+1 to n do 10:

a i,j ← a i,j -l i,k u k,j 11:

end for 12:

end for 13: end for Algorithm 1.2 Forward elimination: Solve Ly=b (left-looking)

1: y ← b 2: for i=1 to n do 3:

y i ← (y i -i-1
j=1 l ij y j )/l ii 4: end for CHAPTER 1. BACKGROUND Algorithm 1.3 Backward substitution: Solve Ux=y (left-looking)

1: x ← y 2: for i=1 to n do 3:

x i ← (x i -n j=i+1 u ij y j )/u ii 4: end for Gaussian elimination is often implemented in-place in practice, and it will also be the case for the variants presented in this thesis. The computed columns of L are successively stored in the lower part of A, overwriting the corresponding coefficients of A that are no longer needed. Similarly, the computed rows of U are stored in the upper part of A. We note that, as Algorithm 1.1 is written, the diagonal coefficients of L have a value 1. Thus, it is possible not to store them explicitly, which leaves just enough space for storing the coefficients of U , including the diagonal. Also, the original matrix A is overwritten when performing the updates. In short, the matrix that is manipulated by actual implementations of the algorithm is, after step k, (L -I k ) + U + A sub , where A sub is obtained from the current matrix A by replacing all unused coefficients by 0, leaving only the k × k coefficients of the bottom-right block. I k is a matrix whose k first diagonal coefficients are 1, the others being 0.

In a less naive approach, one may want to add a pivoting step to Algorithm 1.1. After line 3, we may exchange row k with another row i p ≥ k, and column k with a column j p ≥ k. As a consequence, the result obtained at the end of the algorithm is a factorization P AQ = LU , and the right-hand side of the linear system will have to be adapted as well by performing pivoting operations on b and/or on x.

One of the reasons behind the need to add such a pivoting step is that we cannot perform the division by a k,k if its value is 0. Moreover, we try to avoid dividing by a small value of a k,k for numerical reasons when using an inexact floating-point arithmetic. As a result, a first strategy known as complete pivoting consists of choosing as pivot the largest eligible coefficient, which requires pivoting on both the rows and columns. However, complete pivoting is rather costly, and this is why we will rather consider a cheaper alternative: we only perform interchanges on the rows, taking the largest coefficient of the column as the pivot for instance. This pivoting strategy is known as partial pivoting, and in practice it is often sufficient to obtain numerical stability. As a result, we obtain a decomposition P A = LU of the original matrix.

This algorithm can be adapted to the symmetric case: given a symmetric matrix A, we compute a factorization A = LDL T . The number of coefficients to store and the number of operations to perform on them should both be divided by a factor 2. The pivoting strategies must be adapted in order to keep the matrix symmetric at each step: the pivots are typically chosen on the diagonal instead of the current column.

Adapting Gaussian elimination to sparse matrices

Other versions of the Gaussian elimination handle the case where we want to factorize a matrix A that is structurally sparse, with most of its coefficients being 0. However, with a naive choice of pivots, most of the coefficients of the factors are non-zero: a huge fill-in occurs, and the sparsity is lost. On the other hand, if the pivots are chosen so as to minimize this fill-in, then we may obtain factors L and U that are relatively sparse (although not quite as much as the original matrix). The number of operations performed is greatly reduced as well, because all operations performed on structural zeros are skipped.

The multifrontal method is one of these algorithms. It was first formally introduced in [START_REF] Duff | The multifrontal solution of indefinite sparse symmetric linear systems[END_REF]. It is composed of three main phases:

• During the analysis phase, we choose a main pivoting strategy, choosing pivots on the diagonal in a way that reduces the fill-in of the factorization. This is done by running a certain graph algorithm on the adjacency graph of the matrix. At the end of this phase, we obtain an elimination tree (see Figure 1.1), which is the structure that will be manipulated. It is a representation of the matrix A.

• Factorization phase: The elimination tree of A is modified, by performing the updates of its coefficients. After each of its coefficient has been fully summed, we obtain a tree representing the factors L and U.

• Solution phase: We solve the triangular systems Lx = b and Uy = x.

At the end of the analysis, we obtain an elimination tree of the sparse matrix A. Each of its nodes corresponds to a dense submatrix of A, and is called a front.

We perform a bottom-up traversal of this elimination tree. At each node we perform a partial LU factorization of the current front, as illustrated in Figure 1.2. As a result, we obtain new fully summed parts, that is to say coefficients on which all updates have been performed. They are replaced by L and U , which are submatrices of L and U: they respectively contain all the non-zero coefficients of a set of columns of L, and all the non-zero coefficients of a set of rows of U. The contribution block (CB) corresponds to a certain partial sum of updates, to be added with the right coefficients of the parent front during the assembly step.

Factorization of a frontal matrix

The partial LU factorization of a frontal matrix A 11 A 12 A 21 A 22 can be seen as performing the following steps:

• Compute an LU factorization L 11 U 11 of A 11

• Compute L 21 = A 21 U -1 11 (triangular solve)

• Compute U 12 = L -1 11 A 12 (triangular solve)

• Update the contribution block A 22 , which becomes the Schur complement of the original frontal matrix :

A 22 ← A 22 -L 21 U 12
In order to save some memory space, these four operations are once again performed in-place in practice. As a result, U 11 and L 11 are stored in the same memory space as A 11 (without storing the 1s of the diagonal of L 11 ), and A 21 and A 12 are overwritten by L 21 and U 12 .

We note that all these operations can be performed using well-parallelized linear algebra routines, such as BLAS-3 routines. Therefore, we can expect them to be performed very efficiently in practice.

Assembly

Before performing the factorization of a node, the contributions of all its descendants have to be added to its coefficients. Therefore, the contribution blocks S 1 and S 2 of the children nodes are added through an extend-add operation of the form A ← A ⊕ S 1 ⊕ S 2 : this is the assembly step.

However, the cost for storing the contribution blocks before the assembly can be quite high, and those temporary variables represent a huge part of the memory peak.

Triangular solution

The factorization of A produced a modified elimination tree illustrated in Figure 1.3, containing the factors L and U. We now solve the two triangular systems Lx = b and Uy = x. The nodes are visited from the bottom to the top of the tree for the former, and from the top to the bottom for the latter. On each front, we perform a series of operations on the right-hand side X 1 X 2 using the L factors L 1 L 2 :

• Triangular solve:

X 1 ← L -1 1 X 1 • Update: X 2 ← X 2 -L 2 X 2
Similarly, the operations performed during the backward substitution on the righthand side Y 1 Y 2 using the factors U 1 U 2 are the following :

• Update: Y 1 ← Y 1 -U 2 Y 2 • Triangular solve: Y 1 ← U -1 1 Y 1
Before treating a front, the right-hand side X has to be initialized from the children node, similarly to the assembly step of the factorization.

We may want to solve a linear system with several right-hand sides (RHS) instead of having only one. In this case, x, y, X and Y will not be a 1-dimensional vectors, but matrices with several columns instead. This modification is mostly transparent: for example, matrix-vector products are replaced with matrix-matrix products.

Parallelism

We distinguish two main sources of parallelism in the multifrontal method:

• Different branches of the elimination tree are independent and can be traversed concurrently. This is referred to as tree parallelism.

• The partial LU factorization of large enough fronts may also be performed on several processes/threads: this is node parallelism.

We note that tree parallelism induces an extra memory cost: if the factorizations of two nodes are performed in parallel, then the working space is duplicated. In particular, the contribution blocks are large-size temporary variables. Using tree parallelism will duplicate them, and they will end up taking a huge part of the memory consumption. On the other hand, node parallelism requires a huge amount of data transfers. Thus, node parallelism is mostly suited to the top of the tree, where there is not enough potential for tree parallelism (e.g., the root node is only eligible for node parallelism), while tree parallelism is best suited to the bottom of the tree.

Exploiting data sparsity

In many applications, the matrices manipulated are structurally sparse, with most of their coefficients being zero. As a result, the complexity of the problem can be greatly reduced. However, this aspect fails to fully take into account how sparse the problem really is, which part is significant, and which operations can be performed approximately for numerical reasons. In fact, one may often benefit from data sparsity, given the fact that the matrices may still be compressed much further for a given accuracy. This is also often the case for the frontal matrices taken from an elimination tree of the multifrontal method, especially when the matrix has been obtained from a finite element method. Their off-diagonal blocks correspond to interactions between two groups of variables that are weakly connected. They contain little information, and we could try to compress them.

Low-rank approximations

Definition

Let A be a matrix of size (m, n). A rank-r approximation of A is a matrix T of rank r situated within a certain distance τ of A. In other words, it means that there exist two matrices X ∈ R m×r and Y ∈ R n×r such that ∥A -XY T ∥ ≤ τ . We say that r is the numerical rank of A. If r is sufficiently small, we say that T = XY T is a low-rank approximation of A. Storing X and Y requires fewer coefficients than storing A: such a representation may be used as a matrix compression of A.

We might want to use a threshold τ that is relative to the norm of A instead of a constant. There are several possible choices for the norm ∥.∥. In this thesis, we will normally choose the Frobenius norm:

∀M, ∥M ∥ = ∥M ∥ F = m i=1 n j=1 m 2 ij .

Truncated SVD

One possible choice for computing a low-rank approximation of A is through a trun-

cated SVD. Let U ΣV = min(m,n) i=1 σ i U i V T i be the singular value decomposition of A (SVD).
Given a target accuracy τ , we choose r as the smallest integer such that T = r i=1

σ i U i V T i verifies ∥A -T ∥ ≤ τ .
We know from [START_REF] Eckard | The approximation of one matrix by another of lower rank[END_REF] that using a truncated SVD leads to the optimal low-rank approximation for the Frobenius norm: for a given rank, it leads to the most accurate low-rank approximation. Said otherwise, for a given required accuracy, this technique leads to an approximation of minimal rank. This property of optimality also stands when using the 2-norm instead of the Frobenius norm.

QR decomposition

However, computing a truncated SVD is rather costly in terms of time. In practice one may prefer to perform a non-optimal low-rank approximation, easier to compute. One of those possibilities is to compute a truncated QR factorization with column pivoting. Therefore, we obtain an approximation AP ≈ QR, where P is a permutation matrix, Q is a matrix of dimensions m × r with orthonormal columns, and R an upper triangular matrix of size r × n.

BLR matrices

Let A ∈ R n×n be a dense square matrix partitioned into q ×q blocks. A block low-rank (BLR) representation T of A is a block matrix of the form

T =        T 11 T 12 • • • T 1q T 21 • • • • • • . . . . . . . . . T q1 • • • • • • T qq        , (1.1)
where the blocks A ij of size b × b are approximated by matrices T ij satisfying

∥A ij -T ij ∥ ≤ εβ ij , (1.2)
where β ij > 0. Some of the blocks are left uncompressed, i.e., T ij = A ij . We will refer to them as being "full-rank" (FR). In particular, this is the case for all the diagonal blocks.

The other blocks are approximated as low-rank (LR): each T ij matrix, of rank r ij , is expressed as

T ij =    X ij Y T ij , i > j, Y ij X T ij , i < j, (1.3)
where X ij and Y ij are b × r ij matrices, and where X ij has orthonormal columns.

Even though, in general, each block can be of different dimensions, we assume for simplicity that they are all of the same dimensions b × b, and so n = qb.

Representation (1.3) guarantees the so-called outer orthonormality property [START_REF] Higham | Solving Block Low-Rank linear systems by LU factorization is numerically stable[END_REF][START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance, and scalability[END_REF], which is used in MUMPS. This property will mostly be reused in chapter 2. For simplifications, we will tend to ignore it in chapters 3 and 4, and consider that all low-rank blocks are under the form

T ij = X ij Y T ij .
Importantly, the β ij parameters in Equation 1.2 are used to distinguish two types of BLR compression, local and global, depending on whether block T ij approximates A ij with error ε relative to the local norm β ij = ∥A ij ∥ or relative to the global norm β ij = ∥A∥. In their error analysis of the BLR factorization, [START_REF] Higham | Solving Block Low-Rank linear systems by LU factorization is numerically stable[END_REF] show that global compression achieves a better tradeoff between compression and accuracy, and is therefore to be preferred. Throughout this thesis, we will thus use a global compression, with β ij = β = ∥A∥. On a BLR matrix, it leads to a global error bound that is proportional to ε:

∥A -T ∥ ≤ qε∥A∥.
(1.4)

Adapting the multifrontal method to BLR compression

We observe that, in the multifrontal method, large frontal matrices tend to have a huge potential for BLR compression. For example, Figure 1.4 was obtained from the root node of the elimination tree of a sparse matrix, and most of the blocks are highly compressible. We will now present how to take advantage of the BLR compression described previously in order to reduce the complexity of the multifrontal method. Such modifications were first described during the PhD thesis of Clément Weisbecker (see [START_REF] Weisbecker | Improving multifrontal solvers by means of algebraic Block Low-Rank representations[END_REF], and further improvements were made in the PhD thesis of Théo Mary (see [START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance, and scalability[END_REF]. In order to avoid confusion with the BLR variant of the multifrontal method that we present in this section, we will sometimes refer to the classical multifrontal method presented in section 1.1.2 as the full-rank algorithm (FR).

Factorization of a frontal matrix

Contrary to the full-rank case, where we could apply basic linear algebra operations on large parts of the matrix at once, we now have to follow the block structure of the BLR matrix. Therefore, we must adapt Algorithm 1.1 into a tile algorithm, which manipulates blocks of size b × b as basic elements. The main difference is that the update formula, formerly

a ij ← a ij -a ik a -1
kk a kj , is now generalized to a block:

A ij ← A ij -(A ik U -1 kk )(L -1 kk A kj ).
We note that computing a decomposition A kk = L kk U kk and solving the triangular systems is equivalent to a multiplication by A -1 kk . We now introduce the use of low-rank approximations in order to obtain Algorithm 1.4. The main idea is that we store the computed LU factors as a BLR matrix. Moreover, whenever a low-rank block appears in a formula that involves a matrix product, choosing the right associativity leads to a reduced number of operations.

Algorithm 1.4 LU factorization of a BLR frontal matrix (Left-looking)

1: /* Input: a q × q block frontal matrix A; A = [A i,j ] i=1:q,j=1:q ; q = q fs + p nfs */ 2: for k = 1 to q fs do 3:

for i = k to q do 4: Update (L): A i,k ← A i,k -k-1 l=1 L i,l U l,k 5:
end for 6:

for j = k + 1 to q do 7: Update (U ): A k,j ← A k,j -k-1 l=1 L l,j U k,j 8:
end for 9:

Factor: Compute LU factorization L k,k U k,k = A k,k 10
:

for i = k + 1 to q do 11: Compress (L): compute L i,k ≈ A i,k 12: Compress (U ): compute U k,j ≈ A k,j 13: Solve (L): L i,k ← L i,k U -1 k,k 14: Solve (U ): U k,i ← L -1 k,k U k,j 15:
end for 16: end for 17: for i = q fs + 1 to q do 18:

for j = q fs + 1 to q do 19:

Update (CB):

A i,k ← A i,k -k-1 l=1 L l,j U k,j 20:
if CB compression is activated then 21:

Compress (CB): compute T i,k ≈ A i,k 22:
end if

23:

end for 24: end for When the block (i, j) is in low-rank form, we can reduce the complexity of the Solve step (lines 13 and 14) by performing the following operations:

• If i > j: Solve (L) can be rewritten as Y i,j ← (U T j,j ) -1 Y i,j • If i < j: Solve (U ) can be rewritten as Y i,j ← L -1 i,i Y i,j
As a result, whether the block is considered as low-rank or not, the operations performed for the Solve step require solving a triangular system with several right-hand sides (r ij right-hand sides instead of b if low-rank is used, which systematically reduces the complexity).

Similarly, a product of the form C = L ik × U kj from line 4, 7, or 19 may have its complexity reduced if at least one of the two blocks is in low-rank form. Several cases appear:

• LR×FR product: C = X ik × (Y T ik × U kj ) • FR×LR product: C = (L ik × Y kj ) × X T kj • LR×LR product: C = (X ik × (Y T ik × Y kj )) × X T kj or C = X ik × ((Y T ik × Y kj ) × X T kj )
depending on which side leads to the lower complexity.

Triangular solution on a frontal matrix

We consider the triangular solution phase presented in section 1.1.2, and we adapt it to the case where a front of the factors is stored as a BLR matrix (see Figure 1.5).

Similarly to what have been done for the factorization, it is possible to drastically reduce the complexity because of the BLR compression.

We first rewrite the forward and backward substitutions as tile algorithms, and we obtain Algorithms 1.5 and 1.6. Their main cost is a product between a block of the factor and a block of the RHS (lines 4 and 3 respectively), by calling the function prod.

In the context of the use of BLR compression, the expression prod( M , N ) is considered as a LR×FR product if M is in low-rank form XY T , and evaluated with the following associativity: X × (Y T × N ). On the other hand, if M is a full-rank matrix, then we simply have to perform the matrix product M × N instead.

L X

FS variables: q f s blocks CB variables: 1: for j = 1 to q fs do 2:

q cb blocks
B j ← L -1 jj B j 3:
for i = j + 1 to q do 4:

B i ← B i -prod( L i,j , B j ) 5:
end for 6: end for Algorithm 1.6 Backward elimination (Left-looking algorithm) 1: for i = q fs to 1 do 2:

for j = i + 1 to q do 3:

X i ← B i -prod( U i,j , B j ) 4: end for 5: B i ← U -1
ii B i 6: end for

Low-Rank updates accumulation

Low-rank updates accumulation (LUA) consists in grouping together low-rank updates on the same block rows and/or block columns and applying them with a single matrix multiplication to increase the granularity of the computation.

The LUA technique was originally proposed by Amestoy et al. (2019a) for the outer product operation in the BLR LU factorization, and is used by default in MUMPS. In the BLR factorization, we compute low-rank updates of the form

A ij ← A ij -(X ik Y T ik )(Y kj X T kj ).
In the left-looking variant of the BLR factorization, these low-rank updates can instead be grouped as

A ij ← A ij -Xi,1,j • • • Xi,K,j ×    Ȳ T i,1,j • • • Ȳ T i,K,j   
where either Xikj = X ik (Y T ik Y kj ) and Ȳikj = X T kj or Xikj = X ik and Ȳikj = (Y T ik Y kj )X T kj , depending on the ranks. The product of block matrices in the above expression can be efficiently evaluated as a single matrix multiplication Xij Ȳ T ij .

1.3 Floating-point arithmetic and rounding error analysis

Floating-point arithmetic

Definition

In computer science, it is not possible to represent any real number f ∈ R. Instead, we have to use an inexact representation, for example floating-point arithmetic. The idea is to approximate f in such a way that the relative error of approximation is bounded by a certain constant u. In order to do so, its approximation f has a mantissa with a fixed number of digits, scaled by an exponent allowing cover of a large range of values. Thus, we define a floating-point number as follows:

f = (-1) s × b e-t+1 × m = (-1) s × b e × t-1 i=0 m i b -i (1.5) • b is the base • t is the number of significant digits
• s is the bit of sign (0 or 1)

• e ∈ [[e min ; e max ]] is the exponent • m ∈ [[0; b t [[ is the mantissa, whose representation in base b is m 0 ...m t-1
The constants b, t, e max and e min = 1 -e max define the floating-point format. On the other hand, s, e and m are the parameters corresponding to the encoding of f in this format. We will only focus on binary formats in this thesis, so we will take b = 2.

With this definition, f can have several encodings in a floating-point format. However, such a property is not always acceptable. Therefore, we define the normalized representation of a floating-point number, in which the constraint m ≥ b t-1 has to be respected. As a result, we obtain the property of uniqueness of a normalized representation.

The IEEE-754 standard was first published in 1985 and revised in 2008 (see IEEE Computer Society, 2008). It defines floating-point arithmetics that are portable and reproducible, but also the behavior of their basic operations, their rounding rules, as well as their handling of exceptions. In particular, the standard defines a set of binary interchange formats, having normalized representations whose bit encodings are entirely specified.

In binary interchange formats, the integer e is encoded as a positive unsigned integer, offset by a bias in order to reach negative values as well. We notice that, for a binary format, the first bit of the mantissa is m 0 = 1 for a normalized representation. Therefore, m 0 does not have to be stored explicitly, it can be an implicit bit instead: the mantissa is encoded by the bit string m 1 • • • m t-1 . Also, it is specified that the bit string representing a floating-point number in an interchange format follows the order s, e, m. In particular, the least significant bits of m are stored last. We will make use of this property in chapter 3.

Range of representation

A floating-point format represents a finite set of values, and therefore has a limited range. The result of an operation may be outside this range. If the exponent is beyond e max , we obtain an overflow: the result of the operation is considered as +∞ or -∞.

This may very well lead to an irrelevant result because it causes a propagation of NaN (Not a Number) values. Thus, overflows should generally be avoided at all costs. This is especially the case when handling precision formats such as fp16, whose range is small. As shown in Table 1.1, converting to fp16 a number higher than 10 5 will result in an overflow.

On the other hand, if the value to be stored in a floating-point format is smaller than the minimum value of the format, that is to say if the exponent should be smaller than e min , we say that an underflow occurs. A common choice for handling this case is to round the result to zero. Alternatively, by not imposing the use of a normalized representation, the result can still be approximated as a so-called subnormal number: the first bits of the mantissa are zero. However, the relative accuracy is inferior to the normalized numbers, given the fact that the number of significant bits is reduced. Besides, the use of subnormal numbers is not often implemented in hardware but in software instead, leading to a significant slowdown.

However, in practice the overflows and underflows tend to be rare enough that they usually do not have too much impact on computations. This is especially true for fp64 which has a very large range: [10 -308 ; 10 308 ]. The smaller ranges of fp32 and bfloat16 ([10 -38 ; 10 38 ]) are still generally sufficient for avoiding most issues of overflow and underflow, contrary to fp16. In this thesis, we will consider that the presence of overflows and underflows can be neglected, and that the exponent e can take any integer value.

Unit Roundoff

In addition to the range of its representable values, an important parameter describing the properties of a floating-point format is its unit roundoff u. Let us consider that there is no restriction on the exponent: e ∈ Z. For any real x we note f l(x) its floating-point approximation in the format.

The unit roundoff u of the format is a constant approaching the maximal relative error of approximation of a real number in this format. Such a maximum is obtained when approximating x = 1 + 2 -t by f l(x) = 1 for instance, obtaining a relative error of

|x-f l(x)| |x| ≈ 2 -t
. Therefore, we define u = 2 -t . For every real x within the range, we find that there is a δ such that

f l(x) = x(1 + δ), |δ| < u (1.6)
Equation 1.6 is a rather fundamental property of the rounding operator. It is often used instead of the original definition of floating-point rounding (Equation 1.5). However, it does not contain the same level of information.

In order to carry out a rounding error analysis of an algorithm, we need to make some assumptions about the accuracy of its basic arithmetic operations, namely +, -, × and /. Indeed, the result of such an operation is not always representable in the same format, and the computation cannot be exact in this case. The most common assumption is the standard model of floating-point arithmetic. Given two floating-point numbers x and y, if we note x op y the result of one of the basic operations and f l(x op y) its computed value, this model states that the relative error introduced is at most u:

f l(x op y) = (x op y)(1 + δ), |δ| ≤ u, op ∈ {+, -, ×, /} (1.7)
This property is a consequence of the IEEE-754 standard in which those operations are entirely specified: it is stated that f l(x op y) should give the same result as computing x op y in infinite precision before rounding it to the floating-point format in use. However, this specification gives more information than the standard model. In particular, Equation 1.7 loses track of the fact that some of the operations are performed exactly.

The basics of rounding error analysis

Backward and forward errors

Let us consider a function f : R n → R m . We want to assess the quality of a computed approximation y of y = f (x). For example, if f consists in a series of exact basic operations, such an approximation may be the result of using floating-point arithmetic. There are several possibilities of error indicators available, depending on the context. A first approach is to compute the absolute error committed on the output variable y, i.e., || y -y||. A bound on this term is referred to as the (normwise) forward error.

However, one may prefer to focus on finding for which modified set of data we have actually solved the exact problem, i.e., which perturbation ∆x of the input x is needed in order to obtain y = f (x + ∆x). Among the many ∆x that satisfy this relation we will consider the one having the smallest norm. As a result, a formal definition of the (normwise) backward error of y ≈ f (x) would be:

η(y) = min ||∆x|| ||x|| : ∆x ∈ R n such that y = f (x + ∆x) (1.8)
The concept of backward error is illustrated in Figure 1.6.

Accumulation of rounding errors

We illustrate the definition of the backward error by trying to quantify the impact of rounding errors on the computation of a certain sum s n = n k=1 x k , in the order of its indices. We note s n its value computed using floating-point additions (more generally, we will tend to use the notation var to denote the computed value approximating a variable var in this thesis).

We note δ 1 , • • • , δ n-1 the coefficients obtained for the relative errors of each addition (Equation 1.7). We have the recursive relation

s n = f l( s n-1 + x n ) = (1 + δ n-1 )( s n-1 + x n ),
and s 1 = x 1 . As a result, we give the expression of the computed sum: We choose to express each π k under the form

s n = (π 1 x n + • • • π n-1 x 2 ) + π n-1 x 1 , after having defined each quantity π k = (1 + δ 1 ) • • • (1 + δ k ).
π k = 1 + θ k . Then we obtain |θ k | ≤ γ k , with γ k = ku 1-ku .
The presence of such a weight γ k keeps track of the fact that x n-k+1 appears in k additions.

As a result, we obtain s n = n k=1 (1 + θ k )x k , with a bound on the relative error committed on the input variable

x = (x 1 , • • • , x n ): ∀k, |θ k | ≤ γ k ≤ γ n-1 . Because ||θ|| ≤ γ n-1
, the relative backward error of the computation of the sum is γ n-1 (at most). This constant is quite common in rounding error analysis.

On the other hand, the error committed on the output variable s n is the forward error. For the sum, a rough componentwise bound would be | s n -

s n | = | n k=1 θ k x k | ≤ γ n-1 n k=1 |x k |.
Therefore, the normwise forward error would be ||x||

Matrix-vector product

The rounding error analysis of a matrix-vector product y = Ax is quite similar to that of the sum that we did previously. As explained in detail in [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF], the approximate result y can be expressed as y = (A + ∆A)x, with |∆A| ≤ γ n |A| (backward error). As a consequence, the componentwise forward error bound is | y -y| ≤ γ n |A| |x|. We note that, here and throughout the rest of the thesis, we use the notation A ≤ B to denote the elementwise inequality a ij ≤ b ij on the coefficients of matrices A and B. Moreover, we use the notation |A| to denote the matrix obtained by taking the absolute values of the coefficients of A.

LU factorization

We consider the solution of Ax = b via a Gaussian elimination of A (by combining Algorithms 1.1, 1.2 and1.3). Although we could use a multi-word arithmetic in order to obtain an exact solution, it is much more common to use floating-point arithmetic instead, which is much faster. However, one may wonder whether the impact on the accuracy of the solution would be acceptable or not. Therefore, we would like to obtain its backward error bound. We look for a ∆A such that (A + ∆A) x = b. We know from [START_REF] Rigal | On the compatibility of a given solution with the data of a linear system[END_REF] that such a ∆A of minimal norm verifies:

η(x) = ||∆A|| = ||A x -b|| ||A|| || x|| (1.9)
This result stands for any subordinate matrix norm. η(x) is the relative normwise backward error, defined in Equation 1.8.

We would like to obtain a reasonable bound on ∆A, ideally in O(u)||A||. However, the best obtainable componentwise bound obtained is not as optimistic:

|∆A| ≤ γ 3n | L| | U | (1.10)
Indeed, the coefficients of | L|| U | are generally much larger than the coefficients of A. In order to quantify this increase, we define the growth factor of the factorization:

ρ n = max i,j,k |a (k) ij | max i,j |a ij | (1.11)
where a

(k)
ij is an intermediate result obtained after performing k updates on the original coefficient a ij . James Wilkinson showed in the early 1960s that the numerical stability of Gaussian elimination intrinsically depends on the growth factor. When using partial pivoting, we can find ∆A such that

||∆A|| ∞ ≤ n 2 γ 3n ρ n ||A|| ∞ (1.12)
Therefore, one may wonder how large the growth factor ρ n is. In the worst-case scenario, we may have ρ n = 2 n-1 despite the use of partial pivoting: the backward error obtained is rather catastrophic and the algorithm cannot be considered numerically stable in this case. However, it is highlighted in (Higham, 2002, Section 9.4) that fortunately, "despite the existence of matrices for which ρ n is large with partial pivoting, the growth factor is almost invariably small in practice". Thus, solving linear systems using Gaussian elimination with partial pivoting is stable enough in practice.

More recently, this result of backward stability has been extended to the case of a BLR LU factorization (Algorithm 1.4) in (Higham and Mary, 2021, Theorem 4.3). It is shown that:

A = L U + ∆A, ||∆A|| ≤ (qε + γ q )||A|| + γ c || L|| || U || + O(uε) (1.13)
if the matrix A has q 2 blocks of size b and of rank r, and if c = b + 2r 3/2 + q. If u ≪ ε then the main term of error will be qε||A||.

Mixed-precision algorithms

When an accurate result is needed in scientific computing, the default choice is often to use double precision everywhere, with 64 bits floating-point numbers (fp64). However, lower precision formats like fp32 also exist. They require less storage space and usually run faster: we can typically hope to gain a 2× speedup using fp32 instead of fp64. However, the computations are not as accurate: the relative error of elementary operations (unit roundoff) is higher, and the range is slightly narrower (see Table 1.1).

Mixed-precision algorithms result from the will to do part of the computations in low precision, obtaining time and memory gains, while controlling the accuracy by performing key operations in high precision. In this section, we give some examples of mixed-precision techniques used in the literature in contexts similar to our own. The recent availability of half precision formats in hardware raises the question of how to use them in HPC. In fact, they would allow us to reduce all storage by a factor of 4 compared to double precision (and therefore reduce communication cost), and to possibly have an even better gain regarding computation time.

In some deep learning models where great accuracy is not needed, algorithms tend to use Google's bfloat16 format. The range is the same as in fp32, so the algorithms do not necessarily require many modifications. In linear solvers where we need more precision, fp16's greater relative accuracy would be an advantage. However, fp16 has only a limited range available: positive numbers lie between 6 × 10 -8 and 7 × 10 4 . Standard algorithms can therefore lead to overflow, underflow, or subnormal numbers, all of which are undesirable. This is a serious issue that is almost non-existent for the other main precision formats, and it restricts the use of fp16.

Iterative refinement

One of the best known mixed-precision algorithms is iterative refinement, which was first programmed by Wilkinson in 1948. It aims at improving the accuracy of the solution of a linear system Ax = b.

Algorithm 1.7 Iterative refinement Solve Ax 0 = b in precision u f for i=0:∞ do Compute r i = b -Ax i at precision u r Solve Ad i = r i at precision u s x i+1 ← x i + d i at

working precision u end for

A version of this algorithm consists in computing an LU factorization of the matrix at the beginning and then using it for all the Solve steps. Each task can be performed using a different precision format. For example, with u = u r = fp64 and u f = u s = fp32, the most expensive part of the computation is done entirely in single precision. The computation can therefore be expected to run twice as fast as the double-precision algorithm. This usage was proposed and analyzed by [START_REF] Langou | Exploiting the performance of 32 bit floating point arithmetic in obtaining 64 bit accuracy (revisiting iterative refinement for linear systems)[END_REF]. However, the convergence of this algorithm is not guaranteed if the matrix A is ill conditioned, with a condition number of u -1 f or more. This is why Carson and Higham introduced another version of mixed-precision iterative refinement that does not have this limitation [START_REF] Carson | A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems[END_REF]. In their algorithm, GMRES-IR (GMRES-based Iterative Refinement), each Solve step is computed using GMRES (Generalized Minimal Residual), an iterative method. A lowprecision LU factorization is used as a preconditioner. The reason is that U -1 L -1 A is far better conditioned than A.

In [START_REF] Carson | Accelerating the solution of linear systems by iterative refinement in three precisions[END_REF], the same authors then implemented a version of GMRES-IR with 3 precision formats instead of 2. The Solve step is now handled by GMRES at precision u, except matrix-vector products with A that are computed at precision u r . Their choice of (u f , u, u r ) = (fp16, fp32, fp64) or (fp32, fp64, fp128) allowed them to obtain a much better convergence condition and a smaller forward error than classic iterative refinement using the same working precisions.

Mixed precision on GPUs

GPUs are a type of hardware that is well adapted to matrix computations and to the use of low precisions, which makes them a target of particular interest for mixed-precision algorithms in linear algebra. They handle computations in fp32 at least twice as fast as in fp64. Recent models tend to have support for fp16 and/or bfloat16 arithmetics, with great speedups. This the case for NVIDIA's model A100, available since September 2020, in which fp16 computation offers a 4x speedup compared to fp32 (peak performance).

Moreover, one of the preferred operations on GPUs is block FMA (Block Fuse Multiple-Add), computed in special units referred to as "Tensor Cores" by some vendors. This kind of unit can typically compute operation D = C +AB in one clock cycle, where all matrices have a size of 4 × 4 for instance.

The tensor cores in the NVIDIA Volta, Turing, Ampere and Hopper architectures can perform these FMAs using different combinations of precision formats, including an fp16 input and an fp32 output. The accumulation inside the matrix multiplication is done in fp32. Therefore, mixed precision is implemented at a hardware level, and algorithms may be able to take this fact into account. [START_REF] Haidar | Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers[END_REF] introduced a new class of mixed-precision dense matrix factorization algorithms based on this basic operation. An error analysis of the mixed block FMA as well as the LU factorization based on it was done in [START_REF] Blanchard | Mixed precision block fused multiply-add: Error analysis and application to GPU Tensor Cores[END_REF].

D f p32 = C f p32 + f p32 A f p16 × f p32 B f p16

A scaling algorithm for handling low-precision formats

Using the fp16 format can lead to a huge number of overflows and underflows due to its narrow range. In order to deal with this issue, Higham and Pranesh published in 2019 an algorithm for rounding correctly a matrix to fp16 format (see Higham et al., 2019). First a two-sided diagonal scaling is applied in order to balance the matrix's coefficients. We obtain a factorization D l AD r , with D l and D r diagonal matrices, and every row and column of A has an ∞-norm of 1. Then we multiply A by a scalar to bring the elements of largest magnitude within a factor θ ≤ 1 of the overflow level. We finally round to half precision. The multiplication by a scalar ensures that a large part of the limited range of half precision arithmetic is used, and therefore the number of underflows is reduced.

The authors presented an application to GMRES-IR iterative solver. After rounding in fp16 with the previous algorithm (for negligible cost of O(n 2 ) operations), they compute an LU factorization of the scaled matrix. The parameter θ is chosen carefully to ensure the absence of overflow. With this factorization as a preconditioner, GMRES-IR has a much better convergence rate than using a standard method for rounding to fp16.

A mixed-precision Cholesky factorization

A first method for using mixed precision inside a dense tile Cholesky factorization was proposed in [START_REF] Abdulah | Geostatistical modeling and prediction using mixed precision tile cholesky factorization[END_REF]. The authors chose to address specifically the factorization of a dense covariance matrix in the field of geostatistics, rather than a more general approach. No compression technique is used in this algorithm.

Their mixed-precision approach consists in storing some blocks of the matrix in double precision, and the others in single precision. A simple criterion was used to choose between the two possibilities: every block within a certain distance from the diagonal is kept in double precision, and the others are stored instead in single precision (Figure 1.7). The precision of the computation is based on the precision of the target block being updated. One consequence of this choice is that each computed block is needed at some point both in single precision and in double precision for the updates. As a result, without a careful implementation, this mixed precision algorithm could introduce an extra memory cost (up to 50%).

On the other hand, they succeeded in obtaining a 1.6× performance speedup on massively parallel architectures while maintaining the accuracy necessary for their modeling and prediction.

A mixed-precision representation of H-matrices

The article Ooi et al. (2020a) introduces a mixed-precision algorithm that uses Hmatrices, a low-rank matrix format that is based on hierarchical partitioning of the matrix. No error analysis was considered.

Each low-rank block is separated between a part in double precision and another in single (Figure 1.8). Then we obtain our scaling from diagonal scalings of the 2 factors:

D v = diag(∥v 1 ∥ ∞ , ..., ∥v p ∥ ∞ ) D w = diag(∥w 1 ∥ ∞ , ..., ∥w p ∥ ∞ )
We then scale using diagonal matrix

D = D v × D w .
This type of scaling can guarantee the absence of any overflow and at the same time avoid most underflows. While its usefulness may be limited if we only work with single or double precision, if we consider the use of fp16 format as well, such a scaling method might be helpful due to the huge risk of overflow and underflow associated with the very narrow range of this format.

Using low precision for storing or accessing data

Another strategy for using mixed precision consists in using low precision formats for storage and/or communications, while keeping the computations in high precision formats. This strategy not only reduces the memory footprint and communication volume of the algorithm, but also reduce its runtime, thanks to the reduced data movements.

In fact, there are several examples in the literature where data movement reductions have successfully been turned into time gains. In [START_REF] Mukunoki | Reduced-precision floating-point formats on gpus for high performance and energy efficient computation[END_REF], a set of 7 floating-point formats was used for data storage (the same set as our own, see Table 3.1). On several GPU architectures, they called memory-bound functions from cuBLAS1 , such as GEMV and AXPY, using a custom format for memory accesses and a certain working precision (double or single) for computations. They succeeded in obtaining very decent speedups, up to a factor 2.

Grützmacher et al. ( 2023) developed a similar idea of decoupling the formats for computation and storage, based on Ginkgo's memory accessors2 this time. By doing this, double precision should only be handled at the register level. They applied their method to memory-bound functions in double precision, while accessing single-precision variables. On GPUs as well as on CPUs, they attained a performance comparable to the single-precision functions.

Chapter 2

Dense LU factorization in mixed precision

Introduction

In this chapter, we investigate the potential of combining mixed-precision arithmetic with low-rank approximations. Let A ∈ R m×n and let XΣY T = min(m,n) i=1

x i σ i y T i be its singular value decomposition (SVD). As mentioned in section 1.2.1, given a target accuracy ε, a low-rank approximation T of A satisfying ∥T -A∥ ≤ ε∥A∥ can be built from the truncated SVD T = r i=1 x i σ i y T i , where r is the rank of T . Our starting idea for this work is to ask what precision should be used to store T and to operate on it. In the literature, the truncated SVD (or any other form of lowrank decomposition such as truncated QR with column pivoting), tends to be stored in the lowest possible precision with unit roundoff safely smaller than ε. For example, if ε = 10 -12 and we have access to the floating-point arithmetics defined by the IEEE standard, existing algorithms would use double precision (for which the unit roundoff is u d ≈ 1 × 10 -16 ), because the next lower precision, single precision, has a unit roundoff u s ≈ 6 × 10 -8 that is much larger than the prescribed ε.

However, we explain why and how we can actually exploit much lower precisions, with almost no loss of accuracy. We show that singular vectors associated with sufficiently small singular values can be stored at precisions with unit roundoff larger than ε∥A∥ while maintaining an overall approximation accuracy of order ε∥A∥. For example, with ε = 10 -12 , any singular vector x i whose associated singular value σ i is smaller than ε∥A∥/u s ≈ 2 × 10 -5 ∥A∥ can be stored in single precision. Indeed, the single precision vector x i satisfies ∥ x i -x i ∥ ≤ u s , but the overall error introduced by replacing x i by x i is bounded by ∥(

x i -x i )σ i y T i ∥ ≤ (ε∥A∥/u s )u s = ε∥A∥.
As can be seen from this example, the reason we can afford to convert some singular vectors to lower precision is because the error introduced by this conversion is demagnified by the singular value; hence the error may be safely bounded if σ i is small enough.

In the following, we formalize this intuition with an error analysis considering an arbi-trary number of precisions. Moreover, our analysis applies to any low-rank decomposition of the form T = XY T where X has orthonormal columns. Indeed, the mixed-precision approach presented here is general and can be used for several other low-rank approximations, in particular rank-revealing QR decompositions.

Clearly, the potential of the proposed approach depends on whether the singular values of the matrices to be approximated decay rapidly. In the second part of this chapter, we apply this approach to an important class of matrices: data sparse, rank-structured matrices, whose off-diagonal blocks have low numerical rank [START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF]. We focus in particular on the block low-rank (BLR) format [START_REF] Amestoy | Improving multifrontal methods by means of Block Low-Rank representations[END_REF][START_REF] Amestoy | On the complexity of the Block Low-Rank multifrontal factorization[END_REF], although the approach is also applicable to hierarchical formats. Our numerical experiments demonstrate that the proposed mixed precision low-rank representation presents a very high potential in this context: a large fraction of both the entries needed to represent BLR matrices and the floating-point operations (flops) needed to compute their LU factorization can be switched to lower precisions.

The rest of this chapter is organized as follows. In section 2.2, we describe the proposed mixed precision low-rank representation and show that the loss of accuracy introduced by the use of lower precisions can be rigorously bounded. We then apply this representation to BLR matrices in section 2.3. In section 2.4, we analyze how to compute the LU factorization of a BLR matrix using mixed precision arithmetic. We present numerical experiments on a range of real-life matrices in section 2.5, before concluding in section 2.6.

Throughout the chapter, we define γ

∥A∥ = m i=1 n j=1 a 2 ij 1/2 = min(m,n) i=1 σ 2 i 1/2 .
(2.1)

In our experiments, we will work with IEEE double and single precision arithmetics (denoted as fp64 and fp32, respectively) and bfloat16 arithmetic. The unit roundoffs for these three arithmetics are

u d = 2 -53 ≈ 1 × 10 -16 , u s = 2 -24 ≈ 6 × 10 -8 , and u h = 2 -8 ≈ 4 × 10 -3 , respectively.

Low-rank approximations in mixed precision

Let A ∈ R m×n and let T be a low-rank approximation of A satisfying

∥A -T ∥ ≤ εβ, (2.2)
where ε > 0 is the target accuracy and where β is a scaling parameter chosen by the user: a natural choice is β = ∥A∥, which leads to an accuracy of ε relative to ∥A∥, but other choices are possible, as mentioned in section 1.2.2.

Hereinafter, we refer to the precision that T is stored in as the working precision, and we assume that its unit roundoff u 1 is safely smaller than ε, that is, u 1 ≪ ε.

Given the SVD A = n i=1 x i σ i y T i , it is well known that the approximation of A of lowest rank is given by the truncated SVD

T = XΣY T = r i=1 x i σ i y T i , X ∈ R m×r , Y ∈ R n×r , (2.3)
where the rank r is the smallest integer such that (2.2) is satisfied. Neglecting any noise associated with the working precision, r is given by

r = min k : ∥A - k i=1 x i σ i y i ∥ ≤ εβ = min    k : min(m,n) i=k+1 σ 2 i 1/2 ≤ εβ    .
(2.4)

The goal of this section is to prove that depending on the singular values of T , we can use lower precisions than the working precision (with unit roundoff larger than ε), and still preserve an overall approximation error of order ε. We first carry out our analysis for the SVD, but also provide at the end of this section its extension to other types of low-rank approximation methods, such as rank-revealing QR.

Our analysis assumes that the use of lower precision arithmetic with fewer exponent bits than the working precision does not lead to any overflow or underflow. To ensure that this assumption is satisfied in our experiments, we focus on the use of bfloat16 arithmetic (which has the same range as fp32), rather than fp16 (which has a much narrower range).

Let us assume that p different floating-point arithmetics are available (including the working precision u 1 ), and that their unit roundoffs satisfy

u 1 ≪ ε < u 2 < . . . < u p .
(2.5)

Let us consider a partition of matrix T into p groups

T = XΣY T = X 1 . . . X p     Σ 1 . . . Σ p     Y 1 . . . Y p T , (2.6)
where

T k = X k Σ k Y T k , X k ∈ R m×r k , Σ k ∈ R r k ×r k , Y k ∈ R n×r k (2.7)
is a matrix of rank r k formed by the subset of the singular values and vectors of T assigned to group k.

We now analyze the effect of converting T k to precision u k . We assume that only the singular vectors X k and Y k are converted, whereas the singular values Σ k are kept in precision u 1 . This is because the storage for Σ k is negligible compared with that of X k and Y k . We however note that adapting the analysis to the case where Σ k is also converted to precision u k is straightforward and only slightly increases the constants in the error bounds. We write X k and Y k for the converted vectors, and

T k = X k Σ k Y T k (note that since T 1 is already in precision u 1 , T 1 = T 1 ). The following lemma bounds ∥ T k -T k ∥ for k ≥ 2.
Since X k and Y k are stored in precision u 1 , they only have approximately orthonormal columns. In the following lemma and throughout this chapter, we neglect this loss of orthonormality, which would only introduce a lower order term O(εu 1 ) in the error bounds.

Lemma 2.1. Let T k = X k Σ k Y T
k where X k and Y k have approximately orthonormal columns (stored in precision u 1 ), and let

T k = X k Σ k Y T k be obtained by converting X k and Y k to precision u k . Then ∥ T k -T k ∥ ≤ (2 + √ r k u k )u k ∥Σ k ∥. (2.8) Proof. The converted X k and Y k satisfy, for k ≥ 2, X k = X k + E k , |E k | ≤ u k |X k |, (2.9) Y k = Y k + F k , |F k | ≤ u k |Y k |.
(2.10) Therefore, we have

∥T k -T k ∥ ≤ ∥E k Σ k Y T k ∥ + ∥X k Σ k F T k ∥ + ∥E k Σ k F T k ∥.
(2.11)

For the first term, we observe that

∥E k Σ k Y T k ∥ 2 = ∥E k Σ k ∥ 2 (2.12) = j σ 2 j i e 2 ij (2.13) ≤ j σ 2 j i u 2 k x 2 ij (2.14) = u 2 k j σ 2 j = u 2 k ∥Σ k ∥ 2 , (2.15)
where we have used the fact that the columns of X k have a norm of 1. Therefore

∥E k Σ k Y T k ∥ ≤ u k ∥Σ k ∥.
(2.16)

Similarly, we also have

∥X k Σ k F T k ∥ ≤ u k ∥Σ k ∥.
(2.17)

Finally, for the third term, we have

∥E k Σ k F T k ∥ ≤ ∥E k Σ k ∥∥F k ∥ and so ∥E k Σ k F T k ∥ ≤ √ r k u 2 k ∥Σ k ∥. (2.18)
Reinjecting (2.16), (2.17), and (2.18) into (2.11) yields the result.

Lemma 2.1 shows that converting T k to precision u k introduces an error of order u k ∥Σ k ∥, which is thus proportional to the size of the singular values in Σ k . This fundamental observation is at the foundation of the mixed precision representation that we propose. Indeed, Lemma 2.1 suggests that we can preserve an overall accuracy of order εβ by partitioning the singular values in such a way that ∥Σ k ∥ ≈ εβ/u k .

In order to build such a partitioning where the size of the groups stored in lower precision is as large as possible, it is easy to see that we should start by including the smallest singular values in the last group first, until its norm exceeds εβ/u p ; at this point, we can start building group p-1 with the remaining singular values, and so on. Therefore, the Σ k are formed from consecutive singular values:

Σ k = diag(σ i ), i = i k : i k+1 -1, (2.19)
where the indices i k and i k+1 define which singular values are part of Σ k , and can be computed by the recursive formula:

i k = min i : i k+1 -1 j=i σ 2 j 1/2 ≤ εβ/u k , k ∈ [2 : p], (2.20)
starting with i p+1 = r + 1 and ending with i 1 = 1. We thus end up with a partitioning of the SVD as defined by (2.19)-(2.20). Note that this definition may lead to some empty Σ k , in which case T k is a rank-0 matrix. We note that this partitioning is similar to the Method 3 proposed by Ooi et al. (2020b). Our analysis justifies the use of this partitioning and gives a precise rule to define the p groups depending on the singular values and on the precisions. This partitioning guarantees that ∥Σ k ∥ ≤ εβ/u k for all k ≥ 2 and so, by Lemma 2.1, converting T k to precision u k introduces an error bounded by

∥T k -T k ∥ ≤ (2 + √ r k u k )εβ. (2.21)
By combining (2.21) over k = 2 : p, we readily obtain a bound on the overall error introduced by converting each T k to precision u k .

Theorem 2.1. Let T be a low-rank approximation of

A satisfying ∥A -T ∥ ≤ εβ. If T is partitioned into p groups T k = X k Σ k Y T k as defined by (2.19)-(2.20)
, and the X k and Y k are converted to precision u k , the resulting matrix

T = p k=1 T k = p k=1 X k Σ k Y k satisfies ∥A -T ∥ ≤ 2p -1 + p k=2 √ r k u k εβ. (2.22)
Proof. The triangle inequality ∥A -T ∥ ≤ ∥A -T ∥ + p k=2 ∥T k -T k ∥ together with (2.21) readily yields the result.

Theorem 2.1 shows that the mixed precision low-rank matrix T approximates A with an accuracy of order ε. To first order, the constant in this error bound is 2p -1, instead of 1 for the uniform precision matrix T : the introduction of lower precisions therefore only increases the overall error by a very modest quantity. Moreover, we note that by means of a more sophisticated proof that avoids the use of the triangle inequality, this constant can be reduced to 1 + 2 √ p -1:

∥A -T ∥ ≤ 1 + 2 p -1 + O(u 2 ) εβ.
(2.23) However, we will not use such proofs for the sake of readability, and because the precise value of the constants in the error bounds is unimportant, as long as they are not too large.

Note that while Theorem 2.1 guarantees an approximation error ∥A -T ∥ in O(ε), the errors on the low-rank factors ∥X -X∥ and ∥Y -Y ∥ are in O(u p ), the lowest precision used. Moreover, the loss of orthogonality of X and Y is also in O(u p ). Therefore, one should be careful before using the proposed method in applications that make use of the singular vectors or their orthogonality.

Importantly, the proposed method can be applied to other types of low-rank decompositions, not necessarily based on SVD. For example, Theorem 2.1 can be easily extended to decompositions of the form XBY T , where X and Y have orthonormal columns (the difference with the SVD being that B is not diagonal), bound (2.22) holds with a slightly larger constant; one example of this form is the UTV decomposition [START_REF] Fierro | Low-rank revealing utv decompositions[END_REF]. Our analysis can also be adapted to decompositions of the form XY T , where X has orthonormal columns (but Y does not). This second form is particularly of interest because it applies to rank-revealing QR decompositions. We state the analogue to Lemma 2.1 for XY T decompositions below.

Lemma 2.2. Let T k = X k Y T
k where X k has approximately orthonormal columns (stored in precision u 1 ), and let

T k = X k Y T k be obtained by converting X k and Y k to precision u k . Then ∥T k -T k ∥ ≤ (2 + √ r k u k )u k ∥Y k ∥. (2.24)
Thus, the error introduced by converting group k now depends on ∥Y k ∥, and this means that the XΣY T partitioning (2.19)-(2.20) should be adapted by replacing ∥Σ k ∥ by ∥Y k ∥. Then, it is easy to show that Theorem 2.1 still holds. Lastly, the proposed method can even be applied to decompositions XY T where neither X nor Y have orthonormal columns. In this case the error introduced by the conversion of

T k to precision u k is ∥T k -T k ∥ ≤ (2 + u k )u k ∥X k ∥∥Y k ∥, (2.25)
which shows that there is potential for mixed precision as long as the low-rank components x i y T i are of decreasing norm.

In the rest of this chapter, we will focus on low-rank decompositions of the form XY T , computed by means of a truncated QR factorization with column pivoting.

An important question is under what condition the low-rank compression is beneficial, that is, when does the low-rank approximation T require less storage than the original matrix A ∈ R m×n . In the standard uniform precision case, T = XY T can be represented with r(m + n) entries, and so the condition is

r(m + n) ≤ mn.
(2.26)

With a mixed precision representation, this condition changes due to the fact that entries belonging to groups k ≥ 2 are stored in lower precision. The condition becomes

(m + n) p k=1 c k r k ≤ mn, (2.27)
where c k quantifies the cost of storing a floating-point number in precision u k instead of u 1 . For example, if we use three precisions, fp64, fp32, and bfloat16, (2.27) takes the form

(m + n)(r 1 + 0.5r 2 + 0.25r 3 ) ≤ mn.
Interestingly, the difference between conditions (2.26) and (2.27) means that a matrix that is not "low-rank enough" in uniform precision can become so when using mixed precision arithmetic.

Crucially, the size r k of each group depends on the singular values. Indeed, group k must satisfy ∥Σ k ∥ ≤ εβ/u k , so if the singular values of A decay slowly, most of them must be kept in the first group and little gain can be expected from the use of mixed precision. Conversely, if A possesses many small singular values, the low precision groups will be very large and the use of mixed precision will be very beneficial. Therefore, the potential gains achieved by the proposed approach completely depend on the singular values of the matrix. From now on, we will focus on an important class of matrices that exhibit off-diagonal blocks with rapidly decaying singular values, and therefore present a high potential for the use of mixed precision.

Mixed precision BLR compression

Data sparse matrices are rank-structured matrices most of whose off-diagonal blocks have low numerical rank. In this section, we show how this property can be exploited to represent these matrices in mixed precision. We focus on a specific class of data sparse matrices, called the block low-rank (BLR) format [START_REF] Amestoy | Improving multifrontal methods by means of Block Low-Rank representations[END_REF][START_REF] Amestoy | On the complexity of the Block Low-Rank multifrontal factorization[END_REF](Amestoy et al., , 2019a)). The approach described here could also be extended to other formats, such as hierarchical [START_REF] Hackbusch | Hierarchical Matrices : Algorithms and Analysis[END_REF] or multilevel (Amestoy et al., 2019b) ones.

Background on BLR matrices

As explained in section 1.2.2, a BLR representation T of a dense square matrix A ∈ R n×n has the block q × q form

T =        T 11 T 12 • • • T 1q T 21 • • • • • • . . . . . . • • • • • • . . . T q1 • • • • • • T qq        , (2.28)
where some of the off-diagonal blocks A ij of size b×b have been approximated by matrices T ij of ranks r ij :

T ij =    X ij Y T ij if i > j Y ij X T ij if i < j (2.29)
where X ij and Y ij are b × r ij matrices, and where X ij has orthonormal columns. The other blocks are left uncompressed, with T ij = A ij . In both cases, we ensure that

∥A ij -T ij ∥ ≤ εβ ij , (2.30)
with β ij = ∥A∥ if we use a global threshold.

On a BLR matrix, it leads to the global error bound ∥A -T ∥ ≤ qε∥A∥.

(2.31)

Error analysis of mixed precision BLR compression

We now seek to combine BLR compression with the mixed precision representation proposed in section 2.2. The natural approach is to simply use this mixed precision representation on every low-rank off-diagonal block of the BLR matrix, leaving the fullrank blocks in the working precision u 1 . Then, it is easy to show that (2.22) becomes the resulting BLR matrix T satisfies

∥A ij -T ij ∥ ≤ 2p -1 + p k=2 r (k) ij u k εβ ij , (2.32) where r (k) ij is the rank of the matrix T (k) ij = X (k) ij ( Y (k) ij ) T ,
∥A -T ∥ ≤ q 2p -1 + p k=2 c k u k ε∥A∥,
(2.33)

with c k = max i,j r (k) ij . Proof. Using ∥A -T ∥ 2 = q i=1 q j=1 ∥A ij -T ij ∥ 2 (2.34)
and (2.32), we readily obtain the result.

Compared with the uniform precision bound (2.31), the mixed precision bound (2.33) is thus larger by a modest factor of about 2p-1. Theorem 2.2 therefore shows that we can exploit mixed precision arithmetic in the BLR compression while preserving an accuracy of order ε. First, the full-rank blocks (type 1, dark blue blocks in Figure 2.1, consisting of only the diagonal blocks here) are stored in the working precision (fp64). An example of diagonal block is given in Figure 2.2a, showing that its singular values decay too slowly to benefit from the use of a mixed precision representation. However, there is only a small number of such blocks: the majority of the blocks therefore benefits from the use of lower precisions. Interestingly, the number of low-rank blocks that effectively use all three precisions is quite small (type 2, light blue blocks, example given by Figure 2.2b). Most blocks actually do not need to store any of their entries in fp64. This is a consequence of using global compression: if β ij = ∥A∥ ≫ ∥A ij ∥, fp64 is not needed to represent A ij . In other words, blocks of sufficiently small norm can be stored entirely in lower precision. Among these blocks, we can further distinguish two types: those that are represented in mixed precision using both fp32 and bfloat16 (type 3, green blocks, example given by Figure 2.2c) and those that are stored entirely in bfloat16 (type 4, yellow blocks, example given by Figure 2.2d). The type-4 blocks are those whose norm is smaller than ε∥A∥/u h , where u h = 2 -8 is the unit roundoff of bfloat16. Note that a fifth type of block could arise, those whose norm is smaller than ε∥A∥: these blocks can simply be dropped, that is, replaced by zero (but no such blocks appear in the example of Figure 2.1).

Types of mixed precision blocks

The observation that blocks of small norm can be stored entirely in lower precision is important. It justifies why the simpler approach proposed by [START_REF] Abdulah | Geostatistical modeling and prediction using mixed precision tile cholesky factorization[END_REF] and [START_REF] Doucet | Mixed-precision tomographic reconstructor computations on hardware accelerators[END_REF] can already achieve significant gains. Their approach consists in storing each block in uniform precision, but possibly differing from one block to another. For example, for the matrix in Figure 2.1, all type-2 blocks (light blue) would need to be stored entirely in double precision, but type-3 blocks (green) could be stored entirely in single precision. Moreover, our error analysis provides the criterion that should be used to choose each block's precision: blocks A ij such that ∥A ij ∥ ≤ ε∥A∥/u i can be stored in precision u i .

Mixed precision BLR LU factorization

We now present how to exploit the mixed precision BLR representation described previously in order to accelerate the LU factorization of BLR matrices. There exists several BLR LU factorization algorithms, and here we focus on the so-called UCF (a.k.a. UCFS or FCSU) variant described in Algorithm 2.1. This variant has been successfully used in the literature, for example in the MUMPS (Amestoy et al., 2019a) and PaStiX [START_REF] Pichon | Sparse supernodal solver using Block Low-Rank compression: Design, performance and analysis[END_REF] sparse direct solvers, and its rounding error analysis in the uniform precision case has been carried out by Higham and Mary (2021, sect. 4.2). We note that BLR LU factorization can and usually does incorporate numerical pivoting for stability, but we describe Algorithm 2.1 without pivoting for simplicity. Note that in the UCF variant, the matrix is not compressed from the beginning: instead, we perform the compression on the fly during the LU factorization (step 11 of Algorithm 2.1). To compress the R ik matrices, we use the same truncation criterion as in (2.30)

∥R ik -T ik ∥ ≤ εβ ik , (2.35)
where β ik > 0 is a parameter whose role has been discussed in section 1.2.2 (in practice we will use global compression by setting β ik = ∥A∥).

Algorithm 2.1 BLR LU factorization.

1: /* Input: a q × q block matrix A with blocks A ij of size b × b. */ 2: /* Output: its BLR LDU factors LDU . */ 3: for k = 1 to q do 4:

Update:

5: R kk = A kk -k-1 j=1 L kj D jj U jk .
6:

for i = k + 1 to q do 7: R ik = A ik -k-1 j=1 L ij D jj U jk and R ki = A ki -k-1 j=1 L kj D jj U ji . 8:
end for 9:

Compress:

10:

for i = k + 1 to q do 11:
Compute low-rank approximations T ik ≈ R ik and T ki ≈ R ki .

12:

end for 13:

Factor:

14:

Compute the LU factorization L kk D kk U kk = R kk .

15:

for i = k + 1 to q do 16: Solve L ik D kk U kk = T ik for L ik and L kk D kk U ki = T ki for U ki . 17:
end for 18: end for

The error analysis presented in this section has a double purpose. First, it proves that the numerical stability of the uniform precision BLR LU factorization (proven by [START_REF] Higham | Solving Block Low-Rank linear systems by LU factorization is numerically stable[END_REF]) is preserved in mixed precision arithmetic. Second, for each operation required by Algorithm 2.1, it determines which level of accuracy is needed to maintain the overall error of order ε. Our analysis therefore guides us towards an implementation of mixed precision BLR LU that is both robust and efficient.

One technical difficulty of this analysis is the handling of the scaling factors hidden inside the U factor. To make these details more apparent, we analyze instead the LDU factorization, where L and U are unitriangular matrices (with ones on the diagonal). For the sake of readability, we will not always keep track of lower order error terms; we use the notations ≈ and ≲ to indicate when these terms (of order at most u p ε) have been dropped.

We first analyze each kernel separately. Algorithm 2.1 requires computing products of the form L ij D jj U jk (on line 7), where L ij and U jk may be either uniform precision fullrank blocks or mixed precision low-rank blocks (analysis of sections 2.4.1 and 2.4.2). We also analyze in section 2.4.3 the solution of a triangular system L kk D kk U ki = T ki (needed on line 16), where the right-hand side T ki is a mixed precision low-rank matrix. Finally, we combine the analysis of these kernels to obtain a backward error bound on the mixed precision BLR LU factorization in section 2.4.4.

Low-rank matrix times full-rank matrix

Let us begin with the computation of a product P = BC, where C is a full-rank matrix and B = XY T is a mixed precision low-rank matrix partitioned into p groups

B ℓ = X ℓ Y T ℓ satisfying ∥B ℓ ∥ ≤ εβ/u ℓ for ℓ > 1,
and where the output P is needed under full-rank form.

In which precision should we compute P = BC? The natural approach is to compute each product P ℓ = B ℓ C in precision u ℓ . Then, using Higham and Mary (2021, Lemma 3.2), the computed P ℓ satisfies

P ℓ = B ℓ C + ∆P ℓ , ∥∆P ℓ ∥ ≤ γ (ℓ) c ℓ ∥B ℓ ∥∥C∥ (2.36) with c ℓ = b + r 3/2
ℓ . For ℓ > 1, we thus obtain

P ℓ = B ℓ C + ∆P ℓ , ∥∆P ℓ ∥ ≲ c ℓ εβ∥C∥, (2.37) since γ (ℓ) c ℓ /u ℓ = c ℓ (1 + γ (ℓ) c ℓ ) ≈ c ℓ .
This shows that the partial product P ℓ associated with the part of B stored in precision u ℓ can itself be computed in precision u ℓ , since the introduced error remains of order ε.

The next question is what precision should be used to combine the partial results into P = p ℓ=1 P ℓ . Since for ℓ > 1 ∥P ℓ ∥ ≤ ε/u ℓ β∥C∥, it is easy to see that P i + P j must be computed in precision min(u i , u j ) = u min(i,j) . Knowing this, in order to maximize the performance gains associated with the use of lower precisions, the best approach is to compute p ℓ=1 P ℓ in reverse order. The approach to compute P suggested by our analysis is summarized in Algorithm 2.2.

Algorithm 2.2 Mixed precision low-rank matrix times full-rank matrix. Compute P ℓ = B ℓ C in precision u ℓ .

6:

Update P ← P + P ℓ in precision u ℓ . 7: end for With this algorithm, each component of P ℓ is involved in exactly min(ℓ, p-1) additions, one in each precision u 1 , . . . , u min(ℓ,p-1) . Therefore, the computed P satisfies:

P = p ℓ=1 P ℓ • (J + Θ ℓ ), |Θ ℓ | ≲ u ℓ , (2.38)
where J is the matrix of ones, • denotes the Hadamard product, and the inequality |Θ ℓ | ≲ u ℓ holds componentwise. Overall, we obtain

P = p ℓ=1 (B ℓ C + ∆P ℓ ) • (J + Θ ℓ ) (2.39) = BC + p ℓ=1 B ℓ C • Θ ℓ + ∆P ℓ + ∆P ℓ • Θ ℓ (2.40) = BC + ∆P, ∥∆P ∥ ≲ (c 1 + 1)u 1 ∥B (1) ∥ + p ℓ=2 (c ℓ + 1)εβ ∥C∥ (2.41) = BC + ∆P, ∥∆P ∥ ≲ (pb + r 3/2 + p) max(u 1 ∥B (1) ∥, εβ)∥C∥.
(2.42)

We summarize this analysis in the next theorem.

Theorem 2.3. Let B = p ℓ=1 B ℓ ∈ R b×b be a mixed precision low-rank matrix of rank r such that ∥B ℓ ∥ ≤ εβ/u ℓ for ℓ > 1, and let C ∈ R b×b . If P = BC is computed as described by Algorithm 2.2, then the computed P satisfies

∥ P -BC∥ ≲ c max(u 1 ∥B∥, εβ)∥C∥,
(2.43)

with c = pb + r 3/2 + p.
Theorem 2.3 shows that we can perform many of the flops in Algorithm 2.2 in lower precisions and still maintain an error of order ε. We now prove similar results for the other kernels.

Low-rank matrix times low-rank matrix

Next we analyze the product P = BDC of two mixed precision low-rank matrices

B = X B Y T B and C = Y C X T C
, where we also incorporate a diagonal scaling matrix D, which will be useful for the LU factorization analysis of section 2.4.4.

The product P , which is needed in full-rank form, can be computed in the following three steps:

1. Compute the inner product M = (Y B ) T DY C . 2. Compute the middle product W = X B M (or W = M X T C ).
3. Compute the outer product P = W X T C (or P = X B W ).

A trivial extension of Higham and Mary (2021, Lem. 3.2) to incorporate D shows that if P is computed in uniform precision u, the computed P satisfies

P = BDC + ∆P, ∥∆P ∥ ≲ cu∥B∥∥D∥∥C∥, ( 2
|∆M ℓm | ≲ bu M ℓm |Y Bℓ | T |D||Y Cm |.
(2.45)

By taking norms, we obtain

∥∆M ℓm ∥ ≲ bu M ℓm min(α 1 , α 2 , α 3 ), (2.46)
where

α 1 = ∥B ℓ D∥∥C m ∥, α 1 ≤ εβ B ∥C m ∥/u ℓ if ℓ > 1,
(2.47)

α 2 = ∥B ℓ ∥∥DC m ∥, α 2 ≤ ε∥B ℓ ∥β C /u m if m > 1,
(2.48)

α 3 = ∥B ℓ D∥∥D -1 ∥∥DC m ∥, α 3 ≤ ε 2 β B β C ∥D -1 ∥/(u ℓ u m ) if ℓ, m > 1.
(2.49)

From this we can deduce the optimal choices of precisions u M ℓm that still guarantee an error of order ε.

• If ℓ = m = 1, in general we must take u M 11 = u 1 since ∥B 1 ∥ and ∥C 1 ∥ are not bounded in terms of ε. We obtain

∥∆M 11 ∥ ≲ bu 1 ∥B 1 ∥∥D∥∥C 1 ∥.
(2.50)

• If ℓ = 1 and m > 1, α 2 ≤ ε/u m ∥B∥β C
, and so taking u M 1m = u m yields an error of order ε:

∥∆M 1m ∥ ≲ bε∥B 1 ∥β C .
(2.51)

Similarly, we can take u M ℓ1 = u ℓ and obtain

∥∆M ℓ1 ∥ ≲ bεβ B ∥C 1 ∥. (2.52) • If ℓ, m > 1, we can safely take u M ℓm = max(u ℓ , u m ) = u max(ℓ,m)
. Indeed, if ℓ ≥ m we can use (2.47) and if ℓ < m we can use (2.48), and so, in any case, we have

∥∆M ℓm ∥ ≲ bε max(β B ∥C 1 ∥, ∥B 1 ∥β C ).
(2.53)

Combining (2.50), (2.51), (2.52), and (2.53), we obtain for ℓ, m ≥ 1

∥∆M ℓm ∥ ≲ b max(εβ B ∥C 1 ∥, ε∥B 1 ∥β C , u 1 ∥B 1 ∥∥D∥∥C 1 ∥).
(2.54)

In summary, for any value of ℓ and m, we can compute the product between the part of B stored in precision u ℓ and the part of C stored in precision u m in the lower of the two precisions. This is a crucial observation that allows us to maximize the use of lower precision.

Moreover, in some cases we may actually take u M ℓm > max(u ℓ , u m ) because of (2.49). To see why, let us take an example where ∥D -1 ∥, β B , β C , and ∥A∥ are all approximately equal to 1. In this case, for ℓ, m > 1, (2.46) reduces to

∥∆M ℓm ∥ ≲ bu M ℓm ε 2 /(u ℓ u m ) (2.55)
and so the requirement to obtain an error of order ε is

bu M ℓm ε ≤ u ℓ u m , (2.56)
which thus depends not only on u ℓ and u m , but also on ε. If ε is small enough, (2.56) may be satisfied even for u M ℓm > max(u ℓ , u m ). For example, assume we have three precisions

u 1 = u d = 2 -53 , u 2 = u s = 2 -24
, and u 3 = u h = 2 -8 . We may consider that we also have access to a 0-bit precision format, u 4 = 1, whose only representable value is 0. Then, ignoring the constant b in (2.56):

• The condition u M 22 ε ≤ u 2 2 is satisfied for u M 22 = u 3 if ε ≤ u 2 s /u h ≈ 9 × 10 -13
. Thus, if ε is small enough, M 22 need only be computed in half precision.

• The condition u M 23 ε ≤ u 2 u 3 is satisfied for u M 23 = 1 if ε ≤ u s u h ≈ 2 × 10 -10
. The same holds for u M 32 . Thus, if ε is small enough, the computation of M 23 and M 32 may be skipped altogether. Indeed, by replacing the result of the operation by 0, we would obtain a relative error of 1, which is affordable in this particular case.

• Finally, the condition

u M 33 ε ≤ u 2 3 is satisfied for u M 33 = 1 if ε ≤ u 2 h ≈ 2 × 10 -5
. Again, the computation of M 33 may be skipped in this case.

Going back to the general case, the precise requirement on u M ℓm depends on u ℓ , u m , ε, β B , β C , and ∥D -1 ∥. For global compression (β B = β C = ∥A∥), we obtain

∥∆M ℓm ∥ ≲ bu M ℓm ε 2 ∥A∥ 2 ∥D -1 ∥/(u ℓ u m ).
(2.57)

Middle product W = X B M (or W = M X T C )
We analyze the product m should be summed in reverse order and in increasing precision, since

W = X B M , the case of W = M X T C being analogous. Let W m = p ℓ=1 X Bℓ M ℓm , for m = 1 : p. Assume the product W (ℓ) m = X Bℓ M ℓm is computed in precision u W ℓm , then the computed W (ℓ) m satisfies W (ℓ) m = X Bℓ M ℓm + ∆W (ℓ) m , (2.58) ∥∆W (ℓ) m ∥ ≲ r ℓ u W ℓm ∥X Bℓ ∥∥ M ℓm ∥ ≲ r 3/2 ℓ u W ℓm ∥B ℓ DC m ∥. ( 2
W (ℓ) m +W (ℓ+1) m must be computed in precision u W ℓm . Overall, with u W ℓm = u M ℓm = max(u ℓ , u m ), the computed W m satisfies W m = p ℓ=1 W (ℓ) m • (J + Θ ℓ ), |Θ ℓ | ≲ u ℓ , (2.60) = p ℓ=1 (X Bℓ M ℓm + X Bℓ ∆M ℓm + ∆W (ℓ) m ) • (J + Θ ℓ ), (2.61) = p ℓ=1 W (ℓ) m + ∆ W (ℓ) m = W m + ∆W m , (2.62) with ∥∆ W (ℓ) m ∥ ≲ (b + r 3/2 ℓ + 1) max(εβ B ∥C∥, ε∥B∥β C , u 1 ∥B∥∥D∥∥C∥) (2.63)
and so

∥∆W m ∥ ≲ (pb + r 3/2 + p) max(εβ B ∥C∥, ε∥B∥β C , u 1 ∥B∥∥D∥∥C∥).
(2.64)

Outer product P = W X T C (or P = X B W )
It remains to analyze the final product P = W X T C (or P = X B W , which is analogous). Let P m = W m X T Cm be computed in precision u P m . Then the computed P m satisfies

P m = W m X T Cm + ∆P m , (2.65) ∥∆P m ∥ ≲ r m u P m ∥ W m ∥∥X Cm ∥ ≤ r 3/2 m u P m ∥ W m ∥ ≲ r 3/2 m u P m p ℓ=1 ∥B ℓ DC m ∥.
(2.66)

Since p ℓ=1 ∥B ℓ DC m ∥ is at least as large as ∥B 1 DC m ∥, by (2.48) we must take u P m = u m . Then, (2.66) becomes

∥∆P m ∥ ≲ r 3/2 m max(ε∥B∥β C , u 1 ∥B∥∥D∥∥C∥).
(2.67)

Finally, as previously for W m , we sum P m over m in reverse order and in increasing precision, to obtain a computed P satisfying

P = p m=1 P m • (J + Θ m ), |Θ m | ≲ u m , (2.68) = p m=1 ( W m (X Cm ) T + ∆P m ) • (J + Θ m ), (2.69) = p m=1 P m + (∆W m X T Cm + ∆P m ) • (J + Θ m ),
(2.70)

= P + ∆P, (2.71) with ∥∆P ∥ ≲ (p 2 b + (p + 1)r 3/2 + p 2 + p) max(εβ B ∥C∥, ε∥B∥β C , u 1 ∥B∥∥D∥∥C∥). (2.72)
This concludes the analysis of the product P . We summarize the approach suggested by this analysis in Algorithm 2.3, for which the following theorem holds.

Algorithm 2.3 Mixed precision low-rank matrix times mixed precision low-rank matrix. Initialize W m to zero.

6:

for ℓ = p to 1 do

7:

Compute M ℓm = Y Bℓ DY T Cm in precision max(u ℓ , u m ).

8:

Compute

W (ℓ) m = X Bℓ M ℓm in precision max(u ℓ , u m ). 9: Update W m ← W m + W (ℓ)
m in precision max(u ℓ , u m ).

10:

end for 11:

Compute P m = W m X T Cm in precision u m .

12:

Update P ← P + P m in precision u m . 13: end for Theorem 2.4 (Low-rank times low-rank). Let B = p ℓ=1 B ℓ and C = p m=1 C m be two mixed precision low-rank matrices satisfying

∥B ℓ D∥ ≤ εβ B /u ℓ for ℓ > 1,
(2.73)

∥DC m ∥ ≤ εβ C /u m for m > 1, (2.74)
and D a diagonal matrix, and let P = BDC be computed as described in Algorithm 2.3. Then, the computed P satisfies

P = BDC + ∆P, ∥∆P ∥ ≲ c max εβ B ∥C∥, ε∥B∥β C , u 1 ∥B∥∥D∥∥C∥ , (2.75) with c = p 2 b + (p + 1)r 3/2 + p 2 + p.
2.4.3 Triangular system with low-rank right-hand side

The last kernel that we need to analyze is the solution of a triangular system LDZ = B, where L ∈ R b×b is lower triangular, D is diagonal, and the right-hand side B = Y X T is a mixed precision low-rank matrix (used in Algorithm 2.1, line 16). We analyze the kernel for a lower triangular matrix L, the upper triangular case (ZDU = B) being analogous. For this kernel, the output (the solution Z) is needed under low-rank form.

In the uniform precision case, if the system LDZ = B is solved in uniform precision u, the computed solution Z satisfies (Higham and Mary, 2021, Lemma 3.5)

LD Z = B + ∆B, ∥∆B∥ ≲ bu∥L∥∥D∥∥ Z∥.
(2.76)

Let B ℓ = Y ℓ X T ℓ be the part of B that is stored in precision u ℓ , satisfying ∥B ℓ ∥ = ∥Y ℓ ∥ ≤ εβ/u ℓ for ℓ > 1.
Then, the natural approach to solve LDZ = B in mixed precision is to solve each system LDV ℓ = Y ℓ in precision u ℓ and to define Z ℓ = V ℓ X T ℓ , which yields the mixed precision low-rank solution Z = p ℓ=1 Z ℓ . However, a traditional normwise analysis based on (2.76) does not provide a satisfactory bound here: if we apply (2.76) to LDV ℓ = Y ℓ and use V ℓ ≈ D -1 L -1 Y ℓ , we obtain the bound

LD V ℓ = Y ℓ + ∆Y ℓ , ∥∆Y ℓ ∥ ≲ bεβκ(L)κ(D).
(2.77)

This bound is very weak due to the presence of the normwise condition numbers κ(L)κ(D) = ∥L∥∥L -1 ∥∥D∥∥D -1 ∥.

A stronger bound can be obtained by using a componentwise analysis:

LD V ℓ = Y ℓ + ∆Y ℓ , |∆Y ℓ | ≲ bu ℓ |L||D|| V ℓ |. (2.78) Replacing V ℓ by D -1 L -1 (Y ℓ + ∆Y ℓ ) in the bound on ∆Y ℓ yields |∆Y ℓ | ≲ bu ℓ |L||D||D -1 L -1 Y ℓ | ≤ bu ℓ |L||L -1 ||Y ℓ |.
(2.79)

We can now take norms, obtaining for ℓ > 1

∥∆Y ℓ ∥ ≲ bu ℓ cond(L, Y ℓ )∥Y ℓ ∥ ≤ bεβ cond(L, Y ℓ ) (2.80)
where cond(L, Y ℓ ) is the condition number introduced by Skeel (1979) (Higham, 2002, Eq. (7.13)):

cond(L, Y ℓ ) = ∥|L||L -1 ||Y ℓ |∥ ∥Y ℓ ∥ . (2.81)
Multiplying both sides of (2.78) by X T ℓ on the right yields

LD Z ℓ = B ℓ + ∆Y ℓ X T ℓ .
(2.82) Summing (2.82) over ℓ, we obtain

LD Z = B + p ℓ=1 ∆Y ℓ X T ℓ = B + ∆B, (2.83) ∥∆B∥ ≲ bu 1 ∥L∥∥D∥∥ Z 1 ∥ + pbεβ cond(L), (2.84)
where cond(L) = ∥|L -1 ||L|∥.

The use of intermediate componentwise bounds presents two advantages. First, we obtain a bound with cond(L), which is in general smaller than κ(L) (see Higham (2002, p. 123) for a discussion on the difference between these two quantities). Second and more importantly, we have dropped the matrix D from the term proportional to εβ, which shows that this term is invariant under scaling, and which represents a significant improvement since κ(D) ≈ κ(A) can be arbitrarily large. Importantly, in the case of an LDU factorization with partial pivoting, both L and U are well conditioned, and cond(L) is in practice a small constant. Therefore, in the context of Algorithm 2.1, the mixed precision triangular solution analyzed here is backward stable. However, for a general system LDZ = B, bound (2.84) does not guarantee backward stability, and indeed some examples can be built where the use of mixed precision arithmetic in the solution of the system leads to a large increase of the backward error (we note however that such examples are very hard to find and we were only able to construct one using direct search optimization [START_REF] Higham | Optimization by direct search in matrix computations[END_REF]).

We summarize the proposed approach to solve LDZ = B in Algorithm 2.4 and its error analysis in the following theorem. Solve the triangular system LDV ℓ = Y ℓ in precision u ℓ .

5:

Define Z ℓ = V ℓ X T ℓ (no computation performed: output is low-rank). 6: end for Theorem 2.5. Let L ∈ R b×b be a lower triangular full-rank matrix and let B = p ℓ=1 B ℓ be a mixed precision low-rank matrix satisfying ∥B ℓ ∥ ≤ βε/u ℓ . If the system LDZ = B is solved by Algorithm 2.4, the computed solution Z satisfies

LD Z = B + ∆B, ∥∆B∥ ≲ pbεβ cond(L) + bu 1 ∥L∥∥D∥∥ Z∥.
(2.85)

Putting everything together: error analysis of mixed precision BLR LU factorization

Now that we have analyzed all the kernels of Algorithm 2.1, we are ready to prove the backward stability of the BLR LU factorization in mixed precision arithmetic. We define

λ 1 = max k=1 : q max ∥L -1 kk ∥, ∥U -1 kk ∥, cond(L kk ), cond(U kk ) .
(2.86)

If partial pivoting is performed, λ 1 is almost always small in practice and of order a constant (Higham, 2002, Chapter 8), even though in theory it can only be bounded by Higham, 2002, Lemma 8.6). We also define

2 b -1 (
λ 2 = max i≥j max ∥L ij ∥, ∥U ji ∥ . (2.87) If partial pivoting is performed, λ 2 ≤ b.
Let us bound the error incurred in the computation of some (i, k) block of the L factor, the U factor analysis being similar. For i < k, L ik is obtained by solving

L ik D kk U kk = T ik , (2.88)
where T ik is the compressed form of

R ik = A ik - k-1 j=1 L ij D jj U jk , (2.89)
where L ij and U jk are the LU factors computed at the previous steps, and are represented as mixed precision low-rank matrices, and the product L ij D jj U jk is computed with Algorithm 2.3. Note that if one of L ij or U jk is a full-rank matrix, the analysis is similar and relies on Theorem 2.3; if both are full-rank, the computation is done in uniform precision u 1 ≪ ε and introduces an error term bu 1 ∥ L ij ∥∥D jj ∥∥ U jk ∥. The difficulty is that L ij and U jk are not directly the result of a compression, and so we cannot directly control the norms of

L (ℓ) ij and U (m)
jk , the parts of L ij and U jk stored in precision u ℓ and u m , respectively. Instead, they are given by

L (ℓ) ij ≈ T (ℓ) ij U -1 jj D -1 jj , (2.90) U (m) jk ≈ D -1 jj L -1 jj T (m) jk , (2.91)
where T ij and T jk have been compressed such that

∥T (ℓ) ij ∥ ≤ εβ ij /u ℓ for ℓ > 1, (2.92) ∥T (m) jk ∥ ≤ εβ jk /u m for m > 1.
(2.93)

Therefore, the norms of

L (ℓ) ij and U (m)
jk depend on β ij and β jk , respectively, but also on the scaling factors in D jj . However, one of the two D -1 jj in (2.90)-(2.91) is canceled by the D jj in (2.89) and ∥ L -1 jj ∥ and ∥ U -1 jj ∥ are both bounded by λ 1 . As a result, we can rewrite the product R ik,j = L ij D jj U jk as BDC, where

∥B ℓ D∥ ≲ λ 1 εβ ij /u ℓ for ℓ > 1,
(2.94)

∥DC m ∥ ≲ λ 1 εβ jk /u m for m > 1.
(2.95) By Theorem 2.4, the computed R ik,j satisfies

R ik,j = L ij D jj U jk + ∆R ik,j , (2.96) ∥∆R ik,j ∥ ≲ c max(λ 1 εβ ij ∥ U jk ∥, λ 1 εβ jk ∥ L ij ∥, u 1 ∥ L ij ∥∥D jj ∥∥ U jk ∥).
(2.97)

By (2.89), we obtain a computed R ik R ik = A ik • (J + Θ k ) - k-1 j=1 L ij D jj U jk + ∆R ik,j • (J + Θ j ), (2.98)
where

|Θ j | ≤ γ (1)
j J accounts for the errors in the additions of the products R ik,j to A ik . We thus obtain

R ik = A ik - k-1 j=1 L ij D jj U jk + ∆R ik ,
(2.99)

∥∆R ik ∥ ≲ ku 1 ∥A ik ∥ + k-1 j=1 max λ 1 λ 2 cε max(β ij , β jk ), (λ 2 2 c + j)u 1 ∥D jj ∥ . (2.100)
R ik is then compressed into T ik such that the part of T ik stored in precision u ℓ satisfies

T (ℓ) ik = R (ℓ) ik + E (ℓ) ik , ∥E (ℓ) ik ∥ ≤ εβ ik /u ℓ , (2.101)
and so overall, by Theorem 2.1,

T ik = R ik + E ik , ∥E ik ∥ ≲ (2p -1)εβ ik .
(2.102)

Finally, we solve (2.88) for L ik , and by Theorem 2.5, the computed L ik satisfies

L ik D kk U kk = T ik + F ik , ∥F ik ∥ ≲ pbλ 1 εβ ik + bu 1 λ 2 2 ∥D kk ∥.
(2.103)

Putting together (2.103), (2.102), (2.100), we obtain

L ik D kk U kk = A ik - k-1 j=1 L ij D jj U jk + ∆R ik + E ik + F ik , (2.104)
and so

A ik = k j=1 L ij D jj U jk + ∆A ik , (2.105) ∥∆A ik ∥ ≲ ku 1 ∥A ik ∥ + k j=1 max λ 1 λ 2 cε max(β ij , β jk ), (λ 2 2 c + j)u 1 ∥D jj ∥ . (2.106)
With the choice β ij = β jk = ∥A∥ for j = 1 : k, and since k ≤ q, we obtain

∥∆A ik ∥ ≲ λ 1 λ 2 cqε∥A∥ + qu 1 ∥A ik ∥ + q(λ 2 2 c + q)ρu 1 ∥A∥, (2.107)
where we have used ∥D jj ∥ ≤ ρ∥A∥, where ρ denotes the growth factor. This concludes the case i < k. Bounds analogous to (2.107) hold for i = k and i > k, and so overall we have

A = LD U + ∆A, ∥∆A∥ ≲ q 2 λ 1 λ 2 cε + (λ 2 2 c + q)ρu 1 ∥A∥. (2.108)
We summarize this analysis in the next theorem.

Theorem 2.6 (Mixed precision BLR LU factorization). Let A ∈ R n×n be a BLR matrix partitioned into q 2 blocks of order b. If the BLR LU factorization of A in p precisions described by Algorithms 2.1-2.4 runs to completion, the computed LU factors satisfy

A = LD U + ∆A, ∥∆A∥ ≲ q 2 λ 1 λ 2 cε + (λ 2 2 c + q)ρu 1 ∥A∥, (2.109)
where λ 1 and λ 2 are defined by (2.86)-(2.87), ρ is the growth factor, and c = p 2 b + (p + 1)r 3/2 + p.

Theorem 2.6 therefore proves the backward stability of the mixed precision BLR LU factorization: the computed LU factors give an exact LU decomposition of a perturbed matrix, where the norm of the perturbation ∥∆A∥ is of order ε. The precise value of the constants q 2 λ 1 λ 2 c and q 2 (λ 2 2 c + q)ρ in (2.109) is not of great importance but, as a check, we compare it against the uniform precision bound of Higham and Mary (2021, Thm. 4.3)

A = L U + ∆A, ∥∆A∥ ≲ qε∥A∥ + (b + 2r 3/2 + q)u 1 ∥ L∥∥ U ∥.
(2.110)

Since ∥ L∥∥ U ∥ ≲ n 2 ρ∥A∥ and, with partial pivoting, λ 2 ≤ b, we see that both (2.109) and (2.110) grow as O(n 2 (b + q)ρ).

After Theorem 2.6, not much additional effort is needed to prove the backward stability of the solution of linear systems Ax = v by mixed precision BLR LU factorization. We note that mixed precision arithmetic can also be used in the solution of the triangular systems with the LU factors, but in the interest of space, we omit these details.

Experimental results

Experimental setting

We have written a MATLAB code that implements a mixed precision variant of Algorithm 2.1 that uses Algorithms 2.2-2.4. Our implementation can use any number of arbitrary precisions, where the lower precisions are simulated using the chop function of Higham and Pranesh (2019). To compress the blocks, we use a mixed precision truncated QR decomposition with column pivoting-we omit a detailed description of this algorithm, which we plan to investigate more in depth in future work.

For our experiments, we use the matrices listed in Table 2.1. These matrices are all obtained as the Schur complement of larger sparse matrices (specifically, the root separators of their multifrontal factorization) arising in various applications: P64 comes from the discretization of a Poisson equation, perf009d comes from a structural mechanics problem from EDF (French electricity supplier), the others come from the SuiteSparse collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF].

To confirm experimentally the numerical stability of the algorithms, and to assess the impact of mixed precision arithmetic on their accuracy, we measure backward errors. Rather than measuring the backward error for the LU factorization, which is expensive to compute, we use the computed LU factors to solve a linear system Ax = v, where x is the vector of ones (and v is computed as Ax), and we use the computed solution x to measure the backward error

∥A x -v∥ ∥A∥∥ x∥ (2.111)
given by the Rigal-Gaches theorem (Higham (2002, Thm 7.1), [START_REF] Rigal | On the compatibility of a given solution with the data of a linear system[END_REF]).

To evaluate the potential of using mixed precision for the BLR LU factorization, we will focus on two performance metrics: the storage cost for the factorized matrix and the expected time cost of the factorization. Both these costs depend on the relative performance of each arithmetic. It is easy to measure the storage cost since, for each arithmetic, it is proportional to the number of bits used: an fp32 number requires half the storage of an fp64 one, and a bfloat16 number requires a quarter of the storage. Hence, we have storage cost = (#entries in fp64) + 0.5 × (#entries in fp32) + 0.25 × (#entries in bf16).

(2.112) The time cost will be estimated based on the number of flops performed in each arithmetic. This is a more complex issue because the relative speed of each arithmetic strongly depends on the hardware, the matrix, and several other factors. A practical high-performance implementation of the mixed precision BLR factorization is outside the scope of this chapter but, as a rough indicator, we make the assumption that the speed of each arithmetic is also proportional to the number of bits. Hence, we use the cost model expected time cost = (#flops in fp64) + 0.5 × (#flops in fp32) + 0.25 × (#flops in bf16).

(2.113)

Performance-accuracy tradeoff

The analytical error bounds obtained in section 2.4 show that the use of mixed precision arithmetic should only increase the backward error by a small constant. In this first experiment, we check experimentally (i) that the error increase is indeed small, and (ii) whether the expected time and storage gains obtained by the use of mixed precision justify this error increase, that is, whether the mixed precision variant achieves a better performance-accuracy tradeoff than the uniform precision variant. Indeed, since the mixed precision variant achieves a slightly larger error, to be completely fair, we should compare it against the uniform precision variant with a correspondingly larger ε.

To answer this question, we perform the following experiment in Figure 2.3: for a given matrix, taking several values of ε, we plot the gains in storage and expected time (2.112) and (2.113) as a function of the backward error (2.111). We compare three variants of the BLR factorization: the standard uniform precision variant run entirely with fp64 arithmetic, a two-precision variant using both fp64 and fp32, and a three-precision variant using bfloat16 as well. The figure shows that the two-precision variant achieves a much better performance-accuracy tradeoff than the uniform precision one, and that the threeprecision variant further improves this tradeoff. Indeed, using lower precisions slightly increases the error, but the experiment shows that this increase is largely compensated by the cost reductions. Indeed, the closer a variant is to the top left corner of the plots, the better its tradeoff is: for a given accuracy, it is less costly than the other variants, or, equivalently, for a given cost, it achieves an improved accuracy.

In light of this experiment, and to avoid hand-tuning ε for every variant and every matrix, in the remainder of our experiments we directly compare the variants with the same value of ε.

Results on real-life matrices

In this section we experiment on the Schur complement of real-life matrices listed in Table 2.1. Figure 2.4 compares the backward error (2.111) achieved by the BLR factorization for three values of ε: 10 -12 , 10 -9 , and 10 -6 . For ε = 10 -6 , we compare the uniform fp32 precision BLR factorization with the two-precision one using fp32 and bfloat16; for ε = 10 -9 and 10 -12 , we compare the uniform fp64 precision factorization with both a two-precision variant (using fp64 and fp32) and a three-precision one (also using bfloat16). The figure shows that the use of mixed precision arithmetic does not significantly impact the backward error, leading to an increase of at most an order of magnitude in the worst case (and very often much less than that).

Figure 2.5 shows the associated storage and time gains expected from the use of mixed precision arithmetic. For each matrix, each bar corresponds to a different value of ε. For ε = 10 -12 and 10 -9 , we focus on the gains achieved by the three-precision variant. The y-axis (height of the bars) indicates the number of entries (or number of flops) required by the mixed precision variant as a percentage of the double precision variant. For storage, this percentage would be at least 100% if we could ignore some minor numerical perturbations, and can be larger than that because of the difference between conditions (2.26) and (2.27). Indeed, as explained in section 2.2, there are some blocks that satisfy (2.27) but do not satisfy (2.26): that is, we allow the mixed precision variant to store more entries, because we expect this increase to pay off thanks to the use of lower precisions. The same property stands for the flops, but some operations are also skipped (see section 2.4.2): therefore, the percentage for the flops can also differ from 100%, either being larger or smaller. The colors breakdown in Figure 2.5 shows the proportion of entries that are stored in each precision, and the proportion of flops that are performed in each precision. For ε = 10 -6 , note that the mixed precision algorithm recovers the fact that we do not need any entries or flops in fp64 arithmetic, since ε > u s . For all matrices there is a significant fraction of the entries and flops that can be switched to lower precisions, even for ε = 10 -12 . This is a very positive result that confirms that BLR matrices are amenable to the use of mixed precision arithmetic and that the proposed mixed precision BLR representation can achieve very significant gains with respect to the one in double precision.

The number on top of each bar in Figure 2.5 indicates the resulting gains in storage and expected time achieved by the mixed precision variant compared with the double precision variant. The figure shows very significant reductions, of up to a factor 2.8× in storage and 3.5× in expected time for ε = 10 -9 .

Finally, in Figure 2.6 we perform a similar experiment as in Figure 2.5 for matrices of increasing size belonging to the same Poisson problem class. This experiment highlights an important and valuable property of the mixed precision BLR factorization: the storage and expected time gains increase with the problem size, as a larger and larger fraction of the entries and flops can be safely switched to lower precisions.

Conclusion

We have introduced a novel approach to exploit mixed precision arithmetic for lowrank approximations. Given a prescribed accuracy ε, we have proved in Theorem 2.1 that singular vectors associated with sufficiently small singular values can be stored in precisions with unit roundoff larger than ε while preserving an overall accuracy of order ε. This approach is not only applicable to low-rank matrices built with a singular value decompo- , ε = 10 -9 and ε = 10 -6 respectively. The y-axis shows the number of entries and flops required by the mixed precision variant with respect to the double precision variant (which can be slightly different from 100% because of the different conditions to represent a block under low-rank form (2.26) and (2.27)). The color breakdown gives the proportion of entries and flops in each precision. The number above each bar indicates the resulting factor of gain in storage (2.112) and expected time (2.113) achieved by the mixed precision variant compared with the double precision one.

sition, but also to many other low-rank decompositions, in particular rank-revealing QR (Lemma 2.2).

We have applied this approach to block low-rank (BLR) matrices, for which this new mixed precision low-rank approximation presents a high potential. We have adapted the existing uniform precision BLR LU factorization algorithm (Algorithm 2.1) to exploit the mixed precision representation of the blocks. We carried out the rounding error analysis of this new algorithm and obtained two key results. First, we proved in Theorem 2.6 that the use of mixed precision arithmetic does not compromise the numerical stability of BLR LU factorization recently proven by [START_REF] Higham | Solving Block Low-Rank linear systems by LU factorization is numerically stable[END_REF]. Second, our analysis helps us determine what accuracy is needed for each floating-point operation. The resulting mixed precision BLR algorithms are summarized in Algorithms 2.2, 2.3, and 2.4. Interestingly enough, our analysis suggests that some operations may even be Chapter 3

The multifrontal method in mixed precision

In this chapter we want to improve the performance of the multifrontal method. We adapt it in order to take advantage of mixed precision, based on the algorithms and concepts presented in chapter 2. First, we will consider the use of mixed-precision BLR (MPBLR) as a storage format only. Its aim is to further reduce the cost of storing a BLR matrix. We include it in the multifrontal method in order to reduce its memory consumption. In a second section, we perform computations in mixed precision exploiting the MPBLR compression. By doing this we hope to accelerate the different steps of the multifrontal method.

We implemented our algorithms in the MUMPS multifrontal solver. Some of the solver functionalities presented in this chapter are not public yet, and may be added in a future version. All the experiments were performed on the Olympe supercomputer of the CALMIP center (project P0989). Unless specified otherwise, experiments are performed on 2 computational nodes, using 4 MPI processes, and 18 OpenMP threads per process.

Mixed precision aiming for storage and communication reductions

Using custom precision formats

We first present an alternative MPBLR format aiming at optimizing the storage cost of a matrix. In the next subsection we will integrate it in the multifrontal method.

As seen in chapter 1, a floating-point number is a number of the form

f = (-1) s × t i=0 m i 2 -i × 2 (e-bias)
where s is the sign bit, the integer e is the exponent, and m = m 0 m 1 • • • m t is the mantissa. m 0 = 1 is an implicit bit, and is not actually stored, according to the description of the binary interchange floating-point formats in the IEEE 754 standard (see IEEE Computer Society, 2008). In this case, the unit roundoff is u = 2 -t . If the representation of the exponent stays the same, a format with fewer mantissa bits, less costly in terms of storage, will have a greater unit roundoff. As a consequence, in order to maximize the storage gains, we will try to use a unit roundoff as big as possible.

It is worth mentioning that, according to the description of the interchange floatingpoint formats, the bits are stored in the following order: s, e, m. In particular, the least significant bits of m are stored last.

In chapter 2 we saw that, given a set of formats whose unit roundoffs are u 1 • • • u p , we can find a criterion indicating in which precision format we should store the columns X k and Y k of a low-rank approximation XY T ≈ B. Such a criterion may vary depending on the type of low-rank approximation. For example, we may choose the largest value of u i (i.e., the lowest precision) satisfying ∥X k Y T k ∥ ≤ ε/u i . This condition ensures that the error is bounded by a small multiple of ε, as shown in section 2.3.2.

Another way to consider this is that column k should be stored in any precision format u above a certain target accuracy u

(k) target = ε/∥X k Y T k ∥.
Thus it is understandable that, if we add new precision formats, we will sometimes be able to switch the column to a lower precision format, less costly, while the resulting error bound will be closer to ε (but still below it). The error bound of each column would be equal to ε if we had a large number of available precision formats, enough so that u would vary almost continuously. Such a set of precision formats would minimize the cost for storing an MPBLR matrix (under the constraint of respecting a certain error bound).

We did not try to implement such a continuum of precision formats. However, we understand that having only 2 or 3 precision formats like in chapter 2 may not be entirely satisfactory: we wish to add a few more precision formats to the set, in order to make a trade-off between these two extremes. In this section, MPBLR compression is only used for storage: therefore, in order to define a new precision format, one only needs functions doing the conversions between the working precision and the new format. Such conversions may be implemented in software, which means that it is possible to add virtually any format.

We choose to obtain new custom formats by keeping only the t ′ most significant bits of the mantissa. Conveniently, the discarded bits are all located at the right-hand side of the format, making it easier to truncate them. The remaining bits correspond to a new floating-point format, with the same exponent range but fewer mantissa bits and a higher roundoff u ′ = 2 -t ′ . The rounding mode of this conversion is towards zero.

In order to obtain an implementation that is easy and somewhat efficient, we choose to only allow byte manipulations, instead of bit manipulations: the new formats are obtained by removing the last bytes of the mantissa. The conversion back to high precision consists in padding the missing bits with zeros. As a consequence of these choices, we obtain new formats on 56, 48, 40 and 24 bits respectively, that we name fp56, fp48, fp40 and fp24. We also add a 16-bits format, which turns out to correspond to the definition of bloat16.

We also reuse the standard formats fp64 and fp32. As a result, we obtain a set of seven precision formats, which is described in Table 3.1 As a matter of fact, the very same set of 7 floating-point formats has already been proposed several times in the literature, in particular by [START_REF] Anderson | Efficient multibyte floating point data formats using vectorization[END_REF], [START_REF] Mukunoki | Reduced-precision floating-point formats on gpus for high performance and energy efficient computation[END_REF] and [START_REF] Graillat | Adaptive precision sparse matrix-vector product and its application to Krylov solvers[END_REF].

An example of such a custom floating-point format would be fp24. Its number of exponent bits is the same as fp32: indeed, we deem the range of fp32 to be sufficient, as this hypothesis was already made several times in this thesis. Therefore, when converting a number from fp64 to fp24, it has to be converted to fp32 first. The conversion from fp32 to fp24, shown in Figure 3.1, is performed by discarding the last byte, which corresponds to the 8 least significant bits. The conversion back to fp32 is obtained by padding the last byte with zero. We made a first implementation of these byte manipulations in Fortran.

An interesting perspective would be to increase the number of formats even more, by truncating the mantissa at the bit level instead of the byte level. The sizes of the available formats would take all values between 10 bits (mantissa of size 1), and 64 bits, thus approximating a continuum of floating-point formats. The gain would be especially noticeable for very low precision formats. For example, we may want to convert the last non-discarded columns from bfloat16 to "fp10". For these columns, the gain of switching the truncation operation from the byte level to the bit level would be the most beneficial, reaching an additional storage reduction of 6 16 ≈ 38%. Handling such formats would probably be complex and costly due to data misalignment, but it might be useful when wanting to obtain a near-optimal MPBLR compression, with no regard for the overhead.

When handling very low precision formats, the cost for storing the exponent bits begins to matter. A good illustration of this fact is the format fp10, whose exponent bits occupy 80% of the storage space. For bfloat16, this proportion is already 50%. However, we notice that the range of the exponents in a given format is rather narrow, in particular for the lowest precision formats. This is because we use low precision for the columns having the smallest coefficients: all coefficients of a column stored in precision i tend to be smaller than ε/u i (though this may depend on the kind of low-rank approximation being used). Therefore, another perspective would be to perform a lossless compression of the exponent. For example, we could use a limited number of bits to represent the exponents, by removing the leading bits whose values are zero. In order to know the number of bits to be discarded, we would need to check beforehand the maximal value of the exponent in a column.

Multifrontal method with mixed-precision storage

We want to take advantage of the storage-focused MPBLR format described previously. The main idea is the following: whenever the LU factors are not being used in a computation, we would like them to be stored in mixed precision. Indeed, the storage of the LU factors contributes to a large part of the memory peak of the multifrontal method in practice. However, they are not needed between the end of the factorization of a front and the triangular solution (last step of the multifrontal method).

The main modifications added to Algorithm 1.4 are the following:

• While computing the low-rank compression, we also compute the ranks of the mixedprecision parts r 1 , • • • , r p .

• At the end of the factorization of a front, the low-rank blocks are copied from uniform precision to the mixed-precision low-rank storage format. The original blocks in uniform precision are not needed anymore and may be discarded.

• Before it is used in the triangular solve step, each low-rank block is temporarily converted from mixed precision to uniform precision. The computations remain in uniform precision.

Block-admissibility conditions

In a BLR compression, we use a condition for choosing whether a block should be lowrank or full-rank: this is the block-admissibility condition, already mentioned in section 2.2 (see Equation 2.26). This choice may be based on the "distance" between the groups of variables corresponding to the rows and the columns. However there is also a more pragmatic alternative, which is used in MUMPS: a block is considered as low-rank if and only if the compression leads to a storage reduction. That is to say, the low-rank admissibility condition for a block of size (m, n) is:

r(m + n) ≤ mn. (3.1)
Conveniently, this condition for reducing the storage is also the condition for reducing the number of operations.

However, the situation is a bit different when using BLR compression for computations and MPBLR for storing the factors. The condition for reducing the number of operations stays the same, whereas the condition for reducing storage becomes the following:

(m + n) p k=1 c k r k ≤ mn, (3.2)
where c k corresponds to the cost of storing a floating-point number in precision u k instead of u 1 . For example, if we use two precisions fp64 and fp32, then Equation 3.2 takes the form (m + n)(r 1 + 0.5r 2 ) ≤ mn.

However trivial this aspect may seem, it has proven to be important in practice, as a noticeable fraction of the blocks satisfy the mixed-precision condition (Equation 3.2), but not the uniform-precision condition (Equation 3.1). That is to say, they would be full-rank in BLR but are low-rank in MPBLR.

On the other hand, switching to the block-admissibility condition for mixed precision (Equation 3.2) requires additional operations, which should lead to an overhead. Indeed, some blocks that were considered as full-rank are now compressed, and this compression requires additional operations. Moreover, the maximal values for the ranks have been increased: therefore, more operations may be needed before knowing whether a block is low-rank or not, at least if this choice is done by trying to compress the block. Finally, when using the block-admissibility suited to mixed precision, some blocks are not suited for computations: the ranks are too high, which is not optimal. They should have been considered as full-rank in order to minimize the number operations, as shown in Equation 3.1. By forcefully considering them as low-rank, we increase the number of operations, which should lead once again to an overhead. This is what is done in the current implementation in MUMPS, but it could still be improved in the future.

Implementation in MUMPS

We implemented our algorithm within the sparse solver MUMPS: we can now choose to use an MPBLR compression as a storage format. We added several options in order to customize its use. In particular:

-mpblr_thld: This variable allows the user to modify the thresholds for switching between two precision formats. In chapter 2, we had considered that such a threshold would be ε/u i for precision i. We now use a modified threshold, divided by mpblr_thld. Unless otherwise specified, we set mpblr_thld = 10.

As a result, the theoretical error bounds obtained when converting a low-rank block to mixed precision (with p = 2 formats) is roughly (1+2(p-1)/mpblr_thld)ε = 1.2ε, instead of 3ε obtained in Equation 2.22. This setting was deemed to be conservative enough with regard to the uniform precision, whose error bound is ε: in most applications the numerical error introduced by the use of mixed precision should not be significant. If it is, then an easy correction would be to increase the value of mpblr_thld.

-custom_prec: Allows the user to choose between using only the precision formats supported by the hardware, or a set of custom precision formats. In the current state of our implementation, the first choice results in using fp64 and fp32. Otherwise, we use a subset of our 7 precision formats described in Table 3.1. We should not need any format whose precision is higher than the working precision. Therefore, if MUMPS is compiled in fp32, we will only use 3 precision formats: fp32, fp24 and bf16.

-mixed_adm: Allows to switch between the usual block-admissibility condition suited to uniform precision (see Equation 3.1) and a condition minimizing the storage in mixed precision (see Equation 3.2), at the cost of a few additional operations.

Storage gains and time overhead

In our implementation, we handle two BLR formats at once: one that is used for computations in the working precision, and an MPBLR format that is used for storage. Performing data copies between the two formats, as well as the allocations/deallocations, may cause some overhead. Moreover, handling custom precision formats may have a relatively high cost, due to the fact that all conversions are implemented at the software level instead of the hardware level. Finally, we may want to use the appropriate blockadmissibility condition: this choice results in additional storage reduction, at the cost of an extra overhead. Therefore, choosing of an MPBLR variant results in a trade-off between memory and computation time: despite reducing the memory needed for the factorization even further, the use of more complex storage formats should induce an overhead.

In Table 3.2 we compare the execution time of different steps of MUMPS, as well as the memory peak, when using different levels of compression. The full-rank variant (FR) has no compression, contrary to BLR. MPBLR(2) and MPBLR(7) use sets of 2 and 7 precision formats respectively, and only differ from BLR by adding data conversions (i.e., they keep the same block-admissibility condition as in double precision). We also compare those two variants to their respective optimizations using the block-admissibility conditions suited to mixed precision: MPBLR(2) * and MPBLR(7) * .

As expected, the size of the factors is reduced when increasing the number of precision formats: we obtained reductions of the factor size of up to 24% by using two hardwaresupported precision formats, and up to 32% when using software-implemented custom precision formats. Such a compression of the factors induces a reduction of the memory peak. On the other hand, we confirm that increasing the number of precision formats generally causes an overhead, albeit small. When adding mixed-precision admissibility, the compression becomes a bit more efficient, obtaining storage reductions of up to 27% 38% with sets of 2 and 7 precisions respectively. In general, the time spent in factorization and in solve tends to increase again. In later experiments using MPBLR compressions, we will use this mixed-precision admissibility unless specified otherwise. 7) use sets of 2 and 7 precisions respectively, with the block-admissibility for uniform precision (Equation 3.1). The variants MPBLR(2) * and MPBLR(7) * differ by the use of the block-admissibility condition suited for mixed precision (Equation 3.2). For matrix Geoazur160, the working precision is fp32 and any precision format higher than this will not be used. Therefore, it makes no sense to use variant MPBLR(2), equivalent to uniform precision. When using variant MPBLR(7), only the 3 lowest precision formats of the set are used. †: estimation OOM : out of memory Backward error: ∥A x-b∥∞

∥A∥∞∥ x∥∞

We observed that using a BLR storage format different from the working format causes an overhead with our implementation. However, such overheads may be reduced, or may not even happen at all with an implementation focused on performance. In fact, the data movements have been reduced by doing this: at the very least, a memory-bound algorithm such as the BLR triangular solve should be able to turn such a reduction of the data movements into actual time gains, under the assumption that the conversion functions are sufficiently well tuned (which is not the case for our naive implementation, targeting storage gains only). Indeed, as explained in section 1.4.6 , we could hope to obtain, for a memory-bound algorithm such as this one, a speedup factor equal to the factor of storage gain from mixed precision at most.

A side-effect of using mixed precision is that the backward error could be increased a bit. In fact, we saw in section 2.2 that using mixed precision increased the error bounds of a BLR compression by a certain factor. For a number of precision formats p = 7, we should expect the error bounds to be multiplied by approximately 1 + 2(p -1)/mpblr_thd = 2.2 (a less pessimistic bound would be that, the error bound is multiplied by 1 + 2 √ p -1/mpblr_thd ≈ 1.49, as a consequence of Equation 2.23). Therefore, we could expect to observe a similar increase in the backward error of a sparse linear system. However, we did not observe such a behavior: like in the dense case (see section 2.5), the error stays very similar when adding mixed precision, and we are still far from such an increase.

Compressing contribution blocks in mixed precision

There is another possibility to reduce the memory consumption of the multifrontal factorization. The basic idea is similar to our motivation to use MPBLR as a storage format: when it is kept unused for a long time, data could be compressed as much as possible, in order to reduce the memory peak, even at the cost of some additional operations. Therefore, we allow the compression of the block of non-eliminated variables at the end of a front (the so called Contribution Block, also referred to as CB), by using a BLR compression instead of keeping them uncompressed: in fact, the off-diagonal blocks of the CB have low numerical ranks, similar to the factors. As a consequence, the memory peak will be reduced (as well as the volume of the communication). Such a feature is being implemented in the multifrontal solver MUMPS, and is referred to as option "compressCB".

However, contrary to the use of MPBLR, the time complexity of the compression operations performed here is by no means small, leading in general to a significant time overhead. Moreover, those compressions are used for storage gains only, and they do not contribute to reducing the number of operations like the compressions of the factors. Therefore, such an option should only be used when the memory consumption is critical and a certain overhead is acceptable.

We observed a reduction of the memory peak of up to of 24% when using option Com-pressCB, as illustrated in Table 3.3.

In MUMPS, we added the possibility of switching the BLR compression of the contribution block to mixed precision. Thanks to this, the memory peak is further reduced by up to 15%, which makes a total reduction of up to 35% when considering the combined gains of both options (variants "BLR" and "MPBLR * +MP-CCB * " from Table 3.3). The overhead caused by the copy/conversion operations should be negligible compared to the cost of the extra low-rank compressions that are already added with option CompressCB.

As a bonus, the MPI communication volume has been reduced even more, given the fact that a possibly huge part of the MPI communications now concern mixed-precision data. Table 3.3: Impact of option "compressCB" (CCB), and its mixed-precision variant using 7 precision formats (MP-CCB) on the memory peak and the factorization time. The notation * refers to the use of the block-admissibility condition suited for mixed precision (Equation 3.2).

Reducing the communication volume

Distributed memory parallelism in MUMPS

As explained in section 1.1.3, the multifrontal method can be parallelized by exploiting two kinds of parallelism, referred to as tree parallelism and node parallelism. In MUMPS, distributed-memory parallelism is used in order to take advantage of both sources of parallelism. This is illustrated by Figure 3.2: for the node parallelism, each front is split between several processes. Each process is mapped on a group of rows, on which it will perform all update operations. The master process is associated with the fully summed rows, while the other rows are split between the worker processes.

Sending the LU factors in mixed precision

We consider the case of an LU factorization, whose communication pattern is simpler than an LDL T factorization. We consider a node whose rows are split between several MPI processes (node parallelism), and a worker process that owns a set of rows I. Its main task is to perform update operations on its rows, in right-looking:

A ij ← A ij -L ik U kj , for i ∈ I, for j > k.
The blocks L ik , i ∈ I needed for this operation have been computed locally, and do not need to be sent by another process. However, the blocks U kj have been computed on the master node and need to be sent. Therefore, once the master process has updated a blockrow of U, it sends it to the other processes, so that they may perform the corresponding updates.

When using BLR compression in MUMPS, the volume of communication is reduced by sending the blocks U kj in a low-rank form. We added the possibility of sending these blocks in mixed precision, in order to reduce the volume of communication.

Experimental results

When using this MPBLR compression of the communications on U , as well as option compressCB in mixed precision, most communications are performed in mixed precision. As shown in Table 3.4, the results are quite promising: we can expect a reduction in the communication volume by a factor of 1.9 by adding mixed precision. When combining these gains with the ones from option compressCB, we realize that the reduction of the communication volume is quite substantial, up to a factor 13.

No time gain has been observed yet by switching the MPI communications to mixed precision. However, we can reasonably hope that another implementation, better tuned, would be able to achieve time gains thanks to the communication being in mixed precision, in the case of a large problem treated on a huge number of processes. In fact, Table 3.3 shows that reducing the communication volume may sometimes lead to actual time gains on the appropriate problem: on matrix Poisson200, using option compressCB reduces the factorization time by a factor 3.5, despite performing a greater number of operations.

Mixed precision aiming for time gains

In the previous section we presented an algorithm in which the low-rank blocks are stored using a MPBLR compression. During the triangular solve phase, these blocks were converted back to the working precision (e.g., double precision), and then used in this format for the computations.

However we know from our theoretical analysis in chapter 2 that this might not be optimal: if the precision formats used for storage are also available on hardware (e.g.: fp64 and fp32), then some of the computations can directly be performed on the MPBLR structure, without the need to convert it back to the working precision. This is the case because a low-precision part can often be used in the precision format that was used for storing it. By doing this, the accuracy of the solution is preserved, while taking advantage of the higher speed of low-precision computations.

We will first discuss algorithmic issues when switching computations to mixed precision for the product between a low-rank block and a matrix. Then we will apply it to the triangular solve, in order to obtain time gains. Finally, we will generalize this kernel so that it can be used in the factorization as well.

Main techniques

LR×FR product

One of the main low-rank kernels is the multiplication between a low-rank block XY T and a full-rank matrix M , and we will refer to it as LR×FR kernel. As explained in chapter 1, this operation can be accelerated when using low-rank approximations: if we choose to perform the multiplications in the right order, i.e. X × (Y T × M ), then the complexity is reduced compared to the FR version, from b 3 to 2rb 2 if M is a square matrix of size (b, b).

Algorithm 3.1 LR×FR product: X(Y T M ). Input: a low-rank block XY T and a full-rank matrix M Output:

Z = XY T M 1: W ← Y T M 2: Z ← XW
We now want to further reduce the time complexity by performing the operations in mixed precision, based on the algorithm presented in chapter 2.

If XY T is stored in the mixed-precision form [X 1 • • • X p ][Y 1 • • • Y p ] T , then the operations X k × (Y T k × M
) can be performed in precision k. This is illustrated in Figure 3.4, in the case where the chosen precision formats are fp64 and fp32. The full details are given in Algorithm 3.2, which can handle an arbitrary number p of precision formats (this will prove helpful for supporting half precision in the future). A minor optimization has been taken into account by looping on the precision formats from the lowest one to the highest: by doing this, we can manage to switch some additions to lower precision formats, in case p ≥ 3.

Another minor optimization would be to skip any precision format whose part is empty, making sure that none of these formats require any operation. The implementation is relatively easy: we basically need to loop on the precision formats whose parts are not empty.

× = × + × a copy in fp32 computed in fp64 computed in fp32
Figure 3.4: Illustration of a low-rank product X(Y T M ) performed in mixed precision.

Algorithm 3.2 LR×FR product in mixed precision

Input: a mixed precision low-rank block

XY T = [X 1 • • • X p ][Y 1 • • • Y p ] T and p copies of M in each precision M 1 , . . . , M p . Output: Z = XY T M 1: Z p ← 0 2: for k = p to 1 do 3: W k ← Y T k M k /* performed in precision k */ 4: Z k ← Z k + X k W k /* performed in precision k */ 5: if k > 1 then 6:
Cast Z k to precision k -1 and store the result in Z k-1 .

7:

end if 8: end for 9: Z ← Z 1

Avoid conversions by duplicating a variable in every precision format

Note that Algorithm 3.2 requires on input the copies of the full-rank block M in every precision format. Depending on the context, we may choose to cast this block on the fly, just before entering this kernel. However, most of the time we should prefer to keep those copies of a full-rank block, despite the temporary extra memory cost. Indeed, those copies will be reused many times.

A similar choice is available for the matrix Z on output: it may be duplicated in order to reduce the number of conversions. More specifically, its contributions 1, p] can either be cast on the fly, or be reused in mixed precision. In Algorithm 3.2, we have chosen to do these conversions on the fly, as soon as Z k has been computed. This choice is meaningful if only one update at a time is performed, as in a right-looking algorithm.

Z k = X k Y T k M , for k ∈ [
However, we often want to perform many updates in a row on the block Z, i.e. we may want to perform the following operations: j) . This is typically the case with a left-looking algorithm. In this case, it is still possible to reuse Algorithm 3.2 and do the conversions on the fly, one after the other: we obtain Algorithm 3.3. However, this algorithm is a bit naive, and we may try to better take advantage of this sum.

Z ← Z -q j=1 X (j) Y (j)T M (
Another possibility (see Algorithm 3.4) is to accumulate the contributions of each precision format j) . By doing this, most conversions are avoided. The only ones left are performed on each subsum Z k , when doing the operation Z ← Zp k=1 Z k . Therefore, the number of conversions has been reduced by a factor q overall. The huge majority of the operations performed are now BLAS-like matrix operations, instead of having 1 conversion every 3 matrix operations: we can expect this to be more efficient.

u k in a separate buffer Z k , i.e. perform Z k = q j=1 X (j) k Y (j)T k M (
It is technically possible to implement right-looking algorithms based on Algorithm 3.4: we would still have to compute similar sums, but a large number of them at once, with bad spatial and temporal localities. Also, the size of the temporary buffers Z k would not be small like before. Indeed, we perform a large number of such sums at the same time, and we have to duplicate all the temporary sums. In the worst case (e.g., right-looking factorization), the whole front has to be copied p times. Due to this extra memory cost, Algorithm 3.4 may not be worth using in right-looking.

Algorithm 3.3 Computation of a sum of LR×FR products

Input: a set of mixed precision low-rank blocks j) , Y (j) , M (j) ) /* using Algorithm 3.2 */ 4: end for Algorithm 3.4 Computation of a sum of LR×FR products (accumulation of updates in mixed precision, minimizing the number of conversions)

X (j) Y (j)T =[X (j) 1 • • • X (j) p ][Y (j) 1 • • • Y (j) p ] T , for j ∈ [1; q] and p copies of M (j) in each precision M (j) 1 , . . . , M (j) p , for j ∈ [1; q]. Output: Z = q j=1 X (j) Y (j)T M (j) 1: Z ← 0 2: for all j do 3: Z ← Z + prod(X (
Input: a set of mixed precision low-rank blocks X (j) Y (j)T =[X (j) 1 • • • X (j) p ][Y (j) 1 • • • Y (j) p ] T , for j ∈ [1; q] and p copies of M (j) in each precision M (j) 1 , . . . , M (j) p , for each j ∈ [1; q]]. Output: Z = q j=1 X (j) Y (j)T M (j) 1: Z p ← 0 2: for k = p to 1 do 3:
for all j do 4:

W (j) k ← Y (j)T k M (j) k /* performed in precision k */ 5: Z k ← Z k + X (j) k W (j) k /* performed in precision k */ 6: end for 7:
Cast Z k to precision k -1 and store the result in Z k-1 . 8: end for 9: Z ← Z 1

Combining MPBLR formats for storage and computations

If we want to obtain better speedups, then the set of precision formats used for storage should be the same as the formats used for computations. Indeed, by doing this we remove the need to perform costly copies and conversions on the data, as emphasized at the beginning of the section.

However, in the case where reducing the memory footprint is the main priority, it is in general beneficial to use a custom set of precision formats for storing the factors, and still do the computations using all hardware-supported precision formats. Indeed, we may hope that each computation will be quicker than in double precision. As for casting the custom precision formats to the computation formats, it should not be more costly than casting them to double precision (in both cases, the conversion has to be implemented at the software level). We could also expect to reduce the data movements, which could make the conversion quicker.

As a consequence, we will try to combine the use of the two MPBLR formats at once in order to take advantage of both the optimal storage reduction and computational time gains. In all that follows we will refer to an MPBLR format as MPBLR(sto) if it is used for storage, and we will try to refer to its low-rank blocks as X (sto) Y (sto)T instead of XY T . Similarly, we will refer to an MPBLR format used for computations as MPBLR(calc), and we will call refer to its low-rank blocks as X (calc) Y (calc)T .

When the sets of precision formats used for storage and for computations are identical, we do not need to perform actual data copies between the two formats, because they will refer to the same data structure.

Application to the triangular solve

Algorithms

We apply the use of mixed precision for computations to the triangular solution phase of the multifrontal method. We will consider both the forward elimination and the backward elimination. Both steps are quite similar, and we could use the same techniques for both of them. However, in the current MUMPS implementation, the forward solve is performed in right-looking, and the backward solve is performed in left-looking.

We use Algorithms 1.5 and 1.6 for these two steps. The function prod (at lines 4 and 3 respectively) refers to the product of a block of the factors with a block of the right-hand side (RHS). If the block of the factors is represented as low-rank, then this operation can be seen as a product LR×FR, and it can be performed using Algorithm 3.1.

We now want to switch some of the computations of the forward elimination to mixed precision, in order to obtain Algorithm 3.5. Therefore, we now use mixed precision for the LR×FR products (function prod), as in Algorithm 3.2. However, this kernel needs each block B j to be copied in every precision format: we do so at line 3. We only need to keep copies of one block at a time, which means that the extra memory cost caused by these copies will be negligible.

We also add mixed precision in the left-looking algorithm doing the backward elimination, and we obtain Algorithm 3.6. In left-looking, the updates of a block can be considered as doing the operation B i ← B ij U ij B j . It requires a sum of LR×FR products (via function sum_prods), and it may be computed using one of the algorithms presented previously.

In a first implementation, we used Algorithm 3.3 in order to compute such a sum: the terms are computed independently, with a large number of conversions needed. Another alternative would be to use Algorithm 3.4 for sum_prods, summing together contributions in each precision format.

In both cases we choose to keep copies of each block of the RHS in all precision formats, so that we do not need to cast it on-the-fly before using it. The number of conversions is reduced, but the memory is increased if we are using a large number of right-hand sides.

Algorithm 3.5 Forward elimination (Right-looking algorithm, optionally with mixed precision) 1: for j = 1 to q fs do 2:

B j ← L -1 jj B j 3: Copy B j in precisions u 2 , • • • , u p 4:
for i = j + 1 to q do 5:

B i ← B i -prod(L i,j , B j ) 6:
end for 7: end for Algorithm 3.6 Backward elimination (Left-looking algorithm, with mixed precision)

1: Copy the whole RHS B in precisions u 2 , • • • , u p 2: for i = q fs to 1 do 3: B i ← B i -sum_prods(U i1 , • • • , U iq , B j ) 4: B i ← U -1
ii B i 5: end for

Results

We evaluate the performance of the mixed-precision solve on matrix Queen_4147 from Table 4.3.We report in Tables 3.5 and 3.6 the time performance of the forward elimination and backward elimination respectively for different solve variants, for three values of n rhs .

We begin by analyzing the single RHS case (n rhs = 1). This is rather critical for EDF applications since the solve phase (whose performance can be critical when MUMPS is used as a preconditioner) often involves a unique RHS.

We first consider the MPBLR variants presented in the previous section, referred to as sto2 * and sto7 * here. Having the mixed-precision block-admissibility activated, these two variants minimize the memory cost of the factorization given their respective sets of precisions (2 and 7 precisions), which makes them interesting when reducing the memory is a major concern. However, the time spent in both forward and backward substitutions have been increased drastically, as we had already seen in Table 3.2. Using 7 precisions formats, this time is increased by a factor 2.75.

However, we may want to accelerate the computations of the solve by performing them in mixed precision instead of double precision. This is what we do in the next two solve variants. When using 7 precisions for storage, we successfully accelerated the computations by adding mixed-precision computations, as can be seen in column sto7 * +calc of both tables: the time for forward elimination has been reduced by 20%, and the time for backward elimination by 21%. However, both steps tend to be still slower than doing everything in double precision (column BLR).

If we allow to slightly compromise the storage reductions by using a set of 2 precisions instead of 7, then we obtain rather impressive time gains. The time for forward elimination has been reduced by 56% and the time for backward by 56% when doing computations in mixed precision (sto2 * +calc compared to sto2 * ). Indeed, we do not need to convert and copy the MPBLR matrix from one format to another using this variant: the same data structure is kept for both storage and computation, which proves to be more efficient. In terms of time performance, we even outperform the double-precision variant.

Finally, we analyze the multiple RHS case (n rhs ≫ 1). Tables 3.5 and 3.6 show that the gain achieved with mixed precision is not as good as with a single RHS. In certain cases, it may even be slightly slower than double precision, especially for the backward elimination. The most likely explanation is that in the multiple RHS case, the data movements associated with the RHS are much more costly and become the bottleneck of the BLR triangular solve. Thus, reducing the size of the LU factors does not have an effect as important as in the single RHS case. In fact, this observation is even true when comparing the performance of double precision BLR solve with that of the FR solve: since the RHS is not compressed, the effect of BLR compression on the time performance is also much less significant with multiple RHS. For example, on matrix Queen_4147 the double precision BLR solve is 1.6× faster than the FR one with n rhs = 1, whereas it is only 1.2× faster with n rhs = 250. Table 3.5: Times for forward elimination for the following variants:

-FR: full-rank -BLR: in double precision -sto2 * : using MPBLR with 2 precisions for storage and double precision for computations in solve -sto7 * : using MPBLR with 7 precisions for storage and double precision for computations -sto2 * +calc : using MPBLR with 2 precisions for storage and for computations -sto7 * +calc : using MPBLR with 7 precisions for storage and MPBLR with 2 precisions for computations

The admissibility condition is chosen so that the storage would be minimal. Table 3.6: Times for backward elimination for the same variants as in Table 3.5

In the next chapter, this observation motivates us to rethink the BLR solve algorithm in the context of many right-hand sides.

Block-admissibility condition

In most previous experiments presented earlier in this chapter, we have used the mixed-precision admissibility condition when handling mixed precision, be it for storage or computations (with the exception of Table 3.2). In fact, this condition somewhat optimizes the cost for storing data using a given set of precision formats, while introducing a certain overhead.

We analyze the complexity of a LR×FR product involving a block of the RHS. The matrix sizes involved are m × r and r × n for the low-rank parts, and n × n rhs for the RHS block. Around 2(m + n)r i n rhs operations are performed in precision number i. Having defined c ′ i the time cost for performing an operation in precision format i (compared to doing it in precision 1), and r′ = the complexity of the FR×FR operation, we reach the condition for a block to be worth storing as low-rank in terms of time complexity:

(m + n)r ′ < mn (3.3)
This condition is very similar to the admissibility condition for mixed-precision storage (Equation 3.2), except for the fact that the coefficients c ′ i refer to the gain in speed from format u i compared to u 1 , instead of the gain in storage. However, it is a rather reasonable approximation to consider that ∀i, c i = c ′ i , i.e., the time spent in a computation using a certain precision format is proportional to the space needed to store it. In particular, single precision is supposed to be twice as fast as double precision in our case, which means that c ′ 2 = 0.5c ′ 1 and both mixed-precision admissibility conditions are identical: there is no need to choose.

We now want to assess whether or not these theoretical time complexities are close to the reality. More precisely, we want to confirm whether or not the resulting choice of a mixed-precision block-admissibility condition is the right one in practice. Therefore, we compare the variants MPBLR(calc) (uniform-precision admissibility) and MPBLR(calc) * (mixed-precision admissibility) in Tables 3.7 and 3.8, using a single RHS.

For the forward solve, it is faster to use the admissibility for mixed precision: we obtain a time reduction of 5%. This does not seem to be the case for the backward solve (yet), and further work may be needed.

Towards an application to the factorization

After having successfully accelerated the BLR triangular solution by performing computations in mixed precision, we now want to do so for the factorization as well. We will see how to adapt the algorithms presented in chapter 2 to the implementation of the factorization of a frontal matrix. In particular, we will reuse and generalize the techniques presented in section 3.2.1.

Given an algorithm performing the factorization of a frontal matrix, we want to modify a certain number of its kernels by switching them to mixed precision. Those kernels are rather independent from each other, so it is possible to only switch some of them when modifying a uniform-precision implementation. We can therefore obtain a series of implementations in which mixed precision has been added gradually. One advantage of this method is that it allows for a progressive implementation, in which we can check that everything works as intended after the modification of each kernel. Another advantage is that we can make sure that switching each of those kernels introduces some time gains, and quantify them. Indeed, there is no guarantee that all kernels are worth switching to mixed precision, especially for the most complex ones such as LR×LR, which has the lowest granularity. In the worst-case scenario, it might be beneficial to use an implementation in which some kernels are kept in uniform precision, if it proves to be faster in some cases.

Factorization kernels

We may add mixed precision in the following kernels, independently:

• The LR×FR and FR×LR products • The Solve step (so-called TRSM): as argued in chapter 2, the step X ← L -1 X may be switched to mixed precision: if X = [X 1 • • • X p ], we perform each operation X i ← L -1 X i in precision i.

• The LR×LR products

When performing a product between a LR block and a FR block, we could reuse the algorithms presented in section 3.2.1. However, we may also want the same kernel to be able to perform a FR×LR product since it is very similar. Therefore, we are inclined to describe a more generic approach, where we want to perform an operation A × B, where only one matrix among A and B is low-rank. We decompose this kernel into several steps:

• Do the intermediary product in order to get a mixed-precision low-rank approximation Z × W of the result, via the function computeW . One of these two matrices is the result of the product between the full-rank block and the closest low-rank part of the other block.

• Compute the outer-product Z × W .

Due to implementation details chosen for MUMPS, we sometimes have to deal with the transpose of a block instead of the original block. Therefore, while describing a generic implementation of a mixed-precision product between a LR block and a FR block, we might as well allow any of the two blocks to be transposed.

As for the LR×LR kernel, its implementation details may vary. For example, in the so-called inner product we may compute W ij = Y T Ai X Bj in precision max(i, j). In all cases, we need additional conversions before performing LR×LR products on a block : the part X i should be copied in each lower precision format k > i.

Similarly to the case of a LR×FR product, we notice that a LR×LR product may be split into two steps: (i) a function computeW , which includes the inner-product and the middle-product, returns a low-rank approximation of the result; and then (ii) an outer-product, which is similar in every way to the one performed in the LR×FR kernel.

Right-looking and left-looking factorizations

The factorization may be done in right-looking. In the implementation of MUMPS, this particular case happens when a front is split between several MPI processes: as explained in section 3.1.7, a worker process will perform the updates of its blocks in right-looking, treating a block-row as soon as it is received from the master process.

We may use some of the kernels in mixed precision as explained above. Note that doing so requires the low-rank blocks to have been copied in mixed precision beforehand. If the LR×FR and FR×LR products are performed in mixed precision, then the FR blocks should also be copied in all other precision formats as well. If the TRSM step is performed in mixed precision, then the blocks L kk and U kk of the factors should be copied to all precision formats. In case LR×LR products are performed, additional conversions have to be performed: X i and Y i should be copied to all precision formats lower than i.

It seems worthwhile to perform all these copies and conversions beforehand: indeed, each block will be reused many times.

The same can be applied to a left-looking factorization, but there are a few additional possible optimizations:

• When computing a sum of outer-products in mixed precision, we may once again reuse Algorithm 3.4. This way, all contributions in precision k will be summed together, thus reducing the number of conversions. This is very similar to what was done in the case of the solution phase, except for the fact that we now add the contributions of LR×FR, FR×LR and LR×LR products at once.

• A very promising perspective would be to adapt LUA (low-rank updates accumulation) to the use of mixed precision. As explained in section 1.2.3, it is possible to consider a sum of outer-products j Z (j) W (j) as a product of two larger matrices [Z (1) • • • Z (q) ] × [W (1) ; • • • ; W (q) ]. Although the number of operations performed is left unchanged, the granularity of the matrix operations is improved by doing this, resulting in time reductions. We could adapt the principle behind Algorithm 3.4 in order to compute all the contributions in precision k as one larger matrix product, performed in precision k: j Z

(j) k W (j)
k would be replaced by

[Z (1) k • • • Z (q) k ] × [W (1) k ; • • • ; W (q) k ].
By doing LUA in mixed precision, we would improve the very small granularity induced by the use of mixed precision. In particular, the ranks of the parts in precision u 1 are often small. Some may be empty and are skipped, but there also are many such blocks having a rank r 1 of 1 or 2: the time performance of such matrix products is rather poor. By performing low-rank updates accumulations on each precision format, this aspect should be improved.

Yet another perspective would be to compute the inner-products using LUA as well. Indeed, this kernel is the one whose granularity was the most reduced by the use mixed precision, each matrix product having been split into 4.

Dissociating storage and computation formats

Contrary to the case of the solution phase, combining mixed precision for computations and storage is rather easy, and its overhead should be negligible. We just need to convert the MPBLR(calc) format to an MPBLR(sto) once at the end of the factorization of a frontal matrix. If the precision formats used for computations and for storage are identical, there is no need to perform actual data copies: we can once again reuse the same data structure.

Conclusion

In this chapter, we adapted the multifrontal method in order to benefit from the theoretical performance gains from chapter 2. We first considered the BLR compression in mixed precision (MPBLR) as a storage format only, in order to reduce the size of the factors and the memory peak of the factorization. We succeeded in obtaining reductions of the LU factor size up to 38%, without impacting too much the factorization time nor the error of the solution. We also obtained reductions in the volume of communications up to 50%, which can be a critical aspect when solving very large linear systems on massively parallel architectures.

We also described how to adapt the different steps of the multifrontal method in order to obtain speedups from the use of mixed-precision computations. We implemented the modifications of the triangular solve step. We obtained time reductions up to 12% compared to the standard BLR algorithm, and 20% compared to the mixed-precision version aiming for storage gains only, while keeping the same storage gains as the latter.

In future work we will turn to the implementation of computations in mixed precision during the BLR factorization phase of the multifrontal method, as described in section 3.2, and to obtain actual time gains. A preliminary study of the number of flops eligible for mixed precision seems to indicate that its potential for acceleration is higher than for storage. For example, as illustrated in figure 3.5, the potential factor of acceleration from the use of mixed precision is 1.9× for matrix Poisson224 with ε = 10 -9 , whereas the LU storage reduction is only 1.4×. The green bars correspond to the complexity of the MPBLR variant using 2 precision formats (fp64+fp32), while the blue bars correspond to the complexity added by not using mixed precision. Each group of 3 bars corresponds to values of ε of 10 -6 , 10 -9 and 10 -12 respectively. The percentage below the matrix name indicates which percentage of the entries belong to BLR fronts: the higher the value, the more efficient BLR compression is on this matrix. More importantly, the number above each bar (e.g. ×1.8) indicates the factor of gain expected for the mixed-precision variant, compared to using uniform precision.

Chapter 4

Hybrid algorithm for solve

Introduction

This chapter is concerned with the performance of the triangular solve phase, which is critical in several contexts. Indeed, the BLR factorization is often used as a preconditioner for iterative solvers [START_REF] Amestoy | Combining sparse approximate factorizations with mixed-precision iterative refinement[END_REF][START_REF] Higham | A new preconditioner that exploits low-rank approximations to factorization error[END_REF], which require several iterations and thus solves. Moreover, even in a pure direct solver context, some of the real-life applications where BLR solvers have been the most successful are in the field of geosciences [START_REF] Operto | Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? the Gorgon data FWI case study[END_REF][START_REF] Amestoy | Fast 3D frequency-domain full waveform inversion with a parallel Block Low-Rank multifrontal direct solver: application to OBC data from the North Sea[END_REF][START_REF] Shantsev | Large-scale 3D EM modeling with a Block Low-Rank multifrontal direct solver[END_REF][START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance, and scalability[END_REF], where we need to solve a system AX = B with many right-hand sides (RHS): B ∈ R n×n rhs is a (possibly sparse) matrix with n rhs columns, where n rhs is typically in the thousands or tens of thousands. In this case, the triangular solves LY = B and UX = Y are the bottleneck of the computation. This work started from the observation, made in section 3.2.2, that the time gains from using computations in mixed precision is underwhelming in the case of multiple RHS (see Table 3.5). While trying to understand what we initially thought to be a performance issue specific to mixed precision, we realized that a something similar arises in uniform precision. This is illustrated in Table 4.1, which provides the time for the forward solve LY = B both for the BLR and FR solvers (i.e. with and without compression), depending on the number of RHS, for a few problems of interest. For a single RHS (n rhs = 1), the BLR compression reduces the solve time by significant factors for all problems. However, with multiple RHS (n rhs = 250), the speedup achieved thanks to BLR compression becomes much smaller for all problems (for example, for the Poisson120 problem, the 3.9× speedup with n rhs = 1 becomes only a 1.7× speedup with n rhs = 250).

Therefore, this motivated us to rethink the BLR solve algorithms. The key observation is that with BLR compression, the performance of the triangular solve is memory bound, even for multiple RHS. Therefore, the performance of the BLR solve is mainly determined by its communication costs. Crucially, while BLR compression reduces the size of the LU factors and therefore the cost of accessing them, it does not reduce the size of the RHS, which are uncompressed. The consequence is that when there are many RHS, the cost of accessing them is likely to dominate the cost of accessing the LU factors, thereby reducing (or even canceling) the performance benefits of the BLR compression.

The main contribution of this chapter is to overcome this limitation by proposing new hybrid algorithms that reduce the number of accesses to the RHS and therefore the total volume of communications. We did not consider the use of mixed precision in this chapter, and all our experiments are therefore performed in uniform precision.

The rest of this chapter is organized as follows. After recalling some preliminaries on the existing BLR solve variants in section 4.2, we describe in section 4.3 several novel variants of the BLR solve that reduce its communication costs. We confirm that these new variants are indeed communication-avoiding by performing a theoretical communication volume analysis in section 4.4. To quickly analyze the performance of these new variants and assess their potential, we first develop a simplified prototype code and present our results on synthetic data in section 4.5. Based on these results, we implement a selected subset of the most promising variants in the MUMPS solver and test their performance on a range of real-life applications in section 4.6.

Throughout the chapter, we will discuss the case of the forward solve LY = B without loss of generality. The algorithms and ideas proposed in this chapter also apply to the backward solve UX = Y. For clarity of notation, we rename the forward solve LX = B hereinafter. In our experiments with the MUMPS solver, we will measure the time spent in the computations for the forward solve. The entire triangular solution phase of the solver also consists of the backward solve and some additional non-computational parts (data copies, etc.), whose cost represents a fixed overhead that is independent of the algorithm variants considered in this chapter.

Preliminaries and notations

Notations

With the multifrontal method, the forward solve LX = B amounts to a bottom-up traversal of a tree whose nodes are associated with the frontal matrices. The solution X of the global sparse problem is initialized to the right-hand side B. Then, at each node, a partial forward elimination is performed with the corresponding frontal matrix. To be specific, let L ∈ R m×n be such a frontal matrix, with m ≥ n, and let X ∈ R m×n rhs be the rows of the solution X associated with the row variables of L. We denote as L fs and X fs the top n × n subparts of L and X, respectively, and as L cb and X cb their bottom (m -n) × n subparts. L fs corresponds to the so called "fully summed" (FS) variables of the frontal matrix; these variables are ready to be eliminated by computing X fs ← L -1 fs X fs , which yields the final form of the solution X fs . L cb corresponds to the so-called "contribution block" (CB) variables of the frontal matrix, which are not ready to be eliminated; for these, the solution is merely updated as X cb ← X cb -L cb X fs . X cb is not the final form of the solution, it will be further updated by variables from other fronts, until the fronts that have X cb as fully summed variables are reached. More details about the frontal operations performed and their context can be found in chapter 1.

When using BLR compression, the frontal matrix L is partitioned in q × q fs blocks L ij ∈ R b×b , with q = m/b, q fs = n/b, and where b denotes the block size (we assume for simplicity of notation that it is the same for all blocks). The solution X is also partitioned into q blocks X i ∈ R b×n rhs . L fs is partitioned into q fs ×q fs blocks and L cb is partitioned into q cb × q fs blocks, where q cb = qq fs . The blocks L ij that are low-rank are approximated as

L ij ≈ U ij V T ij , with U ij , V ij ∈ R b×r
, where r denotes the rank of the blocks (we again assume for simplicity of notation that it is the same for all blocks).

We summarize below the notations used for a given front L and its corresponding part of the solution X, which are illustrated in Figure 4.1. Note that these notations differ a bit from those used in chapter 3.

• A, X , B, L, U: the matrix, solution, RHS, and LU factors associated with the global sparse problem;

• A, X, B, L, U , the matrix, solution, RHS, and LU factors associated with a given frontal matrix;

• m, the number of rows of L and X;

• n, the number of columns of L;

• n rhs , the number of columns of X (the number of right-hand sides);

• b, the block size;

• r, the rank of the low-rank blocks;

• L ij ∈ R b×b , the (i, j)th block of L, and U ij V T ij , its low-rank representation;

L X • X i ∈ R b×n rhs , the ith block of X;

• q fs = n/b : the number of block rows in the FS part, also the number of block columns;

• q cb = (m -n)/b: the number of block rows in the CB part;

• q = q fs + q cb = m/b, the total number of block rows.

• We also define q = q fs (q fs -1)/2 + q fs q cb , the total number of off-diagonal blocks in L;

• FS = [1 : q fs ], the set of block indices for FS variables; and

• CB = [q fs + 1 : q], the set of block indices for CB variables.

With the block partitioning defined above, the FS elimination X fs ← L -1 fs X fs leads to the recurrence relation

X i ← L -1 ii (X i - j<i L ij X j ) (4.1)
for i ∈ FS . The CB update X cb ← X cb -L cb X fs takes the form

X i ← X i - j≤q fs L ij X j (4.2)
for i ∈ CB . The BLR representation of L is exploited by computing L ij X j in the above expressions as U ij (V T ij X j ). Computations (4.1) and (4.2) thus involve two types of tasks:

• update(i, j) for i ∈ FS ∪ CB and j ∈ FS verifying i > j: X i ← X i -U ij (V T ij X j
), which can be performed only after trsm(j) has been completed; and

• trsm(i) for i ∈ FS : X i ← L -1

ii X i , which can be performed only after update(i, 1), . . . , update(i, i -1) have all been completed. There are therefore some dependencies between tasks but also some independent tasks which can be performed in any order, possibly concurrently. We next describe two variants using different orders.

Right-looking and left-looking variants

Algorithm 4.1 describes the right-looking (RL) variant, which performs the update(i, j) tasks as soon as they are ready to be performed (eager approach): as soon as trsm(j) has been completed, all the update(i, j) for i > j are immediately performed, as illustrated in Figure 4.2.

Algorithm 4.1 Right-looking variant.

1: for j ∈ FS do 2:

X j ← L -1 jj X j 3: for i > j do 4: X i ← X i -U ij (V T ij X j ) 5:
end for 6: end for Conversely, Algorithm 4.2 describes the left-looking (LL) variant, which performs the update(i, j) tasks as late as possible (lazy approach): for a given i, the update(i, j) are delayed until they are all ready to be performed together, as illustrated in Figure 4.3.

The RL and LL variants perform the same computations in two different orders. In a sequential context, it is not clear whether one variant can be expected to yield better performance than the other. However, in a parallel context, a major difference appears. The RL variant can be parallelized efficiently by executing the loop on block rows (line 3 of Algorithm 4.1) in parallel, since all update(i, j) tasks are independent for a fixed j. This approach is usually efficient because it does not lead to any conflict and the size Algorithm 4.2 Left-looking variant.

1: for i ∈ FS do 2:

for j < i do 3:

X i ← X i -U ij (V T ij X j ) 4:
end for 5:

X i ← L -1
ii X i 6: end for 7: for i ∈ CB do 8:

for j ∈ FS do 9: of the loop, qj, is large enough to expose a high amount of concurrency. In contrast, the parallelization of the LL variant is more difficult. The CB part of the update can be efficiently parallelized by executing the loop on the block-rows (line 7 of Algorithm 4.2) in parallel, since the updates for different block-rows are independent and the size of the loop, q cb , is large enough. However, for the FS part of the computation, the only loop that can be executed in parallel is the loop on the block columns (line 2 of Algorithm 4.2), which presents two difficulties. First, it requires a reduction operation to avoid conflicts since all update(i, j) tasks for a fixed i modify the same block X i . Second, the loop is only of size i -1, with i ≤ q fs , so that there is very little concurrency in the first steps of the loop, and even for the later steps because q fs is typically much smaller than q cb .

X i ← X i -U ij (V T ij X j ) 10

Parallelism in multifrontal solve

In the multifrontal solution, we must carry out several partial solves with frontal matrices following a bottom-up traversal of a tree. As a result, two types of parallelism can be exploited, as mentioned in section 1.1.3.

• Node parallelism consists in processing a given front in parallel, by parallelizing the partial solve as described above for the RL and LL variants. The amount of work required by one partial solve is usually sufficient to be efficiently parallelized only for the largest fronts, which are at the top of the tree.

• Tree parallelism consists in processing multiple fronts on different branches concurrently, using only one process per front. This allows for efficiently parallelizing the bottom of the tree, which consists of many independent fronts of small size.

Since node and tree parallelism are more efficient for the top and bottom of the tree, respectively, the best approach is to combine both types of parallelism. One possibility to do so is to exploit tree parallelism for the bottom layers of the tree, and switch to node parallelism after a given layer (so-called the "L 0 " layer) is reached. This L 0 approach is implemented in MUMPS (L 'Excellent and Sid-Lakhdar, 2014).

With this approach, the frontal triangular solve algorithms described above can be called either in a sequential setting (corresponding to fronts under the L 0 layer) or in a parallel one (corresponding to fronts above the L 0 layer). Both settings are therefore of interest in the following.

An additional source of parallelism lies in the RHS: these can be partitioned into blocks and all blocks can be handled concurrently. This type of parallelism must be used with care because partitioning the RHS in excessively small blocks might degrade performance due to the small granularity of computations. In this work we focus on exploiting parallelism within a single block of RHS.

New hybrid variants of the BLR triangular solve

A novel hybrid variant

As mentioned in section 4.1, the performance of the BLR solve is underwhelming when dealing with many RHS. Based on the description of the RL and LL variants, we can see that one common weakness of both variants is that they require multiple accesses to the entire RHS. Indeed, at each step the RL variant only reads one block of the RHS, but needs to write all the bottom part of the RHS. Conversely, at each step the LL variant only writes one block, but needs to read all the top part of the RHS. We will precisely measure the volume of communications in the next section, but it is clear that when the number of RHS n rhs is large, the accesses to the RHS can represent a much larger volume than the accesses to the BLR factors, and thus constitute the bottleneck of the computation.

Motivated by this observation, we propose in this section a novel BLR solve variant that is based on a hybrid scheme that only needs to read and write one block of the RHS per step. The main idea is to perform the read operations following a right-looking scheme, and the write operations following a left-looking scheme. To do so, we divide the update(i, j) task

X i ← X i -U ij (V T ij X j ) into two separate subtasks: • updateV(i, j): W ij = V T ij X j ; • updateU(i, j): X i ← X i -U ij W ij .
The updateV(i, j) tasks require to read the block X j of the RHS; the updateU(i, j) tasks require to write the block X i of the RHS. Thus, the idea of this hybrid variant is to perform the updateV tasks with a right-looking pattern (accessing X j once and executing immediately all tasks for this j) and the updateU tasks with a left-looking pattern (waiting that all tasks for a given i are ready so as to access X i only once). This is accomplished at the cost of having to store in a temporary workspace the W ij matrices for i > j (W ij is created at step j and consumed at step i). However, these W ij matrices are of small dimension r × n rhs , and so the cost of storing and accessing them should be small. This new hybrid variant is outlined in Algorithm 4.3 and illustrated in Figure 4.4.

Algorithm 4.3 Hybrid variant.

1: for k ∈ FS do 2:

for j < k do 3:

X k ← X k -U kj W kj 4:
end for 5:

X k ← L -1 kk X k 6:
for i > k do 7:

W ik = V T ik X k 8:
end for 9: end for 10: for k ∈ CB do 11:

for j ∈ FS do 12: X k ← X k -U kj W kj 13:
end for 14: end for

Parallelism-driven hybrid variant

One issue with the hybrid variant of Algorithm 4.3 is that it suffers from the same problems as the left-looking variant in a parallel setting: the updateU tasks corresponding to FS blocks of the RHS are performed following a left-looking pattern (loop on line 2 of Algorithm 4.3), which is of small size and requires a reduction.

To overcome this issue, we propose in Algorithm 4.4 a modified hybrid variant that is more amenable to a parallel execution. The idea is to use the hybrid scheme for the CB part of the computation only, which can be efficiently parallelized with a parallel loop on the block rows (line 11 of Algorithm 4.4), and to keep the efficiently parallelizable right-looking scheme for the FS part (with a parallel loop on line 3). Since the FS part is typically much smaller than the CB part, we can still expect to retain most of the benefits associated with the use of the hybrid scheme.

V ik for i > k in right-looking. Each product V T ik X k is stored in W ik . L X W (
: for k ∈ FS do 2: X k ← L -1 kk X k 3: for i > k do 4: if i ∈ FS then 5: X i ← X i -U ik (V T ik X k ) 6: else 7: W ik = V T ik X k 8: end if 9:
end for 10: end for 11: for k ∈ CB do 12:

for j ∈ FS do 13:

X k ← X k -U kj W kj 14:
end for 15: end for

Low-rank updates accumulation

The hybrid variant described above minimizes the number of accesses to the RHS in order to reduce the volume of communications. It also increases arithmetic intensity of the BLR solve, defined as the ratio between the number of flops (which is unchanged for all variants) and the volume of communications. In this section we now propose a further modification of the BLR solve algorithm to further increase its arithmetic intensity. The idea is based on using low-rank updates accumulation (LUA), that is, grouping together low-rank updates on the same block rows and/or block columns and applying them with a single matrix multiplication to increase the granularity of the computation.

As explained in section 1.2.3, LUA can be used to improve the performance of the BLR factorization. We now discuss how to adapt this LUA technique to the BLR solve. The situation is more complicated due to the presence of the RHS, which is not in low-rank form. In fact, we will see that LUA cannot be fully used in either the RL or LL variants of the BLR solve. The hybrid variant allows us to take full advantage of LUA, as we now explain.

updateU accumulation: the updateU(k, j) tasks X k ← X k -U kj W kj can be grouped together for all j as

X k ← X k - J j=1 U kj W kj = X k -U k1 • • • U kJ W T k1 • • • W T kJ T , (4.3)
with J = k -1 for the FS updates and J = q fs for the CB updates. The sum in (4.3) can be efficiently evaluated using only one matrix-matrix product, instead of J smaller ones.

In order to do so, all U kj , j ≤ J, must be allocated contiguously in memory, as well as all W kj , j ≤ J.

updateV accumulation: the updateV(i, k) tasks W ik = V T ik X k can be grouped together for all i as

W T i 0 k • • • W T qk T = V i 0 k • • • V qk T X k (4.4) with i 0 = k + 1.
This operation can also be performed using only one matrix-matrix product instead of qk. In order to do so, all V ik , i > k, must be allocated contiguously, as well as all W ik , i > k. Note that this allocation of the "W " workspaces is not compatible with the one needed for the updateU accumulation: updateU requires each block row to be a contiguous array, whereas updateV requires the same of each block column. As a result, if we wish to accumulate both types of tasks, the "W " arrays need to be transformed (that is, copied) from the updateU allocation scheme to the updateV one during the computation.

Note that the updateU tasks can only be accumulated with a left-looking scheme: thus, the RL variant can only benefit from the updateV accumulation. Conversely, the updateV tasks can only be accumulated with a right-looking scheme, so that the LL variant can only benefit from the updateU accumulation. The hybrid variant uses the right-looking scheme for the updateV tasks and the left-looking scheme for the updateU tasks: it can thus benefit from both LUA strategies.

• 2qrb + qrn rhs + q fs b 2 RO transfers,

• qrn rhs WO transfers,

• qbn rhs RW transfers.

Parallelism-driven hybrid variant Finally, we compute the communication volume of the parallelism-driven hybrid variant (Algorithm 4.4), which is obtained by combining the RL volume for the FS part of the computation with the hybrid volume for the CB part. At each step k ∈ FS , the parallelism-driven hybrid variant requires the following transfers:

• trsm(k) reads L kk and reads/writes X k → b 2 RO transfers and bn rhs RW transfers;

• update(i, k) for i > k and i ∈ FS reads U ik and V ik ; it also reads X k which is already in the fast memory; and finally it also reads/writes X i → 2(q fs -k)br RO transfers and (q fs -k)bn rhs RW transfers;

• updateV(i, k) for i ∈ CB reads V ik and X k , the latter still being in the fast memory; it also writes W ik → q cb br RO transfers and q cb n rhs r WO transfers.

The following additional transfers are then required at each step k ∈ CB :

• updateU(k, j) for j ∈ FS reads U kj and W kj ; it also reads/writes X k , which only needs to be loaded once → q fs br + q fs n rhs r RO transfers and bn rhs RW transfers.

Summing over all steps k ∈ FS ∪ CB , the parallelism-driven hybrid variant therefore requires a total volume of communications of • 2qrb + q fs q cb rn rhs + q fs b 2 RO transfers,

• q fs q cb rn rhs WO transfers,

• (q + q fs (q fs -1)/2)bn rhs RW transfers.

Discussion

First, we note that the RL and LL variants are not completely equivalent in terms of communications: while they require the same overall volume regardless of transfer type, our analysis shows that qbn rhs RW transfers in the RL variant have been replaced with the same volume of RO transfers in the LL variant. Therefore, in a context where RW transfers are more costly than RO ones the LL variant might outperform the RL one, at least in a sequential environment. This could for example occur if a RW transfer requires a first transfer from the slow to the fast memory and then a second transfer in the other direction.

We now seek to determine when the hybrid variant requires less communications than the LL one. To do so, we must make an assumption on the relative cost of RO, WO, and RW transfers. Under the simplifying assumption that RO and WO transfers are equally costly, and neglecting lower order terms in the expression of the communication volume, we obtain a ratio between the LL volume and the hybrid volume approximately equal to 2qrb + 2qrn rhs 2qrb + qbn rhs = 1 + n rhs /b 1 + n rhs /2r . (4.5)

Thus the condition for the hybrid variant to require less communications than the LL one is 2r ≤ b, which we can expect to be always satisfied, since it corresponds to the condition for a block to be low-rank (if r > b/2, the block requires less storage if represented as a full-rank block: see Equation 2.26). We can thus conclude that the hybrid variant should always communicate less than the LL one.

A more important question is when can we expect the hybrid variant to communicate much less than the LL one: that is, when is ratio (4.5) much less than 1? In the regime where n rhs /2r is small (just one or few RHS), the ratio (4.5) is close to 1. The hybrid variant therefore does not significantly reduce the volume of communications for small numbers of RHS. However, in the regime where n rhs /2r ≫ 1, the ratio (4.5) is approximately equal to 2r/n rhs + 2r/b. This shows that (4.5) is much less than 1 if min(n rhs , b) ≫ 2r: that is, the hybrid variant could lead to significant gains in case there is a large number of RHS and the blocks are very low-rank.

As for the parallelism-driven hybrid variant, without surprise it achieves a tradeoff between the hybrid and LL ones. It requires q fs (q fs -1)/2(b -r)n rhs extra transfers corresponding to the FS part of the computation which does not use the hybrid communication pattern. This extra volume can be expected to be small for typical cases where q fs ≪ q cb . We summarize in Table 4.2 the dominant terms in the total communication volume of the different variants. 2qrb + qrn rhs qrn rhs Parallelism-driven hybrid 2qrb + q fs q cb rn rhs q fs q cb rn rhs (q + q 2 fs /2)bn rhs

The above analysis of these hybrid variants also applies to the case where LUA is used. We do not develop a specific analysis for the use of LUA: as mentioned, we can expect LUA to also reduce the volume of communications by reducing the number of cache misses; this is however a more complex phenomenon that our simple communication model used in this section does not capture. In order to assess the potential of the new BLR solve variants, we have first developed a prototype in Fortran, which implements a partial BLR solve of a given frontal matrix. Since we aim to use this prototype for performance analysis only, we use a synthetic random matrix, and we make some further simplifications by forcing both the block size b and the rank of the blocks r to be constant.

The prototype implements the four BLR solve variants: RL, LL, hybrid, and parallelismdriven hybrid. The prototype also allows for the optional use of LUA for the updateU and/or updateV tasks. We observe the actual time spent in the computation, which allow us to compare different variants.

As explained in section 4.2.3, the frontal BLR solve algorithms are needed for large fronts at the top of the tree, where node and tree parallelism are both exploited, but also for smaller fronts at the bottom of tree, where only tree parallelism is exploited. Therefore, our goal is to use the prototype to analyze the performance of the BLR solve in two types of configurations: the first configuration uses tree parallelism only by running 36 instances in parallel (with MPI), each using a single thread; the second configuration uses both node and tree parallelism by running two instances in parallel (with MPI), each of them being parallelized with 18 OpenMP threads.

All our experiments were run on the Olympe supercomputer; each node is equipped with two 18-core Intel Skylake 6140 processors running at 2.3 GHz (for a total of 36 cores per node). We use Intel MKL 2018 for the BLAS libraries.

Performance analysis of hybrid variants

We report our performance results in Figures 4.5 and 4.6 for the tree and node parallelism configurations, respectively. For all experiments we set n rhs = 250 and q fs = q/5. For the node+tree parallelism experiments, we use large fronts (BLR block size b = 500 and q = 100 or 200, leading to a front size of 50, 000 or 100, 000). For the tree parallelism experiments, we use smaller fronts (BLR block size b = 250 and q = 30 or 50, leading to a front size of 7, 500 or 12, 500).

These results show that, as could be hoped, the hybrid variants can be faster than the RL and LL variants in many cases. The gains are the most significant when the rank r is small, and when the problem size is large, with speedups reaching a factor 2.5× in the best case. The parallelism-driven hybrid variant is not as fast as the standard hybrid one in the configuration with tree parallelism only, but is significantly superior, as expected, in the configuration with node parallelism. Indeed, as explained in section 4.3.2, the leftlooking loop at line 2 of Algorithm 4.3 involves a reduction operation on a loop of small size, whose parallelization on 18 threads is not very efficient. This is also confirmed in Figure 4.6 by the observation that for larger values of q the performance of left-looking, hybrid and parallelism-driven hybrid variants gets closer.

Performance analysis of LUA

We finally analyze the time gains obtained when combining the hybrid variant with the LUA approach described in section 4.3.3. We report the performance with and without LUA in Figure 4.7. We use a node+tree parallelism configuration and therefore take the parallelism-driven hybrid variant as baseline. The figure consists of two plots, one where the number of block-rows q is fixed and the rank r varies, and the other where r is fixed and q varies. These two plots illustrate two opposite trends:

• As r increases, the benefits of using LUA diminish, since the granularity of the low-rank updates without accumulation is already large enough to achieve good performance. figure shows that for the larger fronts with node parallelism, the parallelism-driven hybrid variant is the fastest, whereas for the smaller fronts with tree parallelism, the hybrid variant seems to be the best choice for most fronts. This confirms the trends observed with the prototype experiments in the previous section.

The optimal approach therefore seems to be to combine the two types of hybrid variants by using the standard one when only tree parallelism is used and the parallelism-driven when node parallelism is used. To assess the impact of these per-front gains on the total time, we report in Table 4.4 the cumulative time spent in all fronts. This includes the time spent in the full-rank fronts (the very smallest fronts at the bottom of the tree, where BLR compression is not exploited), which are not shown in Figure 4.8. The table compares the standard RL variant to the optimized hybrid variant (which uses the parallelism-driven algorithm only on parallel fronts). The results confirm that the hybrid variant can achieve noticeable time reductions overall.

Conclusion

The performance of BLR sparse triangular solve, which is critical in several applications, is underwhelming when there are many RHS. This is explained by the fact that the computational bottleneck is the memory access to the RHS, which are dense and far heavier than the compressed BLR LU factors. To overcome this limitation, we have proposed novel hybrid algorithms that combine right-looking and left-looking communication patterns to minimize the number of accesses to the RHS. We have carried out a communication volume analysis that proves that these new variants are indeed communicationavoiding. Based on a performance analysis on synthetic data using a simplified prototype of the BLR triangular solve, we have selected a subset of the most promising hybrid variants and implemented them in the widely used MUMPS solver. Using this implementation, we have confirmed the potential of these new variants on several real-life applications, obtaining time reductions up to 20%.

Chapter 5 Conclusion

Summary

In this thesis, we have presented, analyzed, and evaluated several techniques aiming at further improving the performance of dense and sparse direct solvers, on top of using a BLR compression. In particular, we have proposed a new variant of BLR compression in which low-precision floating-point formats are used for storing the least significant columns of a low-rank approximation.

We described such a format in chapter 2 and applied it to reduce both the time and storage complexities of the LU factorization of a BLR matrix. We performed a backward error analysis of the algorithm that demonstrates its stability, in addition to guiding the choice of precision format for each computation. Similarly to the case of introducing BLR compressions, the theoretical gains of adding mixed precision are entirely dependent on the numerical properties of the matrix. We will observe such gains under the assumption that, for most off-diagonal blocks, the singular values decrease rapidly enough. Our simulations showed that our method has a huge potential on sufficiently large BLR matrices : when using 3 precision formats, we observed a theoretical gain of a factor 1.5× to 2× on the storage cost and 2× to 2.5× on the time complexity of the algorithm, while the error stays of the same order of magnitude.

In chapter 3, we adapted the multifrontal method in order to benefit from the theoretical performance gains from chapter 2. Indeed, the multifrontal method relies on a kernel that performs a partial LU factorization of a BLR matrix, sharing many similarities with the factorization of a dense BLR matrix (see Algorithm 2.1). We first considered the BLR compression in mixed precision (MPBLR) as a storage format only, in order to reduce the size of the factors and the memory peak of the factorization. We succeeded in obtaining reductions of the LU factor size up to 38%, without impacting too much the error of the solution. We also obtained reductions of the volume of communications up to 50%, which can be a critical aspect when solving very large linear systems on massively parallel architectures.

We also described how to adapt the different steps of the multifrontal method in to support computations in more floating-point formats than CPUs. Support for fp16 and/or bfloat16 are becoming the norm, and the latest GPU architecture from NVIDIA, Hopper, even supports 8-bit floating-point formats. As a consequence, it would make sense to try adapting our mixed-precision BLR algorithms to GPUs. However, combining BLR computations and GPU may prove to be rather challenging. So far, the granularity of the computations involved was deemed to be too low to try to combine the use of BLR compression with GPU accelerators in MUMPS. Much work would be needed, including rethinking the size of the low-rank blocks, and using low-rank updates accumulation (LUA) as much as possible. We note that, given the high computational efficiency of GPUs, the algorithms tend to be more easily memory-bound, which would be an argument for trying to reduce the data movement by using advanced storage formats such as the MPBLR compression.

Finally, most of the work presented in this thesis is not exclusive to the MUMPS solver, nor even to the multifrontal method. Indeed, the use of MPBLR compression could also be applied, at least partially, to the supernodal method. It should also be interesting to try to apply some of our work to low-rank formats other than BLR, such as HODLR or H-matrices: one might think of switching each low-rank approximation to mixed precision. Finally, we mostly considered the low-rank approximation of a block based on a truncated SVD or truncated QR decomposition in this thesis. However, other compression kernels may be considered as well, as long as we keep the fundamental property of adding a high-accuracy estimation of a block and a small perturbation, the latter being potentially eligible for low precision.

Figure 1 . 1 :

 11 Figure 1.1: Elimination tree of a sparse matrix A

Figure 1 . 2 :

 12 Figure 1.2: Partial LU factorization of a front. The fully summed parts of A are the first rows and the first columns. They become L and U after the partial factorization. The contribution block (CB) corresponds to a certain sum of updates, to be added with the right variables of the parent front during the assembly step.

Figure 1 . 3 :

 13 Figure 1.3: Elimination tree after the factorization of A: it now contains a representation of the factors L and U.

Figure 1 . 4 :

 14 Figure 1.4: Ranks of the blocks T ij in the BLR approximation of matrix P64 (see Table 2.1), for ε = 10 -10 , expressed as a percentage of the block size. Dark blue blocks are full-rank, while light blue blocks have very low ranks.

Figure 1

 1 Figure 1.5: A frontal BLR matrix L and its right-hand side X

Figure 1 . 6 :

 16 Figure 1.6: Illustration of the forward and backward errors for the computation y ≈ f (x).

Figure 1 .

 1 Figure 1.7: Block structure used in Abdulah et al. (2019) for mixed-precision Cholesky

Figure 1 .

 1 Figure 1.8: A mixed-precision representation of a low-rank matrix, used by Ooi et al.

k

  = ku ℓ /(1 -ku ℓ ) for any unit roundoff u ℓ and for any k such that ku ℓ < 1. Given two matrices A and B, we use the notation |A| ≤ |B| to denote the elementwise inequality |a ij | ≤ |b ij |.The unsubscripted norm ∥ • ∥ denotes the Frobenius norm:

Figure 2 . 1 :

 21 Figure 2.1: Precision formats used for each block of a mixed-precision BLR matrix (matrix P64, ε = 10 -10 ).

Figure 2 .

 2 Figure 2.1 shows an example of a mixed precision BLR matrix, plotting for each of its blocks the precisions that are effectively used to represent it. With ε = 10 -10 and with three available precisions (fp64, fp32, and bfloat16), we can distinguish four types of blocks. The singular values of a representative example of each type are plotted in Figure 2.2.

  Diagonal block in position (15,15).

  Near field block in position(15,16).

  Mid field block in position(15,22).

  Far field block in position (15,27).

Figure 2 . 2 :

 22 Figure 2.2: Distribution of the singular values of four blocks of different type in Figure 2.1. The dashed lines indicate the thresholds εβ/u s , εβ/u h , and εβ, with β = ∥A∥ and where u s and u h denote the unit roundoffs of the fp32 and bfloat16 arithmetics, respectively.

1 :

 1 /* Input: a mixed precision low-rank matrix B and a full-rank matrix C. */ 2: /* Output: P = BC. */ 3: Initialize P to zero. 4: for ℓ = p to 1 do 5:

  .44) where c = b + 2r 3/2 . We now consider the case where B and C are partitioned into p groups B ℓ = X Bℓ Y T Bℓ and C m = Y Cm X T Cm , stored in precision u ℓ and u m , respectively. We assume that matrices B and C satisfy ∥B ℓ D∥ ≤ εβ B /u ℓ and ∥DC m ∥ ≤ εβ C /u m for ℓ, m > 1. We analyze each of the three steps separately. Inner product M = Y T B DY C Let us first analyze the computation of the inner product M . Assume M ℓm = Y T Bℓ DY Cm is computed in a given precision denoted as u M ℓm . The computed M ℓm satisfies M ℓm = M ℓm + ∆M ℓm , with

1 :

 1 /* Input: mixed precision low-rank matrices B = X B Y T B and C = Y C X T C and a diagonal matrix D. */ 2: /* Output: P = BDC. */ 3: Initialize P to zero. 4: for m = p to 1 do 5:

Algorithm 2. 4

 4 Solution to LDZ = B (triangular system with low-rank right-hand side).

1 :

 1 /* Input: a mixed precision low-rank matrix B = Y X T , a lower triangular matrix L, and a diagonal matrix D. */ 2: /* Output: a mixed precision low-rank matrix Z, solution to LDZ = B. */ 3: for ℓ = p to 1 do 4:

Figure 2 . 3 :

 23 Figure 2.3: Storage (2.112) and expected time (2.113) for three variants of the BLR LU factorization of matrix perf009d, given as a factor of gain compared with the full-rank factorization, and as a function of the backward error (2.111). Each point corresponds to a run performed with a certain value of ε, which ranges from 10 -16 to 10 -6

  Figure 2.4: Backward error (2.111) for the uniform and mixed precision BLR factorizations, for ε = 10 -12 , 10 -9 , and 10 -6 .

  nb of entries (% of double precision BLR) (a) Storage. nb of operations (% of double precision BLR) (b) Expected time.

Figure 2 . 5 :

 25 Figure2.5: For each matrix, the 3 bars correspond to ε = 10 -12 , ε = 10 -9 and ε = 10 -6 respectively. The y-axis shows the number of entries and flops required by the mixed precision variant with respect to the double precision variant (which can be slightly different from 100% because of the different conditions to represent a block under low-rank form (2.26) and (2.27)). The color breakdown gives the proportion of entries and flops in each precision. The number above each bar indicates the resulting factor of gain in storage (2.112) and expected time (2.113) achieved by the mixed precision variant compared with the double precision one.

Figure 3 . 1 :

 31 Figure 3.1: Conversion from fp32 to fp24, performed by copying the first bytes.

Figure 3

 3 Figure 3.2: Illustration of tree and node parallelism. The shaded part of each front represents its fully summed rows.

Figure 3

 3 Figure 3.3: A node split between several MPI processes. The blue part represents the content of the messages, that is to say the fully summed U factors.

  matrixn rhs time for forward elimination (s) FR BLR sto2 * sto7 * sto2 * +calc sto7 *

  matrixn rhs time for backward elimination (s) FR BLR sto2 * sto7 * sto2 * +calc sto7 *

  i r i , the weighted number of flops needed for the LR×FR product in mixed precision are f ′ = 2(m + n)r ′ n rhs . After comparing it to matrix time for forward elimination (s) MPBLR(calc)MPBLR(calc)

  (a) Size needed for storing the LU factors. (b) Expected time for the factorization. Similarly to section 2.5, such a value is based on the weighted flops, i.e. take the hypothesis that flops in single precision are twice as fast as flops in double precision.

Figure 3 . 5 :

 35 Figure3.5: Estimations of the size of the LU factors and the time spent in the multifrontal factorization, for several matrices obtained from industrial applications. The green bars correspond to the complexity of the MPBLR variant using 2 precision formats (fp64+fp32), while the blue bars correspond to the complexity added by not using mixed precision. Each group of 3 bars corresponds to values of ε of 10 -6 , 10 -9 and 10 -12 respectively. The percentage below the matrix name indicates which percentage of the entries belong to BLR fronts: the higher the value, the more efficient BLR compression is on this matrix. More importantly, the number above each bar (e.g. ×1.8) indicates the factor of gain expected for the mixed-precision variant, compared to using uniform precision.

Figure 4 . 1 :

 41 Figure 4.1: A frontal BLR matrix L and its right-hand side X. For convenience we repeat this figure, which was already present in section 1.2.3.

  Figure 4.2: Step j = 3 of Algorithm 4.1 (right-looking variant).

  Step i = 5 of Algorithm 4.2: Update of a CB block.

Figure 4

 4 Figure 4.3: Left-looking variant.

  updateU tasks for k = 3: access U kj for j < k in left-looking, and overwrite the FS block X k . updateV tasks for k = 3: access

  c) UpdateU for k = 5 : access U kj for j ∈ FS in left-looking, and overwrite the CB block X k .

Figure 4 . 4 :

 44 Figure 4.4: Hybrid variant.

1

 1

Figure 4 . 5 :

 45 Figure 4.5: Time spent in solve with tree parallelism (36 instances with 1 thread per instance), b = 250, n rhs = 250, and q fs = q/5.

Figure 4 . 6 :

 46 Figure 4.6: Time spent in solve with node+tree parallelism (2 instances with 18 threads per instance), b = 500, n rhs = 250, and q fs = q/5.

Figure 4 . 8 :

 48 Figure 4.8: Time spent in each BLR front, in the forward solve of MUMPS.

  that is, the number of columns of X ij and Y ij stored in precision u k . The next result bounds the error introduced by mixed precision BLR compression.

Theorem 2.2 (mixed precision BLR compression). Let T be a BLR approximation of A defined by (2.28)-(2.30) with β ij = ∥A∥ (global compression). If the off-diagonal blocks T ij are represented with the mixed precision representation T ij described in section 2.2,

Table 2 .

 2 1: List of matrices used in the experiments of chapter 2. We use their Schur complement corresponding to the root separator in their multifrontal factorization, whose order n is given in the second column.

	Matrix	n	b Application

  Table3.4: Communication volumes during an LU factorization with the same variants of MUMPS as in Table3.3. The runs are done on 2 computational nodes: 8 MPI processes and 9 OpenMP threads per process. The contributions to the communication volume are the U factors as well as the contribution bloc (CB).

	matrix	variant	Communication volume (GBytes)
			LU CB	Total
	Poisson200	BLR	13.8 203	217
		MPBLR *	6.7 202	209
		BLR+CCB	13.7 21	34
		MPBLR * +MP-CCB * 6.7	10	17
	thmgaz	BLR	3.9 8.6	12.5
		MPBLR	2.0 8.6	10.6
		BLR+CCB	3.9 3.9	7.9
		MPBLR * +MP-CCB * 2.0 2.1	4.2

Table 3 .

 3 7: Impact of the mixed-precision block-admissibility condition on the time spent in the forward solve with 1 RHS. MPBLR(calc) * uses the block-admissibility condition suited to mixed precision, whereas MPBLR(calc) does not. Impact of the mixed-precision block-admissibility condition on the time spent in the backward solve with 1 RHS. MPBLR(calc) * uses the block-admissibility condition suited to mixed precision, whereas MPBLR(calc) does not.

	matrix	time for backward elimination (s)
		MPBLR(calc) MPBLR(calc) * BLR
	Queen_4147 (ε = 10 -10 )	0.336	0.408	0.461
	Table 3.8:

Table 4 .

 4 1: Some motivating examples: forward solve time (s) for the so-called full-rank (FR, with no compression) and BLR solvers on 2 × 18 cores.

		n rhs = 1	n rhs = 250
	Matrix	FR BLR Ratio	FR BLR Ratio
	Poisson120	0.24 0.06 3.9× 1.37 0.78 1.7×
	Geoazur100 0.23 0.10 2.4× 2.31 1.95 1.2×
	atmosmodl	0.14 0.06 2.3× 0.67 0.52 1.3×
	Geo_1438	0.19 0.12 1.6× 1.51 1.35 1.1×
	Queen_4147 1.69 0.38 4.5× 11.13 5.90 1.9×
	Serena	0.20 0.12 1.6× 1.81 1.23 1.5×
	Transport	0.15 0.06 2.6× 0.77 0.73 1.0×

Table 4 .

 4 2: Summary of the communication volume analysis: dominant terms for each variant.

	Variant	Communication volume
		RO	WO	RW
	Right-looking	2qrb		qbn rhs
	Left-looking	2qrb + qbn rhs
	Hybrid

Table 4 .

 4 3: List of sparse matrices used in the experiments of chapter 4.

	Matrix	Order BLR ϵ Arithmetic	BLR compression	source	symmetry
					(factor size, % of FR)
	Poisson120	1.7M	10 -6	d	27%			SPD
	Poisson160	4.1M	10 -6	d	22%			SPD
	Poisson200	8.0M	10 -6	d	19%			SPD
	Geoazur100	1.6M	10 -4	c	52%		S.Operto	UNSYM
	Geoazur140	3.8M	10 -4	c	47%		S.Operto	UNSYM
	atmosmodl	1.5M	10 -6	d	30%		SuiteSparse	UNSYM
	Geo_1438	1.4M	10 -6	d	50%		SuiteSparse	SPD
	Queen_4147 4.1M	10 -6	d	28%		SuiteSparse	SPD
	Serena	1.4M	10 -6	d	40%		SuiteSparse	SPD
	Transport	1.6M	10 -6	d	42%		SuiteSparse	UNSYM
	thmgaz	5.0M	10 -6	d	30%	EDF (code_aster) UNSYM

Table 4 .

 4 4: Time spent in the forward solve in MUMPS; n rhs = 250 and 2 × 18 cores are used for all problems except Poisson200 (4 × 18).

	Matrix	Time (s)		Gain
		Right-looking Optimized hybrid
	Poisson120	0.78	0.69	-12%
	Poisson160	2.15	1.76	-18%
	Poisson200	2.30	1.83	-20%
	Geoazur100	1.95	1.78	-9%
	Geoazur140	5.26	4.70	-11%
	atmosmodl	0.52	0.49	-6%
	Geo_1438	1.35	1.30	-3%
	Queen_4147	5.90	4.80	-20%
	Serena	1.23	1.20	-2%
	Transport	0.73	0.67	-9%
	thmgaz	2.53	2.53	-0%

The cuBLAS library is an implementation of BLAS on top of the NVIDIA CUDA runtime. It allows the user to access the computational resources of NVIDIA GPU.

Ginkgo is a high-performance linear algebra library for manycore systems (including GPUs), with a focus on the solution of sparse linear systems. It is implemented using modern C++.

Acknowledgements

Merci à tous pour votre soutien, votre bienveillance, et vos retours constructifs tout au long de ma thèse. Je tiens à remercier tout particulièrement Théo Mary, qui m'a accompagné et guidé thèse CIFRE à EDF R&D. Merci aux membres de l'ancien groupe I23 d'EDF,

We have evaluated the potential of this mixed precision BLR LU factorization on a range of matrices coming from real-life problems from industrial and academic applications. We have shown that a large fraction of the entries and flops can be safely switched to lower precisions. For ε = 10 -9 , by mixing fp64, fp32, and bfloat16 arithmetics, we obtain reductions in storage of up to 2.8× with respect to double precision BLR. Moreover, assuming fp32 and bfloat16 flops are, respectively, twice and four times faster than fp64 ones, we estimate the expected time gains, predicting reductions of up to 3.5× with respect to double precision BLR. We emphasize that these gains are not achieved at the expense of accuracy: for the same accuracy, the mixed precision variant is less expensive than the uniform precision one, or, equivalently, for a fixed storage or work budget, the mixed precision variant is more accurate (Figure 2.3).

Given the very promising results obtained on large dense matrices with this new mixed precision BLR approach, we will now try to adapt it to the multifrontal method, and to develop its high-performance implementation within the sparse direct solver MUMPS (Amestoy et al., 2019a), which already exploits BLR compression.

Communication volume analysis

In this section we develop a theoretical communication analysis that aims at measuring the total volume of communications required by the different BLR solve variants. In particular we seek to prove that the hybrid variant can significantly reduce the volume of accesses to the RHS, which should provide a benefit in the case where n rhs is large.

To perform the analysis, we use a simple model of a two-level memory hierarchy: a fast but limited memory (such as a cache) and an unlimited but slow memory (such as RAM). To simplify, we assume that we have control over the transfers of data between the two levels of memory, that is, that we can choose which data are discarded from the fast memory to make space for other data that we need to load from the slow memory. We also assume that the fast memory is large enough to accommodate all the data required to perform any given update(i, j) or trsm(j) tasks (essentially one BLR block and two RHS blocks, plus any temporary workspace associated with the computation). Conversely we assume that the fast memory is not large enough to accommodate all the data needed to perform more than one of these tasks, so that after one task is completed all the data required by the next task that were not already used by the previous task need to be loaded from the slow memory.

Under this model, we compute the volume of data that has to be transferred from the slow memory to the fast one for each BLR solve variant. This analysis provides an estimate of the cost of each variant, since the BLR solve tends to be a memory-bound computation in most practical cases.

Throughout the analysis, we distinguish three types of transfers: read-only (RO: the data are used but not modified), write-only (WO), and read/write (RW: existing data are used and modified). As an example, in the operation C ← C -AB, A and B require RO transfers and C requires a RW transfer, whereas in the operation C = AB, C requires a WO transfer.

Analysis

Right-looking variant At each step j ∈ FS , the RL variant (Algorithm 4.1) requires the following transfers:

• trsm(j) reads the diagonal block L jj and reads/writes the RHS block X j → b 2 RO transfers and bn rhs RW transfers;

• each update(i, j) for i > j reads the BLR block U ij V T ij ; X j is also needed but is already in the fast memory → 2(q -j)br RO transfers; each update(i, j) also reads and writes the RHS block X i → (q -j)bn rhs RW transfers.

Summing over all steps j ∈ FS , the RL variant therefore requires a total volume of communications of • 2qrb + q fs b 2 RO transfers,

• qbn rhs + q fs bn rhs RW transfers, where we recall that q = q fs (q fs -1)/2 + q fs q cb denotes the total number of off-diagonal blocks in L.

Left-looking variant At each step i ∈ FS , the LL variant (Algorithm 4.2) requires the following transfers:

• each update(i, j) for j < i reads the BLR block U ij V T ij and the RHS block X j → (i -1)(2br + bn rhs ) RO transfers; each update(i, j) also reads and writes the RHS block X i , but it only needs to be loaded once and can then be kept in the fast memory for all subsequent updates → bn rhs RW transfers.

• trsm(i) reads L ii ; it also reads and writes X i , which is already in the fast memory from the previous updates → b 2 RO transfers.

Then, at each step i ∈ CB , it requires the following transfers:

• each update(i, j) for j = 1 : q fs reads the BLR block U ij V T ij and the RHS block X j → q fs (2br + bn rhs ) RO transfers; each update(i, j) also reads and writes the RHS block X i , but it only needs to be loaded once and can then be kept in the fast memory for all subsequent updates → bn rhs RW transfers.

Summing over all steps i ∈ FS ∪ CB , the LL variant therefore requires a total volume of communications of • 2qrb + qbn rhs + q fs b 2 RO transfers,

• qbn rhs RW transfers.

Hybrid variant At each step k ∈ FS , the hybrid variant (Algorithm 4.3) requires the following transfers:

updateU(k, j) also reads/writes X k , which only needs to be loaded once → bn rhs RW transfers;

• trsm(k) reads L kk ; X k is already in the fast memory → b 2 RO transfers;

• updateV(i, k) for i > k reads V ik and X k , the latter still being in the fast memory; it also writes W ik → (q -k)br RO transfers and (q -k)n rhs r WO transfers.

The following additional transfers are then required at each step k ∈ CB :

• updateU(k, j) for j ∈ FS reads U kj and W kj → q fs (b+n rhs )r RO transfers; updateU(k, j) also reads/writes X k , which only needs to be loaded once → bn rhs RW transfers.

Summing over all steps k ∈ FS ∪ CB , the hybrid variant therefore requires a total volume of communications of • As q increases, the benefits of using LUA increase, since there are more blocks to accumulate together, which leads to a better granularity.

Overall, the use of LUA can lead to noticeable speedups on large fronts with small ranks.

Results on real-life applications with the MUMPS solver

Based on the encouraging results obtained with our prototype in the previous section, we have implemented a subset of variants directly in the MUMPS solver. In addition to the existing RL variant (used so far by MUMPS), we have implemented the LL, hybrid, and parallelism-driven hybrid variants. For now we have not implemented the use of LUA, which we leave for future work.

We now present some experimental results obtained with these new variants of the MUMPS BLR solve on a range of real-life matrices, listed in Table 4.3. All the experiments have been performed on one node of the Olympe supercomputer previously described, using 2 MPI processes and 18 threads per MPI process, except for the tests with the Poisson200 and thmgaz matrices, which do not fit on a single node, and for which 4 × 18 threads were used. The BLR ϵ controls the accuracy of the BLR factorization: the lowrank blocks

We use double precision real arithmetic ("d") for all problems except the Geoazur ones, for which we use single precision complex arithmetic ("c"). First, in Figure 4.8, we plot the time spent in the forward solve of MUMPS for each BLR front. This leads to plots where we can distinguish two distinct groups of points. Those corresponding to the smaller fronts in the lower part of the tree where tree parallelism is used, are each processed using 1 thread. Those corresponding to the larger fronts higher in the tree, where node parallelism is used, are processed using 18 threads. The order to obtain speedups from the use of mixed-precision computations. We implemented the modifications of the triangular solve step and obtained time reductions up to 12% compared to the standard BLR algorithm, while keeping the previous storage reductions.

In the case of multiple right-hand sides, our experiments revealed a performance issue of the BLR triangular solve algorithm. Therefore, in chapter 4 we presented a new variant of this BLR algorithm that has a better data locality, as highlighted by a communication volume analysis. After implementing this variant, we were successful in reducing the computation time of this phase by up to 20%.

Software output

The experiments on new variants of the multifrontal method presented in this thesis were all based on developments made within the MUMPS solver. Some of these developments have been made available as features in advance in the consortium version of MUMPS, available to its industrial partners. None of them have been added to the public version yet.

The options allowing the use of mixed-precision BLR as a storage format were introduced in the consortium version v5.5c in April 2022. The use of computations in mixed precision for accelerating the solve phase without impacting the storage gains was added to the consortium version v5.6c in April 2023. Some of the partners of MUMPS are already taking advantage of these new capabilities. In fact, an option has been added for EDF users of code_aster, allowing them to reduce the memory consumption.

The option allowing to send LU factors in mixed precision in order to further reduce the communication volume (see section 3.1.7) will probably be added in a future release. The hybrid algorithms for accelerating the solve phase in multiple RHS (chapter 4) might one day become an experimental option of MUMPS as well.

Publication output

The contributions presented in chapter 2 have been published in IMA Journal of Numerical Analysis [START_REF] Amestoy | Mixed precision low-rank approximations and their application to Block Low-Rank LU factorization[END_REF]. These contributions and their application to the multifrontal method (chapter 3) were the subject of a talk at the conference CANUM 2022 and a poster at conference Sparse Days 2022.

The MPBLR storage format of MUMPS has also been used to reduce the memory footprint of applications in geophysics. It has contributed to a publication in The Leading Edge [START_REF] Operto | Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? the Gorgon data FWI case study[END_REF]. The reduction of the memory footprint of applications from Florian Faucher has also led to a talk by J.Y.L'Excellent at ECCOMAS Congress 2022.

The contributions from chapter 4 have been submitted to the journal SIAM Journal on Matrix Analysis and Applications (Amestoy et al., 2023a). They were also presented at conference SIAM CSE 2023.

Perspectives

We now briefly discuss a few of the remaining challenges that could be the object of future work.

The most immediate perspective is to implement the computations in mixed precision during the BLR factorization phase of the multifrontal method, as described in section 3.2, and to obtain actual time gains. Such developments are ongoing in the MUMPS solver.

The BLR triangular solve phase could still be accelerated, by trying to combine different methods presented in this thesis. The use of mixed-precision BLR could be combined with the so-called hybrid variant presented in chapter 4, in order to further reduce the data movements in case there are several RHS. Indeed, performing computations in mixed precision generally requires extra accesses to the RHS, which can be mitigated by the use of the hybrid algorithm. Both approaches could also be combined with the so-called LUA variant of the solve. This technique, presented in section 4.3.3, consists in regrouping small matrix-matrix products into bigger ones, thus improving the data locality and making the multicore parallelization easier, by leaving everything to a well-optimized BLAS function. We have already obtained some mitigated time gains on a separate prototype of LUA, although we think it should be possible to obtain better gains with a finely tuned implementation. We believe that mixed-precision variants could greatly benefit from LUA, given the fact that the granularity of the computations involved can be critically low: for some blocks, the parts in the highest precision may even contain only 1 or 2 columns, and it should be beneficial to stack several of them together using LUA. Only half of the operations can benefit from this increase of the granularity, but this aspect can be improved by using our hybrid algorithm.

The use of custom precision formats could still be the object of future research and developments. First, the implementation and integration of the conversion from mixed precision to uniform precision could be improved, by putting an emphasis on the time spent in the conversion, and not only the storage reductions. Indeed, we observed slight time overheads when using mixed precision as a storage format, and larger overheads when using our set of 7 custom precision formats. However, a correctly tuned implementation should be able to remove those overheads. One might even hope that they would turn into speedups for the triangular solve phase, since the data movement is reduced and the algorithm is memory-bound. Indeed, such speedups have been obtained in the literature for other memory-bound problems, as mentioned in section 1.4.6. Second, one may think of using more complex custom precision formats in order to reduce the storage cost of an MPBLR compression even further, as mentioned in section 3.1.1. If we allow operations more complex than byte copies, we can choose the size of the mantissa down to the bit level, thus further saving storage. The number of bits of the exponent could also be slightly reduced in order for the range of representable numbers not to be mostly wasted on unused exponents.

The use of low precision is particularly well suited to GPUs: those architectures tend

Scientific presentations

As speaker:

• Rencontres Arithmétique de l'Informatique Mathématique (RAIM), 28 May 2021

• MUMPS Consortium meeting, 24 June 2021

• Seminar at EDF organized by Olivier Boiteau : Solveurs linéaires HPC pour les études industrielles, 18 November 2021.

• Workshop on mixed precision computing, 30 May 2022, LIP6.

• CANUM, 14 June 2022, Évian-les-Bains.

• Sparse Days, Saint-Girons, 21 June 2022 (poster).

• MUMPS Consortium meeting, 29 June 2022 Subject : Using Block Low-Rank compression in mixed precision for sparse direct linear solvers Abstract: EDF performs numerical simulations in different domains of physics. Several of its software use the library MUMPS in order to perform the costly step of solving sparse linear systems in a way that is generic, robust and efficient. The goal of this work is to develop new techniques for improving the performance gains of an existing functionality of MUMPS, the Block Low-Rank (BLR) compression. By combining several formats of floating-point numbers (mixed precision), it is possible to reduce the time and memory complexities, without compromising the accuracy of the result. Based on an error analysis, we design new variants of the LU factorization of dense matrices. We then adapt this work to the case of a sparse matrix factorization with MUMPS. A first implementation uses our mixed-precision BLR compression as a storage format, thus reducing the memory consumption of MUMPS. A second implementation allows to combine these memory gains with time gains during the resolution of triangular systems. Finally, we study new techniques for improving the data locality of the BLR triangular solve with multiple right-hand sides, and obtain time reductions within MUMPS.