Looked a bit more into this. It was great to be able use Gallium for this by the way. I was able to figure out what the issue is on 0.4. We are not checking error codes returned from ARPACK. The check was added in https://github.com/JuliaLang/julia/commit/955af6b16c94189ec70c467dca445d94d617c989#commitcomment-16399154. Therefore, eigs silently returns the wrong result in your first example even though ARPACK returns a non-zero error code. eigs then happens to error out with the right result for your latest example but that must be pure luck. I think that ARPACK really has a problem here.
On the other hand, looking further into this made me read this remark from dsaupd.f c 3. If M can be factored into a Cholesky factorization M = LL` c then Mode = 2 should not be selected. Instead one should use c Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular c linear systems should be solved with L and L` rather c than computing inverses. After convergence, an approximate c eigenvector z of the original problem is recovered by solving c L`z = x where x is a Ritz vector of OP. and we are not following this advice. Right now, we only handle the generalized problem when B(aka M) can be factored so the wrappers need a fix. I'll file an issue about this. If you know somebody who'd like to work on this, please let us know. On Saturday, August 6, 2016 at 2:03:21 AM UTC-4, Madeleine Udell wrote: > > Actually: I'm not sure we should chalk the error up to ARPACK. Julia 0.4 > is able to produce a (correct, I think) answer to > > eigs(A, C) > > but 0.5 gives an ARPACK error. I don't suppose ARPACK changed between > Julia versions...!? > > On Sat, Aug 6, 2016 at 1:54 AM, Madeleine Udell <[email protected] > <javascript:>> wrote: > >> Andreas, thanks for the investigation. I'll use 0.5 for now, and hope the >> real problems I encounter are within the capabilities of ARPACK. >> >> It's embarrassing to be bested by Matlab... >> >> On Fri, Aug 5, 2016 at 9:23 PM, Andreas Noack <[email protected] >> <javascript:>> wrote: >> >>> I've looked a bit into this. I believe there is a bug in the Julia >>> wrappers on 0.4. The good news is that the bug appears to be fixed on 0.5. >>> The bad news is the ARPACK seems to have a hard time with the problem. I get >>> >>> julia> eigs(A,C,nev = 1, which = :LR)[1] >>> ERROR: Base.LinAlg.ARPACKException("unspecified ARPACK error: -9999") >>> in aupd_wrapper(::Type{T}, >>> ::Base.LinAlg.#matvecA!#67{Array{Float64,2}}, ::Base.LinAlg.##64#71{Array >>> {Float64,2}}, ::Base.LinAlg.##65#72, ::Int64, ::Bool, ::Bool, ::String, >>> ::Int64, ::Int64, ::String, : >>> :Float64, ::Int64, ::Int64, ::Array{Float64,1}) at ./linalg/arpack.jl:53 >>> in #_eigs#60(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void, >>> ::Array{Float64,1}, ::Bool, ::B >>> ase.LinAlg.#_eigs, ::Array{Float64,2}, ::Array{Float64,2}) at >>> ./linalg/arnoldi.jl:271 >>> in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs, >>> ::Array{Float64,2}, ::Array{Floa >>> t64,2}) at ./<missing>:0 >>> in #eigs#54(::Array{Any,1}, ::Function, ::Array{Float64,2}, >>> ::Array{Float64,2}) at ./linalg/arnoldi. >>> jl:80 >>> in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, >>> ::Array{Float64,2}, ::Array{Float6 >>> 4,2}) at ./<missing>:0 >>> >>> and since SciPy ends up with the same conclusion I conclude that the >>> issue is ARPACK. Matlab is doing something else because they are able to >>> handle this problem. >>> >>> Given that 0.5 is almost release, I'll not spend more time on the issue >>> on 0.4. Thought, if anybody is able to figure out what is going on, please >>> let us know. >>> >>> On Friday, August 5, 2016 at 8:47:26 AM UTC-4, Madeleine Udell wrote: >>>> >>>> Setting `which=:LR, nev=1` does not return the generalized eigenvalue >>>> with the largest real parts, and does not give a warning or error: >>>> >>>> n = 10 >>>> C = eye(n) >>>> A = zeros(n,n) >>>> A[1] = 100 >>>> A[end] = -100 >>>> @assert eigs(A, C, nev=1, which=:LR)[1][1] == maximum(eigs(A, C)[1]) >>>> >>>> Am I expected to set nev greater than the number of eigenvalues I truly >>>> desire, based on my intuition as a numerical analyst? Or has eigs broken >>>> its implicit guarantee? >>>> >>>> >> >> >> -- >> Madeleine Udell >> Assistant Professor, Operations Research and Information Engineering >> Cornell University >> https://people.orie.cornell.edu/mru8/ >> (415) 729-4115 >> > > > > -- > Madeleine Udell > Assistant Professor, Operations Research and Information Engineering > Cornell University > https://people.orie.cornell.edu/mru8/ > (415) 729-4115 >
