Set off diagonal elements to zero as described in LAWN3 p18 (#105)

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

Fixes #104

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(σ⁺)

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