# Eigs on huge SparseArrays doesn't use all BLAS threads

**URL:** <https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626>\
**Category:** Performance\
**Created:** [November 23, 2020, 1:25pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626 "2020-11-23T13:25:51Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahmed\_Abouelkomsan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_abouelkomsan/32/19010_2.png) [@Ahmed\_Abouelkomsan](https://discourse.julialang.org/u/Ahmed_Abouelkomsan)\
**Post date:** [November 23, 2020, 1:25pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/1 "2020-11-23T13:25:51Z")

</div>

Hi all,

I am dealing with huge sparse arrays of dimensions greater than 10^6. I am interested in finding only few extremal eigenvalues. To this end, I tried to use the Arpack package but I noticed that if the size of the sparse array get too big, the number of used CPUs in use decreases. For example

```julia
using LinearAlgebra
using Arpack
using SparseArrays
BLAS.set_num_threads(16)

A = sprand(100000,100000,0.001)

eigs(A,nev = 20,which=:SR)

```

Then I get this from the top command

```julia
 PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
38936 ahmed 20 0 19.901g 0.016t 133852 R 1597 2.2 21:09.33 julia

```

Now if I increase the dimension of A

```julia
using LinearAlgebra
using Arpack
using SparseArrays
BLAS.set_num_threads(16)

A = sprand(1000000,1000000,0.001)

eigs(A,nev = 20,which=:SR)

```

Then the top command reads

```julia
PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
38936 ahmed 20 0 34.919g 0.030t 133852 R 155.7 4.1 44:02.02 julia

```

One notices that not all CPUs are being used. The decrease is in fact drastic! Any explanation for this weird behavior or what might be wrong? This has been tested on two machines, one with 12 physical cores and one with 24 physical cores.

I also tested this on similar packages like KrylovKit and it also shows the same behavior which might be related to how Julia itself deals with huge operations on huge sparse arrays. Any help is appreciated

---

<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:** [November 23, 2020, 1:36pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/2 "2020-11-23T13:36:56Z")

</div>

Sparse-matrix operations don’t (can’t) use BLAS. BLAS is only for dense-matrix operations.

---

<div class="post-metadata">

**Author:** ![jagot](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jagot/32/12217_2.png) [@jagot](https://discourse.julialang.org/u/jagot)\
**Post date:** [November 23, 2020, 3:03pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/3 "2020-11-23T15:03:41Z")

</div>

If you’re using a Krylov method (or any other method that relies on matrix–vector products), then perhaps this package: [https://github.com/jagot/ThreadedSparseArrays.jl](https://github.com/jagot/ThreadedSparseArrays.jl) could be useful? It tries to speed-up the multiplication by threading over the output columns, so for matrix–vector products, you have to use `mul!(y, transpose(At), x)` for `At = ThreadedSparseMatrixCSC(A)`, where `A` is your sparse matrix.

The package is mostly a proof-of-concept right now, but seems to work.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 23, 2020, 5:44pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/4 "2020-11-23T17:44:03Z")

</div>

BLAS is used for things like orthogonalization and scalar products in Krylov methods, but direct sparse solvers can’t use it in a significant way ( I think ).

---

<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:** [November 23, 2020, 5:46pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/5 "2020-11-23T17:46:20Z")

</div>

> [@ctkelley](#):
>
> BLAS is used for things like orthogonalization and scalar products in Krylov methods

That’s true, but my impression is that usually the matrix-vector products account for the bulk of the time (unless you are building up a huge subspace without restarting).

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 23, 2020, 6:01pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/6 "2020-11-23T18:01:00Z")

</div>

```julia
@stevengj

 That’s true, but my impression is that usually the matrix-vector products account for the bulk of the time (unless you are building up a huge subspace without restarting).

```

I’ve had cases where we had to do a few hundred Krylovs/Newton because our preconditioner, while scalable, was still pretty bad. The Jacobian-vector product was expensive, but the orthogonalization was expensive enough to be painful. Someone from the Trilinos project pointed me at classical Gram-Schmidt (twice!) for the orthogonalization, and that helped a lot.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 23, 2020, 8:22pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/7 "2020-11-23T20:22:36Z")

</div>

I don’t have much experience with sparse matrices, can you expand on why sparse algorithms (eg spmv) can’t use threading? Because they’re bandwidth-bound?

---

<div class="post-metadata">

**Author:** ![jagot](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jagot/32/12217_2.png) [@jagot](https://discourse.julialang.org/u/jagot)\
**Post date:** [November 23, 2020, 8:27pm UTC](https://discourse.julialang.org/t/eigs-on-huge-sparsearrays-doesnt-use-all-blas-threads/50626/8 "2020-11-23T20:27:29Z")

</div>

I don’t think it’s so much that they can’t, but rather that Julia’s built-in sparse library has not been threaded yet. See e.g. [this open PR](https://github.com/JuliaLang/julia/pull/29525) that the package I mentioned above was based on/inspired by.

Then of course, depending on which storage you use (column- or row-major), some operations are easier to parallelize than others.