diff --git a/base/linalg/arnoldi.jl b/base/linalg/arnoldi.jl index 2faf973164f85..8912f1d2c56ae 100644 --- a/base/linalg/arnoldi.jl +++ b/base/linalg/arnoldi.jl @@ -3,7 +3,7 @@ using .ARPACK ## eigs function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM", - tol=0.0, maxiter::Integer=1000, sigma=0,v0::Vector=zeros((0,)), + tol=0.0, maxiter::Integer=1000, sigma=(), v0::Vector=zeros((0,)), ritzvec::Bool=true, complexOP::Bool=false) n = chksquare(A) (n <= 6) && (nev = min(n-1, nev)) @@ -13,6 +13,7 @@ function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM", T = eltype(A) cmplx = T<:Complex bmat = "I" + if !isempty(v0) length(v0)==n || throw(DimensionMismatch("")) eltype(v0)==T || error("Starting vector must have eltype $T") @@ -20,11 +21,21 @@ function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM", v0=zeros(T,(0,)) end - if sigma == 0 + if isempty(sigma) + # normal iteration mode = 1 + sigma = 0 linop(x) = A * x else - C = lufact(A - sigma*eye(A)) + # always find ev closest in magnitude to sigma ... + which = "LM" + if sigma == 0 + # using invert only + C = lufact(A) + else + # using shift and invert + C = lufact(A - sigma*eye(A)) + end if cmplx mode = 3 linop(x) = C\x diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 849fa20701635..420f75ced82cf 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -389,8 +389,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f * ``ritzvec``: Returns the Ritz vectors ``v`` (eigenvectors) if ``true`` * ``op_part``: which part of linear operator to use for real A (:real, :imag) * ``v0``: starting vector from which to start the Arnoldi iteration - ``eigs`` returns the ``nev`` requested eigenvalues in ``d``, the corresponding Ritz vectors ``v`` (only if ``ritzvec=true``), the number of converged eigenvalues ``nconv``, the number of iterations ``niter`` and the number of matrix vector multiplications ``nmult``, as well as the final residual vector ``resid``. - + ``eigs`` returns the ``nev`` requested eigenvalues in ``d``, the corresponding Ritz vectors ``v`` (only if ``ritzvec=true``), the number of converged eigenvalues ``nconv``, the number of iterations ``niter`` and the number of matrix vector multiplications ``nmult``, as well as the final residual vector ``resid``. If ``which="SM"`` inverse iteration is used to find the eigenvalues with smallest magnitude. .. function:: svds(A; nev=6, which="LA", tol=0.0, maxiter=1000, ritzvec=true) diff --git a/test/arpack.jl b/test/arpack.jl index 8d4602396c2ff..dfba34654f9e3 100644 --- a/test/arpack.jl +++ b/test/arpack.jl @@ -1,4 +1,4 @@ -begin + begin local n,a,asym,d,v n = 10 areal = randn(n,n) @@ -24,6 +24,11 @@ begin (d,v) = eigs(apd, nev=3) @test_approx_eq apd*v[:,3] d[3]*v[:,3] # @test_approx_eq eigs(apd; nev=1, sigma=d[3])[1][1] d[3] + + # test (shift-and-)invert mode + (d,v) = eigs(apd, nev=3, sigma=0) + @test_approx_eq apd*v[:,3] d[3]*v[:,3] + end end