Skip to content

Commit 7f59ef7

Browse files
authored
Integer matrix exponentiation in schurpow (#51992)
This improves type-inference by avoiding recursion, as the `A^p` method calls `schurpow` if `p` is a float. After this, the result of `schurpow` is inferred as a small `Union`: ```julia julia> @inferred Union{Matrix{ComplexF64}, Matrix{Float64}} LinearAlgebra.schurpow([1.0 2; 3 4], 2.0) 2×2 Matrix{Float64}: 7.0 10.0 15.0 22.0 ``` One concern here might be that for large `p`, the `A^Int(p)` computation might be expensive by repeated multiplication, as opposed to diagonalization. However, this may only be the case for really large `p`, which may not be commonly encountered. I've added a test, but I'm unsure if `schurpow` is deemed to be an internal function, and this test is unwise. Unfortunately, the return type of `A^p` isn't concretely inferred yet as there are too many possible types that are returned, so I couldn't test for that.
1 parent f3fe7d2 commit 7f59ef7

File tree

2 files changed

+7
-4
lines changed

2 files changed

+7
-4
lines changed

src/dense.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -490,7 +490,7 @@ end
490490
function schurpow(A::AbstractMatrix, p)
491491
if istriu(A)
492492
# Integer part
493-
retmat = A ^ floor(p)
493+
retmat = A ^ floor(Integer, p)
494494
# Real part
495495
if p - floor(p) == 0.5
496496
# special case: A^0.5 === sqrt(A)
@@ -501,7 +501,7 @@ function schurpow(A::AbstractMatrix, p)
501501
else
502502
S,Q,d = Schur{Complex}(schur(A))
503503
# Integer part
504-
R = S ^ floor(p)
504+
R = S ^ floor(Integer, p)
505505
# Real part
506506
if p - floor(p) == 0.5
507507
# special case: A^0.5 === sqrt(A)

test/dense.jl

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1028,8 +1028,8 @@ end
10281028
@test lyap(1.0+2.0im, 3.0+4.0im) == -1.5 - 2.0im
10291029
end
10301030

1031-
@testset "Matrix to real power" for elty in (Float64, ComplexF64)
1032-
# Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014
1031+
@testset "$elty Matrix to real power" for elty in (Float64, ComplexF64)
1032+
# Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014
10331033
#Aa : only positive real eigenvalues
10341034
Aa = convert(Matrix{elty}, [5 4 2 1; 0 1 -1 -1; -1 -1 3 0; 1 1 -1 2])
10351035

@@ -1065,6 +1065,9 @@ end
10651065
@test (A^(2/3))*(A^(1/3)) A
10661066
@test (A^im)^(-im) A
10671067
end
1068+
1069+
Tschurpow = Union{Matrix{real(elty)}, Matrix{complex(elty)}}
1070+
@test (@inferred Tschurpow LinearAlgebra.schurpow(Aa, 2.0)) Aa^2
10681071
end
10691072

10701073
@testset "diagonal integer matrix to real power" begin

0 commit comments

Comments
 (0)