From 5cf406457cdb7926a48be84e1b1df5eb0e634e4a Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Fri, 3 Jul 2026 08:59:20 -0400
Subject: [PATCH] Fall back to ExpMethodGeneric for immutable matrices in
 exponential!

exponential!(A) with no method previously always selected the in-place
ExpMethodHigham2005, which throws a setindex! error on immutable inputs
such as StaticArrays' SMatrix. ExpMethodGeneric already has a fully
non-mutating code path for immutable matrices, so route them there.

This also makes ForwardDiff.jacobian(exponential!, ::SMatrix) work,
including for Float32 where StaticArrays' own exp overflows its Pade
coefficients through the 3x3 cofactor solve and returns NaN; the
normalized Pade coefficients of ExpMethodGeneric avoid that overflow.

See https://discourse.julialang.org/t/how-to-use-forwarddiff-jl-with-inv-and-matrix-exponential/137880

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
---
 src/exp.jl         | 14 +++++++++++++-
 test/basictests.jl | 15 +++++++++++++++
 2 files changed, 28 insertions(+), 1 deletion(-)

diff --git a/src/exp.jl b/src/exp.jl
index c917dd0e..35c72404 100644
--- a/src/exp.jl
+++ b/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
diff --git a/test/basictests.jl b/test/basictests.jl
index 539ec2cc..6d2bc4ca 100644
--- a/test/basictests.jl
+++ b/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)