Fix typo in estimate_σ⁻ (#88)

Author: andreasnoack

Project.toml

@@ -12,8 +12,9 @@ libblastrampoline_jll = "8e850b90-86db-534c-a0d3-1478176c7d93"
 julia = "1.6"
 
 [extras]
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Quaternions", "Test"]
+test = ["DoubleFloats", "Quaternions", "Test"]

src/svd.jl

@@ -114,7 +114,7 @@ function estimate_σ⁻(dv::AbstractVector, ev::AbstractVector, n1::Integer, n2:
     μ  = abs(dv[n1])
     B1 = μ
     for j in n1:(n2 - 1)
-        μ = abs(dv[j + 1])*(μ*(μ + abs(ev[j])))
+        μ = abs(dv[j + 1])*(μ/(μ + abs(ev[j])))
         B1 = min(B1, μ)
     end
 

test/svd.jl

@@ -1,4 +1,4 @@
-using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
+using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
 
 @testset "Singular value decomposition" begin
     @testset "Problem dimension ($m,$n)" for
@@ -73,7 +73,38 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
     end
 
     @testset "Issue 81" begin
-        m = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
-        @test Float64.(svdvals(big.(m))) ≈ svdvals(m)
+        A = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
+        @test Float64.(svdvals(big.(A))) ≈ svdvals(A)
+
+        A = [
+            0.3   0.0   0.0  0.0  0.0  0.2  0.3   0.0;
+            0.0   0.0   0.0  0.0  0.1  0.0  0.0   0.0;
+            0.0  -0.2   0.0  0.0  0.0  0.0  0.0  -0.2;
+            0.3   0.0   0.0  0.0  0.0  0.2  0.4   0.0;
+            0.0   0.4  -0.2  0.0  0.0  0.0  0.0   0.3;
+            0.2   0.0   0.0  0.0  0.0  0.0  0.2   0.0;
+            0.0   0.0   0.0  0.1  0.0  0.0  0.0   0.0;
+            0.0   0.3  -0.2  0.0  0.0  0.0  0.0   0.3
+        ]
+        @test GenericLinearAlgebra._svdvals!(
+            GenericLinearAlgebra.bidiagonalize!(copy(A)).bidiagonal
+        ) ≈ svdvals(A)
+
+        n = 17
+        A = zeros(Double64, n, n)
+        for j in 1:n, i in 1:n
+            A[i, j] = 1 / Double64(i + j - 1)
+        end
+        @test svdvals(A) ≈ svdvals(Float64.(A))
+
+        # From https://github.com/JuliaMath/DoubleFloats.jl/issues/149
+        n = 64
+        c = Complex{BigFloat}(3//1 + 1im//1)
+        A = diagm(
+            1 => c*ones(BigFloat, n - 1),
+            -1 => c*ones(BigFloat,n - 1),
+            -2 => ones(BigFloat, n - 2)
+        )
+        @test svdvals(A) ≈ svdvals(Complex{Double64}.(A))
     end
 end