Amélioration de

Keywords: linear sparse direct solver, block low-rank compression, ordering

Through the recent improvements toward exascale supercomputer systems, huge computations can be performed in reasonable times by using massively parallelized operations. Unfortunately, the increase of the computational units in these systems does not lead to a rise in the memory available per core. Therefore, this memory limitation forces the scientists/engineers to not only eciently parallelize the operations but also minimize the memory used.

Many scientic and engineering applications have to solve large sparse linear systems of the type Ax = b. Although the direct methods are the most robust solutions for these systems, they are costly in terms of their memory usage and time-to-solution. In this respect, the low-rank representations have been recently introduced into these solvers to reduce the time and memory footprint.

In this work, our goal is to improve the low-rank feature of the block low-rank (BLR) sparse supernodal direct solver PaStiX. For this purpose, we compare some compression methods to determine the fastest kernel, which keeps the representative data with the smallest rank possible. Then, we focus on improving supernodal solver by reducing the number of re-compression during the updates. Firstly, we study the separator reordering strategies to identify the poorly compressible blocks involved in these updates and reduce their occurrences. Secondly, we propose an orthogonal solution to predict the compressibility of the blocks before the numerical factorization. This last approach relies on the use of the level of ll of a symbolic block incomplete factorization. Thanks to these optimizations, the memory usage has been reduced more eectively compared to the state of the art solvers while also improving the time to solution. This thesis is a requested rst step toward a advanced sparse direct solver using hierarchical compression schemes.

Grâce aux récentes améliorations apportées par les nouveaux supercalculateurs exaopiques, des simulations de grande envergure peuvent être eectuées dans des délais raisonnables en utilisant des opérations massivement parallélisées. Malheureusement, l'augmentation du nombre d'unités de calcul dans ces systèmes n'est pas accompagnée d'une augmentation de la mémoire disponible par c÷ur. Par conséquent, cette limitation de la mémoire oblige les scientiques et les ingénieurs à non seulement paralléliser ecacement les opérations, mais aussi à minimiser la mémoire utilisée.

De nombreuses applications scientiques et d'ingénierie doivent résoudre de grands systèmes linéaires creux du type Ax = b. Bien que les méthodes directes soient les solutions les plus robustes numériquement, elles sont coûteuses en termes d'utilisation de la mémoire et de temps de résolution. À cet égard, les techniques de compression de rang faible ont été récemment introduites dans ces solveurs an de réduire leur complexité et leur empreinte mémoire.

Dans ce travail, notre objectif est d'améliorer les représentations de rang faible par bloc dans le solveur direct supernodal creux PaStiX. Pour cela, nous comparons quelques méthodes de compression pour déterminer le noyau le plus rapide, qui conserve les données représentatives avec le plus petit rang possible. Nous nous concentrons sur l'amélioration du solveur supernodal en réduisant le nombre de recompressions lors des opérations de mise à jour. Tout d'abord, nous étudions les stratégies de renumérotation des séparateurs pour identier les blocs peu compressibles impliqués dans ces mises à jour et réduire leurs occurrences. Ensuite, nous proposons une solution orthogonale pour prédire la compressibilité des blocs avant la factorisation numérique. Cette dernière approche s'appuie sur l'utilisation de la notion de niveau de remplissage pour une factorisation incomplète symbolique par blocs. Grâce à ces optimisations, l'utilisation de la mémoire a été fortement réduite, par rapport aux solveurs à l'état de l'art, tout en améliorant le temps de résolution. Cette thèse est un premier pas nécessaire vers un solveur direct creux performant utilisant des schémas de compression hiérarchiques.

Mots-clés : solveur linéaire creux direct, compression à rang faible, numérotation Acknowledgement Writing this thesis was a long journey with a lot of ups and downs for me. I would not succeed it without the help of many people. Firstly, I would like to thank to my director Pierre Ramet for always supporting me both in terms of my work and motivating me psychologically. Secondly, I want to thank to my co-director Mathieu Faverge for always being available for me. I learned a lot from him. Both my directors supported me a lot in terms of work and personal problems. I really feel lucky to work with them. Here, I also want to thank to Gregoire Pichon who is involved in all parts of my work as much as my directors. Anytime I have trouble understanding something, he was my savior :) Among all his responsibilities, he could always nd a time to help me.

I want to express my gratitude to all my jury members. My discussions with my reviewer Alfredo Buttari on my thesis has contributed me a lot and it was very precious for me. I would like to thank to my second reviewer Hatem Ltaief for providing me very nice references and calming me down just before my defense. I thank to both of my jury members Laura Grigori and Francois Pellegrini for their interesting questions which let me to think more broad in my eld.

I feel very lucky to meet/work with all the Inria people (especially the HiePACS members). They contributed me to improve myself both in the eld and personally. I thank to all my colleagues Alena Shilova, Martina Iannacito, Romain Peressoni, Marek Felsoci, Mathieu Verite, Yanfei Xiang, Tony Delarue, Xunyi Zhao, Jean-Francois David, Mahmoud Shouman, Pierre Esterie and Alycia Lisito for all the nice conversations we had and always being supportive.

I really appreciate Romain Peressoni for always supporting me for all my problems and preparing the food for my defense day. In addition, Alena Shilova and Martina Iannacito deserve a special appreciation since they were always by my side whenever I feel stress. And last but not least, I thank to Emmanuel Agullo who taught me how to play fairly in babyfoot :) I am grateful to have my brother Koray Korkmaz as my life-long best friend and supporter. At this point of my life, I am looking forward to continue my journey with Olivier Beaumont, Lionel Eyraud-Dubois and Laercio Lima Pilla who were always very warm and helping since I met them. v

Résumé long

Grâce aux avancées récentes en calcul scientique, nous pouvons désormais remplacer les expériences coûteuses et/ou dangereuses de la vie réelle par des simulations numériques. Les domaines dans lesquels la simulation est utilisée sont vastes et incluent, par exemple, la météorologie, la biologie, la fusion nucléaire ou l'aérodynamique. Bien qu'il soit encore possible d'optimiser ces codes de calcul, les simulations les plus récentes sont déjà capables de fournir des solutions étonnamment précises et rapides.

La simulation de problèmes physiques dans le domaine de la science computationnelle exige que l'espace réel soit discrétisé. Les équations mathématiques continues sont également discrétisées sur ces points d'échantillonnage du domaine entier. L'ensemble des équations qui en résulte peut être représenté comme une matrice. Dans de nombreuses applications, seules les interactions entre les points proches de la discrétisation sont prises en compte. En eet, les interactions entre les points éloignés sont souvent faibles, et les valeurs correspondantes dans la matrice sont donc négligeables. Par conséquent, la matrice résultante comprend de nombreux termes nuls et est appelée une matrice creuse. Plus le nombre de points discrétisés et la précision attendue sont importants, plus le coût de calcul de la solution du système sera élevé. Ainsi, les améliorations dans ce domaine visent à résoudre des systèmes de plus grande taille, tout en réduisant le temps de restitution, pour une précision xée.

Grâce aux récentes améliorations apportées par les nouveaux supercalculateurs exaopiques, des simulations de grande envergure peuvent être eectuées dans des délais raisonnables en utilisant des opérations massivement parallélisées. Malheureusement, l'augmentation du nombre d'unités de calcul dans ces systèmes n'est pas suivie par une augmentation de la mémoire disponible par c÷ur. Par conséquent, cette limitation de la mémoire oblige les scientiques et les ingénieurs à non seulement paralléliser ecacement les opérations, mais aussi à minimiser la mémoire utilisée.

Lors de la résolution de grands systèmes linéaires creux, la structure de la matrice peut être exploitée pour réduire la consommation de mémoire et la complexité des calculs. Parmi toutes les approches permettant de résoudre des systèmes linéaires creux, les méthodes directes représentent les approches les plus robustes au détriment d'exigences élevées en matière de mémoire et de complexité calculatoire. Ces solveurs factorisent la matrice initiale en un produit de matrices triangulaires. Ces systèmes triangulaires résultants peuvent ensuite être résolus pour obtenir la solution du problème.

Selon le type d'applications, il existe plusieurs variantes de ces méthodes directes. Selon les caractéristiques de la matrice, ces méthodes peuvent être basées sur diérents 1 types de factorisation comme LL T , LDL T , ou LU . Ici, L, U et D représentent respectivement une matrice triangulaire inférieure, une matrice triangulaire supérieure et une matrice diagonale. De plus, selon la manière dont les contributions (mises à jour) sont appliquées pendant la factorisation, on peut considérer une approche supernodale ou multifrontale. La première stratégie permet d'appliquer directement toutes les mises à jour et d'éviter les surcoûts mémoires pour leur stockage. La seconde utilise un stockage temporaire pour les contributions an de pouvoir les appliquer plus tard pendant la factorisation.

Il existe de nombreuses études sur la réduction du coût des solveurs directs en considérant divers types d'approximation. Par exemple, on peut considérer des factorisations incomplètes. L'approche consiste alors à abandonner certaines contributions dans la matrice, soit en utilisant des critères numériques, soit en considérant la structure uniquement (avec la notion de niveaux de remplissage). Alternativement, l'utilisation de techniques de compression de rangs faibles permet généralement d'obtenir de meilleures approximations que la factorisation incomplète. Bien qu'ils aient été introduits plus récemment dans le cadre des solveurs creux, de nombreuses études ont déjà prouvé leur ecacité en termes de réduction du temps et de la mémoire pour une précision donnée.

Les techniques de compression de rang faible visent à regrouper les inconnues de la matrice et à approximer les blocs résultants en fonction d'un critère d'admissibilité. Certaines approches regroupent les inconnues de la matrice de manière hiérarchique comme les schémas H, HODLR, H 2 ou HSS, tandis que d'autres adoptent un regroupement non hiérarchique comme l'approche "block low-rank" (BLR). Bien que les schémas hiérarchiques orent de meilleures complexités asymptotiques, un regroupement non hiérarchique est plus simple à mettre en ÷uvre et plus exible dans les environnements parallèles.

Dans cette thèse, nous avons réalisé toutes nos expériences dans le cadre parallèle du solveur PaStiX qui implémente l'approche non hiérarchique BLR. Comme dans la plupart des solveurs directs, la renumérotation des inconnues de la matrice est obtenue en utilisant la méthode de dissection emboîtée. Cette méthode permet de préserver une grande partie de la structure creuse lors de sa factorisation et permet d'appliquer les opérations en parallèle. Dans le cadre des solveurs utilisant une méthode supernodale, le temps de résolution et l'utilisation de la mémoire sont dépendants du moment où les blocs sont compressés. Lors de travaux antérieurs, deux variantes du solveur PaStiX ont été implémentées :

Minimal Memory et Just-In-Time. Dans la première variante, tous les blocs admissibles sont compressés avant de commencer la factorisation. Ainsi, l'utilisation de la mémoire du solveur est réduite de manière optimale. Cependant, l'application des contributions de rangs faibles est coûteuse et peut conduire le solveur à être plus lent que la version sans compression. Dans la seconde variante, les blocs sont compressés seulement après avoir reçu toutes leurs contributions. De cette façon, le temps de résolution est amélioré sensiblement, mais sans réduire de manière signicative la contrainte mémoire. Il est alors nécessaire de trouver une solution intermédiaire entre les deux stratégies Minimal Memory et Just-In-Time.

Les contributions de cette thèse consistent à déterminer le noyau de compression le plus rapide pour obtenir l'approximation la plus ecace des blocs et à fournir une solution hybride entre les stratégies Minimal Memory et Just-In-Time existantes. Grâce à ces améliorations, nous visons à optimiser les besoins en temps et en mémoire des solveurs directs utilisant des approximations de rangs faibles.

Dans le chapitre 1, nous discutons du contexte et des notions nécessaires pour introduire les contributions de cette thèse. Nous commençons le chapitre en expliquant en détail le fonctionnement des solveurs directs creux. Ensuite, nous présentons la condition d'admissibilité des blocs et discutons des représentations hiérarchiques et non-hiérarchiques. Nous mettons l'accent sur l'utilisation du format BLR dans les solveurs creux. Ensuite, nous expliquons le problème de l'application des mises à jour pendant la factorisation dans le contexte d'une méthode supernodale, et plus spéciquement dans le solveur PaStiX, avec les approches Minimal Memory et Just-In-Time. Nous concluons ce chapitre en discutant des approches développées dans d'autres solveurs et de notre positionnement.

Dans le chapitre 2, nous étudions tout d'abord quelques noyaux d'approximation (compression) de blocs an de déterminer le plus ecace en termes de stockage et de temps de compression. Ainsi, la décomposition en valeurs singulières (SVD) permet d'atteindre le stockage optimal à une précision donnée. Cependant, elle est très coûteuse à utiliser en pratique. Nous nous concentrons donc sur des méthodes plus rapides qui fournissent des résultats proches de la SVD. Bien qu'il existe de nombreuses alternatives pour cet objectif, nous ne considérerons que ceux qui répondent à nos besoins et qui sont basés sur la factorisation QR. Ici, nous implémentons toutes ces méthodes de manière homogène pour une meilleure comparaison. En particulier, nous appliquons un critère d'arrêt pour déterminer numériquement la qualité de la compression, ce qui est appliqué de manière originale pour deux des méthodes étudiées dans la littérature. Cette compression basée sur un seuil numérique garantit ainsi une précision prédénie. Dans ce chapitre, nous nous intéressons à la compression de matrices dont la taille cible correspond aux tailles des blocs apparaissant dans notre solveur BLR. De plus, comme PaStiX utilise des noyaux séquentiels pour la compression, nous considérons un environnement également séquentiel. Les résultats numériques montrent que les noyaux de compression utilisant des techniques de projections aléatoires peuvent être compétitifs par rapport à la méthode SVD. En moyenne, toutes les méthodes ont fourni des résultats assez proches de la SVD. Comme la stabilité et la précision sont assurées par la factorisation QR et notre critère d'arrêt, nous déterminons les noyaux les plus appropriés en comparant les temps de compression. Nous avons observé qu'il est nécessaire d'adapter les noyaux de compression en fonction de la taille de la matrice. Ainsi, nous avons proposé d'utiliser la méthode QRCP pour les matrices de petite taille, tandis que la méthode TQRCP est plus appropriée pour les matrices plus grandes.

Dans le chapitre 3, nous avons amélioré une heuristique préexistante, qui réordonne les inconnues de la matrice après le processus de dissection emboîtée an d'améliorer la compressibilité et permet de réduire le nombre total de mises à jour dans lors de la factorisation. Grâce à ce réordonnancement, une identication des blocs peu compressibles est également eectuée. Dans ce chapitre, nous avons optimisé cette heuristique en utilisant une conguration déterminée empiriquement. Plus précisément, nous vérions le besoin de compresser les blocs peu compressibles juste après l'application de toutes les mises à jour (comme avec la stratégie Just-In-Time). Ainsi, cette heuristique conduit à une solution intermédiaire entre les stratégies Just-In-Time et Minimal Memory puisque les blocs hautement compressibles restent compressés avant de commencer la factorisation (comme avec la stratégie Minimal Memory). La compressibilité par rapport à l'heuristique originale est donc augmentée sensiblement. Dans les expériences, nous avons pu mettre en évidence des gains en termes de ops et de temps par rapport à l'heuristique originale.

Dans le chapitre 4, nous avons étudié une autre heuristique pour identier les blocs peu compressibles an de tirer prot des stratégies Just-In-Time et Minimal Memory simultanément. Ici, nous utilisons la notion de niveaux de remplissage utilisée pour les factorisations incomplètes. Dans cette heuristique, nous attribuons une valeur de niveau à chaque bloc. On suppose alors que les blocs présentant des niveaux élevés sont fortement compressibles, tandis que les autres représentent des blocs à forte interaction et sont peu compressibles. Les expériences que nous avons menées sur un large ensemble de matrices issues d'applications réelles ont démontré que notre heuristique parvient à identier ecacement les blocs les plus compressibles. La solution proposée s'exécute jusqu'à 5, 2 fois plus rapidement que Minimal Memory avec seulement une augmentation d'un facteur 1, 38 sur l'utilisation de la mémoire, en considérant la plus haute précision dans les environnements séquentiels et multithreads. Dans l'environnement multithread, elle est même jusqu'à 1, 84 fois plus rapide et surpasse la version originale dans la plupart des cas. Nous avons donc montré que, si le critère sur le niveau de remplissage est correctement réglé, notre heuristique améliore systématiquement les stratégies existantes. Elle améliore la factorisation numérique en termes de mémoire et de temps, et elle améliore la scalabilité pour les environnements parallèles. Nous avons observé qu'un critère utilisant les premiers niveaux de remplissage est susant pour obtenir un gain signicatif, et qu'une tendance claire apparaît pour régler ce critère en fonction de la précision souhaitée.

Introduction

Thanks to the rapid advances in computational science, we can now avoid costly and/or dangerous real-life experiments by simulating them. The elds where simulation is used include but are not limited to meteorology, biology, nuclear fusion and aerodynamics. Despite the need of further improvements in computational science, the up-to-date simulations are already able to provide amazingly accurate and fast solutions.

Simulating real life problems in computational science requires the continuous domain to be discretized into points. Then the continuous mathematical equations are discretized at these sample points of the whole domain. The resulting set of equations can be represented into a matrix form, where the coecients stand for the entries of the matrix. In many applications, like when discretizing a heat equation, only the interactions between close points of the discretization are taken into account. Indeed far point interactions are frequently weak, thus the corresponding values in the matrix are negligible. Therefore, the resulting matrix includes many zero terms and is called a sparse matrix. Then, the larger the number of discretized points and the accuracy expectations, the more expensive the cost of computing the solution of the system. Thus, the improvements in this eld aim at solving larger systems in shorter times at a given accuracy.

Through the recent improvements toward exascale supercomputer systems, huge computations can be performed in reasonable times by using massively parallelized operations. Unfortunately, the increase of the computational units in these systems does not lead to a rise in the memory available per core. Therefore, this memory limitation forces the scientists/engineers to not only eciently parallelize the operations but also minimize the memory footprint.

When solving large sparse systems, the existing zero blocks of the sparse matrix can be exploited to reduce the memory and computational complexities. Among all the solutions to solve sparse linear systems, direct solvers serve as the most robust approaches despite their expensive memory and computational requirements. These solvers factorize the matrix into a multiplication of triangular matrices. Then, the resulting triangular systems can be solved to obtain the solution.

Depending on the application requirements, direct methods can be performed dierently. According to the matrix features, these methods can be based on dierent factorization types like LL T , LDL T , and LU , for example. Here, L, U , and D represent a lower triangular, an upper triangular, and a diagonal matrix, respectively. Additionally, depending on how the contributions (updates) are applied during the factorization, they can adopt a supernodal or a multifrontal approach. The former strategy allows to directly apply all the updates and avoid using extra memory to store those updates. On the other hand, the latter one uses temporary storage for the contributions to be able to apply them later during the factorization. Thus, this approach leads to a larger memory usage compared to the supernodal version.

There are many studies on reducing the cost of direct solvers by approximating them. For instance, incomplete factorization is proposed for this purpose. The approach consists of dropping some data in the matrix, either by using the geometry of the structure (using the ll-in levels concept) or in a numerical way (using a threshold). Alternatively, low-rank compression schemes are oered as better approximation approaches than incomplete factorization. Although they were recently introduced, there are already many studies that have proven their eciency in reducing time and memory at a given precision.

Low-rank schemes aim to cluster the matrix and approximate the resulting admissible blocks according to an admissibility criterion. Some approaches cluster the matrix unknowns in a hierarchical way like H, HODLR, H 2 or HSS schemes, while some of them adopt a non-hierarchical clustering like block low-rank (BLR). Although the hierarchical ones provide better asymptotical complexities, a non-hierarchical clustering is simpler to implement and more exible in parallel environments.

In this thesis, we conduct all our experiments through the parallel framework of the sparse supernodal BLR solver PaStiX. In this solver, the unknowns numbering (ordering) of the matrix is obtained by using the nested dissection method. This method ensures a good sparsity in the factorized matrix and allows parallel operations. We use the existing block low-rank (BLR) features of the solver to compress suitable dense blocks for an approximated solution.

In the context of supernodal solvers, the time to solution and memory usage is aected by the decision on when the blocks are compressed. In this respect, the sparse supernodal direct solver PaStiX oers two solutions: Minimal Memory and Just-In-Time. In the former one, all the admissible blocks are compressed before starting the factorization. Thus, the memory usage of the solver is highly reduced. However, the costly low-rank contributions can lead the solver to be even slower than the classical one. In the latter one, the blocks are compressed just after they have received all their contributions. In this way, time-to-solution is improved, but it does not reduce memory usage.

We aim to compare the hierarchical schemes to an ecient implementation of the BLR version in the context of a supernodal solver. Thus, in this thesis, we focus on improving the BLR scheme in terms of memory usage and/or time-to-solution. This improvement can be seen as a rst step before the hierarchical implementation.

In Chapter 1, we discuss all the necessary background for the remainder of the thesis. Here, we rst explain the sparse direct solvers in detail. Then, we present the low-rank representations with a special emphasis on the use of the BLR format within sparse solvers. Afterward, we explain the problem of applying the updates during the factorization. Here, we specically focus on the solutions within the supernodal PaStiX solver. We conclude the chapter by discussing other solvers and our positioning.

In Chapter 2, we rstly study some block approximation (compression) kernels to determine the most ecient one in terms of storage and compression time. Although there are a lot of dierent candidates for this aim, we compare only a few that meet our needs, which are based on stable QR factorization. Here, we implement all these existing methods similarly for a better comparison. In addition, we apply a stopping criterion to numerically determine the compression amount, which is the rst time applied for two of the studied methods in the literature. This threshold-based compression guarantees a predened precision. In this chapter, we are interested in compressing matrices with a similar size of those arising in our BLR solver. Moreover, as PaStiX uses sequential kernels for the compression purpose, we work in a sequential environment.

In Chapter 3, we improve an existing heuristic, which reorders the matrix entries after the nested dissection process to improve the compressibility and reduce the total number of updates in the solver. Through this reordering, an identication of the poorly compressible blocks is also performed. In this chapter, we rst optimize the original heuristic through an empirically determined conguration. More specically, we verify the necessity of compressing the poorly compressible blocks just after all updates are applied to them (like in Just-In-Time). Thus, this heuristic leads to an intermediate solution between Just-In-Time and Minimal Memory strategies as the highly compressible blocks are compressed before starting the factorization (like in Minimal Memory). For this reason, we propose to further improve the compressibility compared to the original heuristic because we compress both the highly and the poorly compressible blocks. Then, we observe our improvements in terms of ops and time.

In Chapter 4, we focus on another heuristic to identify the poorly compressible blocks to take advantage of Just-In-Time and Minimal Memory strategies at the same time.

