Programming Matrices as Staged Sparse Rows to Generate Efficient Matrix-free Differential Equation Solver

Solving differential equations is a critical task in scientific computing. Domain-specific languages (DSLs) have been a promising direction in achieving performance and productivity, but the current state of the art only supports stencil computation, leaving solvers requiring loop-carried dependencies aside. Alternatively, sparse matrices can represent such equation solvers and are more general than existing DSLs, but the performance is sacrificed. This paper points out that sparse matrices can be represented as programs instead of data, having both the generality from the matrix-based representation and the performance from program optimizations. Based on the idea, we propose the Staged Sparse Row (SSR) sparse matrix representation that can efficiently cover applications on structured grids. With SSR representation, users can intuitively define SSR matrices using generator functions and use SSR matrices through a concise object-oriented interface. SSR matrices can then be chained and applied to construct the algorithm, including those with loop-carried dependences. We then apply a set of dedicated optimizations, and ultimately simplify the SSR matrix-based codes into straightforward matrix-free ones, which are efficient and friendly for further analysis. Implementing BT pseudo application in the NAS Parallel Benchmark, with less than $10\%$ lines of code compared with the matrix-free reference FORTRAN implementation, we achieved up to $92.8\%$ performance. Implementing a matrix-free variant for the High-Performance Conjugate Gradient benchmark, we achieve $3.29\times$ performance compared with the reference implementation, while our implementation shares the same algorithm on the same programming abstraction, which is sparse matrices.