Grégoire Pichon

Eric Darve

Mathieu Faverge

Pierre Ramet

Jean Roman

Sparse supernodal solver using

Keywords: Sparse linear solver, block low-rank compression, PaStiX sparse direct solver, multi-threaded architectures

come

Introduction

Many scientific applications such as electromagnetism, geophysics or computational fluid dynamics use numerical models that require to solve linear systems of the form Ax = b, where the matrix A is sparse and large. In order to solve these problems, a classic approach is to use a sparse direct solver which factorizes the matrix into a product of triangular matrices before solving triangular systems.

Yet, there are still limitations to solve larger and larger systems in a black-box approach without any knowledge of the geometry of the underlying partial differential equation. Memory requirements and time-to-solution limit the use of direct methods for very large matrices. On the other hand, for iterative solvers, general black-box preconditioners that can ensure fast convergence for a wide range of problems are still missing.

In the context of sparse direct solvers, some recent works have investigated the low-rank representations of dense blocks appearing during the sparse matrix factorization, by compressing blocks through many possible compression formats such as Block Low-Rank (BLR), H, H 2 , HSS, HODLR. . . These different approaches reduce the memory requirement and/or the time-to-solution of the solvers. Depending on the compression strategy, these solvers require knowledge of the underlying geometry to tackle the problem or can do it in a purely algebraic fashion.

Hackbusch [START_REF] Hackbusch | A Sparse Matrix Arithmetic Based on H-Matrices. Part I : Introduction to H-Matrices[END_REF] introduced the H-LU factorization for dense matrices. It compresses the matrix into a hierarchical matrix format before applying low-rank operations instead of classic dense operations. The dense solver was extended to consider sparse matrices by using a nested dissection ordering to exhibit the hierarchical structure, as summarized in [START_REF] Hackbusch | Hierarchical Matrices : Algorithms and Analysis[END_REF].

In [START_REF] Grasedyck | Parallel black box H-LU preconditioning for elliptic boundary value problems[END_REF], H-LU factorization is used in an algebraic context. Performance, as well as a comparison of H-LU with some sparse direct solvers is presented in [START_REF] Grasedyck | Performance of H-LU preconditioning for sparse matrices[END_REF]. Kriemann [START_REF] Kriemann | H-LU factorization on many-core systems[END_REF] and Lizé [START_REF] Lizé | Résolution directe rapide pour les éléments finis de frontière en électromagnétisme et acoustique : H-matrices. parallélisme et applications industrielles[END_REF] implemented this algorithm using Directed Acyclic Graphs.

The Hierarchically Off-Diagonal Low-Rank (HODLR) compression technique was used in [START_REF] Aminfar | A fast block low-rank dense solver with applications to finite-element matrices[END_REF] to accelerate the elimination of large dense fronts appearing in the multifrontal method. It was fully integrated in a sparse solver in [START_REF] Aminfar | A fast, memory efficient and robust sparse preconditioner based on a multifrontal approach with applications to finite-element matrices[END_REF], and uses Boundary Distance Low-Rank (BDLR) to allow both time and memory savings, by avoiding the formation of dense fronts.

A supernodal solver using a compression technique similar to HODLR was presented in [START_REF] Chadwick | An Efficient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization[END_REF]. The proposed approach allows memory savings and can be faster than standard preconditioning techniques. However, it is slower than the direct approach in the benchmarks and requires an estimation of the rank to use randomized techniques and to accelerate the solver.

The use of Hierarchically Semi-Separable (HSS) matrices in sparse direct solvers has been investigated in several solvers. In [START_REF] Xia | Superfast Multifrontal Method For Large Structured Linear Systems of Equations[END_REF], Xia et al. presented a solver for 2D geometric problems, where all operations are computed algebraically. In [START_REF] Wang | A Parallel Geometric Multifrontal Solver Using Hierarchically Semis-Separable Structure[END_REF], a geometric solver was developed, but contribution blocks are not compressed, making memory savings limited. [START_REF] Ghysels | An Efficient Multicore Implementation of a Novel HSS-Structured Multifrontal Solver Using Randomized Sampling[END_REF][START_REF] Ghysels | A robust parallel preconditioner for indefinite systems using hierarchical matrices and randomized sampling[END_REF] proposed an algebraic code that uses randomized sampling to manage low-rank blocks and to allow memory savings.

H 2 arithmetic [START_REF] Hackbusch | Data-sparse Approximation by Adaptive H 2 -Matrices[END_REF] has been used in several sparse solvers. In [START_REF] Pouransari | Fast hierarchical solvers for sparse matrices using extended sparsification and low-rank approximation[END_REF], a fast sparse H 2 solver, called LoRaSp, based on extended sparsification was introduced. In [START_REF] Yang | Sparse Hierarchical Solvers with Guaranteed Convergence[END_REF], a variant of LoRaSp, aiming at improving the quality of the solver when used as a preconditioner, was presented, as well as a numerical analysis of the convergence with H 2 preconditioning. In particular, this variant was shown to lead to a bounded number of iterations irrespective of problem size and condition number (under certain assumptions). In [START_REF] Ho | Hierarchical Interpolative Factorization for Elliptic Operators : Differential Equations[END_REF] a fast sparse solver was introduced based on interpolative decomposition and skeletonization. It was optimized for meshes that are perturbations of a structured grid. In [START_REF] Sushnikova | Compress and eliminate" solver for symmetric positive definite sparse matrices[END_REF], an H 2 sparse algorithm was described. It is similar in many respects to [START_REF] Pouransari | Fast hierarchical solvers for sparse matrices using extended sparsification and low-rank approximation[END_REF], and extends the work of [START_REF] Ho | Hierarchical Interpolative Factorization for Elliptic Operators : Differential Equations[END_REF]. All these solvers have a guaranteed linear complexity, for a given error tolerance, and assuming a bounded rank for all well-separated pairs of clusters (the admissibility criterion in Hackbusch et al.'s terminology).

Block Low-Rank compression has been investigated for dense matrices [START_REF] Anton | A Block Low-Rank Multithreaded Factorization for Dense BEM Operators[END_REF][START_REF] Akbudak | Tile Low Rank Cholesky Factorization for Climate/Weather Modeling Applications on Manycore Architectures[END_REF], and for sparse linear systems when using a multifrontal method [START_REF] Amestoy | Improving Multifrontal Methods by Means of Block Low-Rank Representations[END_REF][START_REF] Mary | Block Low-Rank multifrontal solvers : complexity, performance, and scalability[END_REF]. Considering that these approaches are similar to the current study, a detailed comparison will be described in Section 6. One of the differences of our approach with [START_REF] Mary | Block Low-Rank multifrontal solvers : complexity, performance, and scalability[END_REF] is the supernodal context that leads to different low-rank operations, and possibly increases the memory savings. The main contribution of our work is the use of low-rank assemblies to avoid forming dense updates and to maximize memory savings.

The first objective of this work is to combine a generic sparse direct solver with recent work on matrix compression to solve larger problems, overcoming the memory limitations and accelerating the time-to-solution. The second objective is to keep the black-box algebraic approach of sparse direct solvers, by relying on methods that are independent of the underlying problem geometry. In this paper, we consider the multi-threaded sparse direct solver PaStiX [START_REF] Hénon | PaStiX : A High-Performance Parallel Direct Solver for Sparse Symmetric Definite Systems[END_REF] and we introduce a BLR compression strategy to reduce its memory and computational cost. We developed two strategies : Minimal Memory, which focuses on reducing the memory consumption, and Just-In-Time which focuses on reducing the time-to-solution (factorization and solve steps).

During the factorization, the first strategy compresses the sparse matrix before factorizing it, i.e. compresses A factors, and exploits dedicated low-rank numerical operations to keep the memory cost of the factorized matrix as low as possible. The second strategy compresses the information as late as possible, i.e. compresses L factors, to avoid the cost of low-rank update operations. The resulting solver can be used either as a direct solver for low accuracy solutions or as a high-accuracy preconditioner for iterative methods, requiring only a few iterations to reach the machine precision. The main contribution of this work is the introduction of low-rank compression in a supernodal solver with a purely algebraic method. Indeed, contrary to [START_REF] Chadwick | An Efficient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization[END_REF] which uses rank estimations (i.e. a non-algebraic criteria), our solver computes suitable ranks to maintain a prescribed accuracy.

