Type instability for output of matrix exponential

[kshyatt]

julia> using StaticArrays

julia> a = rand(ComplexF64, 4, 4);

julia> sa = SMatrix{4, 4, ComplexF64}(a);

julia> sa^0.24
4×4 Matrix{ComplexF64}:
   0.874908+0.15603im     0.151727+0.0339536im  0.237417+0.0719561im   0.138797-0.157633im
 -0.0531973-0.0687887im   0.840276+0.273785im   0.279529-0.138501im    0.180715-0.042755im
   0.289622+0.114242im   0.0899188-0.226646im    0.85235+0.186083im   0.0758456+0.333684im
   0.139949-0.0276479im  0.0144155-0.0357819im  0.106669+0.0762763im   0.927468+0.233238im

julia> sa^2.0
4×4 SMatrix{4, 4, ComplexF64, 16} with indices SOneTo(4)×SOneTo(4):
 0.229298+3.53602im  0.0232097+2.07716im    -2.15183+2.7979im   -0.885066+2.19474im
 0.629551+2.70017im   0.228516+1.496im      -1.28589+2.23414im  -0.544305+1.95987im
  1.62283+2.20926im   0.673849+1.56577im    0.228108+2.75801im   0.672545+2.01049im
 0.715161+1.50804im   0.143538+0.958167im  -0.334129+1.70432im   0.461188+1.14152im

Presumably the first one is falling back to LinearAlgebra.jl, but not having stable output here is kind of annoying.

[HMegh]

The issue seems to be due to LinearAlgebra.schur, which is used to calculate sa.^p for non-integer p. Here's a workaround based on the ideas in sqrtm.jl:

pow(A::StaticMatrix,p)=_pow(StaticArrays.Size(A),A,p)
_pow(::StaticArrays.Size{S}, A::StaticArray,p) where {S} = SArray{Tuple{S...},ComplexF64,2}(Array(A)^p)

For example:

a = rand(ComplexF64, 4, 4)
sa = SMatrix{4, 4, ComplexF64}(a)


pow(sa,2.4)==sa^2.4 #true
@code_warntype pow(sa,2.4)

Infers the type correctly.

Arguments
  #self#::Core.Const(Main.pow)
  A::SMatrix{4, 4, ComplexF64, 16}
  p::Float64
Body::SMatrix{4, 4, ComplexF64, 16}
1 ─ %1 = Main._pow::Core.Const(Main._pow)
│   %2 = StaticArrays.Size::Core.Const(Size)
│   %3 = (%2)(A)::Core.Const(Size(4, 4))
│   %4 = (%1)(%3, A, p)::SMatrix{4, 4, ComplexF64, 16}
└──      return %4

However, this simplistic approach does not utilize the power of static arrays, and I have not fully tested it.

[kshyatt]

Yeah, I think ideally we would have a fully StaticArray-ified schur, but maybe that is too heavy a lift right now...