Dear all, I am experiencing some problems with eigs, and wondered if anyone had any suggestions for a workaround.
I need to find the smallest positive eigenvalue, lambda, of the generalised eigenvalue problem A*v = lambda*B*v Here, A and B are large, sparse symmetric stiffness matrices from a (fairly involved) finite element model. I have exactly the same implementation in Matlab, and have checked that the matrices A and B I obtain in Matlab match up with the Julia ones. In Matlab, doing [v, lambda] = eigs(A, B, 1, sigma), where 1 specifies the number of eigenvalues, and sigma specifies the shift (I use 1e-12 here, to avoid the zero eigenvalues) returns the "correct" solution. Now, using Julia's eigs, as lambda, v = eigs(A, B, nev=1, sigma=1e-12, ritzvec=true) gives me the error "ERROR: Base.LinAlg.PosDefException(0) in eigs at linalg/arnoldi.jl:38". The issue seems to be that B is positive semi-definite, and is numerically producing some very small negative eigenvalues. Matlab seems able to work around this and obtains sensible results, but Julia is not playing ball with this. Does anyone have an idea for a way around this? Also interesting is that there appears to be some randomness in the Julia implementation. Occasionally eigs doesn't fail and returns a very incorrect answer. Even more occasionally, it returns the correct answer. Many thanks in advance for any insight here, Stef