Here, we use the ll-in levels concept of incomplete factorization. In this heuristic, the blocks with high levels are well compressible, while others represent strong interaction blocks and are poorly compressible. We prove the correctness of our heuristic through the compressibility results of dierent ll-in level blocks. Then, in both parallel and sequential environments, we empirically validate the eciency of our heuristic in terms of memory, ops, and time. Finally, we conclude our work and discuss our perspectives. This work is supported by the Agence Nationale de la Recherche, under grant ANR-18-CE46-0006 (SaSHiMi). Experiments presented in this thesis 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 (https: //www.plafrim.fr/).

Chapter 1

State of the Art

Many scientic applications need to solve large sparse linear systems of the form

Ax = b, (1.1)
where the matrix A is of size n × n. Similarly, the unknown vector (x) and the right-hand-side vector (b) are of size n. Here, we can informally dene a sparse matrix to be "any matrix with enough zeros that it pays to take advantage of them" [START_REF] Hardy | The algebraic eigenvalue problem[END_REF].

Sparse linear systems can be solved through many dierent approaches, which are divided into two main categories [START_REF] Golub | Matrix computations[END_REF]: direct and iterative solvers. Additionally, there are some hybrid methods [GH08, AGG + 13, YL10, RBH12], which benet from both approaches at the same time.

Iterative methods [Saa03] are usually used for very large systems thanks to their lower memory usage and computational cost compared to the direct methods. However, they frequently need a convenient preconditioner to converge to the solution. Therefore, nding an iterative method suitable for arbitrary systems is very challenging. Alternatively, direct methods [DER86, GL81, GHLN88] provide the most robust and reliable solutions. Yet, they are very costly in terms of memory usage and computational complexity. From now on, we only focus on improving the high cost of direct solvers.

All the applications have their own precision requirements and computations at the machine precision is not necessary. As we frequently need to deal with matrices that can be approximated into much lower dimensional blocks, some cheaper (approximated) direct solver alternatives are proposed. For example, incomplete factorization is proposed to neglect some values of the matrix to reduce the solver cost. In this approach, the dropped values can be chosen according to a numerical criterion (threshold) [START_REF] Karypis | Parallel threshold-based ilu factorization[END_REF] or based on the geometrical structure (ll-in levels) [START_REF] Timothy | An unsymmetric-pattern multifrontal method for sparse LU factorization[END_REF]. These solvers are mostly suitable as a preconditioner for the iterative solvers. Recently, low-rank approaches were introduced with better approximation criteria compared to the incomplete factorization. As low-rank approaches are convenient for high precision approximations, they can be used as a cheaper direct solver alternative than the classical version (with some accuracy trade-o). In this thesis, we focus on a sparse low-rank solver, which can be used either as a cheap direct solver, or as a good preconditioner for the iterative methods.

In Section 1.1, we explain the sparse direct solvers. In Section 1.2, we discuss the background on the low-rank schemes. In Section 1.3, we explain the problem of applying the contributions within the low-rank sparse solvers. Then, we provide the solutions of the PaStiX solver [Pic18,PDF + 18] for this problem in Section 1.4. Finally, in Section 1.5, we discuss other solvers to show our positioning compared to them and state our objectives in this thesis.

Sparse Direct Solvers

A straightforward way of nding the x vector in Equation 1.1 is to compute it as x = A -1 b. However, not only does this approach requires A to be invertible, but also inverting the sparse matrix A does not guarantee to conserve its sparsity. In this respect, direct methods can be used to solve these systems in high accuracy.

In this work we focus on Gaussian elimination based direct solvers. These solvers rstly factorize the matrix into a product of triangular matrices depending on the matrix features. For example, symmetric positive denite (SPD) matrices can be factorized through LL T (Cholesky), while symmetric matrices can use LDL T factorization. For the general cases the LU factorization is adopted. Here, L, U and D represent respectively a lower triangular, an upper triangular and a diagonal matrix. Then, following the factorization of the matrix, the resulting triangular systems can be solved to obtain the original solution.

Sparse direct solvers follow four main steps: 1. Find the ordering This step aims to increase the number of zeros in the factorized matrix structure (L and U matrices for the general case) and improve the parallelism.

Perform symbolic factorization

This step allows to allocate the factorized matrix structures before starting the numerical operations.

Factorize

This step factorizes the matrix in-place as the factorized matrix structure is known.

Solve

This step solves the triangular systems. Let us explain each step in more detail in the following sections.

Ordering

After the factorization of a sparse matrix, A, some structural zero values become non-zero in the factorized matrix. These non-zero values that are zeroes in the original matrix are called ll-in. It is trivial to see that reducing the ll-in amount in the sparse matrices improves the storage and computational complexity. For this purpose, an ordering step is performed before the factorization in sparse direct solvers.

Let us introduce some necessary background through Figure 1.1 to understand the ordering step better. In the left side of the gure, we can see the original matrix, A, of order 5 on the left. This matrix can be represented with an adjacency graph, which can be expressed as G(A) = (V, E). Here, V and E represent the graph vertices and edges, respectively. Each variable i in the matrix illustrates the vertex i of the graph. Note that for the unsymmetric matrices, we use the symmetric pattern A + A T for more ecient operations. Then, for each a ij ̸ = 0 or a ji ̸ = 0 for i > j, there is an edge in E, which is represented as a pair (i, j).

On the right side of Figure 1.1, we can see the information about the factorized matrix.

Here, the lled elements on the factorized matrix, L, are represented in red. Similarly, the lled graph G * illustrates both the original and ll-in edges after the factorization, where red color stands for the ll-in edges. The right-most structure, T , illustrates the dependencies when eliminating the unknowns during the factorization. This structure is called elimination tree [START_REF] Schreiber | A new implementation of sparse gaussian elimination[END_REF]. In this tree, the elimination of the unknowns is performed from bottom to top. Here, the siblings are eliminated in parallel since the ll-in occurs only from a node to the ancestors. Now let us explain the importance of the unknown numbering (ordering) through Figure 1.1. As seen at the top row of the gure, the numbering leads to a dense matrix after factorization due to the introduced ll-in. In addition, we can see that each node is eliminated in sequential from bottom to top in the corresponding elimination tree, T . In the bottom row of the gure, the ll-in is avoided through the new numbering of the unknowns. This clever numbering also allows parallelism when eliminating the unknowns, as seen in its elimination tree. That is, the nodes from 1 to 4 can be eliminated in parallel, since there is no dependency between them. On the other hand, the node 5 is eliminated after the other unknowns since it depends on them.

Through Figure 1.1, we can observe the importance of the ordering since the ll-in introduced for large systems can be very crucial and an ecient usage of parallel architectures is necessary. Therefore, we need to order the matrix conveniently, before the factorization. In this respect, the ordering can be represented by a permutation matrix P . Then, instead of solving Equation 1.1, we can equivalently solve (P AP T )P x = P b.

We can order a sparse matrix through some methods like approximate minimum degree (AMD) [START_REF] Amestoy | An approximate minimum degree ordering algorithm[END_REF], minimum local ll (MF) [START_REF] Luce | On the minimum ops problem in the sparse cholesky factorization[END_REF][START_REF] Tinney | Direct solutions of sparse network equations by optimally ordered triangular factorization[END_REF] or nested dissection [START_REF] George | Nested dissection of a regular nite element mesh[END_REF] to successfully reduce the ll-in. The nested dissection method is a classical method in sparse solvers as it is also convenient for providing parallelism. It can be performed through well-known partitioning libraries like Metis [START_REF] Karypis | Metisa software package for partitioning unstructured graphs, partitioning meshes and computing ll-reducing ordering of sparse matrices[END_REF] or Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF]. We adopt this method from now on.

In Figure 1.2, the nested dissection partitioning is illustrated on an adjacency graph. This method recursively divides the graph (and sub-graphs) into two balanced sub-parts through a separator of minimal size. At each recursion, this method separates the sub-parts well to avoid direct interaction between them, so that they are only connected through the separator. As a result, the sub-part interaction blocks become zero and these sub-parts do not depend on each other.

5 4 3 2 1 1 2 3 4 5 1 1 2 1 1 3 1 0 1 4 1 0 0 1 5 1 0 0 0 1 A 3 2 1 4 5 G 1 2 3 4 5 1 1 2 1 1 3 1 0 1 4 1 0 0 1 5 1 0 0 0 1 L 3 2 1 4 5 G * 5 4 3 2 1 T 1 2 3 4 5 1 1 2 1 1 3 1 0 1 4 1 0 0 1 5 1 0 0 0 1 A 3 2 5 4 1 G 1 2 3 4 5 1 1 2 1 1 3 1 0 1 4 1 0 0 1 5 1 0 0 0 1 L 3 2 5 4 1 G * 5 4 3 2 1 T Figure 1.1:
The impact of ordering on a sparse matrix with ve unknowns. On the left, we see the matrices to be factorized (A) with the corresponding graphs (G). On the right, we see the factorized matrix structures (L) with the corresponding lled graph (G * ) and elimination tree (T ). The red color represent the ll-in. The ordering at the top generates no ll-in, while the bottom ordering causes the matrix to be fully lled.

Block symbolic factorization

After the nested dissection partitioning, the symbolic factorization is performed to know the data structure of the factorized matrix before starting any numerical operations. For example, in the right part of Figure 1.2, we can see the data structure of the factorized matrix after the nested dissection ordering. Here, the ll-in which is caused by the separator interactions with the other sub-parts are shown in red and green. The upper red ll-in represents the interaction of the rst and second sub-parts with the third sub-part (separator), while the bottom red ll-in illustrates the interaction of the forth and fth sub-parts with the sixth sub-part (separator). Similarly, the green ll-in represents the interaction of all the sub-parts with the seventh sub-part (largest separator). In practice, this information is obtained through the matrix graph. When a node of the graph is eliminated during the factorization, its neighbors are connected to each other (clique). These new connections are represented as new edges in the lled graph. Then, this information about the factorized matrix allows to allocate the data structures before factorization and collect the data in the blocks for ecient BLAS Level 3 [START_REF] Dongarra | A set of level 3 basic linear algebra subprograms[END_REF] operations.

Each sub-graph, including the separators, obtained after the nested dissection is called supernode. We know that each supernode resulting from the nested dissection process is designed in a way that the unknowns which have similar non-zero structures in the factorized matrix are gathered together.

In Figure 1.2, we can see each supernode and corresponding intra-supernode interactions in each diagonal matrix (in gray). Here, there are seven supernodes resulting Adjacency graph (G). Elimination tree (T).

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Factorized matrix (L). from the nested dissection partition. The red and green dense blocks represent the inter-supernode interactions. Here, each vertical panel consisting of one supernode's diagonal block and its interactions is called column block (or simply panel). Then, both the quotient graph and the supernode elimination tree can be obtained through the symbolic structure of L in Figure 1.2. Here, the quotient graph illustrates the supernode interactions, while the latter one shows the dependencies. Note that the top-most separators of large matrices can cause huge computations. Thus, in practice, they are clustered into smaller sub-graphs (so column blocks) before each panel is assigned to the processors [START_REF] Hénon | PaStiX: A High-Performance Parallel Direct Solver for Sparse Symmetric Denite Systems[END_REF]. In this way the parallelism is increased.

Numerical factorization

The numerical factorization is performed after the preprocessing steps that are explained previously. In order to take advantage of the ecient BLAS Level 3 operations, this factorization works on the blocks of each column block of the sparse matrix. That is, the factorization basically iterates on each panel, where three main steps are applied at each iteration: factorize, solve and update. In the rst step, the diagonal block of the current panel (green in Figure 1.3) is factorized. In the second step, other blocks of the current panel (o-diagonal blocks) are solved. The aected blocks of this step is colored in blue in Figure 1.3. Lastly, the contributions of the current panel to other column blocks are performed. The updated blocks of other panels are colored in red in the gure. Numerical factorization can be performed in two dierent ways: following the multifrontal [START_REF] Du | The multifrontal solution of indenite sparse symmetric linear[END_REF] or the supernodal method. These approaches dier in terms of data locality and parallelism.

In the multifrontal method, the children update only their direct parents in the elimination tree. Therefore, they use some temporary storage (called fronts) for accumulating the contributions to be applied to the red blocks of Figure 1.3.

Opposite to the multifrontal method, supernodal approaches apply the update operations to any level of ancestors in the elimination tree. The supernodal method can be applied in two dierent ways depending on when the updates are performed: left-looking and right-looking approaches. The former approach updates the current supernode blocks just before eliminating it through the contributions from any level of descendants. On the other hand, the right-looking strategy allows to directly apply all the updates to any level of ancestors in the elimination tree, as soon as the contributions are computed. Supernodal methods avoid the extra storage used in the multi-frontal approach for storing the fronts, thus reduce the memory peak compared to it.

In this work, we focus on the right-looking supernodal approach. The right-looking panel-wise Cholesky algorithm can be seen in Algorithm 1. Note that we compute the numerical factorization in-place since memory is a problem for large matrices.

Triangular solves

After the factorization of the matrix into two triangular matrices, we have two steps to solve the resulting system. That is, considering a LU factorization on the matrix A, we end up with LU x = b. Therefore, we rst solve Ly = b (forward substitution) and then we solve U x = y (backward substitution). Note that in this thesis we do not focus on the triangular solve steps.

Low-rank representations

As mentioned before, sparse direct solvers are very costly in terms of memory and time consumption. Nowadays, the introduction of low-rank representations has enabled the use of these solvers for larger systems. These representations basically aim to approximate the matrix. They represent matrices as multiplication of two matrices as A m×n ≈ U m×r V T r×n , where r (called rank) is much smaller than min(m, n). As a result, the memory and computational complexity of the solver can be reduced. The rank can be determined through a predened criterion based on either a xed-rank or a xed-precision to satisfy