A preliminary version of this work appeared in [START_REF] Pichon | Sparse Supernodal Solver Using Block Low-Rank Compression[END_REF]. In this paper, we introduce new orthogonalization methods 2 for the RRQR recompression kernels which leads to a better limitation of the rank growth during the factorization. We present a detailed analysis of the solver with a parallelism study and a comparison of efficiency of the low-rank kernels with respect to the original full-rank kernels. We also evaluate the impact of several parameters on the memory consumption and time to solution such as the blocking sizes, the maximum rank accepted for low-rank forms and the different orthogonalization methods. In Section 2, we go over basic aspects of sparse supernodal direct solvers. The two strategies, introduced in PaStiX, are then presented in Section 3, before detailing low-rank kernels in Section 4. Section 5 compares the two BLR strategies with the original approach, that uses only dense blocks, in terms of memory consumption, time-to-solution and numerical behaviour. We also investigate the efficiency of low-rank kernels, as well as the impact of the BLR solver parameters. Section 6 surveys in more detail related works on BLR for dense and/or sparse direct solvers, highlighting the differences with our approach, before discussing how to extend this work to a hierarchical format (H, HSS, HODLR. . .).

Background on Sparse Linear Algebra

The common approach used by sparse direct solvers is composed of four main steps : 1) ordering of the unknowns, 2) computation of a symbolic block structure, 3) numerical block factorization, and 4) triangular system solves. In the rest of the paper, we consider that all problems have a symmetric pattern given by the pattern of A + A t .

The purpose of the first step is to minimize the fill-in -zero becoming non-zero -that occurs during the numerical factorization to reduce the number of operations as well as the memory requirements to solve the problem. In order to both reduce fill-in and exhibit parallelism, the nested dissection [START_REF] George | Nested dissection of a regular finite element mesh[END_REF] algorithm is widely used through libraries such as Metis [START_REF] Karypis | METIS : A software package for partitioning unstructured graphs[END_REF] or Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide, user's manual[END_REF]. Each set of vertices corresponding to a separator constructed during the nested dissection is called a supernode.

From the resulting supernodal partition, the second step predicts the symbolic block structure of the final factorized matrix (L) and the block elimination tree. This block structure is composed of one block of columns (column block) for each supernode of the partition, with a dense diagonal block and several dense off-diagonal blocks, as presented in Figure 1 for a 3D Laplacian. The goal is to exhibit large block structures to leverage efficient Level 3 BLAS kernels during the numerical factorization. However, one may notice (cf. Figure 1) that the symbolic structure obtained with a general partitioning tool might be composed of many small off-diagonal blocks contributing to larger blocks. These off-diagonal blocks might be grouped together by adding zero to the structure if the BLAS efficiency gain is worthwhile and if the memory overhead induced by the fill-in is limited. Alternatively, it is also possible to reorder supernode unknowns to group off-diagonal blocks together without additional fill-in. A traveling salesman strategy is implemented in PaStiX [START_REF] Pichon | Reordering Strategy for Blocking Optimization in Sparse Linear Solvers[END_REF] and divides by more than two the number of these off-diagonal blocks. Other approaches [START_REF] Ghysels | An Efficient Multicore Implementation of a Novel HSS-Structured Multifrontal Solver Using Randomized Sampling[END_REF][START_REF] Amestoy | Improving Multifrontal Methods by Means of Block Low-Rank Representations[END_REF] perform a k-way ordering of supernodes, starting from a reconnected graph of a separator, to order consecutively vertices belonging to a same local part of the separator's graph. Such reordering technique also allows to reduce ranks of the low-rank blocks as shown in [START_REF] Amestoy | Improving Multifrontal Methods by Means of Block Low-Rank Representations[END_REF]. To introduce more parallelism and data locality, the final structure can then be split into tiles as it is now commonly done in dense linear algebra libraries to fit the computational units granularity. These first two steps of direct solvers are preprocessing stages independent from the numerical values. Note that these steps can be computed once to solve multiple problems similar in structure but with different numerical values.

Finally, the last two steps, numerical factorization and triangular systems solves, perform the numerical operations. We consider here only the first one for the PaStiX solver. During the numerical factorization, the elimination of each supernode (column block) is similar to standard dense algorithms : 1) factorize the dense diagonal block, xxTRF, 2) solve the off-diagonal blocks belonging to this supernode, TRSM, and 3) apply the updates on the trailing submatrix, GEMM. Those steps are then adapted to the low-rank storage format as explained in Section 3.

Block Low-Rank solver

In this section, we describe the main contribution of this paper which is a BLR solver developed within the PaStiX library. First we introduce the notation used in this article, and the basics used to integrate low-rank blocks in the solver. Then, using the newly introduced structure, we describe two different strategies leading to a sparse direct solver that optimizes the memory consumption and/or the time-to-solution. In this work, we use the same static pivoting strategy as the full-rank version of the solver and which was first introduced in a sparse direct solver in [START_REF] Li | SuperLU_DIST : A Scalable Distributed-memory Sparse Direct Solver for Unsymmetric Linear Systems[END_REF].

3.1. Notation ... ... A k ,(1:b k ) A (1: b k ), k A (0), k A ( j ),k A (i),k A (i),k
A (j),( j) A (j),( j) A (i),(j) A (j),(i)

A k ,( j)

A k ,(i) ... ... Let us consider the symbolic block structure of a factorized matrix L, obtained through symbolic block factorization. Initially, we allocate this structure initialized with the entries of A and perform an in-place factorization. We denote initial blocks by A and blocks in their final state by L (or U). The matrix is composed of N cblk column blocks, where each column block is associated with a supernode, or to a subset of unknowns in a supernode when the latter is split to create parallelism. Each column block k is composed of b k + 1 blocks, as presented in Figure 2 where : -

A (0),k (= A k,( 0
)
) is the dense diagonal block ; -A ( j),k is the j th off-diagonal block in the column block with 1 ≤ j ≤ b k , ( j) being a multi-index describing the row interval of each block, and respectively, A k,( j) is the j th off-diagonal block in the row block ; -A (1:b k ),k represents all the off-diagonal blocks of the column block k, and A k,(1:b k ) all the off-diagonal blocks of the symmetric row block ; -A (i),( j) is the rectangular dense block corresponding to the rows of the multi-index (i) and to the columns of the multi-index ( j). In addition, we denote  the compressed representation of a matrix A.

Sparse direct solver using BLR compression

The BLR compression scheme is a flat, non-hierarchical format, unlike others mentioned in the introduction. If we consider the example of a dense matrix, the BLR format clusters the matrix into a set of smaller blocks, as presented in Figure 3. Diagonal blocks are kept dense and off-diagonal blocks are admissible for compression. Thus, these off-diagonal blocks can be represented through a low-rank form uv t , obtained with a compression technique such as Singular Value Decomposition (SVD) or Rank-Revealing QR (RRQR) factorization. Compression techniques are detailed in Section 4.

We propose in this paper to apply this scheme to the symbolic block structure of sparse direct solvers. Firstly, diagonal blocks of the largest supernodes in the block elimination tree can be considered as large dense matrices which are compressible with the BLR approach. In fact, as we have seen previously, it is common to split these supernodes into a set of smaller column blocks in order to increase the level of parallelism. Thus, the block structure resulting from this operation gives the cluster of the BLR compression format. Secondly, interaction blocks from two large supernodes are by definition long distance interactions, and thus can be represented by a low-rank form. It is then natural to store them as low-rank blocks as long as they are large enough. To summarize, if we take the final symbolic block structure (after splitting) used by the PaStiX solver, all diagonal blocks are considered dense, and all off-diagonal blocks might be stored using a low-rank structure. In practice, we limit this compression to blocks of a minimal size, and all blocks with relatively high ranks are kept dense.

From the original block structure, adapting the solver to block low-rank compression mainly relies on the replacement of the dense operations with the equivalent low-rank operations. Still, different variants of the final algorithm can be obtained by changing when and how the low-rank compression is applied. We introduce two scenarios : Minimal Memory, which compresses the blocks before any other operations, and Just-In-Time which compresses the blocks after they received all their contributions. Factorize Â(i),( j) = Â(i),( j) -L(i),k Ûk,( j) LR2LR

A (0),k = L (0),k U k,(0) 7: Solve L(1:b k ),k U k,(0) = Â(1:b k ),k 8: Solve L (0),k Ûk,(1:b k ) = Âk,(1:b k ) 9:
12:

End For

13:

End For 14: End For

Minimal Memory

This scenario, described by Algorithm 1, starts by compressing the original matrix A block by block without allocating the full matrix as it is done in the full-rank and Just-In-Time strategies. Thus, all low-rank blocks that are large enough are compressed directly from the original sparse form to the low-rank representation (lines 1 -4). Note that for a matter of conciseness, loops of compression and solve over all off-diagonal blocks are merged into a single operation. In this scenario, compression kernels and later operations could have been performed using a sparse format, such as CSC for instance, until we get some fill-in. However, for the sake of simplicity we use a low-rank form throughout the entire algorithm to rely on blocks and not just on sets of values. Then, each classic dense operation on a low-rank block is replaced by a similar kernel operating on low-rank forms, even for the usual matrix-matrix multiplication (GEMM) kernel that is replaced by the equivalent LR2LR kernel operating on three low-rank matrices (cf. Section 4).

