# Solving Sparse Linear Systems fast

**URL:** <https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071>\
**Category:** Performance\
**Tags:** sparse, linearsolve\
**Created:** [June 20, 2022, 2:11pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071 "2022-06-20T14:11:41Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![pushkar\_khandare](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pushkar_khandare/32/37323_2.png) [@pushkar\_khandare](https://discourse.julialang.org/u/pushkar_khandare)\
**Post date:** [June 20, 2022, 2:11pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/1 "2022-06-20T14:11:41Z")

</div>

I want to solve **Ax = b** where A is a sparse matrix of size 10^5 x 10^5 (around 18 non-zero entries per row), **x** and **b** are vectors. **A** is constant throughout my simulation, while **b** changes at every step in my loop, and I have recalculate **x** once **b** changes.

I use LU decomposition for **A** and pass this to ldiv!. It takes around 0.05 seconds for the ldiv! step. Strangely (or maybe I am missing something) when I change the values in BLAS.set\_num\_threads I see no difference in the timings reported by @time macro.

I want to parallelize this in hopes of getting a faster solve time. Using CUDA.jl I am unable to use the backslash operator as it does not work for SparseMatrixCSR. The only thing I found was CUSOLVER.csrlsvqr! which takes a staggering 4seconds, which I profiled and noticed that 49% of time is spent on “void csqr\_leftLooking” and 47% on “CUDA memcpy HtoD”. Please note that I had already copied **(A,x,b)** to GPU via CuArray constructs and fed these to CUSOLVER.csrlsvqr!, so the HtoD memcpy is confusing to me.

My question is thus:  
How do I best solve a sparse linear system (on a CPU, or a GPU, whichever is faster for sizes upwards of 10^5 x 10^5) where the factors must be computed ONLY ONCE.

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [June 20, 2022, 2:16pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/2 "2022-06-20T14:16:37Z")

</div>

Hi, welcome @pushkar_khandare  
For // CPU, you may give a try to:

> **[GitHub - JuliaSparse/Pardiso.jl: Calling the PARDISO library from Julia](https://github.com/JuliaSparse/Pardiso.jl)**
>
> Calling the PARDISO library from Julia. Contribute to JuliaSparse/Pardiso.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [June 20, 2022, 2:23pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/3 "2022-06-20T14:23:21Z")

</div>

Did not try myself but you may also use  
[http://juliasmoothoptimizers.github.io/Krylov.jl/dev/gpu/](http://juliasmoothoptimizers.github.io/Krylov.jl/dev/gpu/)  
on GPUs.

I would be interested by your experiments 😉

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 20, 2022, 2:54pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/4 "2022-06-20T14:54:34Z")

</div>

[http://linearsolve.sciml.ai/dev/](http://linearsolve.sciml.ai/dev/)

You can run the whole gambit of methods, along with the documented preconditioners. Whether `KrylovJL_GMRES()`, `KLUFactorization()`, `UMFPACKFactorization()`, the GPU default, `MKLPardisoIterate()` are best is dependent on the sparsity pattern, so it’s hard to guess from this information.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [June 20, 2022, 3:00pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/5 "2022-06-20T15:00:26Z")

</div>

I believe that UMFPACK (facto and solve) is not parallel though.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 20, 2022, 3:10pm UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/6 "2022-06-20T15:10:41Z")

</div>

It uses the BLAS3 kernels, so it’ll multithread if you have the right sparsity pattern.

---

<div class="post-metadata">

**Author:** ![pushkar\_khandare](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pushkar_khandare/32/37323_2.png) [@pushkar\_khandare](https://discourse.julialang.org/u/pushkar_khandare)\
**Post date:** [June 21, 2022, 7:44am UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/7 "2022-06-21T07:44:22Z")

</div>

Thank you for your reply!  
Unfortunately when I applied for an academic license to use Pardiso I got an email stating that academic licences are not being given out starting January 2022.  
Haven’t been able to use Pardiso.jl yet.

---

<div class="post-metadata">

**Author:** ![pushkar\_khandare](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pushkar_khandare/32/37323_2.png) [@pushkar\_khandare](https://discourse.julialang.org/u/pushkar_khandare)\
**Post date:** [June 21, 2022, 8:38am UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/8 "2022-06-21T08:38:13Z")

</div>

Hello Chris,  
Thank you for your help!  
I am attaching the sparsity pattern of my matrix

 ![K_6x6_bc_sparsity](https://global.discourse-cdn.com/julialang/original/3X/0/c/0c3ac390cb49d7b294749370f38fc3238f4a30f8.png)  
This is for a smaller matrix. The first few and last few rows have diagonal entries only, and the rest come from the fact that I am solving a FEM problem on a square mesh.  
I do not know how to see this sparsity pattern and find an efficient preconditioner, or even make an educated guess of which methods might work.  
Any pointers are appreciated! Meanwhile I will go through all solvers from the link you shared.  
The matrices in my code will be of the order 10^5 x 10^5, with the sparsity pattern above.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 21, 2022, 9:01am UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/9 "2022-06-21T09:01:22Z")

</div>

> [@pushkar\_khandare](#):
>
> Unfortunately when I applied for an academic license to use Pardiso I got an email stating that academic licences are not being given out starting January 2022.  
> Haven’t been able to use Pardiso.jl yet.

You don’t need it for MKLPardiso.

MWE with LinearSolve.jl:

```julia
using LinearSolve, LinearSolvePardiso, SparseArrays, Random

A = sparse([1.0 0 -2 3
    0 5 1 2
    -2 1 4 -7
    3 2 -7 5])
b = rand(4)
prob = LinearProblem(A, b)
u = solve(prob, MKLPardisoFactorize(); cache_kwargs...).u

```

should work and automatically install any binaries needed.

> [@pushkar\_khandare](#):
>
> This is for a smaller matrix. The first few and last few rows have diagonal entries only, and the rest come from the fact that I am solving a FEM problem on a square mesh.

This has a regular sparsity pattern, so it should be a case where UMFPACK-type algorithms do well because they should be able to exploit the parallelism of the BLAS3 kernel. I’d try:

```julia
using LinearSolve, LinearSolvePardiso, SparseArrays
A = # your matrix
b = # test b

prob = LinearProblem(A, b)

for alg in (
    MKLPardisoFactorize(),
    MKLPardisoIterate(),
    UMFPACKFactorization(),
    KLUFactorization())

    u = solve(prob, alg).u
end

using AlgebraicMultigrid
ml = ruge_stuben(A) # Construct a Ruge-Stuben solver
pl = aspreconditioner(ml)
solve(prob1, KrylovJL_GMRES(), Pl = pl).u

using IncompleteLU
pl = ilu(A, τ = 0.1) # τ needs to be tuned per problem
solve(prob1, KrylovJL_GMRES(), Pl = pl).u

```

I should make that into a tutorial. You can see this in action in [https://diffeq.sciml.ai/stable/tutorials/advanced\_ode\_example/](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/)

---

<div class="post-metadata">

**Author:** ![pushkar\_khandare](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pushkar_khandare/32/37323_2.png) [@pushkar\_khandare](https://discourse.julialang.org/u/pushkar_khandare)\
**Post date:** [June 21, 2022, 9:24am UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/10 "2022-06-21T09:24:57Z")

</div>

Thanks for such a comprehensive answer Chris!  
I will go through the steps mentioned and post what I find.

---

<div class="post-metadata">

**Author:** ![pushkar\_khandare](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pushkar_khandare/32/37323_2.png) [@pushkar\_khandare](https://discourse.julialang.org/u/pushkar_khandare)\
**Post date:** [June 23, 2022, 9:04am UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/11 "2022-06-23T09:04:44Z")

</div>

Hello Chris,

After running through the examples you gave, and following the excellent caching interface of LinearSolve.jl, I found out that KLU factorisation gives me the faster solve times than UMFPACK (though the factorisation times are quite high for KLU, but since its a one-time operation it does not matter in my case).  
I have a key question. When I change the BLAS.set\_num\_threads value I see no difference in solve times. How can I vary the number of threads that these algorithms use, and choose the optimal number of threads? I want to run my code on my institute cluster, and want to use all resources I can to solve Ax=b in the fastest time possible.

---

<div class="post-metadata">

**Author:** ![rayegun](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rayegun/32/26729_2.png) [@rayegun](https://discourse.julialang.org/u/rayegun)\
**Post date:** [June 23, 2022, 11:22am UTC](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/12 "2022-06-23T11:22:30Z")

</div>

I don’t see a method to change the number of threads used for KLU in the User Guide. I’ll investigate that further, although it may be a compile time option. Off the top of my head I’m not sure which parts are parallelized.

Have you tried setting `OMP_NUM_THREADS` environment variable?