diff --git a/stdlib/LinearAlgebra/src/givens.jl b/stdlib/LinearAlgebra/src/givens.jl index 155d8d6f23ce6..d2ca667212adf 100644 --- a/stdlib/LinearAlgebra/src/givens.jl +++ b/stdlib/LinearAlgebra/src/givens.jl @@ -3,6 +3,9 @@ # givensAlgorithm functions are derived from LAPACK, see below abstract type AbstractRotation{T} end +struct AdjointRotation{T,S<:AbstractRotation{T}} <: AbstractRotation{T} + R::S +end transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R)). Consider using adjoint instead of transpose.") @@ -10,16 +13,19 @@ function (*)(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S} TS = typeof(zero(T)*zero(S) + zero(T)*zero(S)) lmul!(convert(AbstractRotation{TS}, R), copy_similar(A, TS)) end -(*)(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR) -(*)(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR) -function _absvecormat_mul_adjrot(A::AbstractVecOrMat{T}, adjR::Adjoint{<:Any,<:AbstractRotation{S}}) where {T,S} - R = adjR.parent +function (*)(adjR::AdjointRotation{T}, A::AbstractVecOrMat{S}) where {T,S} + TS = typeof(zero(T)*zero(S) + zero(T)*zero(S)) + lmul!(convert(AbstractRotation{TS}, adjR.R)', copy_similar(A, TS)) +end +(*)(A::AbstractVector, adjR::AdjointRotation) = _absvecormat_mul_adjrot(A, adjR) +(*)(A::AbstractMatrix, adjR::AdjointRotation) = _absvecormat_mul_adjrot(A, adjR) +function _absvecormat_mul_adjrot(A::AbstractVecOrMat{T}, adjR::AdjointRotation{S}) where {T,S} TS = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - rmul!(TS.(A), convert(AbstractRotation{TS}, R)') + rmul!(copy_similar(A, TS), convert(AbstractRotation{TS}, adjR.R)') end function(*)(A::AbstractMatrix{T}, R::AbstractRotation{S}) where {T,S} TS = typeof(zero(T)*zero(S) + zero(T)*zero(S)) - rmul!(TS.(A), convert(AbstractRotation{TS}, R)) + rmul!(copy_similar(A, TS), convert(AbstractRotation{TS}, R)) end """ @@ -55,12 +61,11 @@ AbstractRotation{T}(G::Givens) where {T} = Givens{T}(G) AbstractRotation{T}(R::Rotation) where {T} = Rotation{T}(R) adjoint(G::Givens) = Givens(G.i1, G.i2, G.c', -G.s) -adjoint(R::Rotation) = Adjoint(R) -function Base.copy(aG::Adjoint{<:Any,<:Givens}) - G = aG.parent - return Givens(G.i1, G.i2, conj(G.c), -G.s) -end -Base.copy(aR::Adjoint{<:Any,Rotation{T}}) where {T} = Rotation{T}(reverse!([r' for r in aR.parent.rotations])) +adjoint(R::AbstractRotation) = AdjointRotation(R) +adjoint(adjR::AdjointRotation) = adjR.R + +Base.copy(aR::AdjointRotation{T,Rotation{T}}) where {T} = + Rotation{T}([r' for r in Iterators.reverse(aR.R.rotations)]) floatmin2(::Type{Float32}) = reinterpret(Float32, 0x26000000) floatmin2(::Type{Float64}) = reinterpret(Float64, 0x21a0000000000000) @@ -291,7 +296,7 @@ function givens(f::T, g::T, i1::Integer, i2::Integer) where T c, s, r = givensAlgorithm(f, g) if i1 > i2 s = -conj(s) - i1,i2 = i2,i1 + i1, i2 = i2, i1 end Givens(i1, i2, c, s), r end @@ -329,9 +334,7 @@ B[i2] = 0 See also [`LinearAlgebra.Givens`](@ref). """ -givens(x::AbstractVector, i1::Integer, i2::Integer) = - givens(x[i1], x[i2], i1, i2) - +givens(x::AbstractVector, i1::Integer, i2::Integer) = givens(x[i1], x[i2], i1, i2) function getindex(G::Givens, i::Integer, j::Integer) if i == j @@ -386,23 +389,24 @@ function lmul!(R::Rotation, A::AbstractMatrix) end return A end -function rmul!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation}) - R = adjR.parent +function rmul!(A::AbstractMatrix, R::Rotation) + @inbounds for i = 1:length(R.rotations) + rmul!(A, R.rotations[i]) + end + return A +end +function lmul!(adjR::AdjointRotation{<:Any,<:Rotation}, A::AbstractMatrix) + R = adjR.R + @inbounds for i = 1:length(R.rotations) + lmul!(adjoint(R.rotations[i]), A) + end + return A +end +function rmul!(A::AbstractMatrix, adjR::AdjointRotation{<:Any,<:Rotation}) + R = adjR.R @inbounds for i = 1:length(R.rotations) rmul!(A, adjoint(R.rotations[i])) end return A end -*(G1::Givens{T}, G2::Givens{T}) where {T} = Rotation(push!(push!(Givens{T}[], G2), G1)) - -# TODO: None of the following disambiguation methods are great. They should perhaps -# instead be MethodErrors, or revised. -# -# disambiguation methods: *(Adj/Trans of AbsVec or AbsMat, Adj of AbstractRotation) -*(A::Adjoint{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B -*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B -*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B -*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B -# disambiguation methods: *(Diag/AbsTri, Adj of AbstractRotation) -*(A::Diagonal, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B) -*(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B) +*(G1::Givens{T}, G2::Givens{T}) where {T} = Rotation([G2, G1]) diff --git a/stdlib/LinearAlgebra/test/givens.jl b/stdlib/LinearAlgebra/test/givens.jl index c1d0caf7b8883..9f23fe4ffaa61 100644 --- a/stdlib/LinearAlgebra/test/givens.jl +++ b/stdlib/LinearAlgebra/test/givens.jl @@ -3,7 +3,7 @@ module TestGivens using Test, LinearAlgebra, Random -using LinearAlgebra: rmul!, lmul!, Givens +using LinearAlgebra: Givens, Rotation # Test givens rotations @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) @@ -14,7 +14,7 @@ using LinearAlgebra: rmul!, lmul!, Givens end @testset for A in (raw_A, view(raw_A, 1:10, 1:10)) Ac = copy(A) - R = LinearAlgebra.Rotation(LinearAlgebra.Givens{elty}[]) + R = Rotation(Givens{elty}[]) for j = 1:8 for i = j+2:10 G, _ = givens(A, j+1, i, j) @@ -25,14 +25,19 @@ using LinearAlgebra: rmul!, lmul!, Givens @test lmul!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10] @testset "transposes" begin - @test G'*G*Matrix(elty(1)I, 10, 10) ≈ Matrix(I, 10, 10) + @test (@inferred G'*G)*Matrix(elty(1)I, 10, 10) ≈ Matrix(I, 10, 10) @test (G*Matrix(elty(1)I, 10, 10))*G' ≈ Matrix(I, 10, 10) - @test copy(R')*(R*Matrix(elty(1)I, 10, 10)) ≈ Matrix(I, 10, 10) + @test (@inferred copy(R'))*(R*Matrix(elty(1)I, 10, 10)) ≈ Matrix(I, 10, 10) @test_throws ErrorException transpose(G) @test_throws ErrorException transpose(R) end end end + @test (R')' === R + @test R * A ≈ (A' * R')' ≈ lmul!(R, copy(A)) + @test A * R ≈ (R' * A')' ≈ rmul!(copy(A), R) + @test R' * A ≈ lmul!(R', copy(A)) + @test A * R' ≈ rmul!(copy(A), R') @test_throws ArgumentError givens(A, 3, 3, 2) @test_throws ArgumentError givens(one(elty),zero(elty),2,2) G, _ = givens(one(elty),zero(elty),11,12)