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]> 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]> 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
