Repository: JuliaLinearAlgebra/GenericLinearAlgebra.jl (julia)

Issue:
Deflation breaks edge case
It seems https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/108 broke some edge cases.
Consider for instance:
```julia
julia> rank(BigFloat[0 0; 1 -1])
ERROR: MethodError: reducing over an empty collection is not allowed; consider supplying `init` to the reducer
Stacktrace:
  [1] mapreduce_empty(#unused#::typeof(identity), op::Function, T::Type)
    @ Base ./reduce.jl:367
  [2] reduce_empty(op::Base.MappingRF{typeof(identity), typeof(max)}, #unused#::Type{BigFloat})
    @ Base ./reduce.jl:356
  [3] reduce_empty_iter
    @ ./reduce.jl:379 [inlined]
  [4] mapreduce_empty_iter(f::Function, op::Function, itr::SubArray{BigFloat, 1, Vector{BigFloat}, Tuple{UnitRange{Int64}}, true}, ItrEltype::Base.HasEltype)
    @ Base ./reduce.jl:375
  [5] _mapreduce(f::typeof(identity), op::typeof(max), #unused#::IndexLinear, A::SubArray{BigFloat, 1, Vector{BigFloat}, Tuple{UnitRange{Int64}}, true})
    @ Base ./reduce.jl:427
  [6] _mapreduce_dim
    @ ./reducedim.jl:365 [inlined]
  [7] #mapreduce#765
    @ ./reducedim.jl:357 [inlined]
  [8] mapreduce
    @ ./reducedim.jl:357 [inlined]
  [9] #_maximum#783
    @ ./reducedim.jl:999 [inlined]
 [10] _maximum
    @ ./reducedim.jl:999 [inlined]
 [11] #_maximum#782
    @ ./reducedim.jl:998 [inlined]
 [12] _maximum
    @ ./reducedim.jl:998 [inlined]
 [13] #maximum#780
    @ ./reducedim.jl:994 [inlined]
 [14] maximum
    @ ./reducedim.jl:994 [inlined]
 [15] __svd!(B::Bidiagonal{BigFloat, Vector{BigFloat}}, U::Nothing, Vᴴ::Nothing; tol::BigFloat)
    @ GenericLinearAlgebra ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:224
 [16] _svdvals!(B::Bidiagonal{BigFloat, Vector{BigFloat}}; tol::BigFloat)
    @ GenericLinearAlgebra ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:252
 [17] #svdvals!#33
    @ ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:539 [inlined]
 [18] svdvals!(A::Matrix{BigFloat}; tol::BigFloat)
    @ GenericLinearAlgebra ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:567
 [19] svdvals!
    @ ~/.julia/dev/GenericLinearAlgebra/src/svd.jl:564 [inlined]
 [20] svdvals
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/svd.jl:243 [inlined]
 [21] rank(A::Matrix{BigFloat}; atol::Float64, rtol::BigFloat)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/generic.jl:977
 [22] rank(A::Matrix{BigFloat})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/generic.jl:975
 [23] top-level scope
    @ REPL[13]:1
```
In the code `n1 = 1` and `n2 = 2` then because `d = [-1, 0]`, `n2` is decreased to 1 by
https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/8414b44b7ac3325c5e2233728e6b1cf64a416ecf/src/svd.jl#L212
After that
https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/8414b44b7ac3325c5e2233728e6b1cf64a416ecf/src/svd.jl#L224
fails because the interval `n1:n2-1` is now empty and
https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/8414b44b7ac3325c5e2233728e6b1cf64a416ecf/src/svd.jl#L233
fails because `n2 - 1` is now `0`.

Pull request:
Closes #121

Patch:
diff --git a/src/svd.jl b/src/svd.jl
index d8210fd..098a982 100644
--- a/src/svd.jl
+++ b/src/svd.jl
@@ -167,10 +167,7 @@ function __svd!(
     e = B.ev
     iteration = 0
 
-    # See LAWN3 page 6 and 22
-    σ⁻ = estimate_σ⁻!(d, e, n1, n2, tol)
-    fudge = n
-    thresh = tol * σ⁻
+    thresh = tol * zero(first(d))
 
     if B.uplo === 'U'
         while true
@@ -199,6 +196,7 @@ function __svd!(
                     break
                 end
             end
+            @debug "Active submatrix" iteration n1 n2 d[n1] d[n2] e[n1] e[n2 - 1]
 
             # Deflation check. See LAWN3 p21
             # The Demmel-Kahan iteration moves the zero to the end and produces a
@@ -210,6 +208,9 @@ function __svd!(
 
                     # We have now moved a zero to the end so the problem is one smaller
                     n2 -= 1
+
+                    # We'll start over to find the relevant submatrix after deflation
+                    @goto top
                 end
             end
 
@@ -225,7 +226,7 @@ function __svd!(
             fudge = n2 - n1 + 1
             thresh = tol * σ⁻
 
-            @debug "__svd!" iteration n1 n2 d[n1] d[n2] e[n1] e[n2-1] thresh
+            @debug "estimated quantities" σ⁻ σ⁺ fudge thresh
 
             if fudge * tol * σ⁻ <= eps(σ⁺)
                 svdDemmelKahan!(B, n1, n2, U, Vᴴ)


Test patch:
diff --git a/test/svd.jl b/test/svd.jl
index 8fddc86..a091957 100644
--- a/test/svd.jl
+++ b/test/svd.jl
@@ -248,4 +248,9 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
         @test F.U == I
         @test F.Vt == I
     end
+
+    @testset "Issue 121" begin
+        @test svdvals(BigFloat[0 0; 1 -1]) ≈ [sqrt(2), 0]
+        @test svdvals(BigFloat[1 0 0; 0 0 0; 0 1 -1]) ≈ [sqrt(2), 1, 0]
+    end
 end