||A -U V T || ||A|| < ϵ, (1.2)
where ϵ is the precision. At a given precision, the optimal approximation can be computed through the singular value decomposition (SVD). However, this method is too costly to be applied in practice. Therefore, some QR factorization based approaches are oered to reduce this cost, while providing close approximations to SVD. Especially, through the randomized methods, computations on large matrices are avoided by using small representative data [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. As a more detailed explanation on the approximation (compression) kernels is provided in Chapter 2, we skip further information here. Many recent studies have tackled the problem of reducing the memory consumption of linear solvers with low-rank compression. For the sparse solvers, low-rank techniques are more challenging than for dense solvers [CPA + 20, ALS + 18, ALMK17] as the low-rank representations on the full matrix do not exploit the existing sparse structure. That is, after the nested dissection procedure, large zero blocks are introduced. Then, these zero blocks should be ignored to reduce memory and computational complexity. Therefore, in this thesis we focus on a low-rank solver that takes advantage of the existing sparsity of the matrix after the nested dissection and compresses the dense matrix parts independently.

In order to obtain the low-rank representation of a matrix, some hierarchical or non-hierarchical (at) block clustering formats can be used on the dense parts of the matrix. Then, according to a block admissibility criterion, the compressible o-diagonal blocks can be approximated, while the non-admissible ones are kept uncompressed. Let us rst explain the block admissibility condition in 1.2.1. Afterward, we provide the necessary background on the low-rank compression formats in Section 1.2.2, and continue with a more specic case (a at representation, which is called block low-rank) in Section 1.2.3.

Block compression admissibility

In the applications (like the ones with PDE discretization), when the matrix is represented in blocks, some of the blocks can be represented in low-rank form [START_REF] Chandrasekaran | On the numerical rank of the o-diagonal blocks of schur complements of discretized elliptic pdes[END_REF]. For this purpose, some block admissibility criteria can be used to determine the compressibility of the blocks.

To decide which blocks to compress, Hackbusch [START_REF] Hackbusch | Hierarchical Matrices: Algorithms and Analysis[END_REF] dened some admissibility conditions. More specically, he proposed both a weak and a strong admissibility condition.

In order to explain the admissibility criteria, let us rst dene a block with the unknown set σ as row indices and τ as column indices. Then, this block stands for the interaction of the sets σ and τ . Therefore, if σ = τ , this block is a dense diagonal block, which represents the interactions within the same set. If σ ̸ = τ then the block can be represented in low-rank form if the admissibility criterion is respected.

According to the weak admissibility criterion, all the o-diagonal blocks are approximated:

σ × τ is admissible if f σ ̸ = τ.
(1.

3)

The strong admissibility condition relies on the problem geometry and the denition of the diameter of a set of unknowns (diam(σ)), as well as the distance between two sets (dist(σ, τ )). The interaction between two sets of unknowns is then considered admissible if the distance between the sets is suciently larger than both of the diameters of the sets

σ × τ is admissible if f max(diam(σ), diam(τ )) ≤ ηdist(σ, τ ), (1.4)
where η represents a xed value. As we can see, as the distance between the clusters increases, the rank gets smaller since their interaction gets weaker. The least restrictive strong admissibility criterion considers that blocks are admissible only if their distance is larger than 0 (i.e. all blocks except the close neighbors).

In order to exploit the low-rank features, it is necessary to correctly determine the blocks that are worth to be approximated. Thus, in a fully-algebraic context, we need to adopt an algebraic alternative to the admissiblity criteria dened by Hackbusch [START_REF] Hackbusch | Hierarchical Matrices: Algorithms and Analysis[END_REF]. In this respect, we will use the adjacency graph information (for diameter and distance) [START_REF] Grasedyck | Parallel black box h-lu preconditioning for elliptic boundary value problems[END_REF] to determine the block admissibility, instead of using the underlying geometry of the problems. Now let us discuss the clustering approaches to obtain the blocks on which we check the admissibility conditions.

Low-rank compression formats

As mentioned before, in order to approximate the matrices, the matrix unknowns are clustered, so that only the admissible blocks can be compressed. In Figure 1.4, we illustrate an example of the matrix unknown clustering on the left, while we see the corresponding cluster tree on the right. Here, I stands for all matrix unknowns to be clustered. A matrix clustering can be performed by using dierent formats. For example, Hierarchical (H) matrices [START_REF] Hackbusch | A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices[END_REF] are introduced to cluster the unknowns in a hierarchical way. Since this format adopts a strong admissibility, another hierarchical format called HODLR [START_REF] Aminfar | A fast, memory ecient and robust sparse preconditioner based on a multifrontal approach with applications to nite-element matrices[END_REF][START_REF] Chadwick | An Ecient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization[END_REF] is proposed as a weak admissibility criterion based alternative to it. In addition, other hierarchical formats like H 2 [HB02] and Hierarchically Semi-Separable (HSS) [GLR + 16,Xia13a] are oered to take advantage of the nested bases. That is, let us consider a low-rank representation of a block in the form U BV T , where U and V are the orthogonal bases of this block. Then, in these solvers, the bases of a block are based on the bases of the descendants in the corresponding cluster tree. In this way, the storage cost of the solvers is improved. However, as H 2 , which uses a strong admissibility, and HSS, which adopts a weak admissibility, depend on the existence of the nested bases, they are more restrictive compared to H and HODLR formats. As an alternative to the hierarchical formats, some non-hierarchical (at) formats are also proposed. These formats can be seen as a special case of the hierarchical formats (illustrated in Figure 1.4), where all the unknown clusters are the direct children of I as seen in Figure 1.5. Here, we see the cluster tree (on the right) of a simple at block clustering example (on the left).

I 2 1 1 I 2 1 2 I 2 2 1 I 2 2 2 I 2 1 1 I 2 1 2 I 2 2 1 I 2 2 2 I 1 1 I 1 2 I 1 1 I 1 2 I I 1 1 I 2 1 1 I 2 1 2 I 1 2 I 2 2 1 I 2 2 2
Flat representations have larger asymptotical complexities in terms of ops and memory compared to the hierarchical schemes. However, they are simpler to implement and more convenient to parallelize. For example, arranging the block sizes small enough for a shared memory implementation or having good load balance in a distributed environment is easier with these formats compared to the hierarchical ones. Therefore, in this thesis we focus on these at representations, more specically the block low-rank (BLR) format, and explain it in the following section.

Block Low-Rank (BLR) format

Full Rank Low Rank Figure 1.6: Block Low-Rank (BLR) compression on a dense matrix. Brown color represents the dense blocks, blue color represents the compressed matrices and white color stands for the sparsity generated thanks to the compression.

The BLR representation [AAB + 15] of a dense matrix is illustrated in Figure 1.6. Here, we adopt a weak admissibility criterion for the sake of simplicity. In the gure, the dense diagonal blocks are kept in full-rank (in brown) as they represent strong interactions. However, the o-diagonal blocks are compressed into two lower dimensional blocks (in blue) to generate sparsity of the matrix. If strong admissibility criterion was adopted, some of the o-diagonal blocks could be in full-rank. As seen in the gure, BLR format provides similar size blocks.

As mentioned before, to take advantage of the existing zero blocks of the sparse matrices, the BLR representation can be applied only on the dense parts of the matrix, instead of the full matrix. In this way, unnecessary memory usage and computational complexity of the zero blocks are avoided.

We already explained that we obtain the column blocks of the factorized matrix through the nested dissection process, where each column block represents a supernode. Then, the largest separators are also clustered to reduce their sizes as explained in Chapter 3. This procedure results in the BLR representation of the dense parts of the matrix, where the blocks are of similar sizes. In addition, the large supernode interaction blocks are also suitable to be represented in low-rank. In practice, small Figure 1.7: Block Low-Rank (BLR) compression on the dense blocks of a sparse matrix. Brown color represents the dense blocks, blue color represents the compressed matrices and white color stands for the zero parts. blocks and high rank blocks, according to predened criteria, are kept dense to avoid unnecessary compression. We illustrated this kind of BLR representation of a sparse matrix in Figure 1.7. As we see in the gure, the large zero blocks (in white) do not need storage and operations. Now let us specically give the necessary background on the PaStiX solver, which is our focus in this thesis and adopts the BLR scheme. Here, the contributing blocks A and B appear respectively in red and in blue. The updated block C is colored in orange and the impact of the contribution (AB) is in purple.

Block compression in the dense solvers and the sparse multifrontal solvers provides lower update ops and improves the memory footprint. However, sparse supernodal solvers use dierent block sizes in the update operation, where the low-rank cost depends on the largest block. This may result in even more expensive updates than the full-rank one. Nonetheless, supernodal approaches provide more parallelism and have less memory overhead compared to the multifrontal ones.

Figure 1.8 describes two dierent low-rank update (C -= AB) options in a supernodal solver, which will be our focus in this thesis. Here, the updated block (C) is large and the contribution (AB) is small. The A, B and C blocks are respectively in red, blue and orange. The contribution, AB is shown in purple. M and N represent the dimensions of the updated matrix, while m and n stand for the dimensions of the contribution.

The contributing blocks for both strategies in Figure 1.8 are represented in low-rank. Although these blocks can also appear in full-rank within the solver, this is out of our scope. The updated block, C, is represented in full-rank on the left, whereas it is in low-rank on the right.

The non-fully structured strategy (LR2FR) in Figure 1.8a, computes the contribution, AB, cheaply by taking advantage of the low-rank representation of A and B. Additionally, as C is stored in full-rank, a simple addition of the contribution at a cost of order mnr (r being the rank of the contribution) is performed, further improving the computational complexity of the update operation. However, this approach does not reduce the memory usage as all the blocks are rstly allocated in full-rank to take advantage of the cheap addition of order mnr. The fully structured updates (LR2LR) [START_REF] Xia | Ecient structured multifrontal factorization for general large sparse matrices[END_REF] are illustrated in Figure 1.8b. Here, all the blocks are compressed before starting the operations. As a result, the memory usage is reduced. In addition, the contribution (AB) is computed at a low cost similarly to the non-structured version. However, the addition step into C requires a complex low-rank to low-rank update with padding (zeroes are added to match the dimension of C). Let us illustrate the padding through Figure 1.9.

In Figure 1.9, the accumulation on the U and V bases in the structured approach is shown. Here, we observe that after the accumulation, some zero padding is used to align the contribution block (small) and updated block (large) vectors. Then, after the accumulation of several updates, the rank and size of C can highly increase. Thus, an additional re-compression is necessary to maintain a small storage and rank [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF]. This LR2LR operation with re-compression in the context of the supernodal method generates a ops overhead. As a consequence, the factorization time may be highly impacted, and can be even more costly than the full-rank solver, unless the matrix is highly compressible. Now let us introduce the solutions of the sparse supernodal solver PaStiX to the low-rank update problems we discussed.

PaStiX sparse supernodal direct solver

In our work, we conduct all of our experiments through the BLR supernodal sparse direct solver PaStiX [Pic18, PDF + 18]. This solver targets both the real and complex matrices through the LU , LL T , LDL T , LL H and LDL H factorizations. PaStiX allows both static and dynamic shared memory usages, as well as external Parsec [BBD + 13] and StarPU [START_REF] Augonnet | StarPU: a unied platform for task scheduling on heterogeneous multicore architectures[END_REF] run-time supports to exploit accelerators. As this is not our main focus in this thesis, we will not detail it here. Interested readers can nd further information in [Fav09, LFB + 14]. Now let us explain the necessary background for the following chapters.

In PaStiX, two dierent compression strategies that dier depending on the update kernels are implemented to lower either the memory footprint or time-to-solution:

Just-In-Time and Minimal Memory.

Update( A ik , A jk , A ij )
The Just-In-Time strategy in Figure 1.8a aims to reduce time-to-solution through late compression of C. In this strategy, the blocks are compressed only after receiving all the contributions on them, on the y, during the factorization. This operation is represented in blue in Algorithm 2. Then, through the resulting non-structured updates (LR2FR), a simple addition of the contribution is performed, reducing the computational complexity of the updates. However, as all the blocks are allocated in full-rank in the beginning of the factorization, this scenario does not improve the memory peak compared to the full-rank solver.

Update( A ik , A jk , A ij )
The Minimal Memory strategy in Figure 1.8b aims to reduce the memory footprint of the solver by early compression of C. This operation is illustrated in red in Algorithm 3. Indeed, as all admissible blocks are compressed before starting any numerical operation, it greatly improves the memory peak since the factorized matrix structure is never fully allocated. However, the fully structured updates (LR2LR) in the Minimal Memory strategy are complex to perform as explained before. A detailed complexity analysis of the Just-In-Time and Minimal Memory strategies is provided in [PDF + 18]. As we can see, neither Just-In-Time nor Minimal Memory is optimal. In order to take advantage of both the strategies in PaStiX, it is important to come up with a good preselection criterion to distinguish the poorly and highly compressible blocks. In this way, the former ones can be compressed on the y (like in Just-In-Time), while the latter ones are compressed before starting any numerical operation (like in Minimal Memory).

As a result, we can improve both the ops and memory overhead of the solver in an optimal way. Now let us observe our positioning in more detail compared to the other works.

Other solvers

Hierarchical matrices have been rst introduced by Hackbusch [Hac99] for dense matrices. Since this initial work, there has been a lot of interest on this topic from variety of scientic and engineering domains. In this work, the author expressed the full matrix in the H format to approximate it and then proposed fully structured low-rank operations. In the context of sparse solvers, expressing the full matrix in a low-rank format does not take advantage of the existing sparsity. Thus, we approximate only the dense blocks of the sparse matrix.

In [START_REF] Grasedyck | Domain decomposition based h -lu preconditioning[END_REF], an extension of the H matrices to the sparse case is proposed. Here, the geometry of the problem is used to obtain the clustering tree. In this work, the authors adopted the nested dissection partitioning to have large zero blocks in the matrix, so that only the remaining dense parts are approximated. In this way, they reduced the factorization cost of the sparse structured H -LU solver. The drawback of their approach is that they do not eciently exploit the structural zeros: they ignore the zeros of the sub-part to separator interaction. That is, let us consider the rst level separator (green) in the graph of Figure 1.2 and the well-separated sub-graphs generated by this separator. Then, in this gure, the zero blocks in the last ve rows of the corresponding L matrix is not taken into account in [START_REF] Grasedyck | Domain decomposition based h -lu preconditioning[END_REF]. Therefore, they skip some memory and ops reduction. In our work, we perform symbolic factorization after the nested dissection process to determine and exploit all the structural zeros eciently.

A black-box algebraic version of [START_REF] Grasedyck | Domain decomposition based h -lu preconditioning[END_REF] is proposed in [START_REF] Grasedyck | Parallel black box h-lu preconditioning for elliptic boundary value problems[END_REF], where the graph information is used for the clustering without the knowledge of the geometry in the H -LU factorization. Here, they presented the numerical results of both the geometrical and the algebraic H -LU (and H -Cholesky for the SPD case) solver in terms of time and memory. Although the geometrical solver shows better results, the ratio of the algebraic and geometric solvers does not exceed small constants. In addition, the authors also provided the implementation in a parallel context. In the parallel experiments, the geometrical and algebraic solvers demonstrate almost identical speedup results.

In [START_REF] Hackbusch | Hierarchical Matrices: Algorithms and Analysis[END_REF], the adaptation of some hierarchical matrix strategies from dense to sparse matrices are studied. In this work, an approach similar to the Minimal Memory strategy of PaStiX is used, meaning that low-rank updates are performed. Here, all the contributions from the descendants in the tree are considered dense when generating the low-rank representation: structural zeroes are ignored. The assembly is then applied without zero-padding, which allows to reach good performance. In this work, the overhead of low-rank updates does not exist, opposite to our work, as the operations are always performed between blocks of similar sizes. However, it misses an opportunity for memory savings compared to our work since some zeroes are explicitly stored.

In [START_REF] Chadwick | An Ecient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization[END_REF], the authors proposed low-rank format similar to HODLR through the sparse supernodal solver Cholmod [START_REF] Chen | Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate[END_REF] for the SPD matrices. They used xed rank based (depending on the block size) randomized compression kernels both for the o-diagonal blocks of the supernode interactions and the intra-supernode o-diagonal blocks. In practice, determining such a xed rank to keep data with a desired precision is not trivial. In this work, they demonstrated good memory savings compared to the standard direct methods, which makes the low-rank solver suitable for large problems with high memory requirements. Nonetheless, their factorization is slower than the full-rank version and they do not fully take the sparsity of o-diagonal blocks into account. In our work, we allow to determine the rank numerically (with a precision based criterion) and exploit the structural zeroes eciently.

In [START_REF] Xia | Randomized sparse direct solvers[END_REF], HSS approximation is performed on the large fronts of a sparse multifrontal solver through a randomized technique for SPD matrices. The hierarchical low-rank algorithm in this work is extended eciently for general cases in [GLR + 16] through the sparse multifrontal solver StrumPACK. Here, the authors proposed a shared memory implementation. In this work they focused on nding an ecient and robust sparse solver, especially for the discretized PDE systems. They avoided saving the dense matrices explicitly thanks to the adaptive randomized sampling, which in return reduced the memory usage. They showed time improvements compared to the classical multifrontal solver. In [START_REF] Ghysels | A robust parallel preconditioner for indenite systems using hierarchical matrices and randomized sampling[END_REF], the authors improved the same work for distributed memory usage. Here, they focused on using this low-rank solver as a preconditioner.

The algebraic

StrumPACK multifrontal solver provides a fully structured approach, where all the blocks are compressed before the numerical operations. The HSS format in this solver provides better complexity than the BLR format that we adopt. However, opposite to us, the HSS based solver requires the existence of the nested bases.

The block separable (BS) representation in [Gil11] adopts a at block clustering with a weak admissibility. Similarly, the block low-rank (BLR) representation, introduced in [AAB + 15] through the multifrontal sparse direct solver Mumps [START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance, and scalability[END_REF], also uses a non-hierarchical clustering. Here, BLR can be seen as an extended version of BS since it is not restricted to use weak admissibility condition.

In the Mumps solver, the CUFS (Compress, Update, Factor, Solve) technique is a similar approach to the Minimal Memory scenario from PaStiX. Using this technique, low-rank updates are performed within each front, but not between fronts. Thus, as fronts are allocated in full-rank before being compressed, memory savings are less important than the ones of the Minimal Memory strategy. Alternatively, the FCSU (Factor, Compress, Solve, Update) method applied to the fronts is similar to our Just-In-Time scenario.

Here, the updates are applied on the full-rank blocks and the panel compression is done during the factorization. Our Minimal Memory strategy within the supernodal approach enables a larger memory saving by preventing the allocation of the fronts inherent to the multifrontal method.

In [START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance, and scalability[END_REF], the complexity and performance comparison of the multifrontal BLR solver

Mumps and the HSS solver StrumPACK is provided. Here, by using the experimental results, the author oers to use BLR when high accuracy solution is required, while HSS is more convenient for low precisions.

In [START_REF] Marchal | Trading performance for memory in sparse direct solvers using low-rank compression[END_REF], the authors conducted a very similar work to us. Here, a solution to the problem of preselecting the poorly compressible blocks is proposed through the sparse supernodal solver PaStiX. The authors exploit performance models of the update kernel to decide whether or not to delay the compression of some of the blocks. In this work, this decision is taken at run-time during the numerical factorization and requires to generate correct models of the problem. In this initial work, they conducted the work in sequential context, while their nal aim is the parallel environment in the future. In Chapter 3 and in Chapter 4, we focus on nding such a preselection criterion with dierent angles, where we identify the poorly compressible blocks before starting any numerical operations.

In this thesis, we focus on nding algebraic solutions to provide a low accuracy sparse direct solver or high accuracy preconditioner for the widest spectrum of applications. The precision of the solver can be determined by the user. Our aim is to study solutions to reduce the time-to-solution and/or the memory footprint of the BLR solver PaStiX by correctly preselecting the blocks and eciently compressing them.

Chapter 2 Compression Kernels

Many elds like data mining, signal processing and computer vision need to handle large matrices. As the storage and matrix operations are too expensive for these matrices, some approximation (compression) techniques are proposed to reduce the cost. These techniques aim to represent the original matrix, A, in the form

A m×n ≈ U m×r V T r×n (2.1) such that ||A m×n -U m×r V T r×n || < ϵ||A||. (2.2)
Here, the objective is to exhibit an approximation at a required precision (ϵ), where r is much smaller than min(m, n). Unfortunately, determining the most suitable compression technique for arbitrary problems is dicult. Therefore, each application needs to come up with a convenient method according to the requirements like memory, time, stability and accuracy.

As mentioned in Chapter 1, we exploit the existing sparse structure in our sparse solver and split the dense parts of the matrix into blocks (in BLR format). Then, these blocks are approximated by the product of two lower dimensional blocks if they are admissible. In this chapter, our purpose is to compare some interesting compression kernels that can be used on these blocks to improve the solver. More precisely, we want to determine a fast and stable kernel that compresses admissible dense matrix blocks in the form (2.1), with r as small as possible. Here, compressing the data with smaller r values targets to improve both the memory footprint and the op count during the factorization and solve steps.

In Section 2.1, we provide all the necessary background on the compression kernels. In this section, the motivation of dierent kernels, especially the ones on which we will focus, is also explained. The literature review specic to the four compression kernels that we concentrate on is provided in Section 2.2. In Section 2.3, we dene some notations to be used in the algorithms of this chapter. Then, in Section 2.4, four dierent compression methods are studied in detail. In Section 2.5, their complexities are shortly summarized to show the big picture. In Section 2.6, the numerical results of some generated matrices are observed in terms of stability, performance, compression ranks and accuracy, whereas in Section 2.7, the real-life case matrix results are discussed. Finally, in Section 2.8, the conclusion of this chapter is provided.

Background

In order to understand the compression kernels, we rst need to provide the meaning of r in the approximation (2.1). Mathematically, the rank of a matrix (r) stands for the number of linearly independent columns/rows of the matrix. Then, if we represent a matrix as multiplication of two matrices as

A m×n = U m×r V T r×n , (2.3)
r can be determined through the singular value decomposition (SVD), which has the form

A m×n = Ûm×m Σ m×n Vn×n .
(2.4)

Here, Û and V are orthogonal matrices, while Σ is a diagonal matrix with the singular values, σ j , on the diagonal with

σ 1 ≥ σ 2 ≥ ... ≥ σ r > 0 and σ r+1 = σ r+2 = ... = σ n = 0.
Here, the rank is equal to the number of non-zero singular values in Σ.

When we work in nite precision, we are usually interested in only the large values and we omit the small values depending on our precision requirement. In this case, as mentioned in Section 1.2, an approximated rank is determined by using either a xed-rank or a xed-precision based criterion to satisfy the condition (2.2), where r is the rank of the approximation and ϵ stands for the precision. Observe that this rank is numerically computed and it can have dierent values for dierent precisions. From now on, r will only represent approximation ranks in this chapter. Determining a convenient xed-rank to ensure a targeted accuracy is not trivial in practice. Thus, the xed-precision based criterion provides better solutions for abstract problems. We will adopt a xed-precision based criterion in our work.

If the equation (2.4) is written in the block form

A = Û1 Û2 Σ 11 0 0 Σ 22 V1 V2 ,
the approximation (2.1) can be written as

A ≈ U V T = Û1 Σ 11 V1 , (2.5)
where U = Û1 , V T = Σ 11 V1 and Σ 11 is of size r × r.

In this chapter, we will use either

||Σ 22 || 2 ||A|| F = σ r+1 ||A|| F < ϵ (2.6) or || Û2 Σ 22 V T 22 || F ||A|| F = ||Σ 22 || F ||A|| F = n i=r+1 σ 2 i ||A|| F < ϵ (2.7)
to nd the numerical rank r through SVD. The SVD method is not the only way to compress the matrices in the form (2.1). Depending on the application requirements, dierent compression methods can be adopted. For example, for large sparse matrices, some methods that maintain the sparse structure can be used [START_REF] Berry | Algorithm 844: Computing sparse reduced-rank approximations to sparse matrices[END_REF]. These methods are out of our scope. In our work we specically focus on compressing the dense blocks resulting from the BLR representation of sparse matrices as explained in Section 1.2.3. Additionally, we are interested in compressing in a sequential environment since these blocks are not large enough to take advantage of parallel environments. As the precision based numerical stopping criterion (presented in Section 2.1.1) in our work ensures the precision of the compression, our main aim is to choose a stable and fast kernel for our solver, which can keep the representative data with a rank as small as possible.

At a given xed numerical rank, the Singular Value Decomposition (SVD) reaches to the optimal error, through both the spectral and Frobenius norms [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF]. However, the cost of SVD is too high for large matrices, as it requires all the decomposition to complete to determine the rank. Therefore, some alternative methods are proposed to obtain a rank as close as possible to the one obtained using SVD, while being much faster at a given precision. These methods include but are not limited to the adaptive cross approximation (ACA) [START_REF] Bebendorf | Adaptive low-rank approximation of collocation matrices[END_REF], column pivoted QR (QRCP) [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF], rank revealing QR (RRQR) [START_REF] Chan | Rank revealing qr factorizations[END_REF], interpolative decomposition [LWM + 07, MT11], UTV [START_REF] Per-Gunnar Martinsson | randutv: A blocked randomized algorithm for computing a rank-revealing utv factorization[END_REF], and randomized SVD [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF].

The ACA method is proposed as a very ecient method for smooth kernel matrix approximations. However, since the correct construction of the low-rank basis is not ensured, it might fail to converge [BG05]. Therefore, despite its low cost, we prefer more stable methods. Some methods like interpolative decomposition and UTV are powerful in terms of the approximation quality. Nonetheless, the QR factorization with column pivoting (QRCP) methods can produce ranks close to the one of SVD, while being much faster. For example, in [START_REF] Per-Gunnar Martinsson | randutv: A blocked randomized algorithm for computing a rank-revealing utv factorization[END_REF], the authors show that the randomized UTV method they propose is competitive to QRCP, as they have the same asymptotical complexity. However, it can run faster than QRCP only in a parallel environment with matrices larger than the ones we are interested in. In a sequential environment with matrix sizes of our interest (dimensions smaller than 2000), QRCP can run more than 2× faster than UTV. Therefore, QRCP methods are commonly adopted in our context.

The disadvantage of the QRCP methods is that some ill-conditioned matrices, like Kahan matrix, can reach larger ranks than expected [START_REF] Golub | Matrix computations[END_REF]. However, it is worth noting that this problem occurs very rarely in real life. Nevertheless, some rank revealing QR (RRQR) methods are proposed, with an additional cost, to avoid this issue. From now on, we focus on some promising QR methods in our context.

QR based compression kernels as an alternative to SVD

In this section, we briey present the QR based compression methods before providing more specic details on the ones that we will study in this chapter. A QR decomposition has the form

A = QR,
where Q is orthonormal and R is an upper trapezoidal matrix (non-square matrix, where all elements under the diagonals are zero). Opposite to SVD, this decomposition, by its nature, does not provide a rank revealing feature. For this purpose, an orthogonal matrix

Q G (i.e. Q G Q T G = Q T G Q G = I,
where I is the identity matrix) is applied to the original matrix A. In this way, we gather the representative data on the left side to omit the remaining unnecessary ones on the right. Then, considering the block form of the QR compression method as

AQ G = Q 1 Q 2 R 11 R 12 0 R 22 ≈ Q 1 R 11 R 12 , (2.8)
where R 11 is of dimension r × r, we obtain the low-rank form (2.1) with U = Q 1 and

V T = [R 11 R 12 ]Q T G .
Remember that the SVD method is the optimal compression kernel. Therefore, as the QR based compression kernels aim to produce an approximation as close as possible to SVD, the desired approximation criterion of the QR methods can be based on the Frobenius norm of the trailing matrix as

||Q 2 R 22 || = ||R 22 || ≈ ||Σ 22 ||, (2.9)
where ||Q 2 R 22 || is the trailing matrix norm of the block QR representation in the equation (2.8). Then, similar to the criterion (2.7), the compression rank of the QR methods can be determined numerically through the criterion

||R 22 || ||A|| < ϵ.
Note that, opposite to SVD, whenever the stopping criterion is satised (if we can evaluate it), the QR process terminates and Q 2 and R 22 are not computed as they are not needed. That is why, QR compression kernels are cheap alternatives to SVD when they are used to satisfy 2.9.

There are dierent ways to compute the QR decomposition: Gram-Schmidt process, Householder transformation or Givens rotations method [START_REF] Golub | Matrix computations[END_REF][START_REF] Trefethen | Numerical linear algebra[END_REF].

As the Householder transformations are more stable than the Gram-Schmidt and cheaper than the Givens rotations, we adopt the Householder based QR approaches in our work. Note that in this chapter we focus on panel-wise QR factorizations (instead of factorizing the matrix index by index) for exploiting the level-3 BLAS [START_REF] Dongarra | A set of level 3 basic linear algebra subprograms[END_REF] operations to be more ecient on top of modern architectures. The panel-wise QR is illustrated in Figure 2.1 on the left, where b stands for block size. Here, the current panels of b columns and rows, respectively, are in yellow and the trailing matrix is in green. Then, during the factorization, the trailing matrix is updated at the end of each panel iteration through the panel contributions at once, instead of at each column contribution separately.

Although this panel-wise update operation improves the performance, the column pivoted version of the QR factorization requires an additional storage to be able to update the column norms, which is necessary to determine the next column pivot [START_REF] Quintana-Ortí | A blas-3 version of the qr factorization with column pivoting[END_REF]. Now let us focus on more specic background on the compression kernels that we will compare in this chapter.

Dierent features of QR based compression kernels

In this section, we briey discuss the dierences of the QR based compression kernels that we will study in the remainder of the chapter. The rst dierence is the way of gathering the representative data to the left part of the matrix. The second one is whether the randomized sampling is used or not, while the third one is based on when the trailing matrix updates are performed. Let us respectively explain these three dierences in three sub-sections.

Gathering the data to the left of the matrix

The rst dierence of the compression kernels in the following sections is the choice of the Q G matrix in the approximation (2.8). This choice aects how good we can approximate the SVD method. It can simply be chosen as a column permutation matrix to gather the representative data only to the left part of the matrix. For example, the column with maximum 2-norm can be chosen as the pivot at each iteration. Alternatively, the Q G matrix can be selected as a convenient rotation matrix to have more control over the quality of the approximation.

No matter whether a permutation or a general rotation matrix is chosen, the computation of the Q G matrix is inecient and it is a bottleneck for large matrices. Therefore, in the next section we explain a technique to reduce this cost.

Randomization technique

The second dierence of our kernels, randomized sampling, aims to reduce the cost of generating the Q G matrix. This technique projects the original large matrix A m×n on a lower dimensional representative matrix B d×n , where d ≪ m. Through this projection, the cost of memory-bound operations can be reduced.

As an example, we can observe an original matrix and its representative (sample) matrix in a column pivoted QR factorization in Figure 2.1. Here, j stands for an arbitrary index during the factorization and b is the block size. Green color represents the sub-matrices which will be used in the next panel iterations, while the yellow color stands for the panels which are completed in the current iteration. Observe that the row dimension of the sample matrix is chosen as b + p. Here, the oversampling size, p, stands for a small xed number, which is used to improve the quality of the projection. Then, thanks to this projection, we can cheaply compute the panel pivots through the matrix B instead of computing the column pivots through the original matrix. As a result, we can improve the performance for large matrices since column pivoting is a bottleneck in the QR method. Left-looking and right-looking approaches

The last dierence of our methods is whether the decomposition is right-looking or left-looking like the LU factorization as mentioned in Section 1.1.3. This feature is based on when the trailing matrix is updated. Let us explain it through Figure 2.2.

On the left of Figure 2.2, the left-looking approach is presented. Here, j stands for the current index. In this method, the already computed yellow panel contributions are used to update the current green panel. Here, the current panel is updated just before being factorized. In this way, the unnecessary full trailing matrix updates at each panel iteration are avoided. However, in the left-looking approach an additional storage is required compared to the right-looking version. That is, in the classical QRCP algorithm, the column norms should be downdated at each column iteration to be able to determine the next pivot column. When the trailing matrix update is postponed during the QRCP factorization, this downdate operation causes a problem. As explained in [START_REF] Quintana-Ortí | A blas-3 version of the qr factorization with column pivoting[END_REF], although the block-wise update operation requires the earlier row values of A before the block-wise column reductions, these rows should have been already updated for downdating the column norms. Therefore, this problem can be solved at the expense of an extra storage for the previous reector accumulations.

On the right of Figure 2.2, the right-looking approach is illustrated. Here, at the end of each panel iteration, all the trailing matrix (green part) is updated through the contributions from the current panel (yellow parts) in the gure. This feature is costly, especially for large matrices with small ranks, as the trailing matrix after the approximation rank is unused. However, since the panels are small and the updated parts are large for these matrices, this feature is promising for parallel environments.

j j + b j j + b j -b j j -b j W [i] T 2 Y [i] 2
Figure 2.2: Inter-panel left-looking and inter-panel right-looking approaches. The former is illustrated on the left, while the latter is on the right. The contributions of the yellow color parts are used for updating the green parts. Right-looking approach updates all the trailing matrix at the end of each panel iteration, whereas the left-looking one only updates the current panel at the beginning of each panel iteration.

Specications

In this section, we provide the necessary specications to understand the following sections better. For the sake of simplicity, we only study the matrices with real values, where the superscript T represents their transpose. However, all formulas can be easily adapted to complex matrices. In the following sections, the computational complexities of the methods stand for the square matrices of dimension n × n.

In the algorithms, we use indices starting from 0. At the end of the algorithms, R is stored in-place at the upper triangular part of A. Therefore, by applying the inverse permutation or rotation, depending on the method used, we obtain the V T matrix. The U = Q matrix can be obtained by using both the reectors, stored in-place in the lower part of the matrix, and their scalar factors vector. This can be performed easily through the xU N GQR routine of LAPACK [ABB + 07] for example.

The stopping criterion (2.6) allows to reach smaller ranks with SVD compared to the criterion (2.7). However, we use the Frobenius norm of the trailing matrix when applying our stopping criterion within the QR methods. Therefore, for being consistent with the QR methods of this chapter, we will mostly focus on the latter criterion (2.7) with Frobenius norm.

Related work

In this section, we explain the related work on the four QR methods that we study in this chapter, which dier from each other through the features discussed in Section 2.1.2.

In [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF], the authors propose a panel-wise randomized QR with column pivoting (RQRCP) method with a new sample matrix update formula, where they adopt a xed-rank based stopping criterion. The update formula in this work improves the randomization method cost by eliminating the resampling at each panel iteration. Moreover, they also oer a truncated randomized QR with the column pivoting (TQRCP) method to avoid unnecessary trailing matrix updates. Through these methods, they also come up with an ecient truncated SVD algorithm. In this chapter, we will study the RQRCP and TQRCP methods in [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF] respectively in Sections 2.4.2 and 2.4.3.

In [START_REF] Per-Gunnar Martinsson | Householder qr factorization with randomization for column pivoting (hqrrp)[END_REF], the authors independently conduct a similar work to [DG15] on the RQRCP method that is based on Householder reections, with a xed rank criterion. Here, they also contribute to the RQRCP method with a sample matrix update, instead of resampling, which is mathematically equivalent to [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF].

In [XGL17], both [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF] and [START_REF] Per-Gunnar Martinsson | Householder qr factorization with randomization for column pivoting (hqrrp)[END_REF] RQRCP methods are investigated with a xed rank based criterion. The authors analyze the eciency of the sample matrix updates from these two papers, as well as analyzing the reliability of the randomized QR with column pivoting. In addition, they propose a spectrum revealing QR method for the rank revealing problem of the ill-conditioned matrices, thanks to some extra interchanges.

In [START_REF] Martinsson | Blocked rank-revealing qr factorizations: How randomized sampling can be used to avoid single-vector pivoting[END_REF], the author proposes both column pivoted and rotational rank revealing QR algorithms, as well as a power method for a better range convergence, with a xed rank criterion. Depending on the convergence to the singular values, the randomized rotational QR (RQRRT) method can give closer results to SVD by collecting most of the data mass on the diagonals. Although the RQRRT method is more expensive than the QR with column pivoting methods, we will include this method in our study, in Section 2.4.4, to see if a good trade-o between performance and approximation quality can be reached.

In [START_REF] Martinsson | A randomized blocked algorithm for eciently computing rank-revealing factorizations of matrices[END_REF], the authors propose numerically determined rank procedures. Here, the algorithms are more general randomized methods and can be adopted as partial SVD, CUR, interpolative decomposition or RRQR. The oered rank revealing QR can be seen as a modied version of column pivoted Gram-Schmidt. This method is using the numerical stopping criterion similar to us. In this work, they also propose a power method to improve the accuracy when handling slowly decreasing singular values.

In [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF], the authors conduct a detailed study of dierent randomized methods. Here, also a randomized numerical rank decision criterion is proposed.

Note that all the methods we explain in this chapter are existing methods. However, we aim to implement them in a similar fashion for a better experimental comparison. As mentioned before, we determine the compression ranks numerically through our stopping criterion to exploit low-rank features, without any underlying problem knowledge requirement. This criterion is applied to the TQRCP and RQRRT methods for the rst time in the literature. Now let us provide some necessary notations that will be used in the following sections.

Notations

For the sake of simplicity, we use some notations to avoid detailed index specications inside the algorithms of this chapter:

-A [i] represents the sub-matrix part that we use in the i th panel iteration. For example, considering we already factorized 2b columns of A, A [2] represents the matrix A 2b:m;2b:n . The notation can be trivially converted to the vectors.

-At the i th panel iteration, A [i] is represented in the block form

A [i] = A [i] 11 A [i] 12 A [i] 21 A [i] 22
,

where

A [i]
22 stands for the trailing matrix and it is used as A [i+1] . We illustrate the matrix and vector (N ) representations at the i th panel iteration in Figure 2.3.

-A (j) represents the sub-matrix part that we use in the j th column iteration. For example, considering we already factorized 2b + 3 columns of A, A (2b+3) represents the matrix A (2b+3):m;(2b+3):n . The notation can be trivially converted to the vectors.

-At the j th column iteration, A (j) is represented in the block form

A (j) = a (j) 11 A (j) 12 A (j) 21 A (j) 22
, where A (j)

22 stands for the trailing matrix and it is used as A (j+1) . We illustrate the matrix and vector (N ) representations at the j th column iteration in Figure 2.4.

-A (l) represents the l th column of the matrix A. In case of a vector, it stands for the l th index of the vector.

Observe that both A [0] and A (0) represent the original matrix A in our notations. Now let us focus on the compression kernels that we will compare in this chapter.

A [i] 21 A [i] 11 A [i] 12 A [i] 22 = A [i+1] ib (i + 1)b ib (i + 1)b j N [i] 1 N [i] 2
Figure 2.3: Our notations and the sub-parts they stand for during the factorization. On the left we present the panel-wise matrix (A) notations, while on the right we show the vector notations on an arbitrary vector N .

a (j) 11 j + 1 j + 1 j j A (j) 22 = A (j+1) A (j) 12 A (j) 21 j n (j) 1 N (j) 2
Figure 2.4: Our notations and the sub-parts they stand for during the factorization. On the left we present the column-wise matrix (A) notations, while on the right we show the vector notations on an arbitrary vector N .

QR based compression kernels

In this section, we study four dierent QR based compression methods in detail. We explain our most basic method, QR with column pivoting (QRCP), in Section 2.4.1. In Section 2.4.2, the randomized QR with column pivoting (RQRCP) is studied. Then, the left-looking version of the RQRCP method is provided in Section 2.4.3. Finally, a rank revealing QR based compression kernel which adopts a rotational matrix is explained in Section 2.4.4.

QR with column pivoting (QRCP)

QR with column pivoting (QRCP) can be seen as the most basic method in this chapter. The QRCP we adopt is based on the LAPACK xGEQP 3/xLAQP S methods, similar to [START_REF] Drmac | On the failure of rank-revealing qr factorization software a case study[END_REF]. We start with a version of xGEQP 3 which was modied by Alfredo Buttari for the Mumps BLR solver [START_REF] Mary | Block Low-Rank multifrontal solvers: complexity, performance, and scalability[END_REF]. This method uses a panel-wise approach to take advantage of the level-3 BLAS operations and thus reduces the slow memory-cache passes of large matrices.

QRCP has an intra-panel left-looking feature. That is, each column of the panel is updated just before being factorized. On the other hand, it uses a right-looking approach for the inter-panel operations. Therefore, at the end of each panel iteration, the whole trailing matrix is updated with the current panel contributions. In Figure 2.5, we illustrate the intra-panel left-looking feature, while in Figure 2.6 we show the inter-panel right-looking feature of QRCP. In the gures, j stands for the current index and r is the last index of the previous panel. The yellow parts are used for updating the green ones similarly to Figure 2.2.

Algorithm 4 illustrates the QRCP method. Here, C stands for the column norms array (of size n), while P represents the permutation matrix. In practice, P is an array, which keeps all the corresponding pivot indices, p (j) , to swap the columns cheaply. In our code, the column norms are either updated or explicitly computed (only when required) as in [START_REF] Drmac | On the failure of rank-revealing qr factorization software a case study[END_REF]. Then, as we have the 2-norm values of the columns thanks to the column norms array C, whenever our pivot column norm is smaller than the required precision, we compute the trailing matrix norm. More precisely, we compute the Frobenius norm of the trailing matrix by computing the 2-norm of array C. Here, if the trailing matrix norm is also smaller than the required precision (our stopping criterion), we set the numerical rank to the current index value and terminate the function.

As seen in the algorithm, we start each column iteration by checking our stopping criterion. If it is not satised, the column which has the maximum norm-2 value is chosen as the pivot. As the algorithm is left-looking intra-panels, at each column iteration j, we update the current column by using the computed panel reector matrix, Y (j) , and the corresponding part of the panel accumulation matrix, i.e. W (j) T 1 . In Figure 2.5, the yellow parts illustrate these contributing parts to the current green column. After generating the elementary Householder reection for the current column, we compute the accumulations caused by the current row, w (j) T . As we mentioned before, this storage for U pdateN orms(C (j) )

Y (j) 1 Y (j) 2 W (j) T 1 W (j) T 2 j + 1 j + 1 j j r r + b r r + b
▷ Update the column norms 8:

p (j) = argmax l∈[j,...,n] C (l) 9: swap(A (j) , A (p (j) ) ), swap(W (j) T (0) , W (j) T (p (j) -j) ) ▷ Swap pivot column to current index 10:   a (j) 11 A (j) 21   =   a (j) 11 A (j) 21   -Y (j) W (j) T (0)
▷ Update the current column 11:

H (j) = I -y (j) τ (j) y (j) T
▷ Compute the Householder reection 12:

w (j)T = τ (j) (y (j) T   A (j) 12 A (j) 22   -(y (j) T Y (j) )W (j) T ) ▷ Compute the accumulation 13: R (j) 12 = A (j) 12 -Y (j) 1 W (j) T 2
▷ Update current row 14:

A [i+1] = A [i+1] -Y [i] 2 W [i] T 2 ▷ Update the trailing matrix Output : [Q, R, P , r]
the accumulations is needed in the panel-wise column pivoted QR since the trailing matrix is not updated at each iteration, while the column norms should be updated. Note that, the reector and accumulation matrices are the concatenation of each column reection vector, Y = [y (0) , y (1) , ...], and accumulation vector, W T = [w (0) , w (1) , ...] T , respectively. Finally, each column iteration terminates by updating the current row of the matrix.

Let us represent the panel reectors at the end of the i th panel iteration as

Y [i] = Y [i] 1 Y [i] 2
and the accumulation matrix as

W [i] T = W [i] T 1 W [i] T 2 . Then, Y [i]
2 stands for the

A [i] 21 part of Figure 2.3, while W [i] T 2 represents the contributions from A [i]
12 . In other words, the left-most and the upper yellow parts of Figure 2.6 illustrate

Y [i] 2 and W [i] T
2 , respectively. In Algorithm 4, every time the inner loop ends, the computation of R

[i] 11 , R [i]
12 , Y [i] and W [i] is completed. Then, at each panel iteration, the whole trailing matrix is updated by using the contributions from these yellow parts as it has an inter-panels right-looking approach.

Note that all the methods in this chapter use intra-panel left-looking approach.

Therefore, from now on, the algorithms will be called left-looking or right-looking according to their inter-panel update features.

Observe that Line 12 in the panel-wise column pivoted QR algorithm (Algorithm 4) requires large level-2 BLAS operations, i.e. y (j) T A (j)

12

A (j) 22
and y (j) T Y (j) . These level-2 BLAS operations are inecient because of the memory passes for large matrices at each column iteration. Therefore, we will explain a technique to avoid these inecient memory passes in the next section.

The QRCP method reduces the complexity compared to SVD (with the complexity O(n 3 )) to O(n 2 r). Here, the pivot selection operation costs O(n 2 r), while applying the pivot is as cheap as simple column swaps. As QRCP is right-looking inter-panels, the whole trailing matrix updates are performed at a cost of O(n 2 r).

Randomized QR with column pivoting (RQRCP)

As mentioned in the previous section, the need of pivoting leads to many inecient level-2 BLAS operations in the panel-wise QRCP factorization, which reduces the performance for large matrices. In this respect, some randomization methods are proposed to use a lower dimensional representative matrix, instead of the original one, to determine the pivot columns of the panels. For this purpose, a Gaussian matrix with independent and identically distributed random variables (Gaussian i.i.d.), Ω, can be used to project the original matrix onto the sample matrix, B, as

B d×n = Ω d×m A m×n ,
where d ≪ m. Here, we need to get the representative data with small uncertainity to be able to determine the correct pivoting for the current panel. That is, let us consider the j th column of A, a (j) , and the corresponding sample matrix column b (j) .

Then,

||b (j) -b (i) || 2 2 d||a (j) -a (i) || 2 2 -1 ≤ τ ,
with 0 < τ < 0.5 being the threshold and 0 ≤ i, j < n, should be satised with high probability [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF]. This is why we do not choose the row dimension of B as block size b. Instead, we use a pre-dened oversampling size, p. Then, the row dimension is d = b + p. In [XGL17], it is shown that the failure probability of RQRCP is exponentially decreasing for larger oversampling sizes. Although we need this parameter for the quality of the sampling, we should be careful to choose it as small as possible, while respecting the error bounds, to eciently improve the performance of the randomized kernel.

Note that we have to resample the matrix at each panel iteration in the classical approach to be able to determine the panel pivots. However, it is too expensive in practice. Instead, we can update the sample matrix in a cheaper way and use it at each iteration as proposed in [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF][START_REF] Per-Gunnar Martinsson | Householder qr factorization with randomization for column pivoting (hqrrp)[END_REF]. In these works, the authors prove that the performance results of these randomized techniques become close to the unpivoted QR methods, while the accuracy is very close to the one of unrandomized methods. Therefore, in this section and in Section 2.4.3, we adopt the randomized methods through the updated sample versions as in these works.

As explained in [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF], considering the block factorization of the matrix B

BP [i] = ΩAP [i] = Q [i] B S [i] 11 S [i] 12 0 S [i] 22
, the update of the sample matrix can be computed as

B [i] 1 B [i] 2 = Ω [i] A [i] = S [i] 12 -W [i] 11 R [i] 12 S [i] 22 = S [i] 12 -S [i] 11 R [i] -1 11 R [i] 12 S [i] 22
.

Algorithm 5 Randomized QR with Column Pivoting

1: function RQRCP(A, ϵ, b) 2: Ω [0]
▷ Generate the Gaussian i.i.d. matrix 3:

B [0] = Ω [0] A [0]
▷ Compute the sample matrix 4:

for i = 0, 1, ..., ⌜ min(m, n) b ⌝ -1 do 5: [Q [i] B , S [i] , P [i] B , r] = QRCP (B [i] , (b + p)ϵ, b) ▷ Compute permutation 6: A = AP [i] B
▷ Apply the permutation 7:

[

Q [i] , R [i]
11 ] = QR(A [i] ) ▷ Apply QR on the permuted matrix A 8:

  R [i] 12 A [i+1]   = Q [i] T   A [i] 12 A [i] 22   ▷ Update rightmost matrix 9: B [i+1] 1 = S [i] 12 -S [i] 11 R [i] -1 11 R [i]
12 and B

[i+1] 2

= S

[i] 22

▷ Update the sample matrix

Output : [Q, R, P , r]
Now, let us observe the randomized QRCP method, with the sample matrix update formula, in Algorithm 5. Here, the stopping criterion is applied implicitly through the QRCP function, which is used for computing the permutation matrix P B through the sample matrix. As the permutation is computed through the sample matrix, the threshold should also be arranged according to B. As explained in [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF], the column norms of B and A have a constant relative ratio through the Chi-squared distribution [START_REF] Jay | Probability and Statistics for Engineering and Science -8[END_REF]. As a consequence, their residual error ratio is approximately 1 b+p , with b + p being the row dimension of the sample matrix. Therefore, assuming ϵ is the threshold for the matrix A, the threshold for the sample matrix B is computed as

ϵ B = (b + p)ϵ.
As seen in the algorithm, RQRCP is similar to QRCP (Algorithm 4), except for the computation of the permutation matrix and the sample matrix generation/update operations. Then, in this algorithm, the additional cost to QRCP is the cost of generating the sample matrix B (b+p)×n in the beginning, O(n 2 b), as well as its update at all iterations, O(nrb). However, the pivot nding operation using the input matrix A, O(n 2 r), is avoided. The pivoting operation through the sample matrix reduces this cost to O(nrb). Through these computational complexities, we can observe that the randomization method does not aim to reduce the theoretical complexity compared to QRCP. Instead, it targets to reduce the time consuming data movements that occur when computing the column norms through large matrices.

Truncated randomized QR with column pivoting (TQRCP)

As explained in Section 2.1.2, the left-looking feature aims to reduce the complexity by avoiding the whole trailing matrix updates. In [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF], the authors propose the left-looking version of the RQRCP method for this purpose. In this section, we study this left-looking randomized method, which is called truncated randomized QRCP (TQRCP).

Ω [0]
▷ Generate the Gaussian i.i.d. matrix 3:

B [0] = Ω [0] A [0]
▷ Compute the sample matrix 4:

for i = 0, 1, ..., ⌜ min(m, n) b ⌝ -1 do 5: [Q [i]
B , S [i] , P

A = AP [i]
B and W

[i] T G = W [i] T G P [i] B
▷ Apply the permutation 7:

  A [i] 11 A [i] 21   =   A [i] 11 A [i] 21   -   Y [i] G1 Y [i] G2   W [i] T G1 ▷ Update current panel 8: [Q [i] , R [i] 11 , Y [i] ] = QR(   A [i] 11 A [i] 21   )
▷ Apply QR on the current panel 9:

W [i] T = T [i] T (Y [i] T A -(Y [i] T Y [i] G )W [i] T G ) ▷ Compute the accumulations 10: R [i] 12 = A [i] 12 -Y [i] G1 W [i] T G2
▷ Update the row panel 11:

B [i+1] 1 = S [i] 12 -S [i] 11 R [i] -1 11 R [i] 12 and B [i+1] 2 = S [i] 22
▷ Update the sample matrix

Output : [Q, R, P , r]
As explained before, in the right-looking approach, we use panel-wise reectors and panel-wise accumulations to update the whole trailing matrix. In the left-looking approach, we use global reectors, Y G , and global accumulations, W T G , to update only the current panel. In Figure 2.7, we illustrate the Y G and W T G matrices of the left-looking

Y [i] G1 Y [i] G2 W [i] T G1 W [i] T G2 ib (i + 1)b ib (i + 1)b
Figure 2.7: Global reectors, Y G , and global accumulations, W T G , at the beginning of the i th panel iteration. Here, the yellow parts are used to update the current green panel. approach at the beginning of the i th panel iteration. Here, the yellow global reectors and yellow global accumulations are used to update the current green panel.

In Algorithm 6, we provide the code for the left-looking TQRCP method. As explained in [START_REF] Duersch | True blas-3 performance qrcp using random sampling[END_REF], TQRCP follows the same logic as converting the column-wise QRCP to the panel-wise QRCP, but at a panel-wise level. That is, let us consider the composition of all the global reectors, Y

[i+1] G , as Y [i+1] G = Y [i] G Y [i] ,
the block reectors are updated as

(I -Y [i] G T [i] G Y [i] T G )(I -Y [i] T [i] Y [i] T ) = I -Y [i+1] G T [i+1] G Y [i+1] T G , where T [i+1] G is T [i+1] G = T [i] G -T [i] G Y [i] T G Y [i] T [i] 0 T [i] .
Then, the global accumulation matrix, which results in extra storage, for this panel-wise left-looking code is

W [i+1] T G = T [i+1] T G Y [i+1] T G A = W [i] T G W [i] T = T [i] T G Y [i] T G A T [i] T (Y [i] T A -(Y [i] T Y [i] G )W [i] T G ) .
As seen in the algorithm, TQRCP is indeed the panel-wise left-looking version of RQRCP. That is, after each panel iteration, RQRCP updates

A [i] 12 A [i] 22
, while TQRCP stores the contribution accumulations and only updates

A [i]
12 . Observe that TQRCP also applies the stopping criterion implicitly, similarly to RQRCP.

Compared to RQRCP, there is an extra storage cost for the matrix W T G , of size r × n, in the TQRCP method. However, at each panel iteration of the TQRCP method, the update of A As the sample matrix generation, of complexity O(n 2 b), is performed only once at the beginning of the algorithm, it is not necessarily dominant since actually it has a constant 1. Then, thanks to the low cost trailing matrix updates, TQRCP has the lowest computational complexity among the four methods that we study. For large matrices with low ranks, this method highly reduces the update cost of the trailing matrix.

Randomized QR with rotation (RQRRT)

In this section, the rotational QR method in [START_REF] Martinsson | Blocked rank-revealing qr factorizations: How randomized sampling can be used to avoid single-vector pivoting[END_REF] is studied. This method is proposed to x the rank revealing issue of the column pivoted QR method on some matrices.

Remember that the rank revealing QR factorization aims to produce an approximation as close as possible to the optimal approximation (SVD) (see the approximation criterion 2.9). Then as explained in [START_REF] Martinsson | Blocked rank-revealing qr factorizations: How randomized sampling can be used to avoid single-vector pivoting[END_REF], we can apply a convenient orthonormal matrix to A to produce a rank closer to SVD to satisfy this criterion, instead of only permuting the columns.

In the randomized QR factorization with rotation (RQRRT) method, the sample matrix is used similarly to the randomized QRCP. However, this time the original matrix, A, is not permuted through the permutation that is obtained by using QRCP of the sample matrix, B. Instead, an orthonormal matrix is generated through the unpivoted QR factorization of B T . This is indeed an eective method to get the mass to the left of the matrix, since the sample matrix is generated in a way that it really represents the linear dependencies of the original matrix. Therefore, this orthonormal matrix represents the rotation of the input matrix when the QR factorization is applied. When this rotation matrix is applied to matrix A, its representative data moves to the diagonals to have closer results to SVD.

Algorithm 7 stands for the RQRRT compression kernel. Here, the output rotational matrix is equivalent to

Q rot = Q [0] B Q [1] B Q [2] B ...Q [⌜ r b ⌝-1] B
. Let us compare this RQRRT algorithm to the randomized QRCP (Algorithm 5). Their rst dierence is the way of gathering the representative data. In the randomized QRCP algorithm, we nd

AP B = Q A R A ,
where P B is obtained through the QRCP algorithm as

BP B = Q B R B .
Ω [i]
▷ Generate the Gaussian i.i.d. matrix 3:

for i = 0, 1, 2, ..., ⌜ min(m, n) b ⌝ -1 do 4: B [i] = Ω [i] A [i] ▷ Resample B 5: [Q [i] B , S [i] ] = QR(B [i] T , b)
▷ Compute the rotation matrix 6:

A [i] = A [i] Q [i] B
▷ Apply the rotation 7:

[

Q [i] , R [i] 11 ] = QR(A [i] ) ▷ Apply QR on A 8:   R [i] 12 A [i+1]   = Q [i] T   A [i] 12 A [i] 22   ▷ Update the trailing matrix 9: r = StoppingCriterion(R [i] 11 , R [i] 22 , ||A [0] || F , ϵ, b) ▷ Check stopping criterion Output : [Q, R, Q rot , r]
However, in the RQRRT method, the AQ ′ B matrix is used as

AQ ′ B = Q ′ A R ′ A
, where the orthogonal matrix is obtained from the partial QR of the sample matrix as

B T = Q ′ B R ′ B .
This shows the column orthonormalization of B T . Therefore, the leftmost b column spanning the matrix A will have almost the same spanning as the leading b singular vectors of A. Therefore, it gets a closer result to SVD in terms of rank revealing quality [START_REF] Martinsson | Blocked rank-revealing qr factorizations: How randomized sampling can be used to avoid single-vector pivoting[END_REF].

The second dierence of the rotational method is the resampling operation. In the column pivoted method, the sample matrix is cheaply updated at each panel iteration. However, in the rotational method, as its transpose is used, there is no update formula (up to our knowledge).

The third dierence of the rotational method is the explicitly implemented stopping criterion (at Line 9), since we do not call QRCP function anywhere in this code. Here, the input matrix R

[i]
22 stands for A [i+1] . The detailed stopping criterion in RQRRT code is given in Algorithm 8. Everytime this function is called, the trailing matrix norm on Line 2 is explicitly computed. Here, it is important to observe that at the end of the panel iteration, if the criterion (on Line 2) is satised, we check the previous trailing matrix norms by adding each column contribution starting from the last column of the panel (on Line 4). Bear in mind that the value ||(R

[i]
11 ) (l) || 2 does not have to be computed explicitly since it is equal to the magnitude of the rst element of this column. If we check the trailing matrix norms with the accumulations of the left-most columns rst, we can lose stability because of the round-o errors. Note that in the QRCP method, the stability issue is already taken care of and it is out of scope in this thesis.

In Figure 2.8, the stability issue is observed. Here, the result is obtained by using the S-shape matrix (see Section 2.6.1) of size 500 × 500, where the generation rank is 25.

In the gure, RQRRT (right) represents our algorithm and RQRRT (left) stands for the unstable accumulation version. As a reference, the trailing matrix norms at each panel index at generation, n i=K σ 2 i , is also shown (called Reference). The y-axis stands for ||A -U K V T K || F , while x-axis shows the indices, K, of the rst panel iteration. Here, U K is of dimension m × K and represents the K columns of the matrix U . Similarly, V T K stands for the rst K rows of V T .

Figure 2.8 shows that indeed the right-most accumulation gets closer results to Reference. As opposed to our approach, the left-most accumulation diverges after the index 23 supporting our statement.

Algorithm 8 Stopping criterion in RQRRT

1: function StoppingCriterion(R [i] 11 , R [i] 22 , ||A [0] || F , ϵ, b) 2: if ||R [i] 22 || F ≤ ϵ||A [0] || F then 3: for l = b -1, b -2, ..., 1 do ▷ Accumulations starting from right-most column 4: ||R (ib+l) 22 || F = ||R (ib+l+1) 22 || 2 F + ||(R [i] 11 ) (l) || 2 2 5: if ||R (ib+l) 22 || F > ϵ||A [0] || F then return r = ib + l + 1 Output : [r]
For RQRRT, nding the rotation matrix has a complexity of O(nrb), while applying it to the matrix A is of cost O(n 2 r). However, for RQRCP, the pivot nding cost is O(nrb) and the AP B operation is only as cheap as swapping b columns of A. Another cost dierence of RQRRT is the resampling at all iterations, O(n 2 r), compared to the sample matrix update cost in RQRCP, O(nrb). As seen here, the complexity dierences between these two methods can be signicant for large matrix sizes, which makes RQRRT the most costly QR method that we study.

Summary of the complexities

In this section, we aim to summary and comment on the complexities of the compression kernels studied in the previous sections. Table 2.1 illustrates the computation complexities of the compression methods for a generic matrix of size n × n, with a block size b and a compression rank r. Here, only the operations leading to a dierence between the methods are shown for the sake of simplicity. The third column of the table stands for the total cost of each operation, while the last column shows the overall cost of each method. When determining the values in the last column, we assume r > b. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1.00e-14

O(n 2 r) Rotation computing O(nrb) Rotation of A O(n 2 r) Trailing matrix update O(n 2 r)
Table 2.1: Important complexity changes between the compression methods. Here, n, r and b stand for the matrix size, the compression rank and the block size, respectively.

As seen in the table, the truncated method has a total of O(n 2 b) complexity, while all other QR variants have a O(n 2 r) complexity. As in general b < r and TQRCP reduces the trailing matrix update cost, this truncated method is the cheapest procedure for large matrices with low ranks, in a sequential environment. Since SVD has a O(n 3 ) complexity, it is the most costly method, as expected.

For the QRCP, RQRCP and RQRRT methods, the overall complexity is asymptotically the same. However, when each operation complexity is observed, one can see that the usage of the rotational feature in RQRRT introduces large additional costs compared to the other two methods. That is why, it is the most costly QR variant. In RQRCP, the pivoting cost is reduced compared to QRCP. However, there are additional costs for sample matrix generation (in the beginning) and for updating it at each iteration. Thus, this method is advantageous for large matrices, where mostly the communication of the column norms computations are costly. To conclude this table, we expect in our experiments to observe the performance order as (from worst to best, in a sequential environment): SV D ≪ RQRRT < RQRCP < QRCP ≪ T QRCP .

Experiments on generated matrices

In this section, we compare the compression methods on dierent generated matrices, including challenging cases. For the sake of simplicity, the row and column dimensions are chosen to be equal. From now on, the letter r G will represent the generation rank of the matrices, whereas r represents the approximation rank of them.

As a matter of interest, the compression kernels are in sequential within the PaStiX framework. Therefore, the results in this section demonstrates the sequential results, since parallelism is not worth for the matrix sizes we are interested in. For the RQRCP and TQRCP methods, the block size is set as d = b + p, where b = 27 and the oversampling parameter is set to p = 5. For the other methods, the block size is d = b = 32. Let us emphasize that we implement all the compression methods in a similar fashion and we determine the compression ranks numerically. Here, let us also note that setting b = 32 and p = 5 would be a better choice when comparing RQRCP and TQRCP to the other methods, where d = b = 32. However, it will not aect our results observably.

In Section 2.6.1, we rst provide the information about how we generate the experimental matrices. Then, in Sections 2.6.2, 2.6.3, 2.6.4 and 2.6.5, we discuss the stability, performance, numerical ranks and accuracy results on these matrices, respectively.

Generation of the experimental matrices

In this section, ve dierent type of generated matrices are used. These matrices are similar to the ones in [START_REF] Martinsson | Blocked rank-revealing qr factorizations: How randomized sampling can be used to avoid single-vector pivoting[END_REF], and are shown in Figure 2.9, Figure 2.10 and Figure 2.11. In these gures, the x-axis stands for the matrix index, K. The gures on the left illustrate the singular value distribution of a matrix of size 500 × 500, where the generation rank is 100 at the generation precision ϵ = 10 -8 . The gures on the right stand for the relative Frobenius norm of the trailing matrix,

√ n i=K σ 2 i ||A|| F
, of the same matrices. The solid black lines in the gures represent the generation precision.

The rst matrix type is called k-rank and it is observed in Figure 2.9. In order to generate this matrix, a normalized Gaussian i.i.d. matrix is used, where the singular values decrease slowly and suddenly drop to ignorable values at the generation rank. Because of the sudden drop to ignorable values at r G , this matrix type is trivial for the compression methods.

The second and third matrix types are called Z-shape and Z-shape-short, which are shown in Figure 2.10. Both of these matrices are computed as A = U DV T , where U and V are random orthonormal matrices. The singular values (diagonal values of the matrix D) are the same for both cases and decrease rapidly up to the generation rank. However, the Z-shape drops faster to ignorable values at the index 3 2 r G = 150, while Z-shape-short stagnates at 10 -9 . Here, we expect the methods to nd close approximation ranks to r G = 100 for the Z-shape because of the fast drop after the generation rank. Opposite to Z-shape, as seen in the right-most gure, the trailing matrix of the Z-shape-short matrix contains important data at r G with respect to the compression precision 10 -8 . Thus, this case is challenging for our stopping criterion to determine the compression rank close to the generation one.

The last matrix types are called S-shape and S-shape-short, which are shown in Figure 2.11. These matrices are computed as A = U DV T , with random orthonormal matrices U and V . For these matrices, the singular values stagnate at 1 up to the index q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e-21 1e-16

1e-09 1e-08 1e+00 1e+01 0 100 200 300 400 500 K Singular value q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e-21 ) are on the right, while the singular values at each index (σ K ) are on the left. The x-axis stands for the matrix indices, K. The singular values rstly decay slowly, then suddenly drop to very low values. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e-21 1e-16 1e-09 1e-08 1e+00 1e+01 0 100 200 300 400 500 K Singular value q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e-21 ) are on the right, while the singular values at each index (σ K ) are on the left. The x-axis stands for the matrix indices. Both matrices have the same singular values up to the index r G = 100. However, after this point, Z-shape case singular values decrease faster and reach to almost machine precision at the index 3 2 r G , whereas the Z-shape-short case singular values reach to ϵ 10 = 10 -9 at this index. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e-21 1e-16 1e-09 1e-08 1e+00 1e+01 0 100 200 300 400 500 K Singular value q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e-21 ) are on the right, while the singular values at each index (σ K ) are on the left. The x-axis stands for the matrix indices. The singular values of the S-shape and S-shape-short matrices rstly stagnate at the value 1, then they decrease fast and in the end they stagnate again after a heavy tail. Both of the singular values reach the value ϵ = 10 -8 at the index r G = 100. However, the last stagnation point is ϵ 2 = 10 -16 for the S-shape matrix, whereas it is ϵ 10 = 10 -9 for the S-shape-short matrix. r G 2 , and then they drop rapidly to reach ϵ = 10 -8 at index r G . After that point, the S-shape case stagnates at the value ϵ 2 = 10 -16 , while S-shape-short stagnates at ϵ 10 = 10 -9 , which is close to the generation precision. These matrices are challenging experiments for the randomization methods as they have a heavy tail.

In the following sections, we use the compression precision ϵ = 10 -8 , which is equal to the generation precision. We observe the SVD method with two dierent stopping criteria. One of them uses the Frobenius norm of the trailing matrix to determine the compression rank as in 2.7, while the other uses the spectral norm for this purpose as in 2.6. The former, SVDF, provides better comparison to the QR methods, whereas the latter, SVD2, leads to closer ranks to the generation ones.

Stability

In this section, we observe the stability of the compression kernels. In Figure 2.12, the relative residual with respect to the initial matrix,

||A-U K V T K || F ||A|| F
, is observed at each index K. As these values show the approximated values of the relative Frobenius norms that were presented on the right side of Figures 2.9, 2.10 and 2.11, we want to obtain close results to them.

As seen in Figure 2.12, all the methods have good stability. As RQRCP and TQRCP use a relative stopping criterion, which uses an approximation ratio between the original and sample matrix columns, there are more magnitude changes between the panel iterations. However, it will be shown in Section 2.6.5 that these methods still ensure a good accuracy, so this does not indicate a problem. As we can see in the gures, the RQRRT method indeed obtains very close results to SVD.

Another interesting point in these gures is the black horizontal line, which stands for the compression precision ϵ = 10 -8 . That is, observe that our stopping criterion (see Section 2.1) compares the y-axis values in Figure 2.12 to the ϵ value to determine the approximation rank. Therefore, whenever a residual value reaches this black line in the gures, we know that the compression rank is the corresponding index. Therefore, through these stability gures, we observe that for the Z-shape-short case, all the methods get much larger compression ranks than r G = 100. Thus, this case is challenging for our numerical stopping criterion. However, for all the other cases, all the methods reach close approximation ranks to the generation one.

In Figure 2.13, the same information as in Figure 2.12 is shown. However, here the relative residual with respect to the SVDF is presented to better observe the stability compared to the SVD method.

The main observation in Figure 2.13 is that QRCP, RQRCP and TQRCP can have a relative error around 2.5× larger than the one of the SVD around the index r G . However, this relative error is very close to SVDF for the rotational QR variant, which conrms that indeed RQRRT provides closer quality approximations to SVD.

Note that for the k-rank, Z-shape and S-shape test cases, the last values are around the machine precision. Therefore, the extreme relative ratios at large indices are not representative, so we dropped them in the gure. .12: In the stability gures, the relative residual value with respect to the initial matrix,

||A-U K V T K || F ||A|| F , at each index K is illustrated.
||A-U K V T K || F ||A-U K V T K || SV DF F , at each index K is illustrated.
2.6.3 Performance q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Rank = MatSize * is observed for dierent matrix sizes and r G values. In this gure, each row illustrates a dierent test case, whereas each column represents a dierent generation rank. Here, we applied the stopping criterion within the compression methods.

In this section, we observe the performance results of the generated matrices that are explained in Section 2.6.1. In Figure 2.14, each test case is observed with dierent matrix sizes and dierent generation ranks.

In the gure, every grid row shows a dierent test case, while the columns stand for dierent generation ranks. The y-axis of the gures stands for GF lops = M inGF lop t . Here, t (in seconds) represents the time to compress the matrix with the compression rank r for each method, whereas M inGF lop stands for the minimum number of GFlop to compress the matrix of size n × n and rank r. We compute it as

M inGF lop = GF lop GEQRF (n, r) + GF lop U N M QR (n, n -r, r) + GF lop U N GQR (n, r, r).
Here, the GEQRF , U N M QR and U N GQR names are adopted from LAPACK. The GF lop GEQRF (n, r) represents the number of ops to perform a QR factorization on a n × r matrix. GF lop U N M QR (n, n -r, r) stands for the ops to apply the Q to the V part of size r × n of the approximation. Finally, GF lop U N GQR (n, r, r) represents the number of ops required to generate the nal n × r matrix U from the Q matrix in the Householder form.

As seen in Figure 2.14, TQRCP has closer ops to the minimal ops and has always the best performance (up to 20 GFlops) for large matrix sizes, which proves the eectiveness of the left-looking approach in sequential environments.

On average, QRCP reaches the second best performance in the gures, as expected, while it has even the best performance for small matrix sizes. This is an interesting result since the randomization technique in RQRCP was introduced to increase the performance of QRCP by reducing the slow memory passes. However, since the tests run in sequential with matrix dimensions lower than 3500 and RQRCP introduces more computational complexity compared to QRCP, it is not an unexpected result. If the experiments were performed in a parallel environment for larger matrices, RQRCP would have better performance than QRCP.

As expected, RQRRT has the worst performance compared to other QR variants because of its high complexity. However, SVD has the worst performance overall since no matter which stopping criterion is used, this method should apply the full decomposition to be able to determine the compression rank.

The important point in the gures is that the QRCP method has a better performance than TQRCP for small matrix sizes, as mentioned before. However, for the larger sizes TQRCP is the best. Therefore, tuning the methods according to the matrix size is also important since there is not an absolute best method.

Another important point is that in the stability section, we have seen that for the test cases, all the compression ranks are close to the generation rank, except for the Z-shape-short case. This causes a signicant performance reduction for all the methods through this test case. Especially, for the small ranks, the performance loss is very critical for the Z-shape-short case. However, we can also observe that even if the performance is reduced signicantly, the performance order of the methods is still the same as the other test cases, which shows a consistency among all the gures.

Compression ranks

In this section, we observe the compression rank comparison of the methods for the ve generated matrices we introduced. In Figure 2.15 and Figure 2.16, the relative rank q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Rank = MatSize * In the relative rank graphs, the relative rank with respect to SVDF, r r SV DF , for dierent matrix sizes and r G values is observed. In this gure, each row illustrates a dierent test case, whereas each column represents a dierent generation rank. Here, the stopping criterion is applied in the compression methods. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Rank = MatSize * In the relative rank graphs, the relative rank with respect to SVDF, r r SV DF , for dierent matrix sizes and r G values is observed. In this gure, each row illustrates a dierent test case, whereas each column represents a dierent generation rank. Here, we applied the stopping criterion within the compression methods.

with respect to SVDF, r r SV DF , is observed for dierent matrix sizes.

As expected, all the methods reach to the generation rank for the k-rank case. On the other hand, for the Z-shape-short and S-shape-short cases, the rank dierences are higher since these cases have large trailing matrix norms after the generation rank as observed in Figures 2.10 and 2.11, respectively. That is, they are more challenging for our numerical stopping criterion, which uses the Frobenius norm of the trailing matrix. As these test cases result in more rank dierences than other cases, we show them separately in Figure 2. 16. In the gures, the lowest possible ranks are reached by SVD2 since it uses a more strict stopping criterion than the ones that adopt the Frobenius norm of the trailing matrix.

As the rank dierences are higher for the tricky case Z-shape-short, we can easily observe that RQRCP and TQRCP reach closer results to SVD2 compared to all other methods (including SVDF). This is caused by the higher residual norm dierences within the panels for RQRCP and TQRCP as seen in Figure 2.12. The high residual norm dierences can be related with the stopping criterion used in RQRCP and TQRCP, which is based on a sample matrix, whereas other methods use the original matrix norms to terminate. In addition, this dierence can be related to the accumulated round-o eects of small numbers. It is also interesting to observe that the round-o eects of small numbers lead RQRRT to not always reach to the closest results to SVDF. However, as we will observe in Figure 2.17, all the methods provide a good accuracy since we work with xed precision.

Accuracy

In this section, we observe the accuracy of the compression methods through our numerical stopping criterion. In Figure 2.17, the error, ||A-UrV T r || F ϵ||A|| F , for dierent matrix sizes and r G values is observed.

As seen in the gure, since for the k-rank case it is very trivial to reach to the generation rank, there is no observable error dierence between the compression methods.

In the previous section, we observed that RQRCP and TQRCP tend to have smaller ranks than QRCP. Therefore, in Figure 2.17, they tend to have more data in the trailing matrix after the compression. Similarly, since SVD2 nds the smallest ranks in general, the data in the trailing matrix (error) tends to be more than the other methods.

When tackling very small numbers through the stopping criterion within the RQRRT method, the resulting round-o errors can slightly aect the approximation rank and error. As observed in the gure, this does not even cause one digit accuracy changes for this method, since our numerical stopping criterion guarantees a good accuracy. Therefore, we do not observe a noteworthy accuracy dierence between the compression methods. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Rank = MatSize * , for dierent matrix sizes and r G values is observed. Here, we applied the stopping criterion within the compression methods.

Experiments on real-life case matrices

In this section, we compare the compression methods through 810 real-life case matrices. These matrices are downloaded from https://sparse.tamu.edu/ with the criteria that their row and column dimensions are less than 2000. Among them, there are 802 matrices with real data and 8 matrices with complex data. There are both square and rectangular matrices. That is why, in this section, the matrix dimensions are noted as m × n instead of n × n.

Since our stopping criterion ensures a good accuracy, as we observed on the generated matrices, we only present the performance and obtained rank results in this section. In addition, we only provide the SVDF method results (not SVD2) as a reference, since our main point is to compare the QR variants, which use the Frobenius norm of the residuals in the stopping criterion. Note that, in this section, we use the same experimental settings as in Section 2.6. .18: Time proles of the real-life case test matrices. The x-axis ratios are computed as tcurrent t min for each method, while the y-axis stands for each matrix as a number from 1 to 810. Here, the x-axis is in logarithmic scale for the sake of clarity of the gure.

In Figure 2.18, we observe the time ratio, tcurrent t min , of each method for all matrices. In this time prole, a method is better on average, if its curve is closer to the line y = 1.

As expected, SVDF is the slowest method for all the cases, where RQRRT is the slowest QR variant for most of the matrices. For some very small cases, we cannot observe the expected method performances, since advantages of the methods cannot be exploited at these levels. It is interesting to observe that QRCP is the fastest method for the majority of matrices, while TQRCP is the second fastest method. The reason for this is that the real-life case matrices we used are mostly small matrices. In Section 2.6.3, we observed that TQRCP is advantageous for large matrices, so this is an expected result. RQRCP never runs as the fastest method in these experiments. As we do not use large enough matrices to benet from this method, as well as we work in a sequential environment, this method is not promising in our context. Zoomed obtained rank proles of the real-life case test matrices. The x-axis ratios are computed as rcurrent r min for each method, while y-axis represents each matrix as a number from 650 to 810. The y-axis values (from 0 to 649) are omitted to zoom the gure since the methods get very close results at this interval.

In Figure 2.19, we observe the compression rank ratio, rcurrent r min , of each method for some of the 810 matrices. This gure is a zoomed version of the original one in a way that it only stands for 160 matrices, which have observable dierence in the gure. That is, it shows the top-most part (y-axis from 650 to 810) of the original gure. The previous y-axis values (from 0 to 649) are omitted in the gure as the methods obtain similar results at this interval.

In Figure 2.19, the values of RQRCP and TQRCP intersect, so RQRCP does not show up. Both of these methods nd the lowest compression ranks for 734 cases. Supporting the results in Section 2.6.4, even SVD nds larger ranks compared to RQRCP and TQRCP for 140 matrices because of the sample matrix based stopping criterion. Thus, RQRCP and TQRCP seem to be the most promising compression methods to nd the lowest ranks possible. Considering also the performance proles results, TQRCP seems to be better than RQRCP in sequential environments (for the matrices with dimensions less than 2000). That is, TQRCP does not only obtain the lowest ranks possible, but this left-looking method also compresses all the matrices faster than RQRCP (see Figure 2.18). We can compare the quality of using a rotation matrix or a permutation matrix through QRCP and RQRRT. As both use a stopping criterion through the original matrix, the gure proves that RQRRT slightly improves the obtained ranks by getting closer to SVD. However, taking the performance results into account, we do not get good enough trade-o between rank and performance, in our context, through RQRRT. QRCP seems to be the least promising method to nd the minimal rank. However, it is worth noting that the rank dierences between the methods are not important for most of the matrices. Thus, QRCP is still promising for the small matrices, thanks to its performance.

To conclude, as good accuracy is ensured by our stopping criterion, we need to mostly focus on the obtained rank and performance results. Then, QRCP and TQRCP seem to be the most promising compression methods for us. Tuning between these methods according to the matrix size and the compressibility of the matrix (if known) can contribute to minimize the total time of the compression process, as well as the compression rank.

Discussions

In this chapter, we study and compare four dierent stable QR decomposition techniques, which aim to be fast alternatives to SVD when compressing dense matrices.

The most basic QR version is a panel-wise right-looking QRCP method. The second one is the randomized version of the QRCP method, so it is called randomized QRCP (RQRCP). It aims to avoid the costly data movements for large matrices when computing the column norms during the pivoting. Here, we adopt a version, where the resampling is avoided thanks to a cheaper sample matrix update formula. The third method is the left-looking version of the RQRCP method, which is called truncated randomized QRCP (TQRCP). This method reduce the computational complexity by eliminating expensive full trailing matrix updates at each iteration. Therefore, in a sequential environment, this method is the cheapest one especially for large matrices with low ranks. The last method we study is a rotational version of RQRCP. Here, we use a rotation matrix instead of column pivoting, so that we can better approximate SVD.

All the studied methods in this chapter are existing methods. However, we implement them in a similar fashion for a better comparison. In addition, we determine the compression ranks through our numerical stopping criterion at a given precision. This stopping criterion was never used for two of the studied methods, namely TQRCP and RQRRT. We conduct all the experiments in the framework of the PaStiX sparse direct solver, and as a matter of interest, we observe the methods in a sequential environment.

In our experiments, we show that all the methods are accurate thanks to the nature of our numerical stopping criterion. Thus, we mostly focus on the performance and obtained rank results.

We observe that the obtained ranks for the RQRCP and TQRCP methods might be even better than SVDF thanks to the sample matrix based stopping criterion, while still satisfying the accuracy requirement. As TQRCP method also outperforms other methods for larger matrices, it is the best method for these matrix sizes in a sequential environment. We should keep in mind that all the compression methods obtain close ranks for the majority of the matrices. Therefore, performance results determine the best method for small matrices, which is QRCP. Then, tuning between QRCP and TQRCP is necessary according to the block sizes we use in our solver.

The QRCP method has more performance than RQRCP even for the largest matrices we used. Therefore, RQRCP is not an advantageous method for us. Nonetheless, it is still a promising method in other contexts. Here, we could not benet from this compression method as all the experiments are in sequential and our matrix sizes might not be large enough to exploit it.

Chapter 3 Block Low-Rank Clustering

As explained in Section 1.2, direct solvers became more promising for large systems after the introduction of the low-rank representations. Thanks to these representations, the time-to-solution and memory footprint of these robust solvers can be highly reduced. Among the low-rank representations, some consider the full problem and extract the sparsity of the whole matrix from the low-rank representation, while others, like the block low-rank (BLR) representation, compress the dense blocks of the solver independently. In this work, we focus on improving the latter. More specically, this chapter focuses on improving low-rank supernodal direct solvers for sparse systems in two dierent aspects.

Firstly, determining which blocks to compress and when to do it is an important problem for these solvers to be time ecient while beneting from the lower memory consumption oered by the low-rank techniques. Therefore, we want to nd a suitable admissibility criterion to oer an intermediate solution between the existing strategies that we explained in Section 1.4. For this purpose, we gather the poorly compressible blocks together by using the graph information. That is, we do not identify the existing uncompressible blocks. Instead, the identication is done through the reordering of the separator unknowns in a specic way to collect the uncompressible data together. Then, our target is to compress the highly compressible blocks at initialization (like Minimal Memory scenario) to improve the memory footprint, while compressing the poorly compressible ones on the y (like Just-In-Time scenario) to reduce the ops overhead of the solver with a controlled memory overhead compared to Minimal Memory.

Secondly, as explained in 1.1.1, sparse direct solvers need to adopt an ordering strategy which ensures a good sparse structure. However, introduction of the low-rank representations into these solvers made the ordering problem more challenging. In the low-rank solvers, reordering the unknowns to improve the compressiblity can degrade the sparse structure. Therefore, in this chapter we oer an ordering (and clustering) method of the unknowns to improve both the sparse structure and compressibility of the matrix.

Before starting to explain our solutions, we provide the necessary notations and background of this chapter in Section 3.1. In Section 3.2, we mention the related work on the matrix ordering problem. In Section 3.3, we introduce and detail the Projection heuristic [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF], which is our starting point to solve the two problems in this chapter. In 67 Section 3.4, we explain our improvements on this heuristic. In Section 3.5, we discuss on the experimental results and we conclude the chapter in Section 3.6.

Background

In this section, we recall the necessary background information on the matrix ordering problem. Ordering heuristics impact many aspects of the sparse low-rank direct solvers such as the ll-in amount, the sparse structure, the compressibility or the parallelism. In order to improve all these features at the same time, we need to take advantage of several techniques. As mentioned in Section 1.1.1, we adopt the nested dissection method as the coarse level ordering, since this method reduces the ll-in and provides a good parallelism by its nature. Let us explain this recursive method in more detail.

The nested dissection is illustrated on a very simple example in Figure 3.1. Here, two levels of nested dissection procedure are applied on a regular cube (on the left). This method rst splits the cube into two balanced parts through the rst level separator (in grey). Note that the size of the separator should be as small as possible to minimize the ll-in. When this separator is ordered with the largest number in the elimination tree, it ensures that there is no ll-in between both sub-parts. Afterward, the resulting separated parts are also divided into two sub-parts through the second level separators (in green and red), in the same fashion. On the right part of Figure 3.1, we see the rst level separator graph (in grey), where the nodes that interact with the second level separators (called traces) are colored in red and green, similarly to the left-most gure. Let us remind that both the sub-parts and the separators, which are obtained through the nested dissection algorithm, are called supernodes.

In practice, the nested dissection algorithm is applied recursively until the sub-parts are small enough to perform a local ordering heuristic on each of them. This local ordering can be applied for improving both the compressibility and the sparse structure at the same time. However, nding such an ordering is dicult.

In [START_REF] Pichon | Reordering Strategy for Blocking Optimization in Sparse Linear Solvers[END_REF], it has been shown that permutation within each separator can be applied without impacting the ll-in ratio produced by the nested dissection. Indeed, all sparse direct solvers consider the separator blocks originated from this heuristic as dense. We rely on this property to study ordering heuristics within each separator. Before providing more details about the reordering (and clustering) techniques for the separators, let us introduce some important notations for the remainder of the chapter.

As we mentioned in Section 1.1.3, the numerical factorization is performed, in a panel-wise fashion, in three steps: factorize, solve, update. In this respect, let us represent a symmetric matrix to be factorized, A, in a block form as

A = A 11 A 12 A 21 A 22 , (3.1)
where A 22 illustrates the rst level separator, while A 11 stands for all the other unknowns. A 12 and A 21 are the coupling parts, which represent the interactions between Figure 3.1: Two level nested dissection on a regular cube. On the left gure, the rst level grey separator and the second level red and green separators are represented. On the right gure, the graph of the rst level grey separator is shown. The red and green nodes in the right gure represent the second level red and green separator traces on the grey separator. This step factorizes the matrix A 22 In the following sections, we will focus on reducing both the HERK(A 21 , A 22 ) and the POTRF(A 22 ) operation costs through sparse structure and compressibility improvements. For the sake of simplicity, we will explain the reordering methods only on the rst level separator, although it can be applied recursively at each level of the nested dissection. Now let us explain the necessary information on two promising ordering and clustering techniques: K-Way clustering in Section 3.1.1 and Traveling Salesman Problem (TSP) [START_REF] David | The Traveling Salesman Problem: A Computational Study[END_REF] based clustering in Section 3.1.2. Note that the corresponding related work on these two methods are detailed in Section 3.2.

K-Way clustering

As mentioned before, in the sparse supernodal solvers, the clustering techniques can be performed on the separators without changing the ll-in. Performing the K-Way clustering on a separator aims to improve its compressibility, so that the POTRF(A 22 ) cost can be reduced. Let us explain the K-Way clustering through Figure 3.2. In the gure, we see the symbolic factorization (on the left) and the corresponding rst separator graph (on the right) of an 8 × 8 × 8 Laplacian, which is partitioned by Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF]. Here, the K-Way clustering is performed on the rst separator, resulting in four clusters that are represented in four dierent colors. The zoomed part in Figure 3.2a shows the number of external updates on each corresponding block.