A (i),( j) = A (i),( j) -L(i),k Ûk,( j) LR2GE 10:
End For

11:

End For 12: End For

Just-In-Time

This second scenario, described by Algorithm 2, delays the compression of each supernode after all contributions have been accumulated. The algorithm is thus really close to the previous one with the only difference being in the update kernel, LR2GE, at line 9, which accumulates contributions on a dense block, and not on a low-rank form.

This operation, as we describe in Section 4, is much simpler than the LR2LR kernel, and is faster than a classic GEMM. However, by compressing the initial matrix A, and maintaining the low-rank structure throughout the factorization with the LR2LR kernel, Minimal Memory can reduce more drastically the memory footprint of the solver. Indeed, the full-rank structure of the factorized matrix is never allocated, as opposed to Just-In-Time that requires it to accumulate the contributions. The final matrix is compressed with similar sizes in both scenarios.

Summary

For the sake of simplicity, we now compare the three strategies in a dense case to illustrate the potential of different approaches. Figure 4 presents the Directed Acyclic Graph (DAG) of LU operations for the original full-rank version of the solver on a 3-by-3 block matrix. On each node, the couple (i, j) represents respectively the row and column indexes of the block of the matrix being modified.

The Minimal Memory strategy generates the DAG described in Figure 5. In this context, the six off-diagonal blocks are compressed at the beginning. In the update process, dense diagonal blocks as well as low-rank off-diagonal blocks are updated. One can notice that off-diagonal blocks are never used in their dense form, leading to memory footprint reduction. However, as we will describe in Section 4, the LR2LR operation in sparse arithmetic is quite expensive and may lead to an increase of time-to-solution.

The Just-In-Time strategy generates the DAG described in Figure 6. In this case, the six off-diagonal blocks are compressed throughout the factorization. Since those off-diagonal blocks are used in their dense form before being

(1,1) (1,2) (1,3) (2,1) (3,1) (2,2) (3,2) (2,3) (3,3) (2,2) (3,2) (2,3) (3,3) (3,3) Figure 4. Full-rank. (1,2) (1,2) (3,2) (2,2) (1,3) (1,3) (2,3) (3,3) (2,1) (2,1) (3,1) (3,1) (2,3) (2,3) (3,2) (3,2) (1,1) (2,2) (3,3) (3,3)
Figure 5. Minimal Memory.

(1,1)

(1,2) (1,3) (2,1) (3,1) (2,2) (3,2) (2,3) (3,3) (1,2) (1,3) (2,1) (3,1) (2,2) (3,2) (2,3) (3,3) (2,3) (3,2) (3,3) Figure 6. Just-In-Time. Compression GETRF TRSM LR2LR GEMM or LR2GE
compressed, there is less room for memory improvement. However, as we will see later, the LR2GE can still be performed efficiently in the sparse case.

Low-rank kernels

We introduce in this section the low-rank kernels used to replace the dense operations, and we present a complexity study of these kernels. Two families of operations are studied to reveal the rank of a matrix : Singular Value Decomposition (SVD) which leads to smaller ranks, and Rank-Revealing QR (RRQR) which has shorter time-to-solution.

Compression

The goal of low-rank compression is to represent a general dense matrix A of size m A -by-n A by its compressed version  = u A v t A , where u A , and v A , are respectively matrices of size m A -by-r A , and n A -by-r A , with r A being the rank of the block supposed to be small with respect to m A and n A . In order to keep a given numerical accuracy we have to choose r A such that ||A -Â|| ≤ τ||A||, where τ is the prescribed tolerance.

SVD

A is decomposed as UΣV t . The low-rank form of A consists of the first r A singular values and their associated singular vectors such that : Σ r A +1 ≤ τ||A||, u A = U r A , and v t A = Σ 1:r A V t r A with U r A (respectively V r A ) being the first r A columns of U (respectively V). This solution presents multiple drawbacks. It has a large cost of Θ(m 2

A n A + n 2 A m A + n 3 A ) operations to compress the matrix, and all singular values needs to be computed to get the first r A ones. On the other hand, as the singular values of the matrix are explicitly computed, it leads to the lowest rank for a given tolerance.

RRQR

A is decomposed as QRP, where P is a permutation matrix, and QR the QR decomposition of AP -1 . The low-rank form of A of rank r A is then formed by u A = Q (:,r A ) , the first r A columns of Q, and by v t A = R (r A ,:) , the first r A rows of R. The main advantage of this process is that one can stop the factorization as soon as the norm of the trailing submatrix Ã(r A +1:m A ,r A +1:n A ) = A -Q (:,r A ) R (r A ,:) P is lower than τ||A||. Thus, the complexity is lowered to Θ(m A n A r A ) operations. However, it returns an evaluation of the rank which might be a little larger than those returned by SVD and generate more flops during the numerical factorization and solve.

In conclusion, SVD is much more expensive than RRQR. However, for a given tolerance, SVD returns lower ranks. Put another way, for a given rank, SVD will have better numerical accuracy. Thus, there is a trade-off between time-to-solution (RRQR) versus memory consumption (SVD).

Note that for the Minimal Memory scenario, the first compression (of sparse blocks) may be realized using Lanczos's methods, to take advantage of sparsity. However, both SVD and RRQR algorithms inherently take advantage of these zeros. In addition, most of the low-rank compression are applied to blocks stored as dense blocks and represent the main part of the computation.

Solve

The solve operation for a generic lower triangular matrix L is applied to blocks in low-rank forms in our two scenarios :

L x = b ⇔ Lu x v t x = u b v t b .
Then, with v t x = v t b , the operation is equivalent to applying a dense solve only to u b , and the complexity is only Θ(m 2 L r x ), instead of Θ(m 2 L n L ) for the full-rank (dense) representation.

Update

Let us consider the generic update operation, C = C -AB t . Note that the PaStiX solver stores L, and U t if required. Then, the same update is performed for Cholesky and LU factorizations. We break the operation into two steps : the product of two low-rank blocks, and the addition of a low-rank block to either a dense block (LR2GE), or a low-rank block (LR2LR).

Low-rank matrix product

This operation can simply be expressed as two dense matrix products :

 Bt = (u A (v t A v B ))u t B = u A ((v t A v B )u t B ), where u A is kept unchanged if r A ≤ r B (u t
B is kept otherwise) to lower the complexity. However, it has been shown in [START_REF] Anton | A Block Low-Rank Multithreaded Factorization for Dense BEM Operators[END_REF] that the rank r AB of the product of two low-rank matrices of ranks r A and r B is usually smaller than min(r A , r B ). As u A and u B are both orthogonal, the matrix E = (v t A v B ) has the same rank as  Bt . Thus, the complexity can be further reduced by transforming the matrix product to the following series of operations :

E = v t A v B (1)
Ê = v t A v B = u E v t E (2)
u AB = u A u E (3)
v t AB = v t E v t B . (4)

Low-rank matrix addition

Let us consider the next generic operation C = Cu AB v t AB , with m AB ≤ m C and n AB ≤ n C as it generally happens in the supernodal method. This is illustrated for example by the update block A (i),( j) in Figure 2.

If C is not compressed, as it happens in the LR2GE kernel, C will be dense too, and the addition of the two matrices is nothing else than a GEMM kernel. The complexity of this operation grows as Θ(m AB n AB r AB ). If C is compressed as in the LR2LR kernel, C will be compressed too, and

Ĉ = u C v t C -u AB v t AB (5) u C v t C = [u C , u AB ]([v C , -v AB ]) t (6)
where [, ] is the concatenation operator. The update given at Equation ( 5) is the commonly named extend-add operation. Without further optimization, this operation costs only two copies. In the case of the supernodal method, adequate padding is also required to align the vectors coming from the AB and C matrices as shown in Figure 7 for the u vectors. The operation on v is similar. One can notice that in this form, the rank of the updated C is now r C + r AB . When accumulating multiple updates, the rank grows quickly and the storage exceeds the full-rank version. In order to maintain a small rank for C, recompression techniques are used. As for the compression kernel, both SVD and RRQR algorithms can be used.

Recompression using SVD. We start by computing a QR decomposition for both composed matrices :

[u C , u AB ] = Q 1 R 1 and [v C , -v AB ] = Q 2 R 2 . (7)
Then, the temporary matrix T = R 1 R t 2 is compressed using the SVD algorithm described previously. This gives the final Ĉ with :

u C = (Q 1 u T ) and v C = (Q 2 v T ). ( 8
)
The complexity of this operation is decomposed as follows : Θ((m C +n C )(r C +r AB ) 2 ) for the two QR decompositions of equation ( 7), Θ((r C + r AB ) 3 ) for the SVD decomposition, and finally Θ((m

C + n C )(r C + r AB )r C ) for the application of both Q 1 and Q 2 .
Recompression using RRQR. This solution takes advantage of the fact that u C is orthonormal to first orthogonalize u AB with respect to u C . For this operation, we refer to Section 4.3.3 to form an orthonormal basis

