Implement matrix logarithm

Author: eschnett

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{(1,1)}, A::SA) where {SA<:StaticArray}

test/logm.jl

@@ -0,0 +1,20 @@
+using StaticArrays, Test, LinearAlgebra
+
+@testset "Matrix logarithm" begin
+    @test log(SMatrix{0,0,Int}())::SMatrix == SMatrix{0,0,Bool}()
+    @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])
+    @test log(@SMatrix [ -3+1im  -1+2im;-4-5im   5-2im])::SMatrix ≈ log(Complex{Float64}[ -3+1im  -1+2im;-4-5im   5-2im])
+    # 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")