As we can see in Figure 3.2b, the clusters consist of close nodes of the graph and the number of neighbor clusters is minimized. From now on, we will call this kind of clustering as regular. Having a regular clustering reduces the interaction of clusters, so it improves the compressibility of the interaction blocks. However, as the sparse structure of the sparse matrix is not considered in this method, the total number of updates on each block (seen in the zoomed part of Figure 3.2a) is high, increasing the HERK(A 21 , A 22 ) cost.

TSP based clustering

For improving the sparse structure in the solver, a separator clustering method which gathers similar contributions together should be used. In this respect, we adopt the TSP based reordering strategy to cluster the unknowns as in [START_REF] Pichon | Reordering Strategy for Blocking Optimization in Sparse Linear Solvers[END_REF]. From now on, for the sake of simplicity, we abbreviate this clustering as the TSP method. The aim of this method is to reduce the total number of o-diagonal blocks of the coupling matrix by gathering them together, which in return reduces the total number of updates. Note that this refers to the improvement of the HERK(A 21 , A 22 ) step. In Figure 3.3, we see the same information as Figure 3.2. However, in this gure, the separator clusters are generated by the TSP method. In Figure 3.3a, A 21 has a better sparse structure (dense parts stay continuously, not spread into smaller ones), so that the total number of updates on each separator block in the zoomed gure is reduced compared to the K-Way clustering. Therefore, the eciency of the HERK(A 21 , A 22 ) step is improved. Moreover, since fewer and larger blocks in the coupling part might increase the compressibility of A 21 , this method can further improve the memory footprint and computational complexity of the low-rank operations. However, as we can see in Figure 3.3b, compared to K-Way, the separator clustering is more irregular in a way that each color (cluster) is spread over the graph. Thus, clusters have more neighbors and larger diameters. This results in stronger cluster interactions, which reduces the separator compressibility and increases the POTRF(A 22 ) cost.