[u C , u AB ] such that [u C , u AB ] = [u C , u AB ] T. (9)
We follow by applying the RRQR algorithm to :

T [v C , -v AB ] t = QRP. (10)
As for the compression, we keep the k = r C first columns of Q and rows of R to form the final C :

u C = ([u C , u AB ] Q k ) and v t C = R k P. (11)
Note that u C is kept orthonormal for future updates.

When the RRQR algorithm is used, the complexity of the recompression is then composed of : Θ(r C r AB m AB ) to form the intermediate product u t C u AB , Θ(m C r C r AB ) to form the orthonormal basis, Θ(n AB r AB r C ) to generate the temporary matrix used in [START_REF] Xia | Superfast Multifrontal Method For Large Structured Linear Systems of Equations[END_REF], Θ((r C + r AB )n C r C ) to apply the RRQR algorithm, and finally again Θ((r C + r AB )n C r C ) to compute the final u C .

Orthogonalization

Let us consider the orthogonalization of [u C , u AB ] = QT , taking advantage of u C orthogonality. The following recompression of T [v C , -v AB ] t will now be referred as right recompression.

A main issue is that [u C , u AB ] may not be full-rank in exact arithmetic if both matrices share a common spectrum. Thus, if some zero columns can be removed, it will reduce the cost of the right recompression, by ignoring zero rows. In practice, for each zero column, we permute both [u C , u AB ] and [v C , -v AB ] matrices and reduce the rank involved in next computations.

The objective is to reduce the number of operations depending on r C , considering that in many cases r AB < r C as many large off-diagonal blocks receive small contributions. In order to do so, we perform Gram-Schmidt (GS) projections to take advantage of u C orthogonality. Two variants of GS are widely used : Classical Gram-Schmidt (CGS) and Modified Gram-Schmidt (MGS). Both may have stability issues, and several iterations may be performed to ensure a correct (at machine precision) orthogonality. In practice, we use CGS for its locality advantages and perform a second iteration if required. To verify that a second iteration is required, we use a widely used criterion [START_REF] Giraud | A Robust Criterion for the Modified Gram-Schmidt Algorithm with Selective Reorthogonalization[END_REF] and did not experience any orthogonality issue for the set of problems we consider.

We now present several variants for orthogonalization.

Householder QR factorization. This method performs a Householder QR factorization on the full matrix [u C , u AB ] = Q 1 R 1 . Thus, it does not exploit the existing orthogonality in u C , and cannot properly extract zero columns from the final solution to reduce the cost of later operations. However, it is probably the most stable approach and is the mostly tuned kernel in linear algebra libraries such as Intel MKL. The complexity of this operation grows as Θ(m C (r C +r AB ) 2 ).

PartialQR. It performs a projection of u AB into u C :

u AB = u AB -u C (u t C u AB ). ( 12
)
From this projection, we obtain :

[u C , u AB ] = [u C , u AB ] I u t C u AB 0 I . (13)
In practice, we perform two projections to ensure the stability. We now have u C ⊥ u AB and want to orthogonalize u AB to obtain an orthogonal basis. In order to do so, we either perform a Householder QR factorization u AB = u AB R or a rank-revealing QR factorization at the machine precision u AB = Q k R k P to remove zero columns. The first approach can make a good use of Level 3 BLAS operations while the second is less efficient but may construct zero columns and then reduce the right recompression.

Classical Gram-Schmidt. In this variant, we orthogonalize one-by-one each vector of u AB with one or two CGS iterations depending on the criterion presented before. The orthogonalization of the i-th column of u AB to obtain the corresponding column of u AB is performed with :

u AB i = u AB i -u AB 0:i-1 u t AB 0:i-1 u AB i . ( 14
)
Each column is then removed after its orthogonalization if it is a zero column. This reduces the final rank of the update and the cost of the following operations.

In the experimental study, we will compare the three approaches in terms of operation count and performance. For orthogonality, we did not observe any numerical issue in our test problems. Note that for both PartialQR and CGS, only last r AB rows of ([v C , -v AB ]) t are updated. For those two methods, the orthogonalization complexity grows as Θ(m C r C r AB ).

Summary

Table 1 presents the computational complexity for the two low-rank strategies with respect to the original version of the solver. The main factor of complexity is computed under the assumption that m C ≥ m A ≥ m B , r A ≥ r B , m C ≥ n C , and r C ≤ r C . It does not depend on n A but on the ranks r A and r B : there are fewer operations to be performed. On the other hand, the Minimal Memory strategy requires using either the SVD or RRQR recompression, for which the complexity depends on m C and n C , the dimensions of the block C. This explains why this strategy appears as a higher complexity than the original solver.

When considering dense matrices, the low-rank matrix (blocks) are updated by contributions of the same size. In that case, the complexity of the low-rank update becomes asymptotically cheaper than the full-rank updates, and leads to performance gain. This is exploited in dense BLR solvers as [START_REF] Anton | A Block Low-Rank Multithreaded Factorization for Dense BEM Operators[END_REF], and in the CUFS (Compress, Update, Factor, Solve) strategy of the BLR-MUMPS solver, which compresses a dense front before applying operations between low-rank blocks of the same size.

In the supernodal approach, blocks belonging to last level separators receive many small contributions, as illustrated in Figure 1. Thus, r C is often close to or equal to r C and lower than r C + r AB : the rank is often invariant when applying a small contribution, which makes RRQR recompression more efficient than SVD recompression, and especially when an orthogonalization method that reduces the basis on the fly is used.

In terms of complexity, it is important to notice that the LR2LR update depends on the size of the target block C, and not on the contribution size as in full-rank. Thus, low-rank contributions between blocks of similar sizes are asymptotically cheaper than full-rank contributions on these blocks. Dense and multifrontal solvers exploit this property as they work only on regular block sizes. At the opposite, small updates to larger blocks, that regularly appear in supernodal solvers, have a higher cost in low-rank than in full-rank. In summary, we can decompose supernodal solvers in two groups of contributions : 10

Θ(n A r A r B ) (1) : Θ(n A r A r B ) (2) : Θ(r 2 A r B ) (2) : Θ(r A r B r AB ) (2) : Θ(r 2 A r B ) (2) : Θ(r A r B r AB ) (3), ( 4
) : Θ(m A r A r AB ) (3), ( 4
) : Θ(m A r A r AB )
LR matrix addition --

(7) : Θ(m C (r C + r AB ) 2 ) ( 12
) : Θ(m C r C r AB ) (SVD) : Θ((r C + r AB ) 3 ) (10) : Θ(n C (r C + r AB )r C ) (8) : Θ(m C (r C + r AB )r C ) (11) : Θ(m C (r C + r AB )r C ) Dense update Θ(m A m B n A ) Θ(m A m B r AB ) - - Main factor Θ(m A m B n A ) Θ(m A m B r AB ) Θ(m A m B r AB ) Θ(m C (r C + r AB ) 2 ) Θ(m C (r C + r AB )r C )
-the updates that will occur within each supernode at the top of the elimination tree and which are generated by regular split of the supernode to generate parallelism, as well as updates from direct descendants of the sons. These contributions work on regular sizes and will reduce the complexity of the solver as observed in dense and multifrontal solvers ; -the updates from deeper descendants in the elimination tree that generate updates from small blocks to larger nodes of the elimination tree. In these contributions, the complexity is defined by the target block, and thus the cost of the low-rank updates may increase the complexity of the solver with respect to the full-rank solver. The experiments showed that this overhead slows down the solver with the Minimal Memory strategy on small problems, but when the problem size increases this strategy outperforms the full-rank factorization.

The main advantage of the Minimal Memory scenario is that it can drastically reduce the memory footprint of the solver, since it compresses the matrix before the factorization. Thus, the structure of the full-rank factorized matrix is never allocated, and the low-rank structure needs to be maintained throughout the factorization process to lower the memory peak.

Kernel implementation

In practice, as we manage both dense and low-rank blocks in our solver, we adapt the extend-add operation for each basic case to be as efficient as possible. The optimizations are designed for the RRQR version, in practice SVD is only used as a reference for the optimal compression rates.

One of the important criteria for efficient low-rank kernels is the setup of the maximum rank above which the matrix will be considered as non compressible. The setting of this parameter is strategy and application dependent. In practice for an m-by-n matrix, if the main objective is the memory consumption, as in the Minimal Memory strategy, the ranks are limited to mn m+n to reduce as much as possible the final size of the factors without considering the number of flops that are generated. In contrast, if the objective is to reduce the time-to-solution, then the maximum ranks are defined by min(m,n) 4 for which the low-rank operations remains cheaper in number of operations than the full-rank ones. In practice, it also depends on the ranks of the blocks that will be involved in the update and cannot be computed before compression. For a real-life application, the criterion has to be set depending on the number of factorizations (second criterion is more important) and the number of solves (first criterion will reduce the size of the factors and thus the cost of the solve step). The impact of this parameter is studied in 5.6.3.

