Allow exponentiation of Diagonal{*Range}

[ArunS-tack]

Edited diagonal.jl for direct Exponentiation of Diagonal{*Range}. Fixes #40886.

[mcabbott]

Note that this introduces an ambiguity:

julia> dd = Diagonal(rand(100));

julia> @btime $dd .^ 3;
  81.263 ns (3 allocations: 912 bytes)

julia> @btime $dd ^ 3;
  84.845 ns (2 allocations: 1.75 KiB)

julia> @eval LinearAlgebra (^)(D::Diagonal, a::Number) = Diagonal(D.diag .^ a)
^ (generic function with 67 methods)

julia> @btime $dd ^ 3;
ERROR: MethodError: ^(::Diagonal{Float64, Vector{Float64}}, ::Int64) is ambiguous. Candidates:
  ^(A::AbstractMatrix, p::Integer) in LinearAlgebra at /Users/me/.julia/dev/julia/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:445
  ^(A::AbstractMatrix{T}, p::Real) where T in LinearAlgebra at /Users/me/.julia/dev/julia/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:488
  ^(D::Diagonal, a::Number) in LinearAlgebra at REPL[3]:1
Possible fix, define
  ^(::Diagonal, ::Integer)

If you wanted to define something, perhaps it is literal_pow which ought to be forwarded?

julia> @eval LinearAlgebra (^)(D::Diagonal, p::Integer) = Diagonal(D.diag .^ p)
^ (generic function with 68 methods)

julia> @btime $dd ^ 3;
  89.945 ns (1 allocation: 896 bytes)

julia> @eval LinearAlgebra Base.literal_pow(::typeof(^), D::Diagonal, valp::Val) = Diagonal(Base.literal_pow.(^, D.diag, valp))

julia> @btime $dd ^ 3;
  37.726 ns (1 allocation: 896 bytes)

EDIT -- somehow I missed here that the original complaint was about Diagonal(::AbstractRange), and Diagonal(1:3)^3 is indeed an error on master. The signature is (^)(A::AbstractMatrix{T}, p::Integer) where T<:Integer; defining (^)(D::Diagonal, p::Integer) without restricting the eltype of D (as in my @eval) does not seem to introduce ambiguity.

[ArunS-tack]

Note that this introduces an ambiguity:

julia> dd = Diagonal(rand(100));

julia> @btime $dd .^ 3;
  81.263 ns (3 allocations: 912 bytes)

julia> @btime $dd ^ 3;
  84.845 ns (2 allocations: 1.75 KiB)

julia> @eval LinearAlgebra (^)(D::Diagonal, a::Number) = Diagonal(D.diag .^ a)
^ (generic function with 67 methods)

julia> @btime $dd ^ 3;
ERROR: MethodError: ^(::Diagonal{Float64, Vector{Float64}}, ::Int64) is ambiguous. Candidates:
  ^(A::AbstractMatrix, p::Integer) in LinearAlgebra at /Users/me/.julia/dev/julia/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:445
  ^(A::AbstractMatrix{T}, p::Real) where T in LinearAlgebra at /Users/me/.julia/dev/julia/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:488
  ^(D::Diagonal, a::Number) in LinearAlgebra at REPL[3]:1
Possible fix, define
  ^(::Diagonal, ::Integer)

If you wanted to define something, perhaps it is literal_pow which ought to be forwarded?

julia> @eval LinearAlgebra (^)(D::Diagonal, a::Integer) = Diagonal(D.diag .^ a)
^ (generic function with 68 methods)

julia> @btime $dd ^ 3;
  89.945 ns (1 allocation: 896 bytes)

julia> @eval LinearAlgebra Base.literal_pow(::typeof(^), D::Diagonal, valp::Val) = Diagonal(Base.literal_pow.(^, D.diag, valp))

julia> @btime $dd ^ 3;
  37.726 ns (1 allocation: 896 bytes)

Shouldn't it be vectorized dot operation though since its more obvious? Literal_pow multiplies the matrix 3 times but we only need diagonal elements to be cubed.

[dkarrasch]

Literal_pow multiplies the matrix 3 times but we only need diagonal elements to be cubed.

In @mcabbott's suggestion, literal_pow was applied in broadcasted form.

[ArunS-tack]

Literal_pow multiplies the matrix 3 times but we only need diagonal elements to be cubed.

In @mcabbott's suggestion, literal_pow was applied in broadcasted form.

Oh, yes. But currently I m not able to modify the commit as I deleted the repository which was connected to this PR accidentally. Also I can't open another PR for the same as I have another PR pending and I made that directly from the master branch. :/.

[dkarrasch]

You could create a new branch on your new fork, and there make the new changes.

[ArunS-tack]

It shows the commits I have already made in the other PR. If I make another branch and make a new PR, those commits will get added again in the PR.

[dkarrasch]

No worries, I'll take over and prepare a PR based on the comments so far.