As we can see, neither the K-Way nor the TSP method provides an optimal solution. Therefore, beneting both of these methods can lower the POTRF(A 22 ) and HERK(A 21 , A 22 ) costs at the same time.

Related work

The related work on the low-rank solvers can be found in Chapter 1. In this section, we only discuss the literature on the clustering techniques.

In the dense low-rank solvers, the clustering aims to improve the compressibility. However, in the sparse case, in addition to the compressibility, having large enough data for the block operations is also crucial. It allows to exploit the modern architectures for a better performance. Therefore, nding an optimal clustering for the sparse low-rank solvers is a harder problem compared to the dense case.

In [RCL + 18], the authors explore some ordering techniques for hierarchical solver eciency through the StrumPACK solver. Here, they use the geometry of the problems to obtain the cluster trees. Although, they show a huge improvement in the compressibility rate compared to the dense solver, this work is not easily extendable to our context as we work algebraically without depending on the geometry of the problems. In [Beb08], spectral bisection based sparse matrix clustering is proposed, as well as oering the nested dissection for the low-rank clustering. The author also demonstrates the block admissibility in his work. In our context, using graph partitioning techniques on separator graphs is a similar idea to this work since we focus on the separator clustering.

In [START_REF] Chenhan | Geometry-oblivious fmm for compressing dense spd matrices[END_REF], an algebraic dense matrix clustering is proposed. For this purpose, the authors use a result from reproducing kernel Hilbert space theory in [START_REF] Hofmann | Kernel methods in machine learning[END_REF]. According to this result, any SPD matrix corresponds to a Gram matrix of vectors in an unknown Gram space. As a result, they dene their distances by using the matrix entries as inner products. By using these distances they generate a balanced cluster tree, where they avoid computing all the distances through sampling since it is costly for dense matrices.

