From 2c9694f290bc52059a3d6f451615cf77b91e1111 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 9 Apr 2024 18:28:05 +0200 Subject: [PATCH 1/3] overload `Base.literal_pow` for `AbstractQ` --- stdlib/LinearAlgebra/src/abstractq.jl | 5 +++++ stdlib/LinearAlgebra/test/abstractq.jl | 6 ++++++ 2 files changed, 11 insertions(+) diff --git a/stdlib/LinearAlgebra/src/abstractq.jl b/stdlib/LinearAlgebra/src/abstractq.jl index d6e251e620309..3b77d289734db 100644 --- a/stdlib/LinearAlgebra/src/abstractq.jl +++ b/stdlib/LinearAlgebra/src/abstractq.jl @@ -18,6 +18,11 @@ transpose(Q::AbstractQ{<:Real}) = AdjointQ(Q) transpose(Q::AbstractQ) = error("transpose not implemented for $(typeof(Q)). Consider using adjoint instead of transpose.") adjoint(adjQ::AdjointQ) = adjQ.Q +Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{p}) where {p} = + Q*Base.literal_pow(^, Q, Val(p-1)) # for speed +Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{1}) = Q*I +Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{-1}) = Q'*I # for disambiguation + # promotion with AbstractMatrix, at least for equal eltypes promote_rule(::Type{<:AbstractMatrix{T}}, ::Type{<:AbstractQ{T}}) where {T} = (@inline; Union{AbstractMatrix{T},AbstractQ{T}}) diff --git a/stdlib/LinearAlgebra/test/abstractq.jl b/stdlib/LinearAlgebra/test/abstractq.jl index 19b872d685668..f93c70a8ac223 100644 --- a/stdlib/LinearAlgebra/test/abstractq.jl +++ b/stdlib/LinearAlgebra/test/abstractq.jl @@ -39,6 +39,12 @@ n = 5 @test Q'*I ≈ Q.Q'*I rtol=2eps(real(T)) @test I*Q ≈ Q.Q*I rtol=2eps(real(T)) @test I*Q' ≈ I*Q.Q' rtol=2eps(real(T)) + @test Q^3 ≈ Q*Q*Q + @test Q^2 ≈ Q*Q + @test Q^1 ≈ Q*I + @test Q^(-1) ≈ Q'*I + @test (Q')^(-1) ≈ Q*I + @test (Q')^2 ≈ Q'*Q' @test abs(det(Q)) ≈ 1 @test logabsdet(Q)[1] ≈ 0 atol=2n*eps(real(T)) y = rand(T, n) From 124c010b0feabfbf4af344e5817c4a9c92e8f27d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 9 Apr 2024 18:43:52 +0200 Subject: [PATCH 2/3] fix +/-1 exponent --- stdlib/LinearAlgebra/src/abstractq.jl | 4 ++-- stdlib/LinearAlgebra/test/abstractq.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/abstractq.jl b/stdlib/LinearAlgebra/src/abstractq.jl index 3b77d289734db..66f4bd07b7a16 100644 --- a/stdlib/LinearAlgebra/src/abstractq.jl +++ b/stdlib/LinearAlgebra/src/abstractq.jl @@ -20,8 +20,8 @@ adjoint(adjQ::AdjointQ) = adjQ.Q Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{p}) where {p} = Q*Base.literal_pow(^, Q, Val(p-1)) # for speed -Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{1}) = Q*I -Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{-1}) = Q'*I # for disambiguation +Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{1}) = Q +Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{-1}) = inv(Q) # promotion with AbstractMatrix, at least for equal eltypes promote_rule(::Type{<:AbstractMatrix{T}}, ::Type{<:AbstractQ{T}}) where {T} = diff --git a/stdlib/LinearAlgebra/test/abstractq.jl b/stdlib/LinearAlgebra/test/abstractq.jl index f93c70a8ac223..0eb88324e8c20 100644 --- a/stdlib/LinearAlgebra/test/abstractq.jl +++ b/stdlib/LinearAlgebra/test/abstractq.jl @@ -41,9 +41,9 @@ n = 5 @test I*Q' ≈ I*Q.Q' rtol=2eps(real(T)) @test Q^3 ≈ Q*Q*Q @test Q^2 ≈ Q*Q - @test Q^1 ≈ Q*I - @test Q^(-1) ≈ Q'*I - @test (Q')^(-1) ≈ Q*I + @test Q^1 == Q + @test Q^(-1) == Q' + @test (Q')^(-1) == Q @test (Q')^2 ≈ Q'*Q' @test abs(det(Q)) ≈ 1 @test logabsdet(Q)[1] ≈ 0 atol=2n*eps(real(T)) From 55702bedf15ab27dbc82f9d05464ea034108b750 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 9 Apr 2024 19:26:33 +0200 Subject: [PATCH 3/3] follow the matrix case --- stdlib/LinearAlgebra/src/abstractq.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/abstractq.jl b/stdlib/LinearAlgebra/src/abstractq.jl index 66f4bd07b7a16..bf4064c907a2d 100644 --- a/stdlib/LinearAlgebra/src/abstractq.jl +++ b/stdlib/LinearAlgebra/src/abstractq.jl @@ -18,10 +18,9 @@ transpose(Q::AbstractQ{<:Real}) = AdjointQ(Q) transpose(Q::AbstractQ) = error("transpose not implemented for $(typeof(Q)). Consider using adjoint instead of transpose.") adjoint(adjQ::AdjointQ) = adjQ.Q -Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{p}) where {p} = - Q*Base.literal_pow(^, Q, Val(p-1)) # for speed -Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{1}) = Q -Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{-1}) = inv(Q) +(^)(Q::AbstractQ, p::Integer) = p < 0 ? power_by_squaring(inv(Q), -p) : power_by_squaring(Q, p) +@inline Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{1}) = Q +@inline Base.literal_pow(::typeof(^), Q::AbstractQ, ::Val{-1}) = inv(Q) # promotion with AbstractMatrix, at least for equal eltypes promote_rule(::Type{<:AbstractMatrix{T}}, ::Type{<:AbstractQ{T}}) where {T} =