# CUDA eigenvalues of a sparse matrix
**URL:** https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851
**Category:** GPU
**Tags:** question
**Created:** [September 18, 2020, 2:34pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851 "2020-09-18T14:34:23Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@8me](https://discourse.julialang.org/u/8me)
#### Post date: [September 18, 2020, 2:34pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/1 "2020-09-18T14:34:24Z")
I’m trying to optimize the eigenvalue (and eigenvector) calculation for a neutrino oscillation package which is currently under development ([https://github.com/KM3NeT/Neurthino.jl](https://github.com/KM3NeT/Neurthino.jl)).
The problem consists of a lot of small symmetric matrices where the eigenvalues should be determined and I thought maybe I can utilize the GPU for that problem. In order to pass this as a single problem to the GPU via a sparse matrix having all the small matrices on its trace.
For testing I created this matrix:
```julia
A_full = Symmetric(rand(3000, 3000))
```
When doing it in the “classic” way I get:
```julia
@benchmark eigs(A_full)
BenchmarkTools.Trial:
memory estimate: 808.38 KiB
allocs estimate: 2571
--------------
minimum time: 519.974 ms (0.00% GC)
median time: 527.320 ms (0.00% GC)
mean time: 528.424 ms (0.00% GC)
maximum time: 546.705 ms (0.00% GC)
--------------
samples: 10
evals/sample: 1
```
I also refined the example with respect to my initial problem description and created this example:
```julia
dim = 3
n = 1000
A = spzeros(dim*n, dim*n)
for i in 1:n
A[(i-1)*dim+1:i*dim, (i-1)*dim+1:i*dim] = Symmetric(rand(dim, dim))
end
```
and when calculating the eigenvalues using `Arpack.jl` I get:
```julia
@benchmark eigs(A)
BenchmarkTools.Trial:
memory estimate: 810.03 KiB
allocs estimate: 1889
--------------
minimum time: 24.614 ms (0.00% GC)
median time: 27.885 ms (0.00% GC)
mean time: 28.766 ms (0.12% GC)
maximum time: 47.153 ms (0.00% GC)
--------------
samples: 174
evals/sample: 1
```
Using the CUDA package I get:
```julia
A_cuda = CuArray{ComplexF64}(A)
@benchmark CUDA.@sync CUDA.CUSOLVER.heevd!('V','U', A_cuda)
BenchmarkTools.Trial:
memory estimate: 1.23 KiB
allocs estimate: 33
--------------
minimum time: 1.768 s (0.00% GC)
median time: 1.769 s (0.00% GC)
mean time: 1.770 s (0.00% GC)
maximum time: 1.773 s (0.00% GC)
--------------
samples: 3
evals/sample: 1
```
As this is already slower than the first example it leads me to the question whether I have to change my approach?
If I stick to my approach I guess I have to change it to `CUDA.CUSPARSE.CuSparseMatrixCSR(A)` and use specific eigenvalue functionality for sparse arrays, because I have seen there is a function `CUDA.CUSOLVER.csreigs`, but I don’t really get how to use it (with respect to its arguments).
Btw. all the testing was done using a NVIDIA GTX 1080 Ti.
Can you help me out at that point?
Thanks.
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [September 18, 2020, 3:44pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/2 "2020-09-18T15:44:18Z")
I don’t know if that makes a difference but `A` is `Float64` while `A_cuda` is complex right?
---
### Author:  [@8me](https://discourse.julialang.org/u/8me)
#### Post date: [September 18, 2020, 4:31pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/3 "2020-09-18T16:31:58Z")
Yes indeed it is complex, because `heevd` has the following signature:
```julia
help?> CUDA.CUSOLVER.heevd!
No documentation found.
CUDA.CUSOLVER.heevd! is a Function.
# 2 methods for generic function "heevd!":
[1] heevd!(jobz::Char, uplo::Char, A::CuArray{Complex{Float64},2}) in CUDA.CUSOLVER at /home/test/.julia/packages/CUDA/dZvbp/lib/cusolver/wrappers.jl:720
[2] heevd!(jobz::Char, uplo::Char, A::CuArray{Complex{Float32},2}) in CUDA.CUSOLVER at /home/test/.julia/packages/CUDA/dZvbp/lib/cusolver/wrappers.jl:720
```
I chose the comparing examples as I did, because my problem has real numbers. Do you know how to do it in CUDA with real numbers, then I will update my numbers on the benchmark!?
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [September 18, 2020, 6:17pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/4 "2020-09-18T18:17:55Z")
I think the correct cuda function is `syevd`
> **[cuSOLVER :: CUDA Toolkit Documentation](https://docs.nvidia.com/cuda/cusolver/index.html#eig_examples)**
>
> The API reference guide for cuSOLVER, a GPU accelerated library for decompositions and linear system solutions for both dense and sparse matrices.
How to call that from julia though I don’t know.
---
### Author:  [@8me](https://discourse.julialang.org/u/8me)
#### Post date: [September 18, 2020, 10:06pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/5 "2020-09-18T22:06:42Z")
First of all thanks for the help and the hint with `syevd`.
Unfortunately my benchmark gets worse using it:
```julia
A = Symmetric(rand(3000,3000));
A_cuda = CuArray{Float64}(A);
@benchmark CUDA.@sync CUDA.CUSOLVER.syevd!('V','U', A_cuda)
BenchmarkTools.Trial:
memory estimate: 992 bytes
allocs estimate: 20
--------------
minimum time: 15.029 s (0.00% GC)
median time: 15.029 s (0.00% GC)
mean time: 15.029 s (0.00% GC)
maximum time: 15.029 s (0.00% GC)
--------------
samples: 1
evals/sample: 1
```
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [September 21, 2020, 7:05am UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/6 "2020-09-21T07:05:15Z")
The full thing takes 800ms here (also about 500ms using Arpack), so still slower but not that big of a different (while avoiding expensive memory copies).
For the sparse case, isn’t `csreigvsi` the API you need?
```julia
julia> @benchmark CUSOLVER.csreigvsi(dA, rand(T), CUDA.rand(T, 3000), 1e-6, Cint(1000), 'O')
BenchmarkTools.Trial:
memory estimate: 1.47 KiB
allocs estimate: 64
--------------
minimum time: 1.384 ms (0.00% GC)
median time: 2.278 ms (0.00% GC)
mean time: 4.155 ms (0.00% GC)
maximum time: 328.102 ms (0.00% GC)
--------------
samples: 1203
evals/sample: 1
```
(Note that these wrappers are a little rough, and would benefit from a clean-up / higher-level functions.)
---
### Author:  [@8me](https://discourse.julialang.org/u/8me)
#### Post date: [September 21, 2020, 5:28pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/7 "2020-09-21T17:28:36Z")
Hey @maleadt, I think that is what I’m looking for: Thank You!!!
---
### Author:  [@8me](https://discourse.julialang.org/u/8me)
#### Post date: [October 16, 2020, 7:42pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/8 "2020-10-16T19:42:58Z")
Is there a reason why this `CUSOLVER.csreigvsi` just returns one eigenvalue and its associated eigenvector, when I use it exactly your example, @maleadt . I guess there’s some parameter to adjust it? I guess so, because for just getting one eigenvalue there’s also the possibility to use `KrylovKit` and set `howmany` to `1`, which has the following benchmark without GPU usage:
```julia
julia> @benchmark eigsolve(A, 1)
BenchmarkTools.Trial:
memory estimate: 7.92 MiB
allocs estimate: 1432
--------------
minimum time: 15.516 ms (0.00% GC)
median time: 15.844 ms (0.00% GC)
mean time: 16.328 ms (1.43% GC)
maximum time: 21.957 ms (13.21% GC)
--------------
samples: 307
evals/sample: 1
```
---
### Author:  [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)
#### Post date: [November 17, 2021, 1:51pm UTC](https://discourse.julialang.org/t/cuda-eigenvalues-of-a-sparse-matrix/46851/9 "2021-11-17T13:51:20Z")
I need too a few eigenvalues instead than only one. `CUSOLVER.csreigvsi` the exact eigenvalues near my guessing value. However, like `eigs` does, i need the k nearest eigenvalues to the guess, not only one.
Is there a way to do that?
Moreover, what is the rule of `'O'` argument?