The hierarchical interpolative factorization (HIF) offers an efficient way for solving or preconditioning elliptic partial differential equations. By exploiting locality and low-rank properties of the operators, the HIF achieves quasi-linear complexity for factorizing the discrete elliptic operator and linear complexity for solving the associated linear system. In this paper, the distributed-memory HIF (DHIF) is introduced as a parallel and distributed-memory implementation of the HIF. The DHIF organizes the processes in a hierarchical structure and keep the communication as local as possible. The computation complexity is O N log N P and O N P for constructing and applying the DHIF, respectively, where N is the size of the problem and P is the number of processes. The communication complexity is the latency and β is the inverse bandwidth. Extensive numerical examples are performed on the TACC Stampede system with up to 4096 processes. The numerical results agree with the complexity analysis and demonstrate the efficiency and scalability of the DHIF. A great deal of effort in the field of scientific computing has been devoted to the efficient solution of (2). Beyond the O(N 3 ) complexity naïve matrix inversion approach, one can classify the existing fast algorithms into the following groups. The first one consists of the sparse direct algorithms, which take advantage of the sparsity of the discrete problem. The most noticeable example in this group is the nested dissection multifrontal method (MF) method [14, 13, 24] . By carefully exploring the sparsity and the locality of the problem, the multifrontal method factorizes the matrix A (and thus A −1 ) as a product of sparse lower and upper triangular matrices. For 3D problems, the factorization step costs O(N 2 ) operations while the application step takes O(N 4/3 ) operations. Many parallel implementations [3, 4, 27] of the multifrontal method were proposed and they typically work quite well for problem of moderate size. However, as the problem size goes beyond a couple of millions, most implementations, including the distributed-memory ones, hit severe bottlenecks in memory consumption. The second group consists of iterative solvers [9, 29, 30] , including famous algorithms such as the conjugate gradient (CG) method and the multigrid method. Each iteration of these algorithms typically takes O(N ) steps and hence the overall cost for solving (2) is proportional to the number of iterations required for convergence. For problems with smooth coefficient functions a(x) and b(x), the number of iterations typically remains rather small and the optimal linear complexity is achieved. However, if the coefficient functions lack regularity or have high contrast, the iteration number typically grows quite rapidly as the problem size increases. The third group contains the methods based on structured matrices [7, 6, 8, 11] . These methods, for example the H-matrix [16, 18] , the H 2 -matrix [17] , and the hierarchically semi-separable matrix (HSS) [10, 38] , are shown to have efficient approximations of linear or quasi-linear complexity for the matrices A and A −1 . As a result, the algebraic operations of these matrices are of linear or quasi-linear complexities as well. More specifically, the recursive inversion and the rocketstyle inversion [1] are two popular methods for the inverse operation. For distributed-memory implementations, however, the former lacks parallel scalability [22, 23] while the latter demonstrates scalability only for 1D and 2D problems [1] . For 3D problems, these methods typically suffer from large prefactors that make them less efficient for practical large-scale problems. A recent group of methods explore the idea of integrating the MF method with the hierarchical matrix [25, 37, 35, 36, 15, 34, 19] or block low-rank matrix [31, 32, 2] approach in order to leverage the efficiency of both methods. Instead of directly applying the hierarchical matrix structure to the 3D problems, these methods apply it to the representation of the frontal matrices (i.e., the interactions between the lower dimensional fronts). These methods are of linear or quasi-linear complexities in theory with much small prefactors. However, due to the combined complexity, the implementation is highly non-trivial and quite difficult for parallelization. More recently, the hierarchical interpolative factorization (HIF) [20, 21] is proposed as a new way for solving elliptic PDEs and integral equations. As compared to the multifrontal method, the HIF includes an extra step of skeletonizing the fronts in order to reduce the size of the dense frontal matrices. Based on the key observation that the number of skeleton points on each front scales linearly as the one-dimensional fronts, the HIF factorizes the matrix A (and thus A −1 ) as a product of sparse matrices that contains only O(N ) non-zero entries in total. In addition, the factorization and application of the HIF are of complexities O(N log N ) and O(N ), respectively, for N being the total number of degrees of freedom (DOFs) in (2). In practice, the HIF shows significant saving in terms of computational resources required for 3D problems.
The performance of hierarchical interpolative decomposition in a distributed-memory environment is reported in REF .
Distributed-memory Hierarchical Interpolative Factorization
ArXiv
ArXiv
Mathematics Computer Science