eigs for generalized eigenproblems returns zero eigenvectors

The eigs-function returns incorrect eigenpairs in several cases for generalized eigenvalue problems. The returned eigenvectors are almost zero in norm and the returned eigenvalues are wrong.

julia> s=srand(1);
julia> A = sprand(s,100,100,0.3);
julia> B = sprand(s,100,100,0.4);
julia> d,v = eigs(A,B,sigma=0.1,v0=ones(100));
julia> norm(A*v[:,2] - d[2]*B*v[:,2])
2.4912160954221836e-15
julia> norm(v[:,2])
3.3627507194249744e-15
julia> norm(v)
7.970840526504486e-15

Also, the values in d are not eigenvalues:

julia> minimum(svdvals(full(A-d[2]*B)))
0.013397336554077637

I think this is more or less the same problem as the bugs reported in the comments of #24668, i.e. it may be related to requirements on symmetry or positive definiteness of B in ARPACK. ARPACK does indeed try to B-normalize the vectors (as @antoine-levitt points out) and it works in this example when B is posdef:

julia> C=sprandn(s,100,100,0.1);
julia> B=C*C';
julia> d,v = eigs(A,B,sigma=0.1,v0=ones(100));
julia> norm(A*v[:,2] - d[2]*B*v[:,2])
8.777450508400308e-16
julia> norm(v[:,2])
0.40005797301816903
julia> minimum(svdvals(full(A-d[2]*B)))
4.127329501902753e-16
julia> v[:,2]'*B*v[:,2]  # Are eigenvectors B-normalized?
1.0000000000000007 + 0.0im

Computing the norms of the eigenvectors is anyway cheap in comparison to all operations in ARPACK, so one could consider sanity checking that as default behaviour.

It's a standard stable installation:

julia> versioninfo()
Julia Version 0.6.3
Commit d55cadc350 (2018-05-28 20:20 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-4600U CPU @ 2.10GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

JuliaLang / julia

eigs for generalized eigenproblems returns zero eigenvectors #27373