added threshold for very small sv

Author: lostella

src/Arpack.jl

@@ -326,20 +326,26 @@ function _svds(X; nsv::Int = 6, ritzvec::Bool = true, tol::Float64 = 0.0, maxite
     # ind   = [1:2:ncv;]
     # sval  = abs.(ex[1][ind])
 
-    svals = sqrt.(max.(zero(real(otype)), real.(ex[1])))
+    realex1 = real.(ex[1])
+    threshold = max(eps(real(otype))*realex1[1], eps(real(otype)))
+    firstzero = findfirst(v -> v <= threshold, realex1)
+    r = firstzero === nothing ? nsv : firstzero-1 # rank of the decomposition
+    realex1[r+1:end] .= zero(real(otype))
+    svals = sqrt.(realex1)
+
     if ritzvec
         # calculating singular vectors
         # left_sv  = sqrt(2) * ex[2][ 1:size(X,1),     ind ] .* sign.(ex[1][ind]')
         if size(X, 1) >= size(X, 2)
             V = ex[2]
             # We cannot assume that X*V is a Matrix even though V is. This is not
             # the case for e.g. LinearMaps.jl so we convert to Matrix explicitly
-            U = Array(qr(rmul!(convert(Matrix, X*V), Diagonal(inv.(svals)))).Q)
+            U = Array(qr(rmul!(convert(Matrix, X*V), Diagonal([inv.(svals[1:r]); ones(nsv-r)]))).Q)
         else
             U = ex[2]
             # We cannot assume that X'U is a Matrix even though U is. This is not
             # the case for e.g. LinearMaps.jl so we convert to Matrix explicitly
-            V = Array(qr(rmul!(convert(Matrix, X'U), Diagonal(inv.(svals)))).Q)
+            V = Array(qr(rmul!(convert(Matrix, X'U), Diagonal([inv.(svals[1:r]); ones(nsv-r)]))).Q)
         end
 
         # right_sv = sqrt(2) * ex[2][ size(X,1)+1:end, ind ]

test/runtests.jl

@@ -295,6 +295,7 @@ LinearAlgebra.adjoint(A::MyOp) = MyOp(adjoint(A.mat))
 end
 
 @testset "low rank" begin
+    Random.seed!(123)
     @testset "$T coefficients" for T in [Float64, Complex{Float64}]
         @testset "rank $r" for r in [2, 5, 10]
             m, n = 3*r, 4*r
@@ -310,8 +311,8 @@ end
             F = svds(A, nsv=nsv)[1]
 
             @test F.S[1:r] ≈ S
-            @test F.U'*F.U ≈ Matrix{Float64}(I, nsv, nsv)
-            @test F.V'*F.V ≈ Matrix{Float64}(I, nsv, nsv)
+            @test F.U'*F.U ≈ Matrix{T}(I, nsv, nsv)
+            @test F.V'*F.V ≈ Matrix{T}(I, nsv, nsv)
         end
     end
 end