From 351fb25b316400acdbff3d89238f33d8e00fe90a Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Tue, 7 Mar 2023 21:22:22 +0100
Subject: [PATCH] Set off diagonal elements to zero as described in LAWN3 p18

This fixes a convergence issue when the trailing diagonal element
is zero in the bidiagonal input matrix.

Fixes #104
---
 src/svd.jl  | 27 ++++++++++++++++++---------
 test/svd.jl |  8 ++++++++
 2 files changed, 26 insertions(+), 9 deletions(-)

diff --git a/src/svd.jl b/src/svd.jl
index 80c3f3a..b805e69 100644
--- a/src/svd.jl
+++ b/src/svd.jl
@@ -118,21 +118,30 @@ function svdDemmelKahan!(
 end
 
 # Recurrence to estimate smallest singular value from LAWN3 Lemma 1
-function estimate_σ⁻(dv::AbstractVector, ev::AbstractVector, n1::Integer, n2::Integer)
-    λ = abs(dv[n2])
-    B∞ = λ
-    for j = (n2-1):-1:n1
-        λ = abs(dv[j]) * (λ / (λ + abs(ev[j])))
-        B∞ = min(B∞, λ)
-    end
+function estimate_σ⁻!(dv::AbstractVector, ev::AbstractVector, n1::Integer, n2::Integer, tol::Real)
 
+    # (4.3) p 18
     μ = abs(dv[n1])
     B1 = μ
     for j = n1:(n2-1)
         μ = abs(dv[j+1]) * (μ / (μ + abs(ev[j])))
+        if abs(ev[j] / μ) < tol
+            ev[j] = 0
+        end
         B1 = min(B1, μ)
     end
 
+    # (4.4) p 18
+    λ = abs(dv[n2])
+    B∞ = λ
+    for j = (n2-1):-1:n1
+        λ = abs(dv[j]) * (λ / (λ + abs(ev[j])))
+        if abs(ev[j] / λ) < tol
+            ev[j] = 0
+        end
+        B∞ = min(B∞, λ)
+    end
+
     return min(B∞, B1)
 end
 
@@ -154,7 +163,7 @@ function __svd!(
     count = 0
 
     # See LAWN3 page 6 and 22
-    σ⁻ = estimate_σ⁻(d, e, n1, n2)
+    σ⁻ = estimate_σ⁻!(d, e, n1, n2, tol)
     fudge = n
     thresh = tol * σ⁻
 
@@ -208,7 +217,7 @@ function __svd!(
             # current block to determine if it's safe to use shift or if
             # the zero shift algorithm is required to maintain high relative
             # accuracy
-            σ⁻ = estimate_σ⁻(d, e, n1, n2)
+            σ⁻ = estimate_σ⁻!(d, e, n1, n2, tol)
             σ⁺ = max(maximum(view(d, n1:n2)), maximum(view(e, n1:(n2-1))))
 
             if fudge * tol * σ⁻ <= eps(σ⁺)
diff --git a/test/svd.jl b/test/svd.jl
index aff1e9a..b3b6274 100644
--- a/test/svd.jl
+++ b/test/svd.jl
@@ -115,4 +115,12 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
         )
         @test svdvals(A) ≈ svdvals(Complex{Double64}.(A))
     end
+
+    @testset "Issue 104. Trailing zero in bidiagonal." begin
+        dv = [-1.8066303423244812, 0.23626846066361504, -1.8244461746384022, 0.3743075843671794, -1.503025651470883, -0.5273978245088017, 6.194053744695789e-75, -1.4816465601202412e-77, -7.05967042009753e-78, -1.8409609384104132e-78, -3.5760859484965067e-78, 1.7012650564461077e-153, -5.106470534144341e-155, 3.6429789846941095e-155, -3.494481025232055e-232, 0.0]
+        ev = [2.6390728646133144, 1.9554155623906322, 1.9171721320115487, 2.5486042731357257, 1.6188084135207441, -1.2764293576778472, -3.0873284294725004e-77, 1.0815807869027443e-77, -1.0375522364338647e-77, -9.118279619242446e-78, 5.910901980416107e-78, -7.522759136373737e-155, 1.1750163871116314e-154, -2.169544740239464e-155, 2.3352910098001318e-232]
+        B = Bidiagonal(dv, ev, :U)
+        F = GenericLinearAlgebra._svd!(copy(B))
+        @test diag(F.U'*B*F.Vt') ≈ F.S rtol=5e-15
+    end
 end