# Is it possible to parallelize matrix division

**URL:** <https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398>\
**Category:** Performance\
**Created:** [September 29, 2023, 1:44pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398 "2023-09-29T13:44:39Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [September 29, 2023, 1:44pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/1 "2023-09-29T13:44:39Z")

</div>

I have a big (~10000x10000) but sparse matrix A and a vector B. I want to solve the system using backslash

C = A\B ,

and it works, but it is taking too long. Is there a way to parallelize this and utilize multiple cores on the computer?

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [September 29, 2023, 3:24pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/2 "2023-09-29T15:24:56Z")

</div>

It’s possible that the default sparse solver may be parallel already (I’m not sure), or at least to the extent it’s possible/useful for the problem and algorithm.

But you might look into other solvers based on different algorithms. Take a look at some of the packages suggested in [this thread](https://discourse.julialang.org/t/solving-sparse-linear-systems-fast/83071/11). Whether or not they’re parallel, they may offer a speed benefit depending on your problem.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [September 29, 2023, 6:22pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/3 "2023-09-29T18:22:49Z")

</div>

Thanks a lot. UMFPACKFactorization() and KLUFactorization() improve the performance, but I do not see them parallelizing the code. I have checked that BLAS.get\_num\_threads() returns 4, but only one CPU is working on the problem. I have access to a computer with 64 cores, and it would be great if I could utilize it. I am planning to run a calculation with a much bigger size of matrices.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [September 30, 2023, 11:21pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/4 "2023-09-30T23:21:08Z")

</div>

The following code works well and gives a fantastic multithreading speedup. For some reason this code only works for julia version 1.8.2 and not for julia version 1.9.2. I get an error saying “reordering problem”.  
I think there is a bug in julia 1.9.2. Where should I report this?

```julia
using Pardiso
using SparseArrays
using Random
using Printf
using Test
# -----------------
verbose = false
n = 10000
A = sprand(n, n, 0.01)
B = rand(n)

# Initialize the PARDISO internal data structures.
# ps = PardisoSolver()
ps = MKLPardisoSolver()
set_nprocs!(ps, 4)
set_iparm!(ps, 12, 1)
println(get_iparm(ps, 12))
if verbose
    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
end

# If we want, we could just solve the system right now.
# Pardiso.jl will automatically detect the correct matrix type,
# solve the system and free the data
X1 = solve(ps, A, B)

# We also show how to do this in incremental steps.
#println(X1)

```

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [October 1, 2023, 12:31am UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/5 "2023-10-01T00:31:25Z")

</div>

Apologies for the drive-by, but I’m not at my computer right now. You might be interested in [this video](https://www.youtube.com/watch?v=JWI34_w-yYw) by @ChrisRackauckas .

---

<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:** [October 1, 2023, 12:58am UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/6 "2023-10-01T00:58:09Z")

</div>

That should just be `MKLPardisoFactorize()` / `MKLPardisoIterate()` in LinearSolve.jl . Try those and see if it works well.

> [@Shashank](#):
>
> UMFPACKFactorization() and KLUFactorization() improve the performance, but I do not see them parallelizing the code. I have checked that BLAS.get\_num\_threads() returns 4, but only one CPU is working on the problem. I have access to a computer with 64 cores, and it would be great if I could utilize it. I am planning to run a calculation with a much bigger size of matrices.

UMFPACK does parallelize, but OpenBLAS doesn’t use that many threads that well. If you use `using MKL` then UMFPACK may use threads better (and don’t forget to `BLAS.set_num_threads(64)`)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 1, 2023, 9:36am UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/7 "2023-10-01T09:36:09Z")

</div>

[Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl) provides multi-threaded sparse solves.

---

<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:** [October 1, 2023, 1:07pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/8 "2023-10-01T13:07:06Z")

</div>

BTW @Shashank if you can provide the setup for the matrix then we could add it to the SciMLBenchmarks to better track performance of this kind of case. Right now the best I could point to are the two benchmarks we have there:

- [LU Factorization Benchmarks · The SciML Benchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/LinearSolve/LUFactorization/)
- [Finite Difference Sparse PDE Jacobian Factorization Benchmarks · The SciML Benchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/LinearSolve/SparsePDE/)

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [October 1, 2023, 2:01pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/9 "2023-10-01T14:01:15Z")

</div>

Thanks Chris. I am trying to solve a linear integro-pde

```julia
Dr = Differential(r)
Dcth = Differential(cth)
Icth = Integral(cth in DomainSets.OpenInterval(-1,1))

# 2D PDE
eq = cth * Dr(u(r,cth)) + ((1 - cth * cth) / r) * Dcth(u(r,cth)) ~ cemit(r) - d(r) * u(r,cth) + cd(r) * Icth(u(r,cth))

```

where cemit(r), d(r) and cd(r) are functions that are inputs.  
This is just a steady-state solution for a transport equation in spherical geometry. Since this is a linear equation, it is easy to calculate the Jacobian. If we take a grid over r and \cos\theta that has 100 bins each, the jacobian is 10,000 by 10,000. Due to the integral in the equation, the Jacobian is not tridiagonal but broader. The sparsity pattern is shown below,  
 ![sparsity](https://global.discourse-cdn.com/julialang/original/3X/c/5/c54600f5668394f6a3fb32ed2312c06cec81b935.png)

Instead of O(30000) one would expect in a normal pde, I have about a million non-zero entries in the Jacobian. I am working on a project that involves adding nonlinear corrections to this equation, but right now, I am trying to figure out an optimal way of solving the linear part so that I can know which solvers work well in which situation. For now, it seems Pardiso.jl is the best. But I will make plots for all the methods I have tried and post them here so that we have a benchmark for integro-differential equations using single and multiple threads. I was also planning to do the same using a GPU node on the cluster to see whether that gives an advantage in this case.

---

<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:** [October 1, 2023, 2:11pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/10 "2023-10-01T14:11:04Z")

</div>

Is that not a banded matrix? If it’s a banded matrix you should specify it with BandedMatrices.jl and fully specialize on that.

---

<div class="post-metadata">

**Author:** ![Shashank](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashank/32/12323_2.png) [@Shashank](https://discourse.julialang.org/u/Shashank)\
**Post date:** [October 1, 2023, 2:18pm UTC](https://discourse.julialang.org/t/is-it-possible-to-parallelize-matrix-division/104398/11 "2023-10-01T14:18:15Z")

</div>

I did not know this was available. I will try that too. Thanks a lot.