diff --git a/Project.toml b/Project.toml index baffda4..a5ff553 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "GenericLinearAlgebra" uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a" -version = "0.3.17" +version = "0.3.18" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/svd.jl b/src/svd.jl index 8af3949..67259a3 100644 --- a/src/svd.jl +++ b/src/svd.jl @@ -71,7 +71,7 @@ function svdIter!( d[n2] = -G.s * ei1 + G.c * di1 else - throw(ArgumentError("lower bidiagonal version not implemented yet")) + throw(ArgumentError("please convert to upper bidiagonal")) end return B @@ -254,6 +254,7 @@ function __svd!( end function _svdvals!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real} + B = Bidiagonal(B.dv, B.ev, :U) __svd!(B, tol = tol) return sort(abs.(diag(B)), rev = true) end @@ -283,12 +284,31 @@ function _sort_and_adjust!(U, s, Vᴴ) return nothing end +# Convert a lower Bidiagonal to an upper by running a single pass +# of givens rotations from the left. A noop if the Bidiagnoal is +# already upper. +function _lower_to_upper!(B::Bidiagonal, U::Matrix) + if istriu(B) + return B + end + for i in 1:length(B.ev) + G, r = givens(B.dv[i], B.ev[i], i, i + 1) + B.dv[i] = r + B.ev[i] = G.s * B.dv[i + 1] + B.dv[i + 1] = G.c * B.dv[i + 1] + rmul!(U, G') + end + return Bidiagonal(B.dv, B.ev, :U) +end + function _svd!(B::Bidiagonal{T}; tol = eps(T)) where {T<:Real} n = size(B, 1) U = Matrix{T}(I, n, n) Vᴴ = Matrix{T}(I, n, n) + B = _lower_to_upper!(B, U) + __svd!(B, U, Vᴴ, tol = tol) s = diag(B) diff --git a/test/svd.jl b/test/svd.jl index bf211b4..9bcecfe 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -266,4 +266,12 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats @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 + + @testset "Lower Bidiagonal matrices. Issue 157" begin + B = Bidiagonal(randn(4), randn(3), :L) + @test GenericLinearAlgebra._svdvals!(copy(B)) ≈ GenericLinearAlgebra._svdvals!(copy(B')) + U, s, V = GenericLinearAlgebra._svd!(copy(B)) + @test B ≈ U * Diagonal(s) * V' + @test_throws ArgumentError("please convert to upper bidiagonal") GenericLinearAlgebra.svdIter!(B, 1, size(B, 1), 1.0) + end end