On the analysis and evaluation of sparse hybrid linear solvers

A matrix is called sparse if many of its entries are zero. Such types of matrices are generated from discretization problems of many fields like numerical simulations, Fluid Dynamics, Graph theory, Optimization problems, Signal processing, Finance, industry, Linear Programming, Electromagnetics, 2D/3D and many other real-world applications. In general, we call a matrix sparse if we could exploit the sparsity of its elements in terms of memory storage and computation. Different sparse matrix formats are proposed which lead to huge memory and computation savings. Similarly, solving large sparse linear systems becomes increasingly important as a kernel task for such scientific applications. Furthermore, recent years witnessed huge and complex development in modern microprocessor architecture and hardware in general. Besides, the parallel computing methods and languages has shown a kind maturation in solving real-world problems especially, the development of Message Passing Interface (MPI). Consequently, many algorithms and software packages have emerged to exploit the new developments in hardware and software alike. For instance, the development in the hierarchy of memory and computation nodes leads to developing new blocking algorithms which accommodate the small cache memory of modern architecture and increase the floating-point performance. In general, there are two broad categories for solving linear systems mathematically: direct and iterative methods. Direct methods can give high accurate results (10−30) and are more suitable for small problems (current direct methods can solve up to a couple of millions of equations) because of the memory consumption they require (for example, they can not solve large sparse 3D problems which may generate hundreds of millions of unknowns). Besides, they have limited parallel scalability. Preconditioned iterative methods, on the other hand, are more robust, require less memory and easier to parallelize. However, they are problem dependent and can converge faster with good preconditioning methods. These methods are the methods of choice when an approximate solution of the problem is sought. Our focus in this thesis is the parallel software packages used for solving large sparse linear systems. We investigated the various methods and techniques for solving sparse linear systems and concentrate on the open source hybrid approaches as they are more flexible and well-suited the hierarchal structure of modern computers. The word hybrid has different meanings. Hybrid in programming means combining OpenMP and MPI so they work for shared and distributed memory. Hybrid in hardware means working on CPU and GPU. Hybrid in algorithms means combining direct and iterative methods. So, our focus is on the third meaning; combining direct and iterative algorithms in order to get more efficient methods for solving sparse matrices. Therefore, throughout this dissertation, when we mention the word hybrid, we mean direct/ iterative methods. Hybrid solvers are the latest research in developing robust and scalable methods for solving sparse linear systems. These methods combine direct and iterative approaches in certain ways, mostly using Schur complement framework, with the aim of getting the desired features of both direct and iterative methods especially the speed and low memory consumption of the iterative methods and the robustness of the direct methods. Most sparse hybrid solvers use the so called 'Schur complement framework' to combine direct and iterative methods. In this framework, the original matrix is divided into 2×2 block matrices. This approach is also known as sub-structuring domain decomposition method. Graph theory algorithms are responsible for ordering the global matrix and getting such block structure. The first block is a large block diagonal square matrix. Each diagonal block, called internal subdomain, is factorized using direct methods. The second and third block are the interfaces. If the matrix is symmetric, these interface blocks are transpose of each other. The fourth block is a square matrix called Schur complement matrix. This matrix has many desired features well-suited for iterative methods. It has smaller size than the original matrix and it is better conditioned than the original matrix. If the original matrix is symmetric positive definite, this Schur matrix will inherit this property and thus Conjugate Gradient can be used. Although Schur matrix is better conditioned than the original matrix, preconditioning the matrix before applying iterative method is necessary for faster convergence. The factorized internal subdomains are used to get an approximation of the Schur complement matrix and used as a preconditioner. Solving large sparse linear systems can take days even with algorithmic hybrid approaches if performed in sequential. Schur complement framework is more appealing for parallel programming. With increasing development of hardware and supercomputers, these systems can be solved efficiently within few seconds. However, existing algorithms should be modified to accommodate the parallel environment constrains such as scalability and load balancing. PDSLin and Maphys are the best existing public domain Schur-complement based hybrid solvers. They are based on different preconditioning methods which is a crucial ingredient in any iterative solver. PDSLin uses an approximation of the Schur complement as a precontitioner which gives a global view of the domain problem. Maphys uses an approximation of the local assembled Schur matrix which makes the solver more scalable. In our experiments, we thoroughly examined these solvers, test them on different matrix types and compare their results with the state-of-art Superludist direct solver. We also investigated the effect of tuning preconditioning input parameters on PDSLin and Maphys with increasing number of processors. Developed at Lawrence Berkeley National Laboratory (LBNL) by two distinguishing researchers in parallel linear algebra Ichitaro Yamazaki and X. Sherry Li; PDSLin solver is very robust and scales very well with increasing number of processors. However, it is very sensitive to the input parameters especially sparsifying tolerances which is not good. PDSLin is a powerful solver but has a lot of bugs and still needs a lot of work. Our results show that Maphys solver is more stable with the input values than PDSLin and gives better results. PDSLin results are unpredictable and sometimes fails with the slightest change of the input values. Maphys performance is better than both PDSLin and Superlu-dist in terms of time consumption and memory. Besides, the two-level parallelism of PDSLin is more robust than multithreading of Maphys. The serial partitioning time of Maphys is significant in many cases of our experiments and thus it is better to change into a parallel graph partitioner in Maphys than using a sequential partitioner. Sometimes the partitioning time is larger than summation of the other solution steps of Maphys and this is clearly shown in Audik and Freescale cases. Maphys solver also scales very well with increasing number of processors. Our conclusion is that sparse hybrid solvers are more flexible because they have different components and each component can be substituted to accommodate the problem at hand. The developers of the solvers we considered have already aware of this feature and this is clearly seen through the different methods they integrated within their solvers. For example, PDSLin uses either MUMPS or Superlu-dist as a direct method which are different approches of factorizing matrices. Similarly, Maphys uses either MUMPS or PASTIX which also are also different. Using two level parallelism can make hybrid solvers work efficiently and scalable up to thousand number of processors. In this level, Processors are distributed into levels and work concurrently and independently. This method alleviates the problem of increasing Schur complement size. However, load balancing is a challenging problem in this method. Hybrid solvers consume much less amount of memory than either pure direct or iterative methods.