SPAI Preconditioners for HPC Applications

Abstract

Iterative methods to solve systems of linear equations Ax = b usually require preconditioners M to speed convergence. But the calculation of many preconditioners can be notoriously sequential. The sparse approximate inverse preconditioner (SPAI) has particular potential for high performance computing [1]. We have ported the SPAI algorithm to graphical processing units (GPUs) within NVIDIA's CUSP library [2] for sparse linear algebra. GPUs perform well on dense linear algebra problems where data resides for long periods on the device. Since the underlying minimization problems are independent, they are mapped to separate thread-blocks, and an optimized QR algorithm implemented using NVIDIA's CUBLAS library is employed on each.Traditionally the challenge has been to determine a sparsity pattern Sp(M) of the preconditioner dynamically [3], which reduces the condition number of MA to a degree where a preconditioned iterative solver such as GMRES becomes computationally viable. Due to the extremely high performance of the GPU, it is possible to consider initial sparsity patterns much denser than have been previously considered. We therefore consider a fixed sparsity patterns to simplify the GPU implementation.We evaluate the performance of the resulting preconditioner on a standard set of sparse matrices, and compare SPAI to other preconditioners. SPAI OverviewGiven a large, sparse system of equations Ax = b, the sparse approximate inverse technique [3] minimizes the Frobenius norm,which entails solving n decoupled least squares minimization problems, min m k Am k − e k , k = 1, . . . , n. We assume that M has some initial sparsity pattern, e.g., the same sparsity as A. Thus, for each k let J be the set of indices j for which m k (j) = 0, and denote that vector with its compressed representation,m k = m k (J ). Let I be the set of indices i such that A(i, J ) = 0, and = A(I, J ). Finally, defineê k = e k (I). The minimization then becomes the dense problem:The solution to this minimization ism k = R −1 Q Tê k , with = QR, the QR decomposition. The decoupled nature of the minimization problems ensures that all solutionsm k can be found in parallel. Implementation on Graphical Processing UnitsSPAI's underlying numerical kernel is the QR algorithm [4]. Empirical analysis of many sparse matrices indicates that the minimization matrices tend to be "tall and skinny", with rarely more than 200 rows. While there are several GPU QR implementations, e.g., in the MAGMA [5] library, these are optimized for large matrices, not for many independent small matrices. In the former case, the GPU should be applied to the entire matrix. In our case, each small minimization should be mapped to a thread-block.To this end, we have implemented several versions of QR for a single thread-block. Figure 1 clearly indicates that versions which are based on CUSP's own dense linear algebra routines are much less performant than those based on NVIDIA's CUBLAS library, and Householder QR is more effective than Givens QR. Figure 2 indicates that ...

Pipelined Iterative Solvers with Kernel Fusion for Graphics Processing Units

“…A substantial amount of research has been conducted on various preconditioning techniques for iterative solvers on GPUs including algebraic multigrid [Bell et al 2012;Gandham et al 2014;Richter et al 2014;Wagner et al 2012], incomplete factorizations [Li and Saad 2013;Naumov 2012], or sparse approximate inverses [Dehnavi et al 2013;Lukash et al 2012;Sawyer et al 2012]. Nevertheless, hardware-efficient and scalable black-box preconditioners for GPUs are not available, but instead the use of problem-specific information is required [Yokota et al 2011].…”