RFC: selection of shift-and-invert mode in eigs

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,18 +13,29 @@ 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")
     else
         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

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)
 

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