Merge pull request #161 from Godisemo/matrix-exponential-master

Add matrix exponential for 1x1 and 2x2, and fallback to Base

Author: c42f

src/StaticArrays.jl

@@ -8,7 +8,7 @@ import Base: getindex, setindex!, size, similar, vec, show,
              length, convert, promote_op, map, map!, reduce, reducedim, mapreducedim,
              mapreduce, broadcast, broadcast!, conj, transpose, ctranspose,
              hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill,
-             fill!, det, inv, eig, eigvals, trace, vecnorm, norm, dot, diagm,
+             fill!, det, inv, eig, eigvals, expm, trace, vecnorm, norm, dot, diagm,
              sum, diff, prod, count, any, all, sumabs, sumabs2, minimum,
              maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!,
              randexp!, normalize, normalize!, read, read!, write
@@ -90,6 +90,7 @@ include("det.jl")
 include("inv.jl")
 include("solve.jl")
 include("eigen.jl")
+include("expm.jl")
 include("cholesky.jl")
 include("deque.jl")
 include("io.jl")

src/expm.jl

@@ -0,0 +1,46 @@
+@inline expm(A::StaticMatrix) = _expm(Size(A), A)
+
+@inline function _expm(::Size{(1,1)}, A::StaticMatrix)
+    T = typeof(exp(zero(eltype(A))))
+    newtype = similar_type(A,T)
+
+    (newtype)((exp(A[1]), ))
+end
+
+@inline function _expm(::Size{(2,2)}, A::StaticMatrix{<:Any,<:Any,<:Real})
+    T = typeof(exp(zero(eltype(A))))
+    newtype = similar_type(A,T)
+
+    @inbounds a = A[1]
+    @inbounds c = A[2]
+    @inbounds b = A[3]
+    @inbounds d = A[4]
+
+    v = (a-d)^2 + 4*b*c
+
+    if v > 0
+        z = sqrt(v)
+        z1 = cosh(z / 2)
+        z2 = sinh(z / 2) / z
+    elseif v < 0
+        z = sqrt(-v)
+        z1 = cos(z / 2)
+        z2 = sin(z / 2) / z
+    else # if v == 0
+        z1 = T(1.0)
+        z2 = T(0.5)
+    end
+
+    r = exp((a + d) / 2)
+    m11 = r * (z1 + (a - d) * z2)
+    m12 = r * 2b * z2
+    m21 = r * 2c * z2
+    m22 = r * (z1 - (a - d) * z2)
+
+    (newtype)((m11, m21, m12, m22))
+end
+
+# TODO add complex valued expm
+# TODO add special case for 3x3 matrices
+
+@inline _expm(s::Size, A::StaticArray) = s(Base.expm(Array(A)))

test/expm.jl

@@ -0,0 +1,7 @@
+@testset "Matrix exponential" begin
+    @test expm(@SMatrix [2])::SMatrix ≈ SMatrix{1,1}(expm(2))
+    @test expm(@SMatrix [5 2; -2 1])::SMatrix ≈ expm([5 2; -2 1])
+    @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 [1 2 0; 2 1 0; 0 0 1])::SizedArray{Tuple{3,3}} ≈ expm([1 2 0; 2 1 0; 0 0 1])
+end

test/runtests.jl

@@ -25,6 +25,7 @@ using Base.Test
     include("inv.jl")
     include("solve.jl") # Strange inference / world-age error
     include("eigen.jl")
+    include("expm.jl")
     include("chol.jl")
     include("deque.jl")
     include("io.jl")