Merge pull request #237 from mschauer/ms/expm2

Add complex matrix exponential for 2x2

Author: c42f

src/expm.jl

@@ -40,7 +40,28 @@ end
     (newtype)((m11, m21, m12, m22))
 end
 
-# TODO add complex valued expm
+@inline function _expm(::Size{(2,2)}, A::StaticMatrix{<:Any,<:Any,<:Complex})
+    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]
+
+    z = sqrt((a - d)*(a - d) + 4*b*c )
+    e = exp((a + d - z)/2)
+    f = exp((a + d + z)/2)
+    zr = inv(z)
+
+    m11 = (-e*(a - d - z) + f*(a - d + z)) * zr/2  
+    m12 = (f-e) * b * zr
+    m21 = (f-e) * c * zr
+    m22 = (-e*(-a + d - z) + f*(-a + d + z)) * zr/2
+
+    (newtype)((m11, m21, m12, m22))
+end
+
 # TODO add special case for 3x3 matrices
 
 @inline _expm(s::Size, A::StaticArray) = s(Base.expm(Array(A)))

test/expm.jl

@@ -3,5 +3,6 @@
     @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 [ -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])
 end