I would assume scipy uses MKL, and Julia has wrappers for that too: https://github.com/JuliaSparse/MKLSparse.jl. For the GPU CUDA.jl has support for sparse matrices as well.
My only problem too is I just want to write high level mathematical code and have it work. With Python's anaconda distribution, I just download it and the sparse stuff works. I don't want the hassle of dealing with wrappers and code that is not very mature. I know it's a chicken and the egg thing.
Start a julia repl, hit ] for the package manager, paste `add MKLSparse`, write `using SparseArrays, MKLSparse`, run `sprand(100,100, 0.1) * rand(100)`
