Merge 72fff23 into 554d0aa

src/StaticArrays.jl

@@ -6,8 +6,8 @@ import Base: @pure, @propagate_inbounds, getindex, setindex!, size, similar,
              length, convert, promote_op, map, map!, reduce, reducedim,
              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,
-             sum, diff, prod, count, any, all, sumabs, sumabs2, minimum,
+             fill!, det, inv, eig, eigvals, expm, trace, vecnorm, norm, dot,
+             diagm, sum, diff, prod, count, any, all, sumabs, sumabs2, minimum,
              maximum, extrema, mean, copy, read, read!, write
 
 export StaticScalar, StaticArray, StaticVector, StaticMatrix
@@ -49,6 +49,7 @@ include("deque.jl")
 include("det.jl")
 include("inv.jl")
 include("eigen.jl")
+include("expm.jl")
 include("cholesky.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{S<:Real}(::Size{(2,2)}, A::StaticMatrix{S})
+    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{(3,3)} ≈ expm([1 2 0; 2 1 0; 0 0 1])
+end

test/runtests.jl

@@ -22,6 +22,7 @@ using Base.Test
     include("inv.jl")
     include("solve.jl")
     include("eigen.jl")
+    include("expm.jl")
     include("deque.jl")
     include("io.jl")