Eventually, when a low-rank block receives a contribution with a high rank, i.e., (r C + r AB ) ≤ max rank , the C matrix is decompressed to receive the update and then recompressed, instead of directly adding low-rank matrices. This is denoted as Decompression / Recompression of C in the following section.

There is also a trade-off between using a strict compression criterion and a smaller one due to the efficiency of lowrank operations with respect to dense operations. In practice, we add another parameter named rank ratio to control the maximum rank authorized during compression as a percentage of the strict theoretical rank.

Experiments

Experiments were conducted on the Plafrim1 supercomputer, and more precisely on the miriel cluster. Each node is equipped with two Intel Xeon E5-2680v3 12-cores running at 2.50 GHz and 128 GB of memory. The Intel MKL 2017 is used for BLAS and SVD kernels. The RRQR kernel is issued from the BLR-MUMPS solver [START_REF] Amestoy | Improving Multifrontal Methods by Means of Block Low-Rank Representations[END_REF], and is an extension of the block rank-revealing QR factorization subroutines from LAPACK 3.6.0 (xGEQP3) to stop the factorization when the precision is reached.

The PaStiX version used for our experiments is available on the public git repository2 as the tag 6.0.0. The multi-threaded version used is the static scheduling version presented in [START_REF] Lacoste | Scheduling and memory optimizations for sparse direct solver on multi-core/multi-gpu cluster systems[END_REF]. Note that for low-rank strategies, we never perform LL t factorization because compression can destroy the positive-definite property. In the case where the matrix is SPD, we use LDL t factorization for low-rank strategies.

For the initial ordering step, we used Scotch [27] 6.0.4 with the configurable strategy string from PaStiX to set the minimal size of non-separated subgraphs, cmin, to 15. We also set the f rat parameter to 0.08, meaning that column aggregation is allowed by Scotch as long as the fill-in introduced does not exceed 8% of the original matrix.

In experiments, blocks that are larger than 256 are split into blocks of size within the range 128-256 to create more parallelism while keeping sizes large enough. The same 128 criteria is used to define the minimal width of the column blocks that are compressible. An additional limit on the minimal height to compress an off-diagonal block is set to 20. We set the rank ratio to 1 to illustrate the results when the rank is as strict as possible to obtain Flops and/or memory gains. CGS orthogonalization is also used. Note that in Section 5.6 we try different blocking sizes to experiment the impact on the solver, relax the maximum ranks authorized, and compare the orthogonalization strategies.

Experiments were computed on a set of 3D matrices from The SuiteSparse Matrix Collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF], described in Table 2. We also used 3D Laplacian generators (7 point stencils), and defined lap120 as a Laplacian of size 120 3 . Note that when results showing numerical precision are presented, we used the backward error on b : ||Ax-b|| 2 ||b|| 2 .

SVD versus RRQR

The first experiment studies the behaviour of the two compression methods coupled with both Minimal Memory and Just-In-Time scenarios on the matrix Atmosmodj. Table 3 presents the sequential timings of each operation of the numerical factorization with a tolerance of 10 -8 , as well as the memory used to store the final coefficients of the factorized matrix. We can first notice that SVD compression kernels are much more time consuming than the RRQR kernels in both scenarios following the complexity study from Section 4. Indeed, RRQR compression kernels stop the computations as soon as the rank is found which reduces by a large factor the complexity, and this reduction is reflected in the time-to-solution. However, the SVD allows, for a given tolerance, to get better memory reduction in both scenarios.

Comparing the Minimal Memory and the Just-In-Time scenario, the compression time is minimized in the Minimal Memory scenario because the compression occurs on the initial blocks which hold more zero and are lower rank than when they have been updated. The time for the update addition, extend-add operation, becomes dominant in the Minimal Memory scenario, and even explodes when SVD is used. This is expected as the complexity depends on the largest blocks in the addition even for small contributions (see Section 4). Note that this ratio will evolve in favor of the extend-add operation on larger matrices where the ratio of updates of same size becomes dominant with respect to the number of updates from small blocks. For both compression methods, both scenarios compress the final coefficients with similar rates.

The diagonal block factorization time is invariant in the five strategies : the block sizes and kernels are identical. Panel solve, update product, and solve times are reduced in all low-rank configurations compared to the dense factorization and the timings follow the final size of the factors, since this size reflects the final ranks of the blocks.

To conclude, the Minimal Memory scenario is not always able to compete with the original direct factorization on these test cases due to the costly update addition. However, it reduces the memory peak of the solver to the final size of the factors. The Minimal Memory/RRQR offers a 50% memory reduction while doubling the sequential time-to-solution. The Just-In-Time scenario competes with the original direct factorization, and divides by three the time-to-solution with RRQR kernels.

Performance

Figure 8 presents the overall performance achieved by the two low-rank scenarios with respect to the original version of the solver (where lower is better) on the previously introduced set of six matrices. All versions are multithreaded implementations and use the 24 cores of one node. The scheduling used is the PaStiX static scheduler developed for the original version, this might have a negative impact for low-rank strategies by creating a load imbalance. We study only the RRQR kernels as the SVD kernels have shown to be much slower. Three tolerance thresholds are studied for their impact on the time-to-solution and the accuracy of the backward error. The backward errors printed on top of each bar correspond to the use of one refinement step.

τ =10 -4 τ =10 -8 τ =10 -12
(b) Minimal Memory scenario using RRQR. Figure shows that the Just-In-Time/RRQR scenario is able to reduce the time-to-solution in almost all cases of tolerance, and for all matrices which have a large spectrum of numerical properties. These results show that applications which require low accuracy, as the seismic application for instance, can benefit by up to a factor of 3.5. Figure 8(b) shows that it is more difficult for the Minimal Memory/RRQR scenario to be competitive. The performance is often degraded with respect to the original PaStiX performance, with an average loss of around a factor 2, and the tolerance has a much lower impact than for the previous case.

For both scenarios, the backward error of the first solution is close to the entry tolerance. It is a little less accurate in the Minimal Memory scenario, because approximations are made earlier in the computation, and information is lost from the beginning. However, these results show that we are able to catch algebraically the information and forward it throughout the update process. The Minimal Memory scenario is slower than the original solver in most cases, but it is a strategy that efficiently reduces the memory peak of the solver. Figure 9 presents the gain in the memory used to store the factors at the end of the factorization of the set of six matrices with respect to the block dense storage of PaStiX. In this figure, we also compare the memory gain of the SVD and RRQR kernels. We observe that in all cases, SVD provides better compression rate by finding smaller ranks for a given matrix and a given tolerance. The quality of the backward error is in general slightly better with the SVD kernels despite the smaller ranks. The second observation is that the smaller the tolerance (10 -12 ), the larger the ranks and the memory consumption. However, the solver always presents a memory gain larger than 50% with larger tolerance (10 -4 ).

Memory consumption

Figure 10(a) presents the evolution of the size of the factors as well as the full consumption of the solver (factors and management structures) on 3D Laplacians with an increasing size. The memory limit of the system is 128GB. The original version is limited on this system to a 3D Laplacian of 8 million unknowns, and the size of the factors quickly increases for larger number of unknowns. With the Minimal Memory/RRQR scenario, we have now been able to run a 3D problem on up to 36 million unknowns when relaxing the tolerance to 10 -4 . From the same experiment, Figure 10(b) presents the number of operations evolution depending on the Laplacians size. One can note that for a small number of unknowns, Minimal Memory/RRQR scenario performs more operations than the original, full-rank, factorization. However, for a large number of unknowns, the number of operations is reduced by a large factor. For instance, for a 280 3 Laplacian with a 10 -8 tolerance, the number of operations is reduced by a factor larger than 36. It demonstrates the potential of our approach : even if operations are less efficient and may lead to a time-to-solution overhead for relatively small problems, the Minimal Memory strategy enables the computation of larger problem and reduces its time-to-solution. In this experiment, the Minimal Memory strategy becomes faster for Laplacian larger than 150 3 with a 10 -8 tolerance. Figure 10. Scalability of the memory on top and the number of operations on bottom with three tolerance thresholds of the Minimal Memory/RRQR scenario for 3D Laplacians of size n 3 with n ∈ [START_REF] Xia | Superfast Multifrontal Method For Large Structured Linear Systems of Equations[END_REF]330]. The full-rank scenario, in purple, is given as a reference with the numbers computed from the symbolic factorization.