The K-Way clustering is a classical low-rank clustering method to cluster the fronts or separators in sparse solvers. It is an extension from dense low-rank solvers, where it successfully improves the compressibility. The main idea of this clustering is to choose the diameter of the clusters as small as possible, and minimize the total number of cluster neighbors. In this way, the nodes with strong interactions are gathered together and the cluster interactions are minimized through a regular clustering. As a result, the cluster interaction blocks are more compressible, which in return improves the time-to-solution and memory footprint of the solver. In Mumps [AAB + 15, ABLM19b] and StrumPACK [GLR + 16, GLGR17], this clustering is adopted on the fronts. In our context, we perform it on the separators. Note that the K-Way clustering requires a connected graph, but fronts or separators are not always connected. Therefore, it is applied in two steps. Firstly, the vertices of the fronts or separators are reconnected with a halo distance 1 or 2 for obtaining a connected graph. Then, the K-Way clustering can be applied through partitioning libraries like Metis [START_REF] Karypis | Metisa software package for partitioning unstructured graphs, partitioning meshes and computing ll-reducing ordering of sparse matrices[END_REF] or Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF]. Although, this method provides compressible separators, it fails to provide a good sparse structure. In addition, large number of small o-diagonal blocks in this method may degrade the compressibility of the coupling part.

Let us emphasize that having large enough data sizes within the sparse solvers is very important. Methods like Reverse Cuthill-McKee (RCM) [START_REF] George | Computer Solution of Large Sparse Positive Denite Systems[END_REF] are widely used to reduce the number of o-diagonal blocks in the column blocks. However, this approach only considers the intra-supernode interactions, without taking the contributing supernodes into account. Therefore, we adopt the reordering method in [START_REF] Pichon | Reordering Strategy for Blocking Optimization in Sparse Linear Solvers[END_REF], where the authors convert the optimization problem of minimizing o-diagonal blocks number into the traveling salesman problem (TSP). This method provides a good sparse structure in the coupling part by gathering the similar contributions together. In this method, the clustering, following the TSP based unknown reordering, is obtained by using the smart splitting as in [START_REF] Lacoste | Scheduling and memory optimizations for sparse direct solver on multi-core/multi-gpu cluster systems[END_REF]. This smart splitting does not use strictly xed size clusters in order to decrease the total number of o-diagonal blocks. Note that although the TSP method in [START_REF] Pichon | Reordering Strategy for Blocking Optimization in Sparse Linear Solvers[END_REF] is more suitable to be parallelized, in [START_REF] Jacquelin | Fast implementation of the traveling-salesman-problem method for reordering columns within supernodes[END_REF] this method is improved to highly reduce its cost. The authors improve it by reducing the size of TSP and providing a more ecient way of computing the TSP distances. However, even though the TSP reordering provides a good sparse structure, it does not take into account the compressibility of the separator A 22 .

To conclude, the K-Way and TSP methods are promising for the sparse solver separator clustering in dierent perspectives. In order to obtain an optimal solution in terms of both compressibility and sparse structure, we can take advantage of these two methods at the same time. In our work, we propose to improve the Projection heuristic (see Section 3.3), which prots from both of these methods and gathers the uncompressible blocks together.

Projection heuristic

In Sections 3.1.1 and 3.1.2, the existing K-Way and TSP clustering methods and their aims are explained. In this section, we study the Projection heuristic introduced in [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF], which is the main focus in our work. This heuristic targets to improve both compressibility and sparse structure at the same time to outperform the existing clustering methods. In addition, it aims to determine the strong interaction blocks through the adjacency graph, which are poorly compressible. Let us start explaining this heuristic step by step.

Determining the traces

In this step, the vertices of the separator are ordered according to the set of contributions on them, which comes from the closest children in the elimination tree.

In Figure 3.4, we see a generic separator as an example. The gure illustrates the interactions of two levels of the closest children on this separator. These interactions represent the traces on the separator, which are resulted from the nested dissection (see Section 3.1). These traces are colored in red for the rst level children, while they are in green for the second level children. There are two red traces from the rst level children, and four traces from the second level ones. Each pair results in a chromosome shape, as seen in the gure.

In this step, rstly, the trace vertices (preselected vertices that represent the nodes belonging to the traces in the separator graph) are determined and ordered together.

The remaining vertices are in several well-separated groups, called connected components (grey and yellow color in the gure). Each isolated connected component gets the same set of contributions from the two levels of children considered in the elimination tree, so that they have identical sparsity pattern. Note that the connected components which have a smaller size than a given threshold (grey color in Figure 3.4) are gathered together to be large enough with respect to the minimal compressibility size. This component made out of the small ones is called as the remainder component.

Gathering the preselected vertices together and isolating the connected components provides two main benets. Firstly, it leads to a suitable clustering for the connected components, since it increases the distance between them. Thus, the compressibility of their interactions is improved. Secondly, the preselected vertices have strong interactions. Therefore, the blocks that represent their interactions are poorly compressible, and are identied as uncompressible within the Projection heuristic. Let us emphasize that the uncompressible blocks in our context do not represent only the uncompressible blocks but also the poorly compressible blocks for the sake of simplicity.

At the end of this step, representing our separator as

A 22 = A kk A ks A sk A ss (3.2)
the preselected vertices are ordered in A ss and others are ordered in A kk . Here, the only compressible part is A kk , while others are dened as uncompressible. In addition, the sub-diagonal blocks represent the connection of two consecutive column blocks.

Therefore, they are also considered as uncompressible since they have a high probability of having strong interactions. Note that the blocks with small sizes are considered to be non-admissible for low-rank compression, due to their low impact. From now on, all the small blocks and the dense diagonal blocks are automatically dened as non-admissible and never compressed to avoid unnecessary compression.

K-Way and TSP methods within the Projection heuristic

Initially at this step of the heuristic, there are a few large connected components, the traces and the remainder component. In order to reduce the size of the diagonal blocks, a clustering method is applied on each of them.

Each large connected component is clustered through the K-Way method to increase compressibility of the cluster interaction blocks. However, as the traces are considered to have strong interactions, K-Way is not performed on them. They are clustered through the TSP method to improve the sparse structure. The remainder component is also clustered through the irregular TSP method since it consists of small components that are already split in the graph.

Note that each connected component cluster (obtained through the K-Way partitioning) is further reordered through the TSP method. Here, the TSP method only serves as a reordering method to improve the sparse structure at the cluster level, without splitting the unknowns into dierent clusters. Now let us continue with some important implementation details of this heuristic.

Implementation details

As seen in Algorithm 9, at the beginning, the separator graph is isolated from the original one with the function ObtainSubgraph( C ). This routine generates direct connections at a distance 2 from the original graph to obtain the connected separator graph.

In the function ComputeTraces( C, l, d, w ), the traces are determined by checking each vertex of the separator. It is applied in two steps. Firstly, for each separator vertex, it checks the direct connections with the children for each depth level (less than or equal to l). If there is a distance less than or equal to d between them, this vertex is preselected. After preselecting these vertices, the vertices inside the separator, which have distance less than or equal to w to a preselected vertex, are also set as preselected.

Algorithm 9 Clustering and ordering the unknowns within the separator C in the original heuristic.

TSP( remainder component )

In the function IsolateConnectedComponents( C ), the connected components of the separator are isolated. It applies Breadth-First Search (BFS) algorithm before visiting each vertex one time. Every time this algorithm stops before all the vertices are visited, it means that another connected component is generated.

In the connected components loop, the small connected components are all merged together. Therefore, at the end of the loop, the large remainder component can be clustered through the TSP method. On the other hand, each large connected component is rstly clustered through the K-Way procedure, which is followed by a TSP based reordering at the cluster level. Finally, the TSP method is applied on the traces to improve the sparse structure.

In practice three main parameters are used to choose the preselected vertices: Projection depth, Projection distance, Projection width.

Projection depth, l, species how many levels of children should be considered to determine the traces on the separator. If this parameter is set to l, 2 l children in the elimination tree are used to determine the contributions on the vertices. Therefore, this parameter aects the number of traces. If l is too large, there will be a lot of small connected components and it will reduce the compressibility. Thus, it should be chosen small.

Projection distance, d, species the distance of the vertices, from a child to a separator, to be considered as part of the projection. Therefore, it aects the width of the traces. As d is getting larger, the number of preselected vertices will grow because more vertices are accepted as part of the projection.

Projection width, w, is used for widening the traces after their computation. The aim of this parameter is to get a better compressibility of the clusters on each side of the traces by better isolating them.

These parameters should be chosen carefully. That is, the blocks should be compressed as much as possible and there should be large enough connected components for further clustering. Therefore, small parameters should be adopted to reduce the number of preselected vertices and have large connected components. On the other hand, if the traces are not suciently large or if there are not enough number of traces, there will not be large number of well-separated connected components. In this case, the Projection heuristic cannot be exploited suciently.

In this heuristic, the preselected vertices are accumulated from the children, level by level (from the closest children to the farthest), based on the Projection depth (l) parameter. For a separator of size n, if the preselected vertices are more than Θ( √ n), the preselection terminates. Thanks to this stopping criterion, the risk of preselecting too many vertices is eliminated. In addition, let us represent the blocking size by b.

The Projection heuristic is only applied on separators larger than size 16b to suciently benet from it.

It is important to emphasize that the separators, which are obtained with the partitioning tools, are not always connected. That is why, in practice the separator graph is reconnected, using a distance 1 or 2, to allow the K-Way clustering on it. We adopt a distance of 2 for being safe, although in [AAB + 15] it is shown that a distance of 1 is enough for most of the graphs.

Discussions on the Projection heuristic

In the original work [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF], the author shows that the Projection heuristic can reach a 10% time improvement with only a 5% memory overhead compared to the existing K-Way method for the Minimal Memory strategy. However, the original heuristic is not optimal and can be improved. In this respect, we want to both validate the admissibility criterion and increase the performance.

In this chapter, one of our main focus is to empirically validate the admissibility criterion in the Projection heuristic. In Figure 3.5, we have three scenarios to identify the poorly compressible blocks. Here, pink color represents the preselected blocks. These blocks show the interaction with the preselected nodes. Therefore, they should be poorly compressible. Black color stands for the sub-diagonal blocks. As the sub-diagonal blocks represent two consecutive blocks in the matrix, they are possibly poorly compressible because these blocks represent close distances in the graph. Pink and black blocks are our interest in terms of compressibility. Brown color illustrates the preselected node cluster blocks. As seen in the gure, according to the number of preselected nodes and the blocking size, the preselected nodes can be split into dierent clusters. Blue color shows the diagonal blocks of the non-preselected node clusters, while the yellow color stands for all the remaining blocks.

In the left-most scenario, all o-diagonal blocks are considered compressible. In the middle gure, the sub-diagonals are considered as poorly compressible, while in the right-most one, both the sub-diagonals and the preselected blocks are identied as poorly compressible. By comparing these three scenario results, we want to validate the hypothesis on determining the poorly compressible blocks. In addition, we want Figure 3.5: Three dierent scenarios to identify the poorly compressible blocks. Pink color represents the preselected blocks, which are the interaction with the preselected nodes. Brown color illustrates the preselected node cluster blocks. Black color stands for the sub-diagonal blocks. Blue color shows the diagonal blocks of the non-preselected node clusters, while the yellow color is for all the remaining blocks. In the left-most scenario all the o-diagonal blocks are compressed, while in the middle gure the sub-diagonal blocks are kept in full-rank. In the right-most scenario both the sub-diagonal and preselected blocks are represented in full rank.

to experimentally prove that compressing the uncompressible blocks on the y highly improves the computational complexity. After this empirical validation step, we are interested in improving the heuristic in two ways.

Firstly, improving the compressibility of the uncompressible blocks is important when they are compressed on the y. It reduces the computational complexity of the update operations where they contribute in a low-rank form. In addition, in PaStiX, the blocks which are less compressible than a threshold are not compressed, because of their low impact. If we improve the compressibility of the uncompressible blocks, larger number of blocks are compressed as more blocks are eligible according to this threshold. Therefore, the resulting low-rank contributions improve the computational complexity of the operations. Thus, we propose to apply K-Way clustering on the traces for a better ops overhead of the solver. However, note that this improvement does not improve the memory peak since the uncompressible blocks are still allocated in full-rank at initialization. Secondly, as the small components are not necessarily close in the graph, merging them together degrades the separator compressibility, and therefore the computational complexity. Therefore, we propose to eliminate the remainder component and merge them with their neighbors in the graph.

To conclude, our aim is to introduce more regular clusters in the Projection heuristic.

In this way, we primarily expect to improve the computational complexity of factorization.

Improving the Projection heuristic

In the previous section, we discussed our motivation to change the original Projection heuristic in a way to have more regular separator clusters. In this section, we give the We can see two traces from the two children of the rst depth level. details of our changes on this heuristic. However, let us start by a quick comparison of all the clustering techniques that we study in this chapter. In Figure 3.6, we observe clusters on the rst separator graphs of a 120 × 120 × 120 Laplacian matrix with the proposed solutions.

In Figure 3.6a, we observe the TSP method on the whole separator. As we explained before, this method fails to cluster nicely the vertices. All clusters are spread and interleaved in the graph. Although, the outside contributions are gathered together to reduce the total number of updates, the internal clustering of the separator is not suitable to increase its compressibility.

In Figure 3.6b, we see the K-Way clustering on the whole separator, where the separator is clustered regularly to improve its compressibility. Note that in this clustering, each cluster is reordered through TSP at the cluster level to improve the sparse structure of the external contribution. Although the separator compressibility is better compared to Figure 3.6a, the cluster-wise TSP based reordering cannot reach a sparse structure as good as the separator-wise reordering. However, this strategy greatly reduces the cost of the TSP based reordering by applying it to smaller clusters. Last three gures stand for the Projection heuristic. In Figure 3.6c, we see the heuristic before applying clustering on the traces. Here, the Projection distance and width are set to 2 and 1, respectively. We can see the two rst level traces on the separator (in red). The well-separated connected components are regularly clustered through the K-Way method.

In Figure 3.6d, we see the same graph after TSP method is applied on each cluster and trace in the original heuristic. As the connected components are clustered before the reordering method, their regular clustering is not aected. On the other hand, as the trace clusters are obtained through the TSP method, their interaction compressibility is degraded.

Finally, in Figure 3.6e, we see our improved clustering. After Figure 3.6c is obtained similarly to the original heuristic, we merge the small connected components with their neighbors, as well as applying K-Way clustering on the traces. Following this step, the TSP based reordering is still applied on all the regular clusters to rene the sparse structure. Observe that our improved clustering is almost as regular as K-Way method, with an advantage of gathering (identifying) the high interaction nodes together.

Algorithm 10 Clustering and ordering the unknowns within the separator C in the improved heuristic. We start by determining the preselected vertices and isolating the connected components in the same way as in the original heuristic. At this point, we apply the K-Way clustering on the traces to improve the compressibility of the trace interaction blocks. Afterward, we isolate the connected components like in the original version.

In the connected components loop, we merge the small connected components with their neighbor clusters, while they are all merged together through the remainder component in the original heuristic. Therefore, we obtain more regular clusters on the separator, which improves the compressibility of the separator. We do it by checking all the separator vertices. Whenever we nd a neighbor vertex, we connect it to the large component of the determined neighbor vertex. On the other hand, the large connected components are clustered through the K-Way procedure like in the original heuristic.

At this point we apply the TSP based reordering on each cluster to improve the sparse structure at the cluster level similarly to the original version.

Experiments

In this section, we aim to demonstrate the impact of our improvements on the Projection heuristic. All the experiments are performed through the BLR supernodal direct solver PaStiX [PDF + 18], using the miriel nodes of the Plafrim1 supercomputer. Each miriel node is equipped with two Intel Xeon E5-2680v3 12-cores running at 2.50

GHz and 128 GB of memory. For the time results, we use 24 threads, one per core, with the default scheduler of the PaStiX library. The Intel MKL 2020 is used for the BLAS kernels. The minimum block width and height criteria to allow compression are set to 128 and 20, respectively. In the experiments, we use only LDL T and LU factorizations according to the input matrix features. For the SPD matrices, we avoid LL T factorization as the positive denite property can be aected by the compression.

In the experiments, we adopt the real-life matrices, which arise from a real mesh and are given in Table 3.1. They are taken from the SuiteSparse Matrix Collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF]. In addition to these matrices, we also use a Laplacian matrix of size 120 3 .

Let us emphasize that for the sake of simplicity, the TSP name represents the method, where the separator clusters are obtained through the TSP based reordering with smart splitting as in [START_REF] Pichon | Reordering Strategy for Blocking Optimization in Sparse Linear Solvers[END_REF] to gather similar contributions together. Similarly, the K-Way name stands for the K-Way low-rank clustering method followed by a TSP based reordering on each resulting cluster to improve the sparse structure.

We start our experiments by proving the eectiveness of our admissibility criterion and necessity of the late compression of the uncompressible blocks in Section 3.5.1. We then discuss the results of our improvements on the Projection heuristic, through the latest PaStiX version, in Section 3.5.2. Here, we observe the ops and time results on a large set of 31 real-life case matrices. Finally, in Section 3.5.3, we compare our improved heuristic to the existing methods.

Determining the block compression

As a starting point to our experiments, our rst interest is to numerically validate our admissibility criterion, as well as determining whether the poorly compressible blocks should be compressed on the y or not.

In Figure 3.7, the three grid columns stand for the three scenarios in Figure 3.5, respectively. We observe the ops ratio with respect to the full-rank solver for the preselected blocks are poorly compressible, the Projection heuristic, where these blocks are compressed on the y, has the best computational complexity. In Figure 3.7b, we observe the same conguration experiments, with a scenario where the poorly compressible blocks are not compressed at all. As no block is dened as poorly compressible in the left-most gure, the result is the same as Figure 3.7a through this conguration. However, observing Figure 3.7a and Figure 3.7b, we can clearly see that compressing the poorly compressible blocks on the y reduces the computational complexity compared to the scenario where they are not compressed at all. This improvement comes from the cheaper low-rank contributions of the uncompressible blocks compared to their full-rank versions.

To conclude, the Projection heuristic, where both the sub-diagonal and preselected blocks are compressed on the y, reduces the ops overhead on the rst separator computations compared to other clustering methods and congurations. Therefore, from now on, we will adopt this conguration in our experiments and call this optimized version as the OptProj heuristic.

