Implement matrix logarithm

src/StaticArrays.jl

@@ -134,6 +134,7 @@ include("inv.jl")
 include("solve.jl")
 include("eigen.jl")
 include("expm.jl")
+include("logm.jl")
 include("sqrtm.jl")
 include("lyap.jl")
 include("triangular.jl")

src/logm.jl

@@ -0,0 +1,23 @@
+@inline log(A::StaticMatrix) = _log(Size(A), A)
+
+@inline function _log(::Size{(0,0)}, A::StaticMatrix)
+    T = typeof(log(one(eltype(A))))
+    SMT = similar_type(A,T)
+    return SMT()
+end
+
+@inline function _log(::Size{(1,1)}, A::StaticMatrix)
+    T = typeof(log(one(eltype(A))))
+    SMT = similar_type(A,T)
+    return SMT((log(A[1]), ))
+end
+
+function _log(::Size, A::StaticMatrix)
+    eigA = eigen(A)
+    T = typeof(log(one(eltype(A)) * one(eltype(eigA.values))))
+    VT = similar_type(typeof(eigA.values),T)
+    log_values = log.(VT(eigA.values))
+    log_eigA = Eigen(log_values, eigA.vectors)
+    logA = AbstractMatrix(log_eigA)
+    return logA
+end

src/sqrtm.jl

@@ -1,4 +1,3 @@
-
 @inline sqrt(A::StaticMatrix) = _sqrt(Size(A),A)
 
 @inline function _sqrt(::Size{(0,0)}, A::SA) where {SA<:StaticArray}

test/logm.jl

@@ -0,0 +1,25 @@
+using StaticArrays, Test, LinearAlgebra
+
+@testset "Matrix logarithm" begin
+    @test log(SMatrix{0,0,Int}())::SMatrix === SMatrix{0,0,Float64}()
+    @test log(@SMatrix [2])::SMatrix ≈ SMatrix{1,1}(log(2))
+    @test log(@SMatrix [5 2; -2 1])::SMatrix ≈ log([5 2; -2 1])
+    @test log(@SMatrix [4 2; -2 1])::SMatrix ≈ log([4 2; -2 1])
+    @test log(@SMatrix [4 2; 2 1])::SMatrix ≈ log([4 2; 2 1])
+    # log(::Matrix) is inaccurate on at least Julia 1.1...1.2. We need
+    # to record the result explicitly.
+    # @test log(@SMatrix [-3+1im -1+2im; -4-5im 5-2im])::SMatrix ≈ log(Complex{Float64}[-3+1im -1+2im; -4-5im 5-2im])
+    @test log(@SMatrix [-3+1im -1+2im; -4-5im 5-2im])::SMatrix ≈
+        [1.6949262962402925 + 2.56908641247419im 0.4359384105167865 + 0.47910454835485583im
+         -1.7687979183427713 + 0.558514409317816im 1.7199705725159193 + 0.09415380973694587im]
+    # test for identity matrix
+    @test log(SMatrix{2,2}(I))::SMatrix ≈ zeros(2,2)
+    @test log(@SMatrix Complex{Float64}[1 2 0; 2 1 0; 0 0 1]) ≈ log([1 2 0; 2 1 0; 0 0 1])
+
+    for sz in 0:4, typ in (Float64, Complex{Float64})
+        A = rand(SMatrix{sz,sz,typ})
+        expA = exp(A)
+        logexpA = log(expA)
+        @test logexpA ≈ A
+    end
+end

test/runtests.jl

@@ -56,6 +56,7 @@ addtests("inv.jl")
 addtests("solve.jl")
 addtests("eigen.jl")
 addtests("expm.jl")
+addtests("logm.jl")
 addtests("sqrtm.jl")
 addtests("lyap.jl")
 addtests("lu.jl")