The memory of the Just-In-Time scenario has not been studied, because in our supernodal approach, each supernode is fully allocated in a full-rank fashion in order to accumulate the updates before being compressed. Thus, the memory peak corresponds to the totality of the factorized matrix structure without compression and is identical to the original version. To reduce this memory peak, a solution would be to modify the scheduler to a Left-Looking approach that would delay the allocation and the compression of the original blocks. However, it would need to be carefully implemented to keep a certain amount of parallelism in order to save both time and memory. A possible solution is the scheduling strategy presented in [START_REF] Agullo | Robust memory-aware mappings for parallel multifrontal factorizations[END_REF] to keep the memory consumption of the solver under a given limit.

Convergence and numerical stability

Figure 11 presents the convergence of the GMRES iterative solver preconditioned with the low-rank factorization at tolerances of 10 -4 and 10 -8 . The iterative solver is stopped after reaching 20 iterations or a backward error lower than 10 -12 .

With a tolerance of 10 -8 , only a few iterations are required to converge to the solution. Note that on the Audi and Geo1438 matrices, which are difficult to compress, a few more iterations are required to converge. With a larger tolerance 10 -4 , it is difficult to recover all the information lost during the compression, but this is enough to quickly get solutions at 10 -6 or 10 -8 . Note that the GMRES process benefits from the compression, as the solve step.

Parallelism

For the full-rank factorization, supernodes are split between processors depending on the corresponding number of operations. When a target block C receives a contribution, a lock is used to ensure that the block is not modified simultaneously by several threads. For the Just-In-Time strategy, a similar lock is used since the update operations directly apply dense modifications. However, for the Minimal Memory strategy, the update operation is decomposed into two parts : Formation of contribution and Addition of contribution. Because the formation of contribution does not depend on the target block C, the lock is only positioned for the addition of contribution, which may increase the level of parallelism.

Figure 12 presents the speedup of the full-rank factorization and low-rank strategies using tolerances of 10 -4 and 10 -8 for the atmosmodj matrix. The speedup for the full-rank version is above 12, for a relatively small matrix. The speedup of the Minimal Memory strategy is above 11, while the speedup of the Just-In-Time strategy is around 8. As the supernodes distribution cannot predict the ranks and the corresponding number of operations, load balancing may be degraded. Minimal Memory strategy scalability exceeds Just-In-Time strategy because there are less constraints on locks. Note that on recent architectures, the maximal speedup can not be obtained, as the CPU frequency is reduced when all cores are used on the node, while it is increased for single core operation. We thus computed an approximate maximum speedup of 20.7 on this architecture for the Intel MKL BLAS GEMM operation.

Kernels analysis

In this section, we focus on the atmosmodj matrix to illustrate the performance of basic kernels as well as the impact of several parameters. As in practical cases we use RRQR instead of SVD, we will focus on this compression kernel for each strategy.

Performance of basic kernels

We evaluate the performance rate (in GFlops/s) for the full-rank factorization and for the two low-rank strategies, and use a 10 -8 tolerance. Figure 13(a) (respectively Figure 13(b)) presents the runtime distribution among kernels for the full-rank (respectively Just-In-Time) factorization. We can note that in both cases, the Update process is the most time-consuming. For the Just-In-Time strategy, Compression and TRSM are not negligible, because the Update runtime is much reduced with respect to full-rank factorization.

Figure 13(c) presents the runtime distribution among kernels for the Minimal Memory strategy with three different levels from the left to the right : the main steps of the solver, the details for the Update kernels, and the details for the low-rank updates. Note that xx2fr refers to a full-rank update within a compressible supernode (i.e., an update to a dense block which was originally considered compressible) while xx2lr is a low-rank update, xx being one of the four possible matrix products : low-rank/low-rank, low-rank/full-rank, full-rank/low-rank or full-rank/full-rank. Update dense corresponds to blocks that were not considered compressible and managed in dense arithmetic throughout all operations ; as expected the underlying operations represent a small time of computation.

We observe that the low-rank update is the most time consuming part. In practice, the Formation of contribution (cf. Section 4.3.1) is quite cheap, while applying the update is expensive ; as we have seen in Table 1, it depends on the size and the rank of the target C. This update addition is decomposed into three main operations : Orthogonalization of the [u C , u AB ] matrix (cf. Section 4.3.3), Recompression of C + AB (cf. Eq (2)) and Update of the basis (cf. Eq [START_REF] Wang | A Parallel Geometric Multifrontal Solver Using Hierarchically Semis-Separable Structure[END_REF]). If the contribution rank is too high to take advantage of recompression, we perform a Decompression / Recompression of C (cf. Section 4.5).

Table 4 presents the performance of most time-consuming kernels for each type of factorization on a machine where around 32 GFlops/s can be raised for each CPU core when all cores are used. One can note that the performance of kernels for the full-rank factorization is close to the machine peak : the original solver makes good use of Level 3 BLAS even if there are many small blocks. On the other hand, low-rank kernels performing Level 3 BLAS kernels, i.e., Formation of contribution, Update of the basis, TRSM, are running at 1 3 of the peak. It is due to lower granularity generated by the smaller blocks, which reduce the arithmetic intensity of the kernels. For Compression and Recompression of C+AB kernels, the efficiency is even worse due to the behaviour of RRQR. In our implementation, RRQR is a modification of LAPACK xgeqp3 and xlaqps routines. Some stability issues presented in [START_REF] Howell | Cache efficient bidiagonalization using blas 2.5 operators[END_REF] prevent to make efficient use of Level 3 BLAS kernels.

To summarize, low-rank strategies are useful to reduce the overall number of operations. However, due to the poor efficiency of low-rank kernels, the gain in flops does not directly translate into timing reduction. One can expect that for larger problems, the reduction in flops will more easily translate into time-to-solution reduction.

Impact of blocking parameters

In previous experiments, we always used a splitting criterion of 256 with a minimal size of 128 : blocks larger than 256 are split into a subset of blocks larger than 128. For low-rank strategies, we consider only blocks larger than 128 as compressible. In Table 5, we evaluate the impact of using different blocking sizes (between 128 and 256 or between 256 and 512) for the full-rank factorization and for three different tolerances using Minimal Memory and Just-In-Time strategies. The minimum width of compressible supernodes is set to the minimum blocksize : either 128 or 256. All results are performed using 24 threads on the atmosmodj matrix.

Table 5. Impact of the blocking size parameter on the number of operations, the factorization time, and the memory for the atmosmodj on the Full-rank strategy (top) and on both Minimal Memory and Just-In-Time strategies (bottom). Both low-rank strategies are using RRQR kernels.

(a) Full-rank strategy. The first observation concerns the full-rank factorization. The blocking sizes impact both the granularity and the level of parallelism. For a sequential run, it is suitable to use large blocking sizes to increase granularity and thus the efficiency of Level 3 BLAS operations. On the other hand, it may degrade parallel performance if there are many workers and not enough supernodes.

For the Minimal Memory strategy, we have seen in Table 1 that the number of operations depends on the target size and rank and thus increases a lot with the blocking size. This is true even in the case where ranks are not that much impacted by the blocking size. In practice, we observe that for each tolerance, increasing blocking size degrades factorization time. In addition, as less data is considered as compressible, the size of the factors is growing and this will increase the solve cost.

For the Just-In-Time strategy, the impact of blocking size will mostly depend on ranks. If ranks are small, using larger blocks will increase the performance of RRQR and thus reduce the time-to-solution. However if ranks are higher, it will reduce the level of parallelism as in the full-rank factorization. We always observe a gain using a larger blocking size, but the ratio between the use of 256/512 versus 128/256 is decreasing when the tolerance is lower : the granularity gain causes some parallelism issues.

Impact of rank ratio parameter

In previous experiments, we use strict maximum ranks, i.e., the limit ranks to reduce the number of flops or the memory consumption : mn m+n for Minimal Memory strategy and min(m,n)

4
for Just-In-Time strategy. In practice, we have seen that for both strategies, kernel efficiency is poor with respect to classical Level 3 BLAS operations. In Table 6 we evaluate the impact of relaxing the constraint on a strict rank for the Minimal Memory and the Just-In-Time strategies. The ratio parameter corresponds to a percentage of the strict maximum rank (used to obtain smaller ranks) and avoid the overhead of managing low-rank blocks with a high rank by turning back these blocks into full-rank form. All results are performed using 24 threads on the atmosmodj matrix.

For the Minimal Memory strategy, there are two main observations. Firstly, considering that the update complexity depends on the size but also on the rank of the target C, the burden on recompression of high-rank blocks is reduced, and sometimes even the number of operations. Secondly, managing more blocks in full-rank fashion can reduce timeto-solution due to the poor efficiency of low-rank kernels. In practice, we observe that using a relaxed criterion always reduces time-to-solution while having only a little impact on the memory.

