Fall back to ExpMethodGeneric for immutable matrices in exponential!

src/exp.jl

@@ -10,7 +10,15 @@ function alloc_mem(A, method)
     return nothing
 end
 
-exponential!(A) = exponential!(A, ExpMethodHigham2005(A));
+# Immutable matrices (e.g. StaticArrays' SMatrix) cannot use the in-place
+# ExpMethodHigham2005 default, so route them to the non-mutating generic method.
+function exponential!(A)
+    return if ismutable(A)
+        exponential!(A, ExpMethodHigham2005(A))
+    else
+        exponential!(A, ExpMethodGeneric())
+    end
+end
 function exponential!(A::GPUArraysCore.AbstractGPUArray)
     return exponential!(A, ExpMethodHigham2005(false))
 end;
@@ -33,6 +41,10 @@ Computes the matrix exponential with the method specified in `method`. The conte
 needed can be preallocated and provided in the parameter `cache` such that the memory can be recycled when calling `exponential!` several times. The preallocation is done with the command [`alloc_mem`](@ref): `cache=alloc_mem(A,method)`.
 `A` may not be sparse matrix type, since exp(A) is likely to be dense.
 
+If no `method` is given, immutable matrices (e.g. StaticArrays' `SMatrix`) are
+computed out-of-place with [`ExpMethodGeneric`](@ref) and the result is returned
+without modifying `A`.
+
 Example
 
 ```julia-repl

test/basictests.jl

@@ -97,6 +97,21 @@ exp_generic(A) = exponential!(copy(A), ExpMethodGeneric())
     end
 end
 
+@testset "exponential! on immutable (static) matrices" begin
+    for (n, T) in ((3, Float64), (4, Float64), (3, Float32), (4, Float32))
+        A = rand(SMatrix{n, n, T})
+        E = exponential!(A)
+        @test E isa SMatrix{n, n}
+        @test E ≈ exp(Matrix(A))
+        @test !any(isnan, E)
+
+        J = ForwardDiff.jacobian(exponential!, A)
+        Jref = ForwardDiff.jacobian(exp_generic, Matrix(A))
+        @test !any(isnan, J)
+        @test J ≈ Jref
+    end
+end
+
 @testset "exponential! sparse" begin
     A = sparse([1, 2, 1], [2, 1, 1], [1.0, 2.0, 3.0])
     @test_throws ErrorException exponential!(A)