Jacobi or Gauss-Seidel are not really state of the art for solving systems of linear equations. It is more done by preconditioned conjugated gradient (for symmetric positive semi-definite matrices) and preconditioned (F)GMRES (or other Krylov subspace methods) for arbitrary matrices.
The crucial part here is the preconditioner. There was recently (=21th century) huge progress with very efficient preconditioners for Laplacian systems. In Hardness Results for Structured Linear Systems by Rasmus Kyng and Peng Zhang, it is shown
that if the nearly-linear time solvers for Laplacian matrices and their generalizations can be extended to solve just slightly larger families of linear systems, then they can be used to quickly solve all systems of linear equations over the reals.
They go on to say:
This result can be viewed either positively or negatively: either we will develop nearly-linear time algorithms for solving all systems of linear equations over the reals, or progress on the families we can solve in nearly-linear time will soon halt.
But even beyond such deep theory, one can say something about the computational cost of a preconditioner and its efficiency. A matrix is hard to solve, if its condition number is big. The task of the preconditioner is to reduce the condition number, by approximately inverting the matrix. The (asymptotic) computational cost (vs efficiency) can estimated based on separators of the sparse matrix. For example, to separate a N x N x N grid in 3D, you need a separator of size O(N^2). If you would not approximate the system on the separator, you would get a full matrix of size(N^2), i.e. N^4 non-zero entries. So you need to approximate the clique graph on the separator with a sparse (directed) graph, probably an extender graph. One fast preconditioner algorithm for Laplacian systems approximates the clique (among the neighbors of the eliminated node, which has to be added to the graph) in each elimination step by a tree. This is good enough to keep the computational cost under control. The efficiency of this preconditioner in reducing the condition number is also often extremely good in practice, but proving this theoretically would still be a nice achievement.
A recent result on fast solving sparse linear systems replaced the Krylov subspace methods mentioned above by "an efficient, randomized implementation of the block Krylov method": Solving Sparse Linear Systems Faster than Matrix Multiplication by Richard Peng and Santosh Vempala. Even so the paper title says "solving ... faster than matrix multiplication," this does not answer the question asked here, because the paper means dense matrix multiplication, but the question here is about sparse matrix multiplication.
