Sparse matrices generalized eigenvalue problem.

[mseri]

I have two huge sparse matrices A and B generated by FreeFem++.

I'd like to compute what in matlab would be

      [v,lambda]=eigs(A,B,nv,sigma);

(and I need both eigenvalues and relative eigenfunctions)

Sadly with Julia I get:

      ERROR: no method eigs(SparseMatrixCSC{Complex{Float64},Int64}, SparseMatrixCSC{Complex{Float64},Int64}, Int64)
      ERROR: no method eigfact!(SparseMatrixCSC{Complex{Float64},Int64}, SparseMatrixCSC{Complex{Float64},Int64}, Int64, Float64)
       in eigfact at linalg/factorization.jl:543

(and similar with different combinations of arguments).

Is there any way to overcome this problem?

[jiahao]

Looks like we don't currently wrap the ARPACK function to compute the generalized eigenproblem.

[ViralBShah]

I have been planning to do it for a while. Shouldn't be too difficult.

[ViralBShah]

That said, I won't be able to get to it until next week - so anyone who wants to take a crack should probably do it.

[rodrperezi]

Hi, any news on this? I'm having the same problem
Thanks!

[ViralBShah]

Initial support is now in place. I suspect that this is going to be slow, because of allocation of a vector during the matvec, but I haven't actually benchmarked it yet. Do try it out...

julia> a = sprand(100,100,0.3);

julia> b = sprand(100,100,0.3);

julia> d,v = eigs(a,b,sigma=0.1);

julia> norm(a*v[:,2] - d[2]*b*v[:,2])
8.156335076785677e-8

[timholy]

I timed a few iterations of your demo above. Most of the time it completes in 0.2-0.3 seconds, but occasionally it takes 7-8 seconds (without changing a and b). Very curious. I haven't checked your algorithm, is there an element of randomness to it?

[jarlebring]

The initial support test of @ViralBShah didn't check if eigenpairs were computed, and in many cases they are not eigenpairs, unfortunately. Although the residual is almost zero, so are the eigenvectors in many cases:

julia> a = sprand(100,100,0.3);
julia> b = sprand(100,100,0.3);
julia> d,v = eigs(a,b,sigma=0.1);
julia> norm(a*v[:,2] - d[2]*b*v[:,2])
5.417169902738541e-15
julia> norm(v[:,2])
9.350489717877918e-15

Also, the values in d are not eigenvalues:

julia> minimum(svdvals(full(a-d[2]*b)))
0.02991546632563478

I think the origin of this is the same as origin of the bugs reported in the comments of #488. Credits to @meleg for identifying this.

[ViralBShah]

Maybe open a new issue?