From fea18b121b304054591dfe7c37369b4a5e30082a Mon Sep 17 00:00:00 2001
From: Erik Schnetter <schnetter@gmail.com>
Date: Thu, 10 Feb 2022 12:14:50 -0500
Subject: [PATCH 1/3] Implement matrix logarithm

---
 src/StaticArrays.jl |  1 +
 src/logm.jl         | 23 +++++++++++++++++++++++
 src/sqrtm.jl        |  1 -
 test/logm.jl        | 20 ++++++++++++++++++++
 test/runtests.jl    |  1 +
 5 files changed, 45 insertions(+), 1 deletion(-)
 create mode 100644 src/logm.jl
 create mode 100644 test/logm.jl

diff --git a/src/StaticArrays.jl b/src/StaticArrays.jl
index d13ca558..96eaedbd 100644
--- a/src/StaticArrays.jl
+++ b/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")
diff --git a/src/logm.jl b/src/logm.jl
new file mode 100644
index 00000000..2671eaa4
--- /dev/null
+++ b/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
diff --git a/src/sqrtm.jl b/src/sqrtm.jl
index 17078fde..1c7d3bd9 100644
--- a/src/sqrtm.jl
+++ b/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}
diff --git a/test/logm.jl b/test/logm.jl
new file mode 100644
index 00000000..d9230bf1
--- /dev/null
+++ b/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
diff --git a/test/runtests.jl b/test/runtests.jl
index a0b27b0a..38f1e5b1 100644
--- a/test/runtests.jl
+++ b/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")

From 2e6690e9078536caf226ba63a3989b158d14215a Mon Sep 17 00:00:00 2001
From: Erik Schnetter <schnetter@gmail.com>
Date: Mon, 14 Feb 2022 14:21:24 -0500
Subject: [PATCH 2/3] Record test result for logm

---
 test/logm.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/test/logm.jl b/test/logm.jl
index d9230bf1..5befa41d 100644
--- a/test/logm.jl
+++ b/test/logm.jl
@@ -6,7 +6,12 @@ using StaticArrays, Test, LinearAlgebra
     @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])
+    # 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])

From 6f09c13ba7238b3dbec1dc92ec8daf1caad6600c Mon Sep 17 00:00:00 2001
From: Erik Schnetter <schnetter@gmail.com>
Date: Mon, 14 Feb 2022 14:52:57 -0500
Subject: [PATCH 3/3] Test return type of log(SMatrix{0,0})

---
 test/logm.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/logm.jl b/test/logm.jl
index 5befa41d..6701ccf7 100644
--- a/test/logm.jl
+++ b/test/logm.jl
@@ -1,7 +1,7 @@
 using StaticArrays, Test, LinearAlgebra
 
 @testset "Matrix logarithm" begin
-    @test log(SMatrix{0,0,Int}())::SMatrix == SMatrix{0,0,Bool}()
+    @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])