Fix subproblem detection. (#82)

* Fix subproblem detection.

Closes #81

* Some debug info

* Fix computation of first Givens rotation in SVD when n1>1.

Author: andreasnoack

src/svd.jl

@@ -24,7 +24,7 @@ function svdIter!(B::Bidiagonal{T}, n1, n2, shift, U = nothing, Vᴴ = nothing)
         d = B.dv
         e = B.ev
 
-        G, r = givens(d[n1] - abs2(shift)/d[n1], e[n1], 1, 2)
+        G, r = givens(d[n1] - abs2(shift)/d[n1], e[n1], n1, n1 + 1)
         lmul!(G, Vᴴ)
 
         ditmp       = d[n1]
@@ -155,16 +155,16 @@ function __svd!(B::Bidiagonal{T}, U = nothing, Vᴴ = nothing; tol = 100eps(T),
 
             # Search for largest sub-bidiagonal matrix ending at n2
             for _n1 = (n2 - 1):-1:1
-                if abs(e[_n1]) < thresh
+                if _n1 == 1
                     n1 = _n1
                     break
-                elseif _n1 == 1
+                elseif abs(e[_n1 - 1]) < thresh
                     n1 = _n1
                     break
                 end
             end
 
-            debug && println("n1=", n1, ", n2=", n2, ", d[n1]=", d[n1], ", d[n2]=", d[n2], ", e[n1]=", e[n1], " e[n2]=", e[n2-1], " thresh=", thresh)
+            debug && println("n1=", n1, ", n2=", n2, ", d[n1]=", d[n1], ", d[n2]=", d[n2], ", e[n1]=", e[n1], " e[n2-1]=", e[n2-1], " thresh=", thresh)
 
 
             # See LAWN 3 p 18 and

test/svd.jl

@@ -19,6 +19,7 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
 
         F = svd(A)
         @test vals ≈ F.S
+        @show norm(F.U'*A*F.V - Diagonal(F.S), Inf)
         @test F.U'*A*F.V ≈ Diagonal(F.S)
     end
 
@@ -70,4 +71,9 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions
         @test abs.(FA.V'*Float64.(FAb.V))  ≈ I
         @test abs.(FA.V'*Float64.(FAtb.U)) ≈ I
     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)
+    end
 end