RFC: Introduce AdjointRotation to avoid subtyping AbsMat (#46233)

Author: dkarrasch

stdlib/LinearAlgebra/src/givens.jl

@@ -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
 
 """
@@ -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])

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)