Experiments on the OptProj heuristic

In this section, we compare the OptProj heuristic to the OptProj versions with our improvements. We do so by using a large set of real-life matrices through the last version of PaStiX. We set the Projection depth, distance and width parameters to 3, 1 and 1, respectively.

In Figure 3.8, the left-most gures are the factorization ops proles, whereas the right-most gures stand for the factorization time proles. The x-axis shows the ops and time overhead with respect to the minimum for each matrix, respectively. A method is better on average if its curve is closer to x = 1.

In the gure, the yellow curve stands for our OptProj heuristic, where we merged the small connected components with their neighbors. This OptProj version with Eliminated Remainder Component is abbreviated as OptProj+ERC. As we can see in the gures, the ops overhead of our contribution is very close to the basic adapted heuristic (in blue). However, in the time proles, we can clearly see the improvement of our new heuristic version. The reason of more time improvement compared to the ops dierence is because of the operation eciency dierences of Just-In-Time and Minimal Memory. That is, the small components are merged with their high interaction neighbor trace clusters, so that the uncompressible data is increased. As more data is compressed on the y and the Just-In-Time operations are more ecient than the Minimal Memory operations [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF], the factorization time of our heuristic is improved.

In the gure, the black curve stands for our OptProj heuristic, where we performed K-Way on the Traces (abbreviated as OptProj+KWT). As we can see in the gures, since our new heuristic degrades the sparse structure, it can even have worse computational complexities compared to the basic adapted heuristic, especially at lower precisions. This is caused by the larger number of smaller coupling matrix blocks, which are not eligible for compression. On the other hand, this method separates the trace nodes well. Thus, it avoids the interaction between these trace clusters. Then, the resulting null blocks of the separator through the OptProj+KWT strategy allows a time improvement in the gures, despite its computational complexities in total. In the time proles, we can see that OptProj+KWT is the fastest method for one third of the matrices.

3.5.3

Comparing OptProj+ERC heuristic to the existing solutions q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Factorization flops Factorization time Figure 3.9: The ops/time proles for real-life case matrices. The y-axis represents the number of matrices. The x-axis stands for the results with respect to the optimal for each matrix. The yellow curve stands for our OptProj heuristic with Eliminated Remainder Component, where the small components are merged with their neighbors.

In this section, we want to observe our contribution improvements compared to the existing methods. We have chosen the OptProj+ERC version as the most promising method through the Figure 3.8. Therefore, we will compare this version to the existing methods.

Observing the ops proles of the Figure 3.9, our OptProj+ERC method has a better computational complexity than the existing methods on average, as we expected. As TSP method does not provide separator compressibility, there are more operations in this method. On the other hand, although K-Way has a good separator compressibility, our heuristic gathers the strong interaction blocks together and does not compress them before the factorization. In this way, our heuristic provides more compressibility ratio in the compressible blocks, while avoiding the costly low-rank updates of the strong interaction blocks. As a result, the OptProj+ERC method has the best ops prole on average.

On the other hand, when we observe the time proles, we can see the eect of the sparse structure. The TSP method becomes more competitive in terms of factorization time, while it still cannot be the best average method because of its high ops overhead. It is interesting to observe that the OptProj+ERC heuristic is worse than the K-Way method on average, although it has a better ops prole. The reason is that determining the trace nodes is not trivial in practice. Choosing large enough traces to isolate the connected components and gather all the strong interaction unknowns together is crucial. However, when we gather the preselected vertices together, we can split the similar contribution blocks even more compared to K-Way. As a result, the OptProj+ERC heuristic cannot improve the average time overhead compared to the K-Way method at lower precisions.

Discussions

In this chapter, we improved the Projection heuristic [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF] to outperform the existing K-Way and TSP methods. The K-Way strategy provides regular clusters, where the close nodes of the graph are gathered together. This kind of clustering improves the compressibility of the separator cluster interaction blocks, which reduces the POTRF(A 22 ) cost. However, it does not provide large enough data sizes to take advantage of the underlying architectures. On the other hand, the TSP method can obtain fewer and larger o-diagonal blocks by gathering similar contributions together. As a result, it reduces the HERK(A 21 , A 22 ) cost. Although this method provides a good sparse structure, it does not improve the separator compressibility.

The Projection heuristic is proposed to get advantage of these two methods in an optimal way and provide an ecient admissibility criterion for the solver. In this heuristic rstly the strong interaction nodes are gathered (preselected) together to be able to determine the poorly compressible blocks. As they are uncompressible, the corresponding nodes are clustered through the TSP method to only improve the sparse structure. On the other hand, the non-preselected node groups are clustered through K-Way to increase the compressibility of the cluster interaction blocks. However, as this method is not optimal yet, we proposed some improvements for this heuristic in our work.

The experiments on the Geo_1438 matrix in Section 3.5.1 validated the eectiveness of the admissibility criterion of this heuristic. Moreover, these experiments showed the necessity of the late compression on the poorly compressible blocks. As a result, knowing that the uncompressible blocks should be compressed on the y, we proposed to improve the compressibility of these blocks. As the poorly compressible blocks are the interaction blocks of the trace nodes, we proposed to perform a regular clustering (K-Way) on them.

In addition, we proposed a new heuristic for the small connected components that were originally merged together. Indeed, as they are spread over the graph, their clustering is irregular. Therefore, we proposed to merge them with their closest neighbor for improving the regularity of the clusters, and thus their compressibility.

We provided our optimized Projection heuristic (OptProj) performance results

through the last version of the PaStiX solver. In the experiments of our contributions with regular clusters, we could observe an improvement compared to the original version of the heuristic. However, it has a lower impact than expected compared to the existing methods through the OptProj heuristic versions. It is because choosing the preselected vertices is a dicult problem in practice. Choosing a lot of vertices can reduce the compressibility, while not choosing enough vertices can lead the connected components to not be well-separated.

As another compressibility criterion, we can also change our perspective to use the ll-in levels concept of the incomplete factorization. In the next chapter, we explain this ll-in levels based heuristic and prove the eectiveness of this method.

Chapter 4 ILU Levels Based Preselection Heuristic

Direct solvers are widely used for their robustness, despite their high cost in terms of memory and time. As mentioned in Section 1.2, low-rank representations are introduced into direct solvers to tackle these problems. Through the block low-rank (BLR) representation within sparse solvers, dense blocks of the sparse matrix are approximated into lower dimensional matrices. This approximation successfully reduces the memory and/or computational complexity problems with a controlled precision loss. A detailed study on low-rank approximations of dense matrices using various methods can be found in Chapter 2.

As explained in Section 1.3, in sparse supernodal solvers, reducing both the time and the memory footprint through the low-rank representations is challenging. In these solvers, various block sizes are involved in the block updates. The cost of the low-rank update depends on the largest block involved, contrary to the cost of full-rank update that only depends on the contribution size. That is why, the update operation might be even more expensive than the full-rank one.

The sparse supernodal direct solver PaStiX oers two opposite strategies: favor memory peak reduction over time-to-solution (Minimal Memory), or prefer time-to-solution by delaying the data compression (Just-In-Time). As explained in Section 1.4, the former suers from the costly update operations, while the latter reduces it without improving the memory usage compared to the full-rank solver. Therefore, nding a compromise between time and memory footprint reduction becomes crucial.

Identifying the potential compressibility of each block is a key problem to benet from both strategies in PaStiX. By compressing the most compressible blocks as soon as possible, we can signicantly reduce the memory footprint of the solver. On the other hand, by delaying the compression of the poorly compressible blocks, we can get advantage of cheap update operations on them, without degrading the memory footprint much.

In Chapter 3, we study a heuristic to tackle this problem, where the unknowns are further ordered after the nested dissection algorithm. This heuristic aims at three improvements: 1) the uncompressible blocks are gathered together, 2) the remaining blocks have a better compressibility and 3) the sparse structure of the matrix is improved. The heuristic has two steps. Firstly, the data that represents strong interaction is gathered together. The resulting strong interaction blocks represent the poorly compressible blocks (named preselected blocks), for which Just-In-Time is more adapted, avoiding the overhead of Minimal Memory with expensive low-rank updates. In addition, the remaining blocks represent weak interactions, so that they are more compressible. At the second step of the heuristic, the data is further ordered to improve the compressibility and sparse structure of the matrix. For large matrices, the sparse structure improvement is necessary to reduce the number of memory accesses, and thus the time-to-solution. Unfortunately, this heuristic does not improve the solver as much as we expected and needs to be further improved.

In this chapter, we propose a new heuristic based on incomplete factorization to dene levels of admissibility (compressibility) for the blocks (see Section 1.2.1). Through this heuristic, we can successfully identify the compressible blocks and exploit the low-rank features of the solver. We prove the eectiveness of our heuristic with numerical experiments in both sequential and parallel environments. Note that the related work of this chapter can be found in Chapter 1.

In Section 4.1, we start by providing background details on the general idea behind incomplete LU factorization. In Section 4.2, we detail the new heuristic, which denes the non-compressible blocks to exploit the existing compression strategies in PaStiX. We present our numerical results in Section 4.3. Here, we rst prove the applicability of our heuristic as an admissibility criterion in Section 4.3.1. In Section 4.3.2, we analyse the experiments on a large set of matrices in sequential, while in the remainder of Section 4.3, the multi-threaded experiments are discussed. In Section 4.4, we conclude on the results of this work and its perspectives.

Background

The incomplete LU (ILU) factorization is a well-known method to get a general preconditioner for the iterative solvers [START_REF] Hénon | On nding approximate supernodes for an ecient block-ILU(k) factorization[END_REF]. It is an approximated version of the LU factorization, where part of the information is dropped [Saa03] and has the form A ≈ LU = LU + R. Here, the matrix R carries the negative values of the dropped elements.

Multiple dropping heuristics have been studied through the years targeting either parallelism and eciency of the implementation or a good numerical accuracy of the preconditioner. Among the dropping heuristics, there exist non-numerical solutions: using the position, or the ll-in levels [START_REF] Timothy | An unsymmetric-pattern multifrontal method for sparse LU factorization[END_REF], ILU(k), and numerical solutions: using a threshold [START_REF] Karypis | Parallel threshold-based ilu factorization[END_REF] value, ILU(τ ). Through the years, their block-based variants have been studied to improve their eciency [Gup17,BSV21,ARF + 19]. While we may consider low-rank solvers as a solution similar to ILU(τ ), we exploit in this study the block ILU(k) approach to enrich the block low-rank solver. The ILU method with the ll-in levels denition, as we use, is rst suggested in 1981 [START_REF] Watts | A conjugate gradient truncated direct method for the iterative solution of the reservoir simulation pressure equation[END_REF] and improved by a graph-based denition through ll-path theorem [START_REF] Rose | Algorithmic aspects of vertex elimination on directed graphs[END_REF] in [START_REF] Hysom | A scalable parallel algorithm for incomplete factor preconditioning[END_REF]. Figure 4.1 illustrates the idea of the ll-in levels that is strongly related to the ordering of the unknowns. On the left, the Fill-in entries that may occur during the numerical factorization are represented in red (level 1) and black (level 3). matrix non-zeroes pattern is represented, where the blue crosses are the original entries. On the right, the graph associated with this matrix is shown. During the numerical factorization, some entries may become non-zeroes (ll-in). On Figure 4.1, these entries are represented in red and black, both on the matrix and as new edges on the graph. We can dene the ll-in level as the length of the path connecting the two unknowns in the original graph. The path connecting 3 and 5 (and 3 and 4) in red goes only through 1 (resp.

2). Thus, the level of the ll-in between these two unknowns is 1. We can also see that 4 and 5 are connected at a level 3 (the path goes through 1, 3 and 2). As the ll-in level gets higher, the value of the new entry becomes smaller as it represents far interactions in the graph. That is, the ll-in levels can represent the generalized algebraic distance in the low-rank strong admissibility condition (see Section 1.2.1), without any knowledge of the geometry. Therefore, the ILU factorization can be implemented by dropping the values which have higher ll-in levels than a predened maximum level. Similarly, this procedure can be applied in a block-wise fashion. Thus, considering the block ll-in levels as an admissibility criterion for low-rank compression, we can algebraically decide on which blocks the compression should be delayed to reduce the overhead of the fully structured updates.

The new ll-in level based compressibility heuristic

As mentioned in Section 1.4, the Minimal Memory strategy suers from the complexity overhead of the LR2LR update kernel. We propose to exploit ILU ll-in levels to identify, at low cost, the blocks with large ranks. These blocks will increase the cost of the update step while providing only a small memory reduction. Once identied, it is possible to postpone the compression of these blocks as late as possible to replace the LR2LR kernels by LR2FR. This comes at the cost of a controlled memory overhead if the identication is correct. In Section 4.1, we mentioned that as the ILU ll-in levels get larger, the magnitude of the entries gets smaller. Thus, blocks with large level values should have smaller ranks and can be kept compressed to save memory, while blocks with small level values may have high ranks and are better candidates for delayed compression. Algorithm 11 presents the main steps to compute the ll-in levels of the blocks. This algorithm performs the same loops as the numerical factorization focusing only on the ll-in level information and the symbolic structure of the factorized matrix L. Initially, all blocks are considered with level 0, if they are part of the original matrix A, or ∞ if they are created by ll-in (lines 1-2). Then, the main loop updates the levels according to the formula given in [Saa03], which is adapted to the block-wise algorithm (Line 6). Note that as PaStiX uses the symmetric pattern structure of A+A T even for LU factorization, Cholesky-based algorithms are presented. For strongly non-symmetric matrices, the ll-in levels of the blocks in L and U are computed separately for better identication. This very cheap algorithm is in fact fully integrated within the parallel matrix initialization to avoid an extra loop over the structure, and such that its cost is completely hidden to the user.

( L ik , L jk , A ij )
Now that the ll-in levels are computed, the numerical factorization can be adapted to exploit this information. Algorithm 12 presents the proposed algorithm with a generic parameter maxlevel, which allows to set the new admissibility criterion (in orange color). First, lines 1 -3 compress the admissible blocks, which have a ll-in level larger than maxlevel. These blocks have the smallest ranks and are compressed before starting the numerical factorization loop (like Minimal Memory in red). Thus, they will be involved in LR2LR updates and are the most important ones to compress to reduce the memory footprint. On the other hand, the non-admissible blocks will be involved in LR2FR updates to reduce the ops count overhead while inducing a small memory overhead. These blocks are still compressed (lines 6 -8), just after the factorization of the diagonal block (like Just-In-Time in blue), to reduce the cost of the following operations: solve and updates. In the remainder of the chapter, we will refer to this BLR sparse factorization as ILU (k), with k the maximum level of the admissibility criterion. Note that choosing the right k value for a given problem is important. As a matter of fact, the larger the number of non-admissible blocks, the higher the memory overhead.

It is important to observe that ILU (-1) is the Minimal Memory scenario as all the admissible blocks are compressed during the initialization. On the opposite, ILU (∞) corresponds to the Just-In-Time scenario as all blocks are compressed only after all the cumulative updates are performed.

As a consequence, in this chapter, we propose to improve supernodal methods by trading a small memory overhead for lower ops count and a better time-to-solution.

For that purpose, we implement our solution in the sparse direct solver PaStiX [Pic18, PDF + 18], which supports the BLR compression scheme.

Experiments

In the following sections, the experiments are run for a set of 31 real case matrices taken in the SuiteSparse Matrix Collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF]. The experimental setup is similar to Section 3.5. The times shown are the average of 3 runs on each matrix. Data reported in the graphs are only related to the numerical factorization. The solve step is never considered as it is not impacted by this new algorithm. Interested readers can nd an experimental backward error study on the Minimal Memory and Just-In-Time strategies in [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF] that shows the numerical stability of the methods. Here, the author shows that both methods fulll the accuracy requirement given by users, despite the fact that Just-In-Time is more accurate thanks to the late compression scheme in the factorization. Note that as our new heuristic is an intermediate solution between Minimal Memory and Just-In-Time, it is also stable.

Compressibility statistics

We rst want to validate the hypothesis that using the ll-in levels is a good heuristic to classify the blocks based on their compressibility ratios. Figure 4.2 reports the cumulative memory consumption for all the 31 studied matrices in our experiments with three dierent tolerance criteria and by ll-in levels. The purple bars show the memory consumption of the structurally non-admissible blocks (too small to be compressed) and is identical for all precisions. The red part corresponds to the admissible blocks. The dark red is the memory footprint when they are compressed. It naturally increases with The bars report the cumulative memory of the 31 matrices. Purple represents the memory consumed by blocks below the size criteria, and red the memory of the admissible blocks for compression. Light red is the portion that can be saved when compressed. a higher precision. The light red shows the amount of memory that can be saved by compressing the blocks. One can observe that the lower the precision requirement, the higher the gain. These results show that the compression ratio of the admissible blocks increases with the levels. It conrms the original hypothesis that ll-in levels can help to better tune the admissibility criterion in order to save ops for a controlled memory overhead. Furthermore, this parameter needs tuning to adapt to the tolerance. One can see that for a tolerance of 1e -12 , only levels greater than 2 oer more than 40% memory savings, while all levels at 1e -4 provide more than 40% memory reduction. As a consequence, the ll-in level used to dene the admissibility criterion will need to be adapted to both the tolerance and the maximum memory overhead dened by the user.

Impact of the ll-in level heuristic on the sequential version

This section discusses the sequential experiments. Figure 4.3 shows the memory peak, factorization ops and factorization time proles obtained for dierent precisions. We study the impact of the rst ll-in levels (0 to 4) with respect to Minimal Memory (-1) and Just-In-Time (∞). Each curve represents the number of matrices within a percentage overhead of the best solution for each metric and matrix.

First, as expected, the lower the ll-in level chosen for admissibility, the lower the memory peak of the solver. One can observe that the impact of the ll-in level increases as the precision decreases. This conrms the trend already observed on shows percentages with respect to the best method for each metric and matrix, while y-axis represents the matrix count in cumulative way.

consumes up to 6.6 times more memory at 1e -4 , while it drops to 3.4 times at 1e -12 . Additionally, at this high precision, high levels of ll-in are able to reach the best memory peak. This means that potential ops reduction is possible without negatively impacting the memory. Second, when observing the ops count evolution, the results are naturally reversed. The higher levels of ll-in are better to generate less ops. One can observe that for low precision (1e -4 ), a level of 0 is enough to reach the same ops count as the Just-In-Time scenario (ILU (∞)). When increasing the precision, more levels need to be considered non-admissible to lower the ops count to its minimal value. Except some corner cases, levels 1 or 2 are enough for 1e -8 , and respectively 3 or 4 for 1e -12 .

Finally, the time proles follow the same trend as the ops proles, with larger dierences between the ILU (k) methods. This can be explained by the disparity of the LR2LR and LR2FR eciency, as well as their variation in number that may increase the phenomenon already observed on ops. Thus, one can observe that ILU (0) is the ordered by families and by increasing the memory ratio of ILU (-1) at a tolerance of 1e -12 . As we can observe, the compression ratio of all matrices varies a lot with the matrices and the tolerances. It also reects the fact that increasing the ll-in level does not necessarily mean an extra memory usage in high precision problems as the dots are merged together. On the right gure, this trend is summarized with boxplots representing the average gain.

One can observe that the Just-In-Time (ILU (∞)) strategy has the same memory peak as the full-rank version, and that the lower the level, the lower the memory peak. One can also observe that increasing the level of admissibility at low precision seems to greatly increase the memory footprint with respect to the best one. However, the impact remains moderate for high precision and it can be aorded to highly reduce the time-to-solution, as observed in Figure 4.3.

The summary on the right of Figure 4.4 conrms the fact that the memory consumption increases with the higher ll-in levels and with the higher precision. The results at 1e -12 for the new heuristic at levels 1 and 2 show interesting results. As a matter of fact, they provide in average 25% memory improvement compared to the full-rank version. Note that at this precision, even ILU (∞) does not manage to eciently accelerate the full-rank version on the ve right-most cases, which are poorly compressible. On the other matrices, in addition to this memory saving, levels 1 and 2 are also faster than the full-rank version.

To conclude, it is dicult to give a single level as the optimal solution. However, depending on the problem, as well as the precision and memory restrictions, the level can be tuned to provide a solution that outperforms the Minimal Memory strategy in terms of time, for a small controlled memory overhead. Moreover, it can even have a speedup compared to the Just-In-Time strategy, while reducing the memory footprint.

Impact of the ll-in level heuristic on the multi-threaded version

This section presents the results of the previous experiments in a multi-threaded environment with 24 threads. Figure 4.5 shows the time proles of the multi-threaded numerical factorization on the set of 31 matrices. Memory peak and ops are not reported as they are identical to the sequential ones.

We can observe that the impact of the new heuristic is even greater in the multi-threaded environment. The shift in the memory usage induced by the heuristic improves the memory bandwidth in the multi-threaded context and allows getting better performance. The ILU (k) heuristic performs better respectively with a level of 0, 2 and 4, for tolerances of 1e -4 , 1e -8 , and 1e -12 . The proposed solution outperforms the Just-In-Time strategy as it can be especially seen at 1e -4 . In this parallel context, the large memory reduction improves the memory contention of the threads, while it merely degrades the ops count with respect to the ILU (∞). ILU (0), which initially stores only the blocks of the original matrix A, clearly outperforms the other versions at low precisions. When increasing the precision (1e -8 or 1e -12 ), the extra ops count for small values of k degrades the performance. However, one can observe that ILU (2) and ILU (4) (at 1e -8 and 1e -12 , respectively) are good competitors with ILU (∞) in terms of time, while they still induce memory savings as opposed to ILU (∞). the ll-in level between 0 and 4, at 1e -8 and using 24 threads.

To better highlight the high gain obtained with our new heuristic, we compare it to the full-rank solver and to the previous low-rank strategies existing in the PaStiX solver.

Figure 4.6 presents the time proles of the multi-threaded numerical factorization with 24 threads on the left, as well as the memory peak at 1e -8 precision on the right. On shows that the new heuristics is up to 45 times faster than the full-rank solver, and is up to almost 2 times faster than the previous fastest implementation ILU (∞) except for the two CurlCurl matrices which would require higher ll-in levels to get better performance due to the high ranks of their contributions.

Overall, the proposed heuristic, if correctly tuned, is the fastest method for all matrices and outperforms the ILU (∞) solution while providing an eective memory footprint reduction. Furthermore, the memory footprint reduction is similar to the one obtained with the ILU (-1) strategy in two thirds of the cases (using a ll-in level of 0 in most cases), and reaches at most a 39% overhead in the last set of matrices. One can notice that in two cases (the two CurlCurl matrices), the new heuristic outperforms the ILU (-1) strategy thanks to a slight change in the order of allocation of the blocks. In these specic cases, the best memory peak is even obtained with a ll-in level of 1 that better delayed the workspace allocation to compress the blocks of the matrix. In this section we focus the study on one of the large matrices from the collection: the Serena matrix. The size of the matrix is N = 1 391 349 for 32 961 525 of non zeroes and it requires 28.6 TFlops to be factorized in full-rank, which makes it a good average test case of the collection. Figure 4.7 presents a performance study of dierent compression tolerances with 24 threads. The ILU(tuned) corresponds to the use of levels chosen for Tolerance ILU (tuned) 1e -4 ILU (0) 1e -6 ILU (2) 1e -8 ILU (2) 1e -10 ILU (4) 1e -12 ILU (4) each tolerance as reported in Table 4.1. The full-rank result is shown in the gure as a reference for the speedup observation. In conclusion, the new heuristic outperforms the former two solutions on all precisions with a small ll-in level of 0 to 4, relatively to the precision. As shown previously, at these levels, this solution even provides an important memory gain as opposed to the fastest existing solution ILU (∞). Table 4.2 reports the detailed factorization times of the Serena matrix with dierent numbers of threads and tolerances. The values inside parentheses present the speedup compared to the sequential run of each method at the corresponding precision. The ILU (k) heuristic levels are chosen as in Table 4.1. The level selection has been limited to the range 0 to 4 to ensure a useful memory saving compared to ILU (∞). That is why, with a small number of threads at 1e -12 , our method cannot run faster than ILU (∞). However, the ILU (k) heuristic benets from a higher scalability, which allows it to quickly outperform other solutions as the number of threads increases. These results, while not reported in this chapter, are similar to the ones observed on the set of matrices used for the experiments. To conclude, the ILU (k) heuristic, through the correct tuning of the levels, solves the original issue of the ops overhead of ILU (-1). Furthermore, thanks to its better memory footprint, it provides a better scalability in parallel environments and outperforms the fastest original solutions.

