expm implementation for sizes >= 3, relies on LU-based ldiv (\)

Author: schmrlng

src/expm.jl

@@ -66,9 +66,58 @@ end
     (newtype)((m11, m21, m12, m22))
 end
 
-# TODO add special case for 3x3 matrices
-if VERSION < v"0.7-"
-    @inline _exp(s::Size, A::StaticArray) = s(Base.expm(Array(A)))
-else
-    @inline _exp(s::Size, A::StaticArray) = s(Base.exp(Array(A)))
+# Adapted from implementation in Base; algorithm from
+# Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008
+function _exp(::Size, A::StaticMatrix{<:Any,<:Any,T}) where T
+    # omitted: matrix balancing, i.e., LAPACK.gebal!
+    nA = maximum(sum(abs.(A), Val{1}))    # marginally more performant than norm(A, 1)
+    ## For sufficiently small nA, use lower order Padé-Approximations
+    if (nA <= 2.1)
+        A2 = A*A
+        if nA > 0.95
+            U = @evalpoly(A2, T(8821612800)*I, T(302702400)*I, T(2162160)*I, T(3960)*I, T(1)*I)
+            U = A*U
+            V = @evalpoly(A2, T(17643225600)*I, T(2075673600)*I, T(30270240)*I, T(110880)*I, T(90)*I)
+        elseif nA > 0.25
+            U = @evalpoly(A2, T(8648640)*I, T(277200)*I, T(1512)*I, T(1)*I)
+            U = A*U
+            V = @evalpoly(A2, T(17297280)*I, T(1995840)*I, T(25200)*I, T(56)*I)
+        elseif nA > 0.015
+            U = @evalpoly(A2, T(15120)*I, T(420)*I, T(1)*I)
+            U = A*U
+            V = @evalpoly(A2, T(30240)*I, T(3360)*I, T(30)*I)
+        else
+            U = @evalpoly(A2, T(60)*I, T(1)*I)
+            U = A*U
+            V = @evalpoly(A2, T(120)*I, T(12)*I)
+        end
+        expA = (V - U) \ (V + U)
+    else
+        s  = log2(nA/5.4)               # power of 2 later reversed by squaring
+        if s > 0
+            si = ceil(Int,s)
+            A = A / T(2^si)
+        end
+
+        A2 = A*A
+        A4 = A2*A2
+        A6 = A2*A4
+
+        U = A6*(T(1)*A6 + T(16380)*A4 + T(40840800)*A2) +
+            (T(33522128640)*A6 + T(10559470521600)*A4 + T(1187353796428800)*A2) +
+            T(32382376266240000)*I
+        U = A*U
+        V = A6*(T(182)*A6 + T(960960)*A4 + T(1323241920)*A2) +
+            (T(670442572800)*A6 + T(129060195264000)*A4 + T(7771770303897600)*A2) +
+            T(64764752532480000)*I
+        expA = (V - U) \ (V + U)
+
+        if s > 0            # squaring to reverse dividing by power of 2
+            for t=1:si
+                expA = expA*expA
+            end
+        end
+    end
+
+    expA
 end

test/expm.jl

@@ -4,5 +4,15 @@
     @test expm(@SMatrix [4 2; -2 1])::SMatrix ≈ expm([4 2; -2 1])
     @test expm(@SMatrix [4 2; 2 1])::SMatrix ≈ expm([4 2; 2 1])
     @test expm(@SMatrix [ -3+1im  -1+2im;-4-5im   5-2im])::SMatrix ≈ expm(Complex{Float64}[ -3+1im  -1+2im;-4-5im   5-2im])
-    @test expm(@SMatrix [1 2 0; 2 1 0; 0 0 1])::SizedArray{Tuple{3,3}} ≈ expm([1 2 0; 2 1 0; 0 0 1])
+    @test expm(@SMatrix [1 2 0; 2 1 0; 0 0 1]) ≈ expm([1 2 0; 2 1 0; 0 0 1])
+
+    for sz in (3,4), typ in (Float64, Complex{Float64})
+        A = rand(typ, sz, sz)
+        nA = norm(A, 1)
+        for nB in (0.005, 0.1, 0.5, 3.0, 20.0)
+            B = A*nB/nA
+            SB = SMatrix{sz,sz,typ}(B)
+            @test expm(B) ≈ expm(SB)
+        end
+    end
 end