Skip to content

RFC: Introduce AdjointRotation to avoid subtyping AbsMat #46233

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Sep 6, 2022
Merged

RFC: Introduce AdjointRotation to avoid subtyping AbsMat #46233

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
34e61ed
Introduce `AdjointRotation` to avoid subtyping `AbsMat`
dkarrasch Jul 27, 2022
3ef7e72
clean up
dkarrasch Jul 31, 2022
fd7b2e7
adjustments
dkarrasch Aug 2, 2022
fce8190
extend tests
dkarrasch Aug 22, 2022
e82be7b
Merge branch 'master' into dk/rotation
dkarrasch Aug 29, 2022
3c7df81
Merge branch 'master' into dk/rotation
dkarrasch Aug 31, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 35 additions & 31 deletions stdlib/LinearAlgebra/src/givens.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,29 @@
# 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.")

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

"""
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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])
13 changes: 9 additions & 4 deletions stdlib/LinearAlgebra/test/givens.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down