For the Just-In-Time strategy, the maximum rank criterion cannot be set theoretically, because it depends on all ranks within a same supernode. However, contrary to the Minimal Memory approach, the number of operations being really reduced with respect to full-rank factorization, there is a time-to-solution gain and it seems suitable to compress as many blocks as possible. In practice, relaxing the max-rank criteria is only interesting for 10 -4 tolerance, for which the granularity is really small.

Note that in both cases, relaxing the burden on large ranks increases the size of the factors, and in the same way the cost of the solve. Thus, selecting a suitable criteria will depend on the application and the ratio between the number of factorizations and the number of solves.

Orthogonalization cost

In previous Minimal Memory experiments, we used CGS as orthogonalization process. As presented in Section 4.3.3, some other approaches can be investigated. In Table 7, we present the impact of using CGS, Householder QR or PartialQR on the number of operations as well as on the efficiency of the solver. It only impacts intermediate ranks and not the final size of the factors. As predicted by its complexity, Householder QR factorization is more expensive than both CGS and PartialQR. The difference between CGS and PartialQR is related to the number of columns of zero. While CGS can remove those columns during computations, PartialQR can only deal with those zero after all operations : the number of operations increases for each tolerance.

However, in terms of time, PartialQR outperforms CGS especially for small tolerances (10 -12 ) due to the efficiency of Level 3 BLAS operations. One can expect that for a larger tolerance (10 -4 ), CGS may be faster than PartialQR, because the smaller granularity will degrade the assets of Level 3 BLAS operations.

Discussion

In this section, we discuss the positioning of our low-rank strategies with the closest related works and we give some perspectives in extending this work to a hierarchical format.

Sparse H solver. In [START_REF] Hackbusch | Hierarchical Matrices : Algorithms and Analysis[END_REF], different approaches using H-matrices are summarized, including the extension to the sparse case. The authors use a nested dissection ordering as in our solver, and thus there is no fill-in between distinct branches of the elimination tree. However, contributions of a supernode to its ancestors are considered as full, in the sense that all structural zero are included to generate the low-rank representation. Thus, they do not have extend-add (LR2LR) operations between low-rank blocks of different sizes : assembly is performed as in our Minimal Memory scenario, but without zero padding. The memory consumption is higher with their approach because some structural zeros are not managed : unlike them, we perform a block symbolic factorization to consider sparsity in all contributions between supernodes.

Dense BLR solver. A dense BLR solver was designed by Livermore Software Technology Corporation (LSTC) [START_REF] Anton | A Block Low-Rank Multithreaded Factorization for Dense BEM Operators[END_REF]. In this work, the full matrix is compressed at the beginning and operations between low-rank blocks are performed. This approach performs low-rank assembly as in our Minimal Memory scenario. As it handles dense matrices, the extend-add process is performed between low-rank matrices of the same size and zero padding is not required. Thus, the LR2LR operation is less costly than the full-rank update in this context.

Sparse Multifrontal BLR solver. A BLR multifrontal sparse direct solver was designed for the Mumps solver. The different strategies are summarized in [START_REF] Mary | Block Low-Rank multifrontal solvers : complexity, performance, and scalability[END_REF] and a theoretical study of the complexity of the solver for regular meshes is presented in [START_REF] Amestoy | On the complexity of the block low-rank multifrontal factorization[END_REF]. In the current implementation, the fronts are always assembled in a full-rank form before being compressed. The authors suggest as an alternative solution to assemble the fronts panel by panel to avoid allocating the full front in a dense fashion, but this was not implemented and it will impact the potential parallelism. Mumps is a multifrontal solver, and BLR is considered when eliminating a dense front and not between fronts while our supernodal approach has a more global view of all supernodes.

Our scenario Just-In-Time is similar to the FCSU (Factor, Compress, Solve, Update) strategy from [START_REF] Mary | Block Low-Rank multifrontal solvers : complexity, performance, and scalability[END_REF] where the fronts are compressed panel by panel during their factorization and where only full-rank blocks are updated (LR2GE kernel). The LUAR (Low-Rank Update Accumulation with Recompression) groups together multiple lowrank products to exploit the memory locality during the product recompression process. This could be similarly used in the Just-In-Time strategy, but would imply larger ranks in the extend-add operations of the Minimal Memory strategy.

The CUFS (Compress, Update, Factor, Solve) [START_REF] Mary | Block Low-Rank multifrontal solvers : complexity, performance, and scalability[END_REF] strategy is the closest to the Minimal Memory scenario, as it performs low-rank assembly (LR2LR kernel) within a front. However, the fronts are still first allocated in full-rank to assemble the contributions from the children. Within a front, the CUFS strategy is similar to what was proposed by LSTC for dense matrices. In that case, since blocks are dense and low-rank operations are performed between blocks of equal sizes, zeros padding is not necessary as in the Minimal Memory strategy.

We believe that there is more room for memory savings using the Minimal Memory strategy for two reasons : 1) we avoid the extra cost of fronts inherent to the multifrontal method and 2) the low-rank assembly avoids forming dense blocks before compressing them. We demonstrated that the Minimal Memory strategy reduces the number of operations performed on 3D Laplacians and other matrices coming from various applications, which is the main contribution of this paper. A theoretical complexity study has to be performed to investigate if such an approach performs asymptotically the same number of operations as Just-In-Time or Mumps strategies with a constant factor overhead. Low-rank assembly without randomization techniques has only been studied in [START_REF] Xia | Superfast Multifrontal Method For Large Structured Linear Systems of Equations[END_REF] for 2D problems with the knowledge of the underlying problem, and computing the complexity bounds for a more general case is still an open research problem, we are currently investigating.

