Dense matrix algebra has the benefit that besides the matrix dimensions, the data itself has no impact on performance: the access pattern is entirely predictable and known as soon as you know the matrix sizes (side note: which is why just-in-time libraries such as libxsmm exist). For most cases, you can rely on MKL or Eigen to produce near-optimal performance for your dense linear algebra operations.

Sparse matrices do not have that luxury. Sparse matrix ops have access patterns that entirely depend on your data. You simply cannot make a library that covers every possible distribution of non-zeroes in an optimal way. At best you have a toolbox of a large number of representations and algorithms that you choose for each different use-case. Engineering simulations and solvers can have block or cluster patterns of non-zeroes, which can be exploited by doing dense operations for the dense blocks. Large real-world graphs often have an exponential non-zero distribution and are often extremely sparse, which might need an hash-based join algorithm for an axpy-kernel. Are you running a tight CG-solver loop? Then you might want to bundle multiple spmv-kernels to operate multiple vectors on the same matrix (for similar reasons as states in the original article). You are distributing with Spark? Unless your non-zero distribution is uniform, your matrix is really hard to partition evenly. If your matrix never changes, you may consider preprocessing using a partitioner (METIS, Mondriaan).

In short, it is highly unlikely that the library you are using is anywhere near optimal for your use case. Do explore alternatives based on your domain-insights, it is not uncommon to net double digit performance improvements.