Study of the Serena matrix in multi-threaded environment

Discussions

The behavior of sparse supernodal direct solvers using low-rank compression highly depends on when the compression is performed. On one hand, all admissible blocks can be compressed before the factorization (as it happens with the Minimal Memory / ILU (-1) strategy). It allows high memory savings, but in the specic case of supernodal methods, it induces an expensive ops overhead during the low-rank updates. On the other hand, admissible blocks can be compressed after they have received all their updates (as for the Just-In-Time / ILU (∞) strategy). It reduces signicantly the ops count as in dense solvers, and thus the execution time. However, the memory peak of allocating all the L matrix blocks in full-rank, before the numerical factorization, still exists. In order to reduce this memory peak, the allocation of the blocks should be carefully delayed to their rst access.

In this chapter, we proposed a new heuristic to estimate the compressibility of each block and constructed an algorithm that is a compromise between the two strategies mentioned previously. The new heuristic, named ILU (k), identies poorly compressible blocks similarly to the ILU method, which identies the most important data. It relies on the block ILU ll-in levels to dene an algebraic distance to compute low-rank admissibility of the blocks. The purpose of dening the admissibility is to propose an intermediate solution that accelerates the Minimal Memory solution, while it slightly increases the memory consumption. Moreover, it gives the chance to tune the levels according to the precision, the matrix properties, and the characteristics of the machine to better exploit the advantages of both the Minimal Memory and the Just-In-Time strategies.

The experiments that we conducted on a large set of 31 real matrices demonstrated that the ILU (k) heuristic manages to identify eciently the most compressible blocks.

The solution proposed runs up to 5.2 times faster than Minimal Memory with only a 1.38 times increase of memory usage for high precision in both sequential and multi-threaded environments. Moreover, due to the elimination of the null blocks before the numerical factorization, the ILU (k) heuristic is also able to run 1.4 times faster than the Just-In-Time strategy in sequential, with a much lower memory consumption (0.23 times less). In the multi-threaded environment, it even goes up to 1.84 times faster and outperforms it for most of the cases. We showed that through correctly tuned ll-in levels, the ILU (k) heuristic can be used as an improved version of the existing strategies. It improves the numerical factorization in terms of both memory and time, and it improves the scalability for parallel environments.

Despite the high gain of the ILU (k) heuristic, it remains dicult to know beforehand which level will provide the best improvement. Section 4.3 has shown that only the rst levels are worth consideration, and that a clear trend appears on the level to use depending on the tolerance. implementation.

In Chapter 2, we proved the eciency of the QRCP methods on obtaining small enough ranks. Thus, we do not need to adopt more expensive approaches (like RQRRT) to reach closer results to SVD. For this reason, our further studies should mainly target to improve the performance, rather than the obtained ranks, compared to the QR methods we observed.

Block low-rank clustering In Chapter 3, we studied the Projection heuristic [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF]. This heuristic reorders the separator nodes to increase the separator compressibility and reduce the total number of updates within the solver. The former is applied through regular clusters, while the latter is done by gathering similar contributions together. It additionally aims to identify the poorly compressible (uncompressible) blocks, which represent the interaction blocks of the trace nodes.

In this work, we improved the Projection heuristic by conguring it correctly and clustering the separators in a regular way. We proposed two changes to obtain these regular clusters to increase the separator compressibility. Firstly, we merged the small node groups to a large neighbor cluster in the graph. In this way, the remainder component that includes the small node groups, which are spread along the graph, is eliminated. Secondly, we performed the K-Way (regular) clustering on the traces. Indeed, as we decided to compress the trace interaction blocks on-the-y, improving their compressibility is also crucial. In our experiments, we showed an improvement compared to the original heuristic. However, we could not obtain the results we expected.

In this chapter, we implemented the Projection heuristic, where only the small components are merged to their neighbors (OptProj+ERC). Additionally, we tested the heuristic, where only the traces are clustered through the K-Way partitioning (OptProj+KWT). However, we did not conduct experiments on the version, which implements both improvements at the same time. This should be our perspective at this point.

This chapter is devoted to improve the clustering after the coarse level nested dissection partitioning. Alternatively, we can focus on improving the nested dissection method itself. For example, in order to have similar contributions from children to the separators, the nested dissection procedure can be changed in a way to align the traces on the separators better [START_REF] Pichon | On the use of low-rank arithmetic to reduce the complexity of parallel sparse linear solvers based on direct factorization technique[END_REF]. As a result, we can increase the compressibility through well-separated large connected components, by getting advantage of the symmetric structure.

ILU levels based preselection heuristic In Chapter 4, our purpose is to identify the poorly compressible blocks with a dierent angle than Chapter 3. In this respect, we used the ll-in levels concept of the incomplete factorization in a block-wise manner. Then, according to our heuristic, a block is highly compressible if it has a high ll-in level.

Otherwise, it is uncompressible. In both sequential and multi-threaded environments, we proved the eciency of our heuristic.

In the experiments, we adopted predened values (k) to determine how many levels are considered poorly compressible. Through the experiments with dierent k values, we concluded that this value should be tuned according to memory and time restrictions of the applications. That is, smaller values provide reduced memory usage, while larger ones lead to a faster solver. However, determining the most advantageous level value is dicult. Therefore, as a perspective, we would like to introduce tuning techniques to automatically infer the best value of k depending on the given tolerance, as well as the properties of the machine and the matrix used.

In the future, we can use our ILU ll-in levels based compressibility information within the performance model of [START_REF] Marchal | Trading performance for memory in sparse direct solvers using low-rank compression[END_REF] as a parameter. In this way, we can take advantage of both heuristics at the same time.

As another future work, we can also identify the uncompressible blocks with another admissibility criterion. For example, we can set our criterion based on the k-way cluster distances. Then, the smaller the distance, the stronger the interaction (so uncompressible).

What is next? As mentioned before, thanks to the non-hierarchical clustering with similar size blocks, the BLR representation is simpler to implement and more convenient in parallel environments. For example, it is easier to arrange suitable block sizes to t into the memory of a shared-memory implementation. Similarly, we can reach a good load balance in a distributed environment through the blocks of similar size in the BLR representation. However, hierarchical schemes provide better asymptotical complexities compared to at ones. Therefore, some intermediate solutions between these techniques, like BLR-H [START_REF] Ida | Lattice h-matrices on distributed-memory systems[END_REF] and Multilevel-BLR [START_REF] Amestoy | Bridging the gap between at and hierarchical low-rank matrix formats: the multilevel BLR format[END_REF], are proposed to take advantage of both formats at the same time.

In BLR-H, the BLR scheme is used at the coarse level to obtain the blocks to be distributed. In this way, an ecient load-balanced distribution of the matrix is targeted. Then, the distributed parts are represented in a hierarchical way to improve the memory and computational complexity of the solver. In Multilevel-BLR, the BLR scheme is applied recursively to the full-rank diagonal blocks of the previous BLR representation. As a result, both BLR-H and Multilevel-BLR can provide complexities close enough to the hierarchical solvers (although not as good as them) and they are more promising to be ecient in the parallel context.

We already mentioned that our next step is to implement hierarchical schemes within the PaStiX solver to improve the time and the memory complexities. After this step, our long-term perspective should be to apply a hybrid scheme, like BLR-H and Multilevel-BLR. In this way, we can exploit both hierarchical and at formats at the same time. Here, we aim to oer a sparse direct solver which provides the best intermediate solution between the most ecient in terms of time/memory usage and the most exible one in parallel environments, compared to the other alternatives. Through these improvements, we target larger linear systems arising from dierent applications that need fast, accurate and memory ecient solutions in parallel environments. As a short term perspective, we should conduct a detailed study of our solver in mixed precision. We know that it is possible to control the accuracy, time and memory of direct solvers through dierent oating-point arithmetic. For example, adopting low precision arithmetics (like single or half precision) can highly speed up the solver and reduce the storage. However, low precision does not lead to high accuracy solutions. Alternatively, using higher precisions (like double precision) provides high accuracy despite the expensive time and storage requirements. We can use a mixed-precision solver to both ensure the required accuracy and reduce the solver cost. In this respect, we need a detailed comparison of our low-rank solver to the mixed-precision solver. Additionally, we need a clever implementation of the mixed-precision strategy within our low-rank solver to further improve it in terms of time and storage, while ensuring a given accuracy.

As another perspective, we can improve the low-rank feature of our solver to better adapt to the domain decomposition application needs. Domain decomposition methods [QV99,DW94,Mat08] can be used to reduce the time-to-solution of sparse linear systems. They are applied by using the underlying geometry and the partial dierential equation before the discretization step. Here, the main idea is to split the domain into subdomains, each of which has a new boundary condition on the interface. Then, the subdomains are solved separately in a way that the solution continuity is ensured through the interface. The solution of each discretized subdomain is performed through some algebraic solvers like direct methods, iterative methods or hybrid methods. When we reorder a subdomain matrix into interior and interface subparts, the former subpart is sparse and the latter is dense. Then, in the interior part, it is more convenient to adopt a robust solver with high accuracy. Thus, a low-rank solver with high accuracy threshold can be used for these parts. On the dense interface parts, we prefer an iterative solver to improve the time-to-solution, where a preconditioner is needed for the fast convergence. For the preconditioning aim, a low-rank direct solver with low accuracy threshold can be applied. As we can observe, we can use a low-rank direct solver in both the interior and interface parts within a domain decomposition application. However, the parameters of the solver should be adapted depending on the matrix parts that it is used on. In this respect, we should automatically tune the low-rank parameters of the PaStiX solver within the domain decomposition applications.

Figure 1 . 2 :

 12 Figure 1.2: Nested dissection applied on the adjacency graph. The factorized matrix structure is represented on the right. The corresponding quotient graph and elimination trees are shown on the left bottom.

Figure 1 . 3 :

 13 Figure 1.3: The three steps that are performed during the factorization of an arbitrary column block.

Figure 1 . 4 :

 14 Figure 1.4: Clusters obtained (on the left) and the corresponding cluster tree (on the right).

Figure 1

 1 Figure 1.5: Four clusters obtained by a at clustering (on the left) and the corresponding cluster tree (on the right).

1. 3

 3 Figure 1.8: Representation of two dierent low-rank update options (C -= AB). Here, the contributing blocks A and B appear respectively in red and in blue. The updated block C is colored in orange and the impact of the contribution (AB) is in purple.

Figure 1

 1 Figure 1.9: Accumulation of the U and V bases during the accumulation of two dierent size low-rank blocks. Here, the large C (brown) matrix gets the AB (purple) contribution of smaller size. Therefore, there is zero padding.

Figure 2

 2 Figure2.1: Panel-wise QR with column pivoting on the original matrix A m×n (on the left) and on the sample matrix B (b+p)×n (on the right). Here, b + p ≪ m, where b stands for the block size and p is the oversampling size, which is used to obtain more representative data on the sample matrix. The yellow panels represent the parts which are completed after the current panel iteration. The green parts illustrate the trailing matrices which will be used in the next panel iterations. The upper triangular yellow matrices stand for the factorized partial R matrices.

Figure 2 Figure 2 ▷

 22 Figure 2.5: Left-looking intra-panel update of QRCP. The updated column is in green, whereas the contributing parts are in yellow.

Algorithm 6

 6 Truncated Randomized Block QR with Column Pivoting 1: function TQRCP(A, ϵ, b) 2:

B

  , r] = QRCP (B [i] , (b + p)ϵ, b) ▷ Compute permutation 6:

  size (n -r -b) × (n -r -b), is avoided. This reduces the total trailing matrix update cost to O(nr 2 ) in TQRCP compared to O(n 2 r) in RQRCP. The complexity of computing W T G in TQRCP is O(n 2 b).

Algorithm 7

 7 Randomized Block QR with Rotation 1: function RQRRT(A, ϵ, b) 2:

Figure 2 . 8 :

 28 Figure2.8: The round-o error impact when the trailing matrix norm columns are accumulated starting from right or left. The y-axis stands for ||A-U K V T K || F in logarithmic scale, while x-axis stands for each index in the rst panel iteration. The result represents the S-shape matrix (see Section 2.6.1) of size 500 × 500, where the generation rank is 25.

Figure 2 . 9 :

 29 Figure 2.9: The k-rank matrix of size 500 × 500, where the generation rank is r G = 20n 100 = 100 at the generation precision ϵ = 10 -8 . The Frobenius norms of the trailing matrices at each index (

Figure 2

 2 Figure 2.10: The Z-shape and Z-shape-short matrices of size 500 × 500, where the generation rank is r G = 20n 100 = 100 at the generation precision ϵ = 10 -8 . The Frobenius norms of the trailing matrices at each index ( √ n i=K σ 2 i

Figure 2

 2 Figure 2.11: The S-shape and S-shape-short matrices of size 500 × 500, where the generation rank is r G = 20n 100 = 100 at the generation precision ϵ = 10 -8 . The Frobenius norms of the trailing matrices at each index ( √ n i=K σ 2 i

  Figure2.12: In the stability gures, the relative residual value with respect to the initial matrix,

Figure 2 . 13 :

 213 Figure 2.13: In the relative stability graphs, the relative residual value with respect to the SVDF,

Figure 2 . 14 :

 214 Figure 2.14: In the performance gures, GF lops = M inGF lop t

  Figure2.15: In the relative rank graphs, the relative rank with respect to SVDF, r r SV DF , for dierent matrix sizes and r G values is observed. In this gure, each row illustrates a dierent test case, whereas each column represents a dierent generation rank. Here, the stopping criterion is applied in the compression methods.

  Figure2.16: In the relative rank graphs, the relative rank with respect to SVDF, r r SV DF , for dierent matrix sizes and r G values is observed. In this gure, each row illustrates a dierent test case, whereas each column represents a dierent generation rank. Here, we applied the stopping criterion within the compression methods.

Figure 2

 2 Figure 2.17: In the accuracy gures, the error, ||A-UrV T r || F ϵ||A|| F

  Figure2.18: Time proles of the real-life case test matrices. The x-axis ratios are computed as tcurrent t min for each method, while the y-axis stands for each matrix as a number from 1 to 810. Here, the x-axis is in logarithmic scale for the sake of clarity of the gure.

  Figure2.19: Zoomed obtained rank proles of the real-life case test matrices. The x-axis ratios are computed as rcurrent r min for each method, while y-axis represents each matrix as a number from 650 to 810. The y-axis values (from 0 to 649) are omitted to zoom the gure since the methods get very close results at this interval.

A

  11 and A 22 . Using this representation, we can dene the factorization steps as 1. POTRF(A 11 ) This step factorizes the matrix A 11 2. TRSM(A 11 , A 21 ) This step solves the blocks in A 12 and A 21 3. HERK(A 21 , A 22 ) This step updates A 22 4. POTRF(A 22 )

  Figure 3.2: 8 × 8 × 8 Laplacian partitioned using Scotch and K-Way clustering on the rst separator. The gure on the left stands for the symbolic factorization structure. In the zoomed gure, the total update count on each block of the rst separator (A 22 ) is represented. The graph on the right shows the clustering of the unknowns inside the rst separator.

  Figure 3.3: 8 × 8 × 8 Laplacian partitioned using Scotch and reordering with smart splitting on the rst separator. The gure on the left stands for the symbolic factorization structure. On the zoomed gure, the total update count on each block of the rst separator (A 22 ) is represented. The graph on the right shows the clustering of the unknowns inside the rst separator.

Figure 3

 3 Figure 3.4: A generic separator with two levels of traces represented on it. First level separators are in red, while the second level separators are in green. There are two red traces from the rst level children, and four traces from the second level ones. Each pair separator leads to a chromosome shape. Small node groups (in grey) appear between the separators.

1

  : function OriginalProjection(C, l, d, w)

  Figure 3.6: The rst separator graphs of a 120 × 120 × 120 Laplacian matrix. For the Projection heuristic, the distance and width parameters are set to 2 and 1, respectively. We can see two traces from the two children of the rst depth level.

1

  : function Projection(C, l, d, w) ( blocksize, C i ) 11: for each cluster K i do 12: TSP(K i ) Algorithm 10 stands for our version of the Projection heuristic. Our improvements compared to the version [Pic18] is written in blue. Let us explain this algorithm, comparing to the original heuristic.

Figure 3 . 8 :

 38 Figure3.8: The ops and time proles for real-life case matrices. The y-axis represents the number of matrices. At the left-most gures, the x-axis stands for the ratio f lops method f lops min for each matrix. Similarly, right-most gures show the relative time ratios with respect to the minimum. OptProj+ERC stands for the OptProj method with Eliminated Remainder Cluster. Similarly, OptProj+KWT represents the OptProj version, where K-Way is applied on the Traces.

Figure 4 . 1 :

 41 Figure 4.1: An adjacency matrix (on the left) and its associated graph (on the right).Fill-in entries that may occur during the numerical factorization are represented in red (level 1) and black (level 3).

Figure 4

 4 Figure 4.2: Potential memory saving based on the tolerance criterion and ll-in levels.The bars report the cumulative memory of the 31 matrices. Purple represents the memory consumed by blocks below the size criteria, and red the memory of the admissible blocks for compression. Light red is the portion that can be saved when compressed.

Figure 4 . 3 :

 43 Figure 4.3: Memory peak, factorization ops and factorization time proles with dierent precisions for sequential runs. Each color stands for a dierent ILU (k) level. The x-axisshows percentages with respect to the best method for each metric and matrix, while y-axis represents the matrix count in cumulative way.

Figure 4

 4 Figure 4.5: Time proles of the multi-threaded runs for dierent precisions. Each color shows a dierent ILU (k) level result.

Figure 4 . 6 :

 46 Figure 4.6: Best results achievable for time and memory with ILU (k) heuristic, tuning the ll-in level between 0 and 4, at 1e -8 and using 24 threads.

Figure 4 .

 4 Figure4.6a, ILU (best time) refers to the best factorization time obtained with the new heuristic, where the ll-in level is taken in the range ILU (0) to ILU (4) included. Similarly, in Figure4.6b the ILU (best memory) stands for the best memory consumption reachable with a ll-in level in the same range. Note that on Figure4.6b, the memory consumption of the full-rank solver and the ILU (∞) are mingled, as they are identical. The left gure shows that the new heuristics is up to 45 times faster than the full-rank solver, and is up to almost 2 times faster than the previous fastest implementation ILU (∞) except for the two CurlCurl matrices which would require higher ll-in levels to get better performance

Figure 4 . 7 :

 47 Figure 4.7: Serena matrix time results of dierent precisions with 24 threads. The ILU (tuned) corresponds to the smallest level, which runs faster than ILU (∞).

  ij ) = min(lvl(L ij ), lvl(L ik ) + lvl(L jk ) + 1)

	Algorithm 11 Cholesky-based ILU ll-in levels initialization
	1: for all block L ij in L do	▷ Initialize the block ll-in levels in L
	2:	lvl(L ij ) = (A ij ̸ = 0) ? 0 : ∞
	3: for all column block L * k in L do	▷ Set the block ll-in levels
	4: 5:	for all block L ik in L * k do for all block L jk in L * k (with j > i) do
	6:	lvl(L

Table 4 .

 4 Table 4.1: Corresponding levels for ILU (tuned) at each precision. 2: Factorization times for the Serena matrix. Speedup with respect to the sequential runs are written inside parentheses.

	NbThreads	1e -4	1e -8	1e -12	Full-Rank

1.2. LOW-RANK REPRESENTATIONS

1.5. OTHER SOLVERS

2.7. EXPERIMENTS ON REAL-LIFE CASE MATRICES

https://www.plafrim.fr

Geo_1438 matrix at 1e -8 tolerance. Data reported on the graphs are only related to the rst separator. Each color represents dierent operation cost, while each bar stands for a dierent clustering technique. The right-most conguration only changes the Projection heuristic results as the preselected blocks, which represent the trace node interactions, only exist in this method. For the Projection heuristic, the depth, distance and width parameters are set to 3, 1 and 1, respectively. Figure 3.7a stands for the scenario where the poorly compressible blocks are compressed on the y. Here, the left-most scenario, where all o-diagonal blocks are compressed before any operation, is the most costly one for all the clustering techniques. This is because of the costly LR2LR update kernels that explained in Section 1.3. However, in the middle gure, when the poorly compressible sub-diagonal blocks are compressed on the y, the ops overhead of the low-rank operations is reduced. In the right-most gure, proving the hypothesis that both the sub-diagonal blocks and the Figure 4.4: Memory peak ratio of the ILU (k) heuristic with respect to full-rank on the 31 test matrices. On the left, the detailed information is presented for each matrix and precision. On the right, the information is summarized with boxplots showing minimum, maximum, median, rst quartile and third quartile for each ll-in level.

best average solution at 1e -4 . It even outperforms the Just-In-Time strategy by a factor up to 1.4×. To explain this performance, we recall that in the Just-In-Time strategy, many null-rank blocks are allocated and later compressed, while in the new ILU (k) heuristic they may never be allocated. Indeed, they are originally null blocks and at low precision they may receive only null contributions. Thanks to these savings, ILU (0) at this tolerance almost doubles the memory footprint in the worst case with respect to the Minimal Memory strategy, but it remains 0.23× the memory consumption of the ILU (∞) and the full-rank versions.

The observations in higher precisions are similar to the lower precision, but with higher ll-in levels. At 1e -8 , only the levels above 1 compete with the Just-In-Time strategy in terms of time, as well as providing a controlled memory overhead with respect to the best solution. At 1e -12 , the levels higher than 3 are required to get the best factorization time, which reduces the gain one can obtain on the memory footprint. However, solutions with level 1 is a good compromise at this precision. It is up to 8.5× faster than ILU (-1), with only up to 1.49× more memory usage.

Figure 4.4 presents the detailed ratios of the memory peak of the dierent levels of admissibility with respect to the full-rank solution. On the left gure, the matrices are

Conclusion and Future Work

Recently, the introduction of low-rank representations into sparse direct solvers allowed to successfully reduce their memory and/or time complexities. Block low-rank (BLR) format, compared to other alternatives, provides more exibility in parallel environments and it is simpler to implement. We believe BLR format is more advantageous in supernodal solvers in terms of memory since the fronts inherent to the multifrontal solvers are not allocated in these solvers.

In this work, we focused on improving the sparse supernodal BLR solver PaStiX. This solver adopts the nested dissection method to obtain the coarse level matrix ordering to generate more sparsity in the factorized matrix and provide better parallelism. Following the nested dissection procedure, it uses the BLR approximation scheme on the resulting large dense blocks of the factorized matrix structure, exploiting the existing sparsity. We contributed PaStiX in three ways, each of which is explained in a devoted chapter. Now let us provide the conclusions and perspectives of each chapter (so each improvement). Then, we explain the long-term future work.

Compression kernels In Chapter 2, we studied and compared four dierent QR based compression methods as an alternative to the costly SVD kernel. We implemented all the methods in a similar fashion with similar numerical stopping criteria. Let us emphasize that determining the rank through a numerical criterion is new for the TQRCP and RQRRT methods that we studied. As the blocks of BLR representation are not large enough to take advantage of parallelism when compressing them, PaStiX adopts sequential compression kernels. Therefore, we only observed the useful matrix sizes in our context and in a sequential environment.

We tested the QR based compression kernels by using some generated matrices. Here, we saw that each method ensures a good stability and accuracy. Thus, we mainly concentrated on the obtained ranks and performance results through both the generated and the real-life case matrices. These results showed that the QRCP method is the most promising QR alternative for the BLR format blocks.

In the current PaStiX solver, compression is performed on the blocks of size less than 300. In the future, the blocks could be much larger after the hierarchical scheme implementation. Thus, although QRCP seems to be the best kernel within the up-to-date PaStiX version, the TQRCP can outperform it within the future solver. Therefore, comparing these two kernels should be our perspective after the hierarchical scheme 103