Extension to hierarchical formats. With the aim of extending our solver to hierarchical compression schemes, such as H, HSS, or HODLR, we consider graphs coming from real-life simulations of 3D physical problems. From a theoretical point of view, the majority of these graphs have a bounded degree or are bounded density graphs [START_REF] Miller | Density graphs and separators[END_REF], and thus good separators [START_REF] Lipton | A separator theorem for planar graphs[END_REF] can be built. For a n-vertices mesh, time complexity of a direct solver is in Θ(n 2 ), and we expect to build a low-rank solver requiring Θ(n Let us consider the last separator of size Θ(n 2 3 ) for a 3D mesh, and one of the largest low-rank blocks of this separator in a hierarchical clustering. They have asymptotically the same size. Previous studies have shown that such a block may have a rank of order Θ(n 1 3 ). For the Just-In-Time scenario, maintaining such a block in a dense form before compressing it requires Θ(n 4 3 ) memory. Thus, we will still have the same memory peak, but we might encounter a large overhead when compressing off-diagonal blocks with current RRQR and SVD kernels. For the Minimal Memory scenario, we have seen that the cost of the solver can be split in two stages. The low-level one from the elimination tree that generates small updates to large contribution blocks and that might become more expensive with the hierarchical compression, and the high-level where blocks fit the hierarchical structure and generate flops savings. To overcome the issue of the low-level contributions, new ordering techniques need to be investigated to minimize the number of updates on larger off-diagonal blocks. We will also investigate the use of randomization techniques to lower the complexity of the updates from the bottom of the elimination tree.

Conclusion

We presented a new Block Low-Rank sparse solver that combines an existing sparse supernodal direct solver PaStiX and low-rank compression kernels. This solver reduces the memory consumption and/or the time-to-solution depending on the scenario. Two scenarios were developed. For the set of real-life problems studied, Minimal Memory saves memory up to a factor of 4 using RRQR kernels, with a time overhead that is limited to 2.8. Large problems that could not fit into memory when the original solver was used can now be solved thanks to the lower memory requirements, especially when low accuracy solutions. For larger problems, one can expect that the reduction of the number of operations will translate into a time-to-solution reduction. We experienced this behaviour with large Laplacians : we are now able to solve a 330 3 unknowns Laplacian while the original solver was limited to a 200 3 unknowns Laplacian, the time-to-solution being reduced over 150 3 unknowns.

Just-In-Time reduces both the time-to-solution by a factor up to 3.5 and the memory requirements of the final factorized matrix with similar factors to Minimal Memory. However, with the current scheduling strategy, this gain is not reflected in the memory peak.

Two compression kernels, SVD and RRQR, were studied and compared. We have shown that, for a given tolerance, both approaches provide correct solutions with the expected accuracy, and that RRQR, despite larger ranks, provides faster kernels. In addition, we demonstrated that the solver can be used either as a low-tolerance direct solver or as a good preconditioner for iterative methods, that normally require only a few iterations before reaching the machine precision. A comparison with other preconditioners (AMG, ILU(k)) will be performed in future work to measure the impact of using a low-rank factorization as a preconditioner.

In the future, new kernel families, such as RRQR with randomization techniques, will be studied in terms of accuracy and stability in the context of a supernodal solver. To further improve the performance of Minimal Memory and close up the gap with the original solver, aggregation techniques on small contributions will also be studied. This will lead to the extension of this work to hierarchical compression in large supernodes that could further reduce the memory footprint and the solver complexity.

Regarding Just-In-Time, future work is focused on studying smart scheduling strategies that combine Right-Looking and Left-Looking approaches in order to find a good compromise between memory and parallelism for the target architecture. Sergent et al. [START_REF] Sergent | Controlling the Memory Subscription of Distributed Applications with a Task-Based Runtime System[END_REF] studied scheduling with memory constraints for dense factorization using BLR, and we expect it can be adapted for sparse direct solvers. This will follow up recent work on applying parallel runtime systems [START_REF] Lacoste | Scheduling and memory optimizations for sparse direct solver on multi-core/multi-gpu cluster systems[END_REF] to the PaStiX solver.

Figure 1 .

 1 Figure 1. Symbolic factorization of a 10 × 10 × 10 Laplacian matrix partitioned using Scotch.

Figure 2 .

 2 Figure 2. Symbolic block structure and notations used for the algorithms for the column block of index k, and its associated blocks. The larger gray boxes highlight the fact that the contributions modify only a part of the target blocks.

Figure 3 .

 3 Figure 3. Block Low-Rank compression.

Algorithm 1

 1 Right looking block sequential LU factorization with Minimal Memory scenario. /* Initialize A (L structure) compressed */ 1: For k = 1 to N cblk Do 2: Â(1:b k ),k = Compress( A (1:b k ),k ) 3: Âk,(1:b k ) = Compress( A k,(1:b k ) ) 4: End For 5: For k = 1 to N cblk Do 6:

For j = 1 to b k Do 10 :

 10 For i = 1 to b k Do /* LR to LR updates (extend-add) */ 11:

Algorithm 2 1 : 7 : 8 :

 2178 Right looking block sequential LU factorization with Just-In-Time scenario. For k = 1 to N cblk Do 2: Factorize A (0),k = L (0),k U k,(0) /* Compress L and U off-diagonal blocks */ 3: Â(1:b k ),k = Compress( A (1:b k ),k ) 4: Âk,(1:b k ) = Compress( A k,(1:b k ) ) 5: Solve L(1:b k ),k U k,(0) = Â(1:b k ),k 6: Solve L (0),k Ûk,(1:b k ) = Âk,(1:b k ) For j = 1 to b k Do For i = 1 to b k Do /* LR to dense updates */ 9:

Figure 7 .

 7 Figure 7. Accumulation of two low-rank matrices when sizes do not match.

( a )

 a Just-In-Time scenario using RRQR.

Figure 8 .

 8 Figure 8. Performance of both low-rank strategies with three tolerance thresholds. The backward error of the solution is printed on top of each bar.

Figure 9 .

 9 Figure 9. Memory peak for the Minimal Memory scenario with three tolerance thresholds and both SVD and RRQR kernels. The backward error of the solution is printed on top of each bar.

Figure 11 .

 11 Figure 11. Convergence speed for the Minimal Memory/RRQR scenario with two tolerance thresholds.

4 Figure 12 .

 412 Figure 12. Speedup of the factorization for the atmosmodj matrix with 1 to 24 threads.

Figure 13 .

 13 Figure 13. Breakdown of the most time-consuming kernels for the full-rank strategy (top-left), the Just-In-Time strategy (top-right), and the Minimal Memory strategy (bottom) on the atmosmodj case using the RRQR kernels with τ = 10 -8 . The analysis of the Minimal Memory strategy is shown at three different levels from the left to the right : the global operations, the partition of the updates between low-rank and full-rank, and the details of the low-rank updates.

Table 1 .

 1 Summary of the operation complexities when computing C = C -AB t .

		GEMM (Dense)	LR2GE (Just-In-Time)	LR2LR (Minimal Memory)
			SVD	RRQR	SVD	RRQR
			(1) :
	LR matrix product	-

Table 2 .

 2 Real-life matrices used in experiments.

	Matrix	Precision Method	Size	Field	Ops (TFlops) Memory (GB)
	Atmosmodj Audi Hook Serena Geo1438	d d d d d	LU LL t LDL t LDL t LL t	1 270 432 943 695 1 498 023 1 391 349 1 437 960 geomechanical model of earth atmospheric model structural problem model of a steel hook gas reservoir simulation	12.1 5.5 8.6 28.6 18.0	16.3 9.5 12.7 21.7 20.1

Table 3 .

 3 Cost distribution on the Atmosmodj matrix with τ = 10 -8 .

		Full-rank	Just-In-Time SVD RRQR	Minimal Memory SVD RRQR
	Factorization time (s)
	Compression	-	4.1e+02 3.4e+01	1.8e+02 5.6+00
	Block factorization (GETRF)	7.2e-01	7.4e-01	7.3e-01	7.8e-01	7.6e-01
	Panel solve (TRSM)	1.7e+01	6.9e+00 7.4e+00	7.6e+00 7.9e+00
	Update
	Formation of contribution	-	-	-	9.9e+01 4.2e+01
	Addition of contribution	-	-	-	3.0e+03 7.3e+02
	Dense udpate (GEMM)	4.6e+02	1.3e+02 9.7e+01	2.8e+01 2.4e+01
	Total	4.7e+02	5.5e+02 1.4e+02	3.6e+03 8.1e+02
	Solve time (s)	6.3e+00	1.9e+00 3.0e+00	1.5e+00 3.2e+00
	Factor final size (GB)	16.3	6.95	7.49	6.85	7.31
	Memory peak for the factors (GB)	16.3	16.3	16.3	6.85	7.31

Table 4 .

 4 Kernels efficiency per core for full-rank, Just-In-Time and Minimal Memory strategies, on atmosmodj with τ = 10 -8 .

	Strategy	Name	Nb calls	Time (s)	GFlops Perf (GFlops/s)
	-	GETRF	74665	981.6	6.4	6.5
	Full-rank	TRSM Update	1572464 10993468 446003.1 12869.6 17862.0 421.2	23.6 28.9
		Compression	44544	32831.8	146.6	4.5
	Just-In-Time	TRSM	1572464	9312.8	111.1	11.9
		Update	10993468 126681.3	1464.1	11.6
		TRSM	1572464	7355.9	103.8	14.1
		Update
		dense	5782698	4838.5	24.6	5.1
	Minimal Memory	xx2fr xx2lr	2207642	31280.6	285.0	9.1
		Formation of contribution	3036024	39319.2	493.7	12.6
		Orthogonalization	2688031 115957.9	630.2	5.4
		Recompression of C + AB	2499058 301204.3	755.4	2.5
		Update of the basis	2499058	54132.7	742.5	13.7
		Decompression / Recompression of C	556650 256821.7	2236.3	8.7

Table 6 .

 6 Impact of the maximum rank ratio on the number of operations, the factorization time, and the memory for the atmosmodj case with RRQR kernels for both Minimal Memory and Just-In-Time strategies.

		Precision		τ = 10 -4			τ = 10 -8			τ = 10 -12
		Ratio	1	0.5 0.25	1	0.5 0.25	1	0.5	0.25
		Operations (TFlops) 1.8	1.9	2.3	4.8	4.8	4.5	8.3	7.2	7.4
	Minimal Memory	Fact. time (s)	48.7 49.7 43.7 100.0 85.4 56.0 146.0 105.1 65.6
		Memory (GB)	4.7	4.7	5.4	6.53	7.0 8.94 8.31	9.26 11.71
		Operations (TFlops) 0.5	0.6	0.9	0.8	1.4	2.6	1.3	2.8	5.1
	Just-In-Time	Fact. time (s)	11.9 11.6 10.9 14.2	15.4 17.8 19.0	20.8	25.3
		Memory (GB)	4.9	5.0	5.9	6.71 7.23 9.12 8.43	9.36 11.73

Table 7 .

 7 Impact of the orthogonalization method on the number of operations and the factorization time for Minimal Memory strategy on atmosmodj.

	Precision		τ = 10 -4			τ = 10 -8		τ = 10 -12
	Orth. Method	CGS QR PartialQR CGS	QR	PartialQR CGS	QR	PartialQR
	Operations (TFlops)	1.8	2.8	2.0	4.8	8.0	5.1	8.3	13.2	8.7
	Fact. time (s)	50.5 60.3	48.5	99.1 128.9	96.3	145.2 191.8	138.6

https://www.plafrim.fr

https://gitlab.inria.fr/solverstack/pastix

) operations. For memory requirements, the direct approach leads to an overall storage of Θ(n

), while we target a Θ(n log(n)) complexity.

Acknowledgments

This material is based upon work supported by the DGA under a DGA/Inria grant. Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https ://www.plafrim.fr/).