A Direct Elliptic Solver Based on Hierarchically Low-rank Schur Complements

A parallel fast direct solver for rank-compressible block tridiagonal linear systems is presented. Algorithmic synergies between Cyclic Reduction and Hierarchical matrix arithmetic operations result in a solver with $O(N \log^2 N)$ arithmetic complexity and $O(N \log N)$ memory footprint. We provide a baseline for performance and applicability by comparing with well known implementations of the $\mathcal{H}$-LU factorization and algebraic multigrid with a parallel implementation that leverages the concurrency features of the method. Numerical experiments reveal that this method is comparable with other fast direct solvers based on Hierarchical Matrices such as $\mathcal{H}$-LU and that it can tackle problems where algebraic multigrid fails to converge

Hierarchical interpolative factorization for elliptic operators: differential equations

This paper introduces the hierarchical interpolative factorization for elliptic partial differential equations (HIF-DE) in two (2D) and three dimensions (3D). This factorization takes the form of an approximate generalized LU/LDL decomposition that facilitates the efficient inversion of the discretized operator. HIF-DE is based on the multifrontal method but uses skeletonization on the separator fronts to sparsify the dense frontal matrices and thus reduce the cost. We conjecture that this strategy yields linear complexity in 2D and quasilinear complexity in 3D. Estimated linear complexity in 3D can be achieved by skeletonizing the compressed fronts themselves, which amounts geometrically to a recursive dimensional reduction scheme. Numerical experiments support our claims and further demonstrate the performance of our algorithm as a fast direct solver and preconditioner. MATLAB codes are freely available.Comment: 37 pages, 13 figures, 12 tables; to appear, Comm. Pure Appl. Math. arXiv admin note: substantial text overlap with arXiv:1307.266

Optimal-complexity and robust multigrid methods for high-order FEM

The numerical solution of elliptic PDEs is often the most computationally intensive task in large-scale continuum mechanics simulations. High-order finite element methods can efficiently exploit modern parallel hardware while offering very rapid convergence properties. As the polynomial degree is increased, the efficient solution of such PDEs becomes difficult. This thesis develops preconditioners for high-order discretizations. We build upon the pioneering work of Pavarino, who proved in 1993 that the additive Schwarz method with vertex patches and a low-order coarse space gives a solver for symmetric and coercive problems that is robust to the polynomial degree. However, for very high polynomial degrees it is not feasible to assemble or factorize the matrices for each vertex patch, as the patch matrices contain dense blocks, which couple together all degrees of freedom within a cell. The central novelty of the preconditioners we develop is that they have optimal time and space complexity on unstructured meshes of tensor-product cells. Our solver relies on new finite elements for the de Rham complex that enable the blocks in the stiffness matrix corresponding to the cell interiors to become diagonal for scalar PDEs or block diagonal for vector-valued PDEs. With these new elements, the patch problems are as sparse as a low-order finite difference discretization, while having a sparser Cholesky factorization. In the non-separable case, the method can be applied as a preconditioner by approximating the problem with a separable surrogate. Through the careful use of incomplete factorizations and choice of space decomposition we achieve optimal fill-in in the patch factors, ultimately allowing for optimal-complexity storage and computational cost across the setup and solution stages. We demonstrate the approach by solving a variety of symmetric and coercive problems, including the Poisson equation, the Riesz maps of H(curl) and H(div), and a H(div)-conforming interior penalty discretization of linear elasticity in three dimensions at p = 15