Repository: JuliaLinearAlgebra/GenericLinearAlgebra.jl (julia)

Issue:
svd stuck in infinite loop for certain matrices
Using svd on certain matrices leads to an infinite loop:
```julia
using LinearAlgebra, GenericLinearAlgebra
m = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
svd(m) # ok
svd(BigFloat.(m)) # stuck in infinite loop
svdvals!(Float64.(m), debug=true) # also stuck, so it's not BigFloat's fault.
```
The svd of some permutations of the matrix works:
```julia
perm = [2,3,1,4]
svd(BigFloat.(m[perm,perm])) # ok
```
For this particular matrix, any small perturbation (e.g. `m + 1e-10*rand(4,4)`) seem to work fine but there are other matrices, such as 
```julia
m2 = [0.4 0 0 0 0 0 0 0; 0 0.8 0 0 0 0 1 0; 0 0 -0.3 0 0 0 0 0; 0 0 0 -0.4 0 0 0 0; 0 0 0 0 0.4 0 0 0; 0 0 0 0 0 -0.3 0 0; 0 1 0 0 0 0 0.8 0; 0 0 0 0 0 0 0 -0.4]
```
where `svdvals!(m2+1e-2*rand(8,8), debug=true)` often gets stuck.

The problem may be related to the code linked below. If I replace `shift^2 < eps(B.dv[n1]^2)` with `true`, the infinite loop is avoided, but I don't know if this is a proper thing to do.
https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/9f9b46f05df3eb13e3c889c64c9ed2ca72f450bb/src/svd.jl#L181-L188

Pull request:
Found by @rabbott999 in https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/81#issuecomment-1125082565

Closes https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/81

Patch:
diff --git a/Project.toml b/Project.toml
index 0024b1b..8fb2356 100644
--- a/Project.toml
+++ b/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"]
diff --git a/src/svd.jl b/src/svd.jl
index 1a89e04..30718c3 100644
--- a/src/svd.jl
+++ b/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 patch:
diff --git a/test/svd.jl b/test/svd.jl
index 8ebc5f7..6652ca0 100644
--- a/test/svd.